Typicality of Prethermalization

TL;DR

This paper analytically demonstrates that prethermalization, a prolonged non-thermal state before eventual thermal equilibrium, is a common feature in a broad class of weakly perturbed integrable many-body systems.

Contribution

It extends Deutsch's framework to analytically predict the typicality of prethermalization in general weakly perturbed integrable systems.

Findings

Prethermalization occurs in a wide range of weakly perturbed integrable systems.

Analytical predictions confirm the typicality of prethermalization.

The framework generalizes previous understanding of relaxation dynamics.

Abstract

Prethermalization refers to the remarkable relaxation behavior which an integrable many-body system in the presence of a weak integrability-breaking perturbation may exhibit: After initial transients have died out, it stays for a long time close to some non-thermal steady state, but on even much larger time scales it ultimately switches over to the proper thermal equilibrium behavior. By extending Deutsch's conceptual framework from Phys. Rev. A 43, 2046 (1991), we analytically predict that prethermalization is a typical feature for a very general class of such weakly perturbed systems.

Equations

$$\displaystyle\!\!\!\!\!\!\!\!{\cal A}_{0}(t):=\mbox{Tr}\{\rho_{0}(t)A\}=\sum_{mn}\rho^{0}_{mn}(0)\,A^{0}_{nm}\,e^{i\frac{E^{0}_{n}-E^{0}_{m}}{\hbar}t}\,,\ \ \$$

$${\rho }_{mn}^0(t):={}_0⟨m|{\rho }_0(t)|n{⟩}_0,A_{nm}^0:={}_0⟨n|A|m{⟩}_0,$$

$$H=H_0+V$$

$$U_{mn}:=⟨m|n{⟩}_0,$$

$$\mathcal{A}(t)=\sum _{mn}\sum _{\mu \nu \sigma \tau }{U}_{m\mu }{U}_{n\nu }^{\ast }{U}_{n\sigma }{U}_{m\tau }^{\ast }{\rho }_{\mu \nu }^0(0)A_{\sigma \tau }^0{e}^{i\frac{E_n-E_m}{\hslash }t}.$$

$$\Delta (t):=\mathcal{A}(t)-{\mathcal{A}}_0(t)$$

$$\Delta (t)$$

$${\gamma }_{\tau \mu }(t)$$

$$\overline{|\Delta (t)|}$$

$$f(t)$$

$$Y(t)$$

$$W(t)$$

$$Z(t)$$

$$\displaystyle{\rm{Prob}}\big{(}\,|\Delta(t)|\leq\epsilon\,\Delta_{\!A}\big{)}\geq 1-f(t)/2\epsilon\ ,$$

$$\overline{{\gamma }_{\tau \mu }(t)}$$

$$g_{\mu }(t)$$

$$g(t):=\overline{g_{\mu }(t)}$$

$$g(t)=[g(t){]}^{\ast }.$$

$$f(t)=4\sqrt{1-[g(t){]}^2}.$$

$$g_{\mu }(t)$$

$$h_{\mu }(E)$$

$$h_{\mathrm{BW}}(E):=\frac{1}{2\pi }\frac{\Gamma }{E^2+{\Gamma }^2/4}$$

$$f(t)=4\sqrt{1-e^{-\Gamma |t|/\hslash }}\le 4\sqrt{\Gamma |t|/\hslash }.$$

$$\Gamma :=2\pi {\sigma }_v^2/D$$

$$|{\tilde{\mathcal{A}}}_0(t)-{\mathcal{A}}_0(t)|\le {\Delta }_A|t|\underset{1\le n\le N}{\max }|{ϵ}_n|/\hslash .$$

$${\sigma }_v^2\ll \frac{1}{32\pi }\frac{\hslash D}{t_{\mathrm{rel}}}.$$

$${\sigma }_v\gg D.$$

$$\Delta$$

$$\alpha$$

$$\beta$$

$${\Delta }_A:=a_{\max }-a_{\min }$$

$$\Vert A\Vert ={\Delta }_A/2,$$

$$\alpha$$

$$\beta$$

$$\eta :=\gamma -\bar{\gamma }.$$

$$\alpha$$

Full text

Isolated many-body quantum systems are known to equilibrate, i.e., expectation values exhibit an initial relaxation and then spend most of their time close to a constant value, provided some rather weak preconditions are fulfilled. Furthermore, thermalization is expected for so-called non-integrable systems, i.e., the long-time behavior is well approximated by a microcanonical ensemble. (Possible exceptions, e.g., due to many-body localization, are tacitly ignored here.) In contrast, integrable systems usually exhibit quite significant deviations from such a thermal long-time behavior. All these issues have been extensively explored in the literature, as reviewed, among others, in Refs. gog16 ; ale16 ; bor16 ; mor18 . They are not the subject of our present work but rather will be taken for granted.

Our main issue is the question of how the temporal relaxation of an integrable system changes in response to a weak integrability-breaking perturbation. More specifically, we will derive a rigorous bound for the difference between unperturbed and perturbed expectation values, implying that those changes remain over a long period of time negligibly small for a very large class of weak perturbations. Our approach is conceptually akin to Deutsch’s seminal work on thermalization, treating the perturbations along the lines of random matrix theory deu91 . In particular, we will exploit Deutsch’s result concerning the ultimate thermalization of the perturbed systems. With respect to the unperturbed (integrable) system, we will moreover take for granted that its initial relaxation is not extremely slow, and that it exhibits clearly observable deviations from a thermal long-time behavior. Altogether, we are thus left with a very large class of perturbations with the following quite remarkable property, henceforth named prethermalization: Initially, the perturbed system closely follows the unperturbed relaxation towards a non-thermal steady state, but on even much larger time scales, there must be a clearly visible transition to the ultimate thermal behavior.

Originally, the term prethermalization was introduced by Berges, Borsányi, and Wetterich ber04 for matter under extreme conditions in a quasi-steady state far from equilibrium, which nevertheless exhibits some genuine thermal properties, however without any reference to the concept of integrability. Our present, somewhat different notion of prethermalization has been independently established by Moeckel and Kehrein in Ref. moe08 . During recent years, these and further slightly differing guises of prethermalization have been explored in numerous theoretical div1 ; div2a ; div2b ; ber16 as well as experimental div3 investigations, see also the recent reviews lan16 ; mor18 and further references therein.

Incidentally, the particular examples in Moeckel and Kehrein’s original work moe08 and also in some subsequent studies div2a are beyond the above mentioned realm of our present approach: If the unperturbed system is initially at thermal equilibrium or in the energy ground state, as it is the case in moe08 ; div2a , then the unperturbed dynamics is trivial, and also the signatures of prethermalization after adding a weak perturbation remain too small for our purposes.

Against our treating the perturbations as random matrices (in the unperturbed energy basis), one might object that the “true” perturbation in any concrete physical model is not a random matrix. In particular, the true matrix is often banded deu91 ; fei89 ; fyo96 ; gen12 , i.e., the typical magnitude of its entries decreases with increasing distance from the matrix diagonal. Furthermore, for non-interacting systems perturbed by few-body interactions, the matrix will be very sparse, i.e., only a small fraction of its entries is non-zero bor16 ; bro81 ; fyo96 ; fla97 .

To overcome these concerns, we will consider ensembles of random matrices which can be tailored to emulate the basic features of many concrete models, such as sparsity, bandedness, and other statistical characteristics bor16 ; deu91 ; fei89 ; fyo96 ; gen12 ; bro81 ; fla97 . The true perturbation is thus expected to be contained as one specific matrix in such a properly tailored ensemble as well. (For simplicity, one may imagine matrices of large but finite dimension, whose entries can assume only a finite number of different possible values, as it is the case in any numerical investigation. If each possible value has non-zero probability, there is a finite chance to sample the true matrix from the ensemble.) Hence, if one could prove that some property applies to all members of the ensemble, the property would also apply to the true model. Our main result consists in a slightly weaker statement, namely that the property “prethermalization” at least applies with overwhelming probability when randomly sampling perturbations from the ensemble (“typicality of prethermalization”). It is therefore still very reasonable to expect that the true model is not one of the extremely unlikely exceptions. An illustrative example (spin chain model) is provided in the Supplemental Material sup . Analogous arguments are routinely adopted in random matrix theory, which is well known to be extremely successful in practice bor16 ; bro81 , though its applicability has to our knowledge not been rigorously justified in any concrete physical example. Similar considerations also apply to many other “non-systematic” but practically very well established approximations, such as density functional theory or Boltzmann equations beyond the validity limits of their derivation.

We will demonstrate typicality of prethermalization for a great variety of different ensembles. The resulting total set of all admitted perturbations is therefore extremely large. This seems to us a quite noteworthy finding in itself, independent of the question whether some particular model is covered or not. Moreover, to actually exclude some particular model, it would have to be untypical with respect to every one of those various ensembles. We finally remark that most applications of random matrix theory focus on the ensemble-averaged behavior and take for granted that most individual matrices behave very similarly to the average bor16 ; bro81 . In our present approach, no such extra assumption will be needed.

The unperturbed system is described by a Hamiltonian $H_{0}$ with eigenvalues $E^{0}_{n}$ and eigenvectors $|n\rangle_{\!0}$. The unperturbed evolution of an arbitrary initial state $\rho(0)$ can thus be written as $\rho_{0}(t)=e^{-iH_{0}t/\hbar}\rho(0)e^{iH_{0}t/\hbar}$ and the expectation value of any given observable $A$ as

[TABLE]

where, depending on the specific system under consideration, the indices $m$ and $n$ run from $1$ to infinity or to some finite upper limit.

Likewise, the perturbed system

[TABLE]

exhibits eigenvalues $E_{n}$ and eigenvectors $|n\rangle$. Focusing on the same initial state $\rho(0)$ as before, the expectation value ${\cal A}(t)$ under the perturbed dynamics is then given by the same formulas as in (1) and (2), except that all indices “[math]” must be omitted. In terms of the unitary basis transformation matrix

[TABLE]

this expectation value can be further rewritten as

[TABLE]

The quantity of foremost interest is the difference

[TABLE]

between the perturbed and the unperturbed expectation values. Taking into account $\rho^{0}_{mn}(0)e^{i(E^{0}_{n}-E^{0}_{m})t/\hbar}=\rho^{0}_{mn}(t)$ (see above (1)), it follows with (1) and (5) that

[TABLE]

where $\delta_{mn}$ is the Kronecker delta.

Finally, instead of one particular perturbation $V$ in (3), we consider a statistical ensemble of different $V$’s, and we indicate averages over the ensemble by an overline. This randomization of $V$ is inherited by the Hamiltonian $H$ in (3) and thus by the eigenvalues $E_{n}$, the eigenvectors $|n\rangle$, the $U_{mn}$ in (4), and the $\gamma_{\tau\mu}(t)$ in (8). On the other hand, $H_{0}$, $\rho(0)$, and $A$ are considered as arbitrary but fixed (non-random), hence the same must apply to $E_{n}^{0}$, $|n\rangle_{\!0}$, $\rho_{0}(t)$, and to the matrix elements in (2).

The first main result of our Letter consists in the general rigorous bound

[TABLE]

where $\Delta_{\!A}$ is the measurement range of $A$ (largest minus smallest eigenvalue). The quite tedious derivation has been relegated to the Supplemental Material sup .

Applying Markov’s inequality to (9), it follows for any $\epsilon>0$ that

[TABLE]

where the left hand side denotes the probability that $|\Delta(t)|\leq\epsilon\,\Delta_{\!A}$ when randomly sampling perturbations $V$. For sufficiently small $f(t)$, the difference $\Delta(t)$ in (6) will thus be negligible for the vast majority of all $V$’s.

Our first assumption regarding the so far arbitrary ensemble of $V$’s is as follows: Multiplying the unperturbed energy eigenvectors $|n\rangle_{\!0}$ by arbitrary factors $\sigma_{n}\in\{\pm 1\}$ leaves the $V$ ensemble invariant. Hence also the statistical properties of $\gamma_{\mu\nu}$ in (18) remain unchanged if all the matrix elements $U_{mn}$ in (4) are multiplied by arbitrary factors $\sigma_{n}\in\{\pm 1\}$. As a consequence (see also sup ), the ensemble average of (8) must vanish unless $\tau=\mu$,

[TABLE]

To justify this assumption we note that randomly flipping the signs of the $|n\rangle_{\!0}$ leaves all physical properties unchanged but randomizes the signs of the true perturbation matrix elements $V_{mn}^{0}$. Hence, it is appropriate to adopt a random matrix model with the above invariance property.

As stated in the introduction, the unperturbed system is assumed to exhibit equilibration but not thermalization. Implicitly, this requires a macroscopically well defined system energy; i.e., there must exist a microcanonical energy interval $I_{\mathrm{mc}}:=[E-\delta E,E]$ so that only energies $E^{0}_{n}\in I_{\mathrm{mc}}$ exhibit non-negligible level populations $\rho^{0}_{nn}(0)$. The number of energies $E_{n}^{0}$ contained in $I_{\mathrm{mc}}$ is denoted by $N$ and, without loss of generality, we assume that $n\in\{1,...,N\}$ for all those $E_{n}^{0}$’s. Furthermore, whenever $E_{n}^{0}\not\in I_{\mathrm{mc}}$, we adopt the idealization that $\rho^{0}_{nn}(0)$ is strictly zero f1 . The Cauchy-Schwarz inequality $|\rho^{0}_{mn}(0)|^{2}\leq\rho^{0}_{mm}(0)\rho^{0}_{nn}(0)$ then implies that in (1) only summands with $m,n\in\{1,...,N\}$ actually contribute. As usual, we take for granted that $N$ is huge (exponentially large in the system’s degrees of freedom gol10a ), while the local level density remains close to $D:=\delta E/N$ throughout the interval $I_{\mathrm{mc}}$.

Given that only indices $m,n\in\{1,...,N\}$ actually matter in (1), we can and will assume that their range is extended to arbitrary integer values and that the energies $E_{n}^{0}$ and the matrix elements $V^{0}_{mn}:=\,_{0}\!\langle m|V|n\rangle_{\!0}$ are (re-)defined for arbitrary integers $m,n\not\in\{1,...,N\}$ by way of “extrapolating” in a physically natural way their properties for $m,n\in\{1,...,N\}$.

As a first example, we consider the particularly simple case that $E^{0}_{n+1}-E^{0}_{n}=D$ for all $n$, and that the statistical properties of the matrix elements $V^{0}_{mn}$ do not depend separately on $m$ and $n$, but only on the difference $m-n$. As a consequence (see also sup ), the statistical properties of (16) remain invariant when simultaneously adding an arbitrary integer $\nu$ to all indices on the right hand side (but not on the left hand side). Upon averaging, one can thus infer that $\overline{g_{\mu}(t)}=\overline{g_{\mu+\nu}(t)}$, hence

[TABLE]

is a well-defined ($\mu$-independent) function.

Under the additional assumption that all statistical properties of the diagonal matrix elements $V^{0}_{nn}$ are identical to those of $-V^{0}_{nn}$, one can finally show sup that

[TABLE]

To justify this assumption we note that the diagonal elements of the true $V$ in (3) can always be re-adjusted to vanish on the average. A symmetrization procedure for the remaining distribution will be provided later.

Introducing (15)-(18) into (11)-(13), and taking into account that $\sum_{\nu}\rho^{0}_{\nu\nu}=1$ yields $Y(t)=[g(t)]^{2}$, $Z(t)=g(t)$, and $W(t)=0$, hence (10) takes the form

[TABLE]

One readily infers from (4), (16), and (17) that $g(0)=1$ and that $|g(t)|\leq 1$ for all $t$. Furthermore, it is convenient to rewrite (16) as

[TABLE]

The quantity $F_{\mu}(E):=h_{\mu}(E-E_{\mu}^{0})$ plays a key role in random matrix theory under the name strength function or local spectral density of states fyo96 ; bor16 . Specifically, one finds that the ensemble average $\overline{h_{\mu}(E)}$ is very well approximated by the Breit-Wigner distribution

[TABLE]

under conditions which, together with the concomitant definition of $\Gamma$, will be discussed in more detail shortly. Introducing this result into (20) yields $\overline{g_{\mu}(t)}=e^{-\Gamma\,|t|/2\hbar}$, and with (17), (19) we obtain

[TABLE]

Eqs. (14) and (23) represent our main results. In the remainder of the Letter we focus – as usual in random matrix theory fyo96 ; bor16 ; bro81 – on the case that all $V_{mn}^{0}$ with $m\geq n$ are statistically independent of each other (those with $m<n$ follow from $V^{0}_{nm}=[V^{0}_{mn}]^{\ast}$), that the statistics only depends on $m-n$ (see above (17)), and that $V_{mn}^{0}$ and $-V_{mn}^{0}$ are equally likely (see above Eqs. (15), (18) and sup ).

If all $V_{mn}^{0}$ are furthermore real and Gaussian distributed with variance $\sigma_{\!v}^{2}$, the result (22) with

[TABLE]

was obtained by Deutsch deu91 .e been worked out by Fyodorov et al. in Refs. fyo96 , including distributions with a pronounced delta peak at zero, corresponding to sparse random matrices $V_{mn}^{0}$. In addition, they also admitted the possibility of banded matrices f5 . We have further extended their analytical supersymmetry approach, and moreover performed extensive numerical explorations, showing that the key results (22)-(24) remain valid also for complex $V_{mn}^{0}$’s and under still considerably weaker assumptions regarding their statistics. A few illustrative examples are provided in the Supplemental Material sup .

Multiplying $V$ in (3) by an extra factor $\lambda$ (coupling strength) entails a factor $\lambda^{2}$ in (24), hence the characteristic time scale in (23) decreases as $\lambda^{-2}$, in agreement with previous findings for the persistence of the prethermalized state mor18 ; ber16 . However, we note that our inequality (9) admits strictly speaking no conclusions regarding the actual appearance of non-small differences in (6).

In order to abandon the requirement of equally spaced energies $E_{n}^{0}$ (see above Eq. (17)), let us consider an unperturbed Hamiltonian $\tilde{H}_{0}$ with the same eigenvectors $|n\rangle_{\!0}$ as the original $H_{0}$, but with modified energies $\tilde{E}_{n}^{0}=E_{n}^{0}+\epsilon_{n}$. In view of (1) one anticipates that the corresponding expectation value $\tilde{\cal A}_{0}(t)$ still remains close to ${\cal A}_{0}(t)$ for sufficiently small $\epsilon_{n}$ and not too large $t$. Indeed, it can be rigorously shown sup that

[TABLE]

Taking for granted that the unperturbed Hamiltonian $\tilde{H}_{0}$ exhibits equilibration but not thermalization (see beginning of the Letter), we denote by $t_{\mathrm{rel}}$ its relaxation time; i.e., $\tilde{\cal A}_{0}(t)$ remains very close to some (non-thermal) equilibrium value ${\cal A}_{\!\!\;\rm{eq}}$ for (almost) all $t\geq t_{\mathrm{rel}}$. It follows with (25) that also ${\cal A}_{0}(t)$ exhibits practically the same initial relaxation behavior and then remains close to ${\cal A}_{\!\!\;\rm{eq}}$ for quite some time, provided $|\epsilon_{n}|\ll\hbar/t_{\mathrm{rel}}$ for all $n=1,...,N$. Recalling that the energy level density is exponentially large in the degrees of freedom gol10a , these conclusions must actually apply to rather general non-equidistant energies $\tilde{E}_{n}^{0}$.

Returning to our perturbed systems of the form (3), where the considered ensemble of $V$’s satisfies the rather weak assumptions mentioned above, we can thus conclude from (6), (14), (23), and (24) that also the perturbed expectation values ${\cal A}(t)$ exhibit an initial relaxation and then remain close to ${\cal A}_{\!\!\;\rm{eq}}$ for quite some time f3 , at least for the vast majority of perturbations $V$, and provided they are sufficiently weak so that

[TABLE]

On the other hand, ultimate thermalization for most such $H$’s in (3) has been established in Refs. deu91 ; rei15b . Recalling the considerations at the beginning of our paper, all those “typical” $H$’s thus exhibit prethermalization.

Finally, upon defining modified perturbations $\tilde{V}$ via $\tilde{V}^{0}_{mn}:=V^{0}_{mn}-\delta_{mn}\epsilon_{n}$, we can conclude with Eq. (3) that $\tilde{H}_{0}+\tilde{V}=H$. Hence, the vast majority of those perturbations $\tilde{V}$ of $\tilde{H}_{0}$ must, again, entail prethermalization. In doing so, the $\epsilon_{n}$’s are often expected to be so small that the resulting ensemble of $\tilde{V}$’s is almost identical to the original ensemble of $V$’s (see below). More generally, since the modified energies $\tilde{E}_{n}^{0}$ need no longer be ordered by magnitude, even substantially more general ensembles of $\tilde{V}$’s than of $V$’s are actually admitted, see also sup . Along similar lines, also a possible asymmetry of the $V_{nn}^{0}$ distribution can be removed, as announced below (18).

Next we turn to the question of how far perturbations which satisfy (26) are “weak” in some physically meaningful sense. Quite obviously, such considerations are only possible in terms of non-rigorous arguments and rough estimates.

First of all, typicality of thermalization, as invoked below (26), trivially fails for vanishing perturbations and hence may possibly still fail for extremely weak perturbations ale16 ; bra15 . Yet, a closer inspection of the non-perturbative approach from Refs. deu91 ; rei15b suggests f4 that typicality of thermalization generally does apply provided $\Gamma\gg D$ (cf. Eq. (24)) and thus

[TABLE]

In particular, the diagonal matrix elements $V^{0}_{nn}$ are then typically much larger than the level spacing $D$, thus corroborating the claim below (26) that the ensembles of $\tilde{V}$’s and of $V$’s are often quite similar.

Second, while the unperturbed system $\tilde{H}_{0}$ is assumed not to thermalize for the given initial condition $\rho(0)$, one still expects that it exhibits the usual thermodynamic properties when the system state happens to be the microcanonical ensemble $\rho^{0}_{\mathrm{mc}}:=N^{-1}\sum_{n=1}^{N}|n\rangle_{\!0\,0}\!\langle n|$ corresponding to the energy window $I_{\mathrm{mc}}:=[E-\delta E,E]$ introduced below Eq. (16). Denoting by $\Omega(E)$ the number of energy levels $\tilde{E}_{n}^{0}$ below $E$, by $k_{B}$ Boltzmann’s constant, and by $S(E):=k_{B}\ln[\Omega(E)]$ the entropy, the temperature is thus given by $T(E):=1/S^{\prime}(E)$. Moreover, $\delta E$ must not exceed $k_{B}T(E)$, otherwise the level density would no longer be (approximately) constant throughout $I_{\mathrm{mc}}$ (as assumed below Eq. (16)). It now seems reasonable to say that a perturbation is weak if it does not notably change the thermal equilibrium properties ($S(E)$, $T(E)$, heat capacity, state of matter, etc.) of the unperturbed system. Closer inspection of the approach from Refs. deu91 ; rei15b implies that the perturbations are weak in this sense as long as $\Gamma\ll k_{B}T(E)$ (otherwise, regions with different level densities start to “interact” via the perturbation). According to (24) this amounts to $\sigma_{v}^{2}\ll k_{B}T(E)\,D$. Focusing on the special (largest possible) choice $\delta E=k_{B}T(E)$ and exploiting $D:=\delta E/N$, we arrive at $\sigma_{v}\ll\sqrt{N}D$ and $\sigma_{v}^{2}\ll\delta E\,D$. The first relation complements the lower bound from (27). Since $N$ is exponentially large in the degrees of freedom, the range of admitted $\sigma_{v}$ values is thus still very large. The second relation agrees with (26) if $t_{\mathrm{rel}}$ is comparable to $\hbar/\delta E$. As shown in Ref. bal17 , this is indeed the case for a quite large class of Hamiltonians $\tilde{H}_{0}$, observables $A$, and initial conditions $\rho(0)$.

Alternatively, a perturbation may be considered as weak if the perturbed expectation value ${\cal A}(t)$ remains for (almost) all sufficiently large times $t$ close to the expectation value $\mbox{Tr}\{\rho^{0}_{\mathrm{mc}}A\}$, which the unperturbed system *would * assume in thermal equilibrium. By similar arguments as above, one can see that this alternative weak perturbation criterion is essentially equivalent to the one from the previous paragraph *and * the condition (27).

Altogether, Eq. (26) thus seems to be a physically very natural weak perturbation condition, and it appears reasonable to conjecture that prethermalization will in general be ruled out if (26) is violated. We plan to further pursue this issue in our future work.

In summary, prethermalization has been established for a very large class of integrable (non-thermalizing) Hamiltonians $\tilde{H}_{0}$ and weak perturbations $\tilde{V}$, which closely imitate the essential features of many particular examples of interest in this context. Adopting the common lore of random matrix theory fyo96 ; bor16 ; bro81 ; deu91 , the same conclusion is thus expected to apply “typically” or “with overwhelming likelihood” also to any given such example, unless there is some a priori reason (inappropriate choice of the ensemble, another non-thermalizing system very near-by etc.) why the specific example at hand must be one of the very rare exceptions with respect to every admitted $\tilde{V}$ ensemble gol10a ; gol10b ; rei15b ; rei15 ; bal17 . Remarkably, the same predictions also apply to any other $\tilde{H}_{0}$ which exhibits equilibration but not thermalization, for instance due to many-body localization effects gog16 ; nan15 ; gol17 .

I Derivation of Eqs. (m9)-(m13)

In this section, Eqs. (m9)-(m13) will be derived. Moreover, it will be shown that $Y(t)\in[0,1]$ and $1-Y(t)+W(t)\geq 0$, hence both roots on the right hand side of (m10) are non-negative real numbers.

To begin with, we recall the definitions in Eqs. (m1)-(m8). We also recall that the average over the random perturbations $V$ in Eq. (m3) is indicated by an overline. Accordingly, $H_{0}$, $\rho(0)$, and $A$ are considered as arbitrary but fixed (non-random), hence the same follows for $E_{n}^{0}$, $|n\rangle_{\!0}$, $\rho_{0}(t)$, $\rho^{0}_{mn}(t)$, and $A_{nm}^{0}$. In contrast, $H$, $E_{n}$, $|n\rangle$, $U_{mn}$, and $\gamma_{\tau\mu}(t)$ are random quantities.

Next, we rewrite (m7) as

[TABLE]

where, for notational convenience, all arguments “$t$” and indices “[math]” are temporarily omitted. Note that all factors on the right hand side of (3) and hence $\beta$ on the left hand side are non-random quantities. In contrast, $\alpha$ in (2) is a random quantity due to the random $\gamma$’s on the right hand side. Moreover, its average $\bar{\alpha}$ gives rise to a correlation-type quantity on the right hand side of (2). Roughly speaking, the main idea behind (1) is to later split things into such a correlation-type part and a “rest”.

Denoting by $a_{\mathrm{max}}$ and $a_{\mathrm{min}}$ the largest and smallest eigenvalues of $A$, and considering that the possible outcomes of a measurement are the eigenvalues of the observable $A$, it follows that

[TABLE]

represents the measurement range of $A$. Since all possible outcomes of any realistic measurement must be finite, we can take for granted that $a_{\mathrm{max}}$, $a_{\mathrm{min}}$, and $\Delta_{\!A}$ are finite. Our next observation is that we can add an arbitrary constant to $A$ without changing the main quantity of interest, namely the difference $\Delta(t)$ in (m6). Without loss of generality we thus can and will assume that $a_{\mathrm{min}}=-a_{\mathrm{max}}$, and hence

[TABLE]

where $\|A\|$ denotes the operator norm of $A$ (largest eigenvalue in modulus). While $\Delta$ on the left hand side of (1) thus remains unchanged when adding some constant to $A$, the same does not apply separately to the two terms on the right hand side, as can be seen by closer inspection of (2) and (3). Hence, while our special choice in (5) is irrelevant with respect to $\Delta$, it does matter (and actually is optimized) with respect to the splitting of $\Delta$ into $\alpha$ and $\beta$.

We recall that the $\rho_{\mu\nu}$ in (2) are – according to Eq. (m2) – the matrix elements of $\rho_{0}(0)=\rho(0)$ in the unperturbed basis $|n\rangle_{\!0}$, and likewise for $A_{\sigma\tau}$. In the same vein, $\gamma_{\tau\mu}$ in (m8) may be viewed as the matrix elements ${}_{0}\!\langle\tau|\gamma|\mu\rangle_{\!0}$ of an operator $\gamma$ (put differently, $\gamma$ is uniquely defined via those matrix elements). Note that – in view of (m8) – this operator $\gamma$ may in general be non-Hermitian. Likewise, the Kronecker delta $\delta_{\tau\mu}$ in (3) may be viewed as the matrix elements of the identity operator. As a consequence, one can rewrite (2) and (3) as

[TABLE]

I.1 Evaluation of $\alpha$

In order to further evaluate $\alpha$ in (6), we define the auxiliary (non-Hermitian, random) operator

[TABLE]

Replacing $\gamma$ in (6) by $\eta+\bar{\gamma}$, a straightforward calculation yields

[TABLE]

By rewriting the right hand side of (10) as $\mbox{Tr}\{A(\eta\rho\eta^{\dagger})\}$, observing that $\eta\rho\eta^{\dagger}$ is a non-negative Hermitian operator, and evaluating the trace by means of the eigenbasis of $A$, one finds that $|\alpha_{1}|\leq\|A\|\mbox{Tr}\{\eta\rho\eta^{\dagger}\}$. With (5) we thus arrive at

[TABLE]

where

[TABLE]

is a real-valued, non-negative random variable.

Turning to $\alpha_{2}$ in (11), we note that $\rho$ is a non-negative Hermitian operator, hence there exists a Hermitian operator, which we denote as $\sqrt{\rho}$, with the property that $(\sqrt{\rho})^{2}=\rho$. Exploiting the cyclic invariance of the trace, the right hand side of (11) can thus be rewritten as $\mbox{Tr}\{(\sqrt{\rho}\bar{\gamma}^{\dagger}A)(\eta\sqrt{\rho})\}$. Considering $\mbox{Tr}\{B^{\dagger}C\}$ as a scalar product between two (not necessarily Hermitian) operators $B$ and $C$, the Cauchy-Schwarz inequality takes the form $|\mbox{Tr}\{B^{\dagger}C\}|^{2}\leq\mbox{Tr}\{B^{\dagger}B\}\mbox{Tr}\{C^{\dagger}C\}$. Choosing $B=A\bar{\gamma}\sqrt{\rho}$ and $C=\eta\sqrt{\rho}$ we can thus infer from (11) that

[TABLE]

Exploiting the cyclic invariance of the trace once more, the second trace in (14) can be identified with $X$ in (13). Likewise, the first trace can be identified with $\mbox{Tr}\{A^{2}\bar{\gamma}\rho\bar{\gamma}^{\dagger}\}$. Similarly as above (12), the latter trace can be upper bounded by $\|A\|^{2}Y$, where

[TABLE]

Together with (5), we thus can conclude that

[TABLE]

Next, we introduce $\eta$ from (8) into (13) to obtain

[TABLE]

The last step follows along similar lines as above (6). It follows that

[TABLE]

From (m8), we can infer that

[TABLE]

hence (20) can be rewritten as

[TABLE]

Taking into account the definition (m4), it follows that $\theta_{mn}=\delta_{mn}$, hence (22) takes the form

[TABLE]

Again, (m4) implies that the last sum equals $\delta_{\mu\nu}$, hence $Q_{\mu\nu}=\delta_{\mu\nu}$, and (19) yields

[TABLE]

Technically speaking, the possibility to exactly evaluate the products of four $U$-matrix elements appearing in (19) via (21) is one of the key points of our present approach. In particular, this is the only place where such products (without averaging over the $V$ ensemble) actually appear.

Introducing (25) into (17) and averaging on both sides implies

[TABLE]

Upon comparison with (15), we can conclude that

[TABLE]

Since both $X$ in (13) and $Y$ in (15) must be non-negative real numbers, it follows that

[TABLE]

Observing that (9) implies $|\alpha|\leq|\alpha_{1}|+2|\alpha_{2}|$, we thus can infer from (12), (16), and (28) that

[TABLE]

I.2 Evaluation of $\beta$

By means of the definition

[TABLE]

one readily verifies that

[TABLE]

and with (7) that

[TABLE]

Rewriting (30) as

[TABLE]

it follows, similarly as above (14), that

[TABLE]

and hence

[TABLE]

Similarly as above (15), the first factor on the right hand side of (35) can be upper bounded by $\|A\|^{2}S_{+}$, yielding with (5) the result

[TABLE]

From (32) we can conclude that $|\beta|\leq|R|$ and hence

[TABLE]

Next, we rewrite (36) as

[TABLE]

In the second step in (39), we exploited $Z^{\ast}=\mbox{Tr}\{\rho\bar{\gamma}^{\dagger}\}$, $\mbox{Tr}\{\rho\}=1$, and Eq. (15). It follows that $S_{+}S_{-}=(1+Y)^{2}-(Z+Z^{\ast})^{2}=(1-Y)^{2}+W$, where

[TABLE]

and hence with (27) and (38) that

[TABLE]

We finally remark that $S_{+}S_{-}=\bar{X}^{2}+W$ is a real and non-negative number according to (36). Due to (28) it follows that also $W$ is a real number and satisfies $\bar{X}^{2}+W\geq 0$.

I.3 Evaluation of $\overline{|\Delta|}$

Eq. (1) implies $|\Delta|\leq|\alpha|+|\beta|$ and hence $\overline{|\Delta|}\leq\overline{|\alpha|}+\overline{|\beta|}$. Recalling that $\beta$ is a non-random quantity (see below (3)) yields $\overline{|\beta|}=|\beta|$. With (29) and (42) it then follows that

[TABLE]

Note that since $X$ is a real-valued, non-negative random variable (see below (13)), the same applies to $\sqrt{X}$. Observing that $F(x):=\sqrt{x}$ is a concave function for all $x>0$, we can exploit Jensen’s inequality to obtain $\overline{\sqrt{X}}\leq\sqrt{\bar{X}}$. Moreover, we can conclude from (28) that $\bar{X}\leq\sqrt{\bar{X}}$. Altogether, (43) thus implies

[TABLE]

Finally, we can recast $Y$ from (15) by means of similar arguments as above (6) into the form

[TABLE]

Likewise, $Z$ from (40) takes the form

[TABLE]

Reintroducing the omitted indices “[math]” and arguments “$t$” (see below (3)), we recover from (27), (40), (44)-(46) our final result

[TABLE]

These results are identical to Eqs. (m9)-(m13). We also recall that $Y(t)$ and $W(t)$ are known to be real numbers with $Y(t)\in[0,1]$ (see (28)) and $1-Y(t)+W(t)\geq 0$ (see below (42) and above (44)). Altogether, we thus recover the findings announced at the beginning of this section.

II Random matrix description of a spin chain model

In this section, we illustrate the approximation of a concrete physical model system by the random matrix approach from the main paper. As a particularly simple and common example we consider an (integrable) Heisenberg spin-1/2 chain in the presence of a weak integrability breaking perturbation.

Similarly as in the main paper, the perturbed Hamiltonian is written in two alternative forms,

[TABLE]

where $\tilde{H}_{0}$ and $\tilde{V}$ are the “bare” operators, while $H_{0}$ and $V$ are their “dressed” counterparts. The detailed definition of those operators will be provided below, and the essential idea behind those definitions will be discussed in Sec. II.3.

II.1 The model

As announced, we focus on a perturbed Heisenberg spin-1/2 chain with “bare” operators (cf. Eq. (52))

[TABLE]

where $\sigma^{\alpha}_{l}$ (with $\alpha=x,y,z$ and $l=1,...,L$) are Pauli matrices acting on site $l$ and $\vec{\sigma}_{l}:=(\sigma_{l}^{x},\sigma_{l}^{y},\sigma_{l}^{z})$. The unperturbed Hamiltonian (53) amounts to $L$ spins with nearest neighbor interactions and is well known to be integrable, while the perturbation (54) consists of integrability breaking next-to-nearest neighbor interactions. The units are chosen so that

[TABLE]

and the (in these units) dimensionless parameter $\lambda$ quantifies the perturbation strength.

Following the main paper, the eigenvalues and eigenvectors of $\tilde{H}_{0}$ in (53) are denoted as $\tilde{E}^{0}_{n}$ and $|n\rangle_{\!0}$, and the matrix elements of the perturbation (54) are abbreviated as $\tilde{V}^{0}_{mn}:=\,_{0}\!\langle m|\tilde{V}|n\rangle_{\!0}$. For a spin chain of length $L=13$ and perturbation strength $\lambda=0.2$, the unperturbed eigenvalues $\tilde{E}_{n}^{0}$ and the perturbation matrix $\tilde{V}^{0}_{mn}$ are illustrated in Fig. 1. In this example, the indices have been chosen so that the $\tilde{E}_{n}^{0}$’s are ordered by magnitude, where $n=1,...,2^{L}$ and $2^{L}$ is the Hilbert space dimension of the model (52). Closer inspections shows that it is always possible to choose the eigenvectors $|n\rangle_{\!0}$ so that all matrix elements $\tilde{V}^{0}_{mn}$ are real numbers. Without loss of generality, the eigenvectors have been chosen in this way in all numerical results of this section.

Fig. 1 depicts the so obtained matrix $\tilde{V}^{0}_{mn}$, the diagonal matrix elements $\tilde{V}^{0}_{nn}$, and the energy levels $\tilde{E}^{0}_{n}$ in the first, second, and third rows, respectively. In each case, we show a view of the full Hilbert space in the left panels and a close-up of the middle of the spectrum in the right panels.

The coarse structure in Fig. 1(a) already appears to be qualitatively “random”. Looking at the individual matrix elements in detail in Fig. 1(b), however, reveals correlations that manifest themselves very roughly speaking as blocks of vanishing and non-vanishing entries. Comparing off-diagonal and diagonal matrix elements, we observe that the former are typically one to two orders of magnitude smaller than the latter. For the rest, at least in the middle of the spectrum shown in Fig. 1(d), the diagonal matrix elements essentially look like random numbers.

In Fig. 1(e) and (f), we display the distribution of energy levels $\tilde{E}^{0}_{n}$ by means of the function $\Omega(E)$ from the main paper, which counts the number of $\tilde{E}_{n}^{0}$’s with the property $\tilde{E}_{n}^{0}\leq E$. Put differently,

[TABLE]

where $\Theta(x):=\int_{-\infty}^{x}\delta(y)\,dy$ is the Heaviside step function. The example in Fig. 1(e) nicely illustrates a property, which is generic for many-body systems, and which is taken for granted in the main paper, namely that the energy spectrum gives rise to a well-defined local level density, which may change notably only on scales much larger than the corresponding mean level spacing. (Otherwise, the Boltzmann entropy $S(E):=k_{\mathrm{B}}\ln\Omega(E)$ would not lead to reasonable thermodynamic properties of the system, see end of the main paper.) Here and in the following, we focus on the middle of the spectrum, where the level density is almost constant over a particularly large energy interval, but similar results would be obtained for any other energy interval which is not too large and not too close to the upper and lower ends of the spectrum. (With increasing $L$, those restrictions are expected to become weaker and weaker).

We numerically explored the above model along the same lines as for the dressed model in Sec. II.2 and found that it can not be satisfactorily approximated by our present random matrix approach. The main reason seems to be (see also Sec. III.1 below) that the diagonal matrix elements $\tilde{V}^{0}_{nn}$ in Fig. 1(c) exhibit rather large fluctuations (upon variation of $n$), which are not reflected in the considered random matrix model fyo96 .

II.2 Dressed operators

The main purpose of the “dressed” operators $H_{0}$ and $V$ in (52) is to overcome the above mentioned mismatch between the bare operators and our present random matrix approach. To this end, we first combine the unperturbed levels $\tilde{E}^{0}_{n}$ with the diagonal matrix elements $\tilde{V}^{0}_{nn}$, defining

[TABLE]

Very roughly speaking, the idea behind this definition and the following considerations is to effectively absorb into the unperturbed Hamiltonian those parts of the perturbation $\tilde{V}$ that do not break integrability, thus isolating the integrability-breaking non-diagonal contributions.

Given that the fluctuations of the $\tilde{V}^{0}_{nn}$’s are still much smaller than the energy scale over which the level density exhibits notable variations, and that those fluctuations are essentially unbiased random numbers (see above), one expects that also the $q_{n}$’s in (57) exhibit some well-defined local density, which is moreover practically equal to the local density of the $\tilde{E}_{n}^{0}$’s. On the other hand, since the fluctuations of the $\tilde{V}^{0}_{nn}$ are significantly larger than the mean level spacing $D$, the so-defined $q_{n}$’s are no longer ordered by magnitude. To restore this order, we relabel the eigenvectors $|n\rangle_{\!0}$ so that $q_{n}\leq q_{n+1}$ for all $n$. In all that follows, this modified ordering of the labels $n$ is implicitly understood. (Obviously, the physical properties of the problem at hand do not depend on how we order the $n$’s.)

As in the main paper, we finally require that $H_{0}$ in (52) should exhibit the same eigenvectors $|n\rangle_{\!0}$ as $\tilde{H}_{0}$, whereas the eigenvalues $E_{n}^{0}$ must be equally spaced, with a level spacing $D$ which reproduces the mean level spacing of the $\tilde{E}_{n}^{0}$’s (or, equivalently, the $q_{n}$’s) within the microcanonical energy window $I_{\mathrm{mc}}:=[E-\delta E,E]$ of actual interest. Recalling that $N$ denotes the (large but finite) number of energies $\tilde{E}_{n}^{0}$ within $I_{\mathrm{mc}}$, it follows that $D=\delta E/N$. Accordingly, we set

[TABLE]

with a still arbitrary constant $c$. In order to satisfy the second identity in Eq. (52), the diagonal matrix elements of the dressed perturbation must thus be defined as

[TABLE]

while the off-diagonal matrix elements must remain unchanged, $V^{0}_{mn}:=\tilde{V}^{0}_{mn}$, apart from the above mentioned reordering of the labels $n$.

By way of their above construction, we expect that the dressed operators $H_{0}$ and $V$ now meet the characteristics of generic random matrix ensembles such as those from Ref. fyo96 . This is verified numerically in Fig. 2 for the same model and data as in Fig. 1, but presented in terms of the dressed matrix $V$ instead of the bare $\tilde{V}$. (The distribution of energies $E^{0}_{n}$ as defined in Eq. (58) is rather boring to look at and therefore not plotted.) As the microcanonical energy window, we selected the central $20\%$ of states, and the constant $c$ in Eq. (58) was chosen such that the average of the $V^{0}_{nn}$’s vanishes within this window.

While the coarse view in Fig. 2(a) looks similar to the one from Fig. 1(a), the close-up in Fig. 2(b) shows that also the fine structure of the perturbation now indeed looks very much like a “true” (possibly sparse) random matrix. The diagonal elements $V_{nn}^{0}$ are separately reproduced in Fig. 2(c) and (d). In particular, Fig. 2(d) confirms that the local fluctuations are now much smaller than those of $\tilde{V}_{nn}^{0}$ in Fig. 1(e).

Similarly as when comparing “true” random numbers with numerically generated pseudo-random numbers, it is in general very difficult to quantitatively determine how closely the $V_{mn}^{0}$’s in Fig. 2 emulate a “true” random matrix. We therefore focus on the specific quantity which is at the heart of our present work, namely the (shifted) local density of states (see also (m21))

[TABLE]

where $E_{m}$ and $|m\rangle$ denote (as in the main paper) the eigenvalues and eigenvectors of the perturbed system $H$ in (52).

As detailed in the main paper, the ensemble average of $h_{n}(E)$ is theoretically predicted to be given by

[TABLE]

where $\sigma_{v}^{2}$ is the ensemble averaged variance of the $V^{0}_{mn}$’s (with $m\not=n$ and $|m-n|$ not too large). In our numerics, this prediction is compared with

[TABLE]

where $I_{N}$ represents the set of indices with $E_{n}^{0}\in I_{\mathrm{mc}}$ and $N$ is their number. Furthermore, $\sigma_{v}^{2}$ in (64) is now the numerically determined variance of the $V_{mn}^{0}$’s. Finally, the delta function appearing in (62) is replaced by a pre-delta function

[TABLE]

essentially amounting to a binning of the numerical data with bin width $\epsilon$.

Estimating the ensemble average of $h_{n}(E)$ from (60) by the column average $\hat{h}(E)$ of a single realization in Eq. (64) is common practice in random matrix theory (sometimes referred to as self-averaging or ergodicity property), and is justified by the large number of states $N$ within the microcanonical energy window and the translational invariance of the ensemble of perturbations, meaning that $V^{0}_{mn}$ and $V^{0}_{m+k,n+k}$ exhibit similar statistical properties for all $k$ and thus may be viewed as different samples of the random matrix ensemble. In particular, the column average in Eq. (64) is thus expected to converge for large $N$ towards the ensemble average of $h_{n}(E)$.

While the theoretical predictions strictly speaking apply to infinitely large matrices, the numerically considered matrices are of large but finite dimension $2^{L}=8192$. In order to minimize finite size artifacts, which are known to be particularly pronounced near the borders of the matrices, we choose the interval $I_{N}$ to comprise the central $20\%$ of states as before, such that $N=1638$.

In Fig. 3, we compare the so obtained numerical results for the function $\hat{h}(E)$ with the theory from Eqs. (62) and (63) for three different perturbation strengths $\lambda$ in (54). We emphasize that there are no fit parameters: The width $\Gamma$ of the Breit-Wigner distribution (62) is obtained from the values of $\sigma_{v}^{2}$ (see figure) and $D=0.002$ computed from the actual matrix elements $V^{0}_{mn}$ and energy levels $\tilde{E}^{0}_{n}$. The plots reveal good agreement between numerics and theory, even for the moderate perturbation strength $\lambda=0.3$, where $\Gamma/\delta E\approx 0.14$.

II.3 Concluding remarks

Finally, it may be worthwhile to recall the following key ideas from the main paper regarding the connection of the above considerations with the main objectives of the paper: In general, the two alternative unperturbed Hamiltonians $H_{0}$ and $\tilde{H}_{0}$ in (52) give rise to different time-dependent expectation values. However, the difference of those expectation values can be bounded with the help of our inequality (m25). Closer inspection of how the right hand sides of (m24) and of (m25) numerically scale with the number of spins $L$ in (53) indicates that for sufficiently large $L$ the left hand side of (m25) in fact becomes arbitrarily small (compared to $\Delta_{\!A}$) for all times $t$ with $|t|\leq 1/\Gamma$. Since we are mainly interested in large $L$ and $|t|\leq 1/\Gamma$ in our present work, the two unperturbed Hamiltonians $H_{0}$ and $\tilde{H}_{0}$ thus give rise to expectation values which can be considered as practically indistinguishable. Therefore it is justified to work with the dressed perturbation $V$ in (52) instead of the bare $\tilde{V}$ in our above considerations.

Besides those times with $|t|\leq 1/\Gamma$, also the long-time limit plays an essential role. Specifically, the long-time average of (m1) takes the form

[TABLE]

where $\langle{\cal A}_{0}(t)\rangle_{t}$ indicates the long-time average of ${\cal A}_{0}(t)$, and where we exploited that the spectrum of $H_{0}$ is non-degenerate (see (58)). In particular, this long-time average refers to the dynamics governed by $H_{0}$. But since the eigenvectors of $\tilde{H}_{0}$ and $H_{0}$ are identical, the above result must also be equal to the long-time average when the dynamics is governed by $\tilde{H}_{0}$. (In contrast to $H_{0}$, the spectrum of $\tilde{H}_{0}$ may exhibit degeneracies; in such a case, the eigenvectors $|n\rangle_{\!0}$ are tacitly chosen so that $\rho(0)$ is diagonal within every energy eigenspace.) One thus can conclude that also with respect to the long-time properties it is justified to work with the dressed perturbation $V$ in (52) instead of the bare $\tilde{V}$.

III Numerical illustration of generalized

random matrix ensembles

In this section we provide three illustrative numerical examples, indicating that the analytical predictions from Ref. fyo96 remain valid for considerably more general random matrix ensembles than those actually admitted therein (see also main paper below (m24)). In particular, while in random matrix theory it is usually assumed that all matrix elements of some given ensemble are statistically independent of each other, we will demonstrate that many results, e.g. Eqs. (62) and (63), are still applicable in the presence of various, quite notable correlations. Since such correlations naturally arise in real systems (see also Sec. II), this is an important observation with regard to their modeling by random matrices.

The general framework and the notation is as in the previous section (and in the main paper). In particular, we again compare the analytical prediction (62) with the numerically obtained counterpart (64), i.e., we will numerically estimate the ensemble average of $h_{n}(E)$ from (60) by the column average $\hat{h}(E)$ of a single realization in Eq. (64). While the theory strictly speaking applies to infinitely large matrices, our numerical examples will be of large but finite dimension $M$. In order to minimize the corresponding finite size artifacts, which are known to be particularly pronounced near the borders of the matrices, the interval $I_{N}$ in (64) is chosen so that it remains sufficiently far away from the lower and upper limits at $n=1$ and $n=M$, respectively. As before, the delta function, which enters into (64) via (60), is approximated by the pre-delta function (67), effectively amounting to a binning of the numerical data.

III.1 First example

As our first example, we again start out from the spin model (53)–(54) with $L=13$, hence $M=2^{L}=8192$. But now, we modify the actual perturbation $\tilde{V}$ by manually setting all diagonal elements $\tilde{V}^{0}_{nn}$ equal to zero, or equivalently, we consider dressed operators which no longer satisfy (52) but rather are defined via

[TABLE]

While the resulting Hamiltonian $H_{0}+V$ is no longer expected to faithfully model some “natural” physical system, it still amounts to an instructive “artificial” example: On the one hand, the large fluctuations of the diagonal elements $\tilde{V}^{0}_{nn}$ (see Fig. 1(c)) are eliminated. On the other hand, the remaining off-diagonal elements still exhibit a considerable amount of “structures” or “correlations” (see Figs. 1(b) and 4(a)). Such a perturbation matrix $V^{0}_{mn}$ may thus be considered either as an extremely “unlikely/exceptional/untypical” member of one of the random matrix ensembles admitted in the analytical explorations from fyo96 , or as a “typical” member of an ensemble which is not admitted in those analytical explorations.

Focusing, as in the previous section, on a microcanonical energy window consisting of the central $20\%$ of states, we numerically evaluated $\hat{h}(E)$ from (64), the local density of states $1/D$, as well as the sample variance of perturbation elements to estimate $\sigma_{v}^{2}$. Hence there are again no free parameters in the theory. Comparing $\hat{h}(E)$ to the Breit-Wigner formula (62) with width (63) in Fig. 5, we observe that the theory still describes the numerical estimate quite well. In particular, this substantiates the assertion at the end of Sec. II.1.

III.2 Further examples

In the following two examples, the total Hilbert space dimension is $M=2^{12}=4096$ and the unperturbed system is chosen to exhibit Poissonian level statistics with a mean spacing $D=1/512\approx 0.002$, reflecting the distribution of gaps of typical integrable systems. (Equally spaced levels as in (60) or Wigner-Dyson distributed spacings would lead to practically the same results.) As before, the considered perturbations are quite artificial and are not meant to reflect actual physical systems. Rather, their purpose is to illustrate the generality of the result (62), (63).

In our first artificial ensemble of perturbations, only every $r$th entry starting from the first minor diagonal is nonvanishing. Hence the resulting matrices are sparse with a very regular entry pattern. Moreover, the nonvanishing entries can only take values $+\lambda$ and $-\lambda$, each with probability $\frac{1}{2}$. An example of such a perturbation matrix is shown in Fig. 4(b). The resulting numerical estimate $\hat{h}(E)$ along with the theory prediction is shown in Fig. 6 for $r=4$ and two different choices of $\lambda$. As before, there are no fit parameters involved, since $D$ was fixed a priori and $\sigma_{v}^{2}=\lambda^{2}/r$ by construction. Even though the dimension of the considered matrices is rather small, at least vastly smaller than the size of the Hilbert space for typical many-body systems, theory and numerics agree very well.

For our second artificial example, we uniformly scatter $M_{\mathrm{isl}}$ circular islands of diameter $s_{\mathrm{isl}}$ across the perturbation matrix. Within each island, we alternate entries $+\lambda$ and $-\lambda$ in the radial direction. More precisely, if $(\hat{m}_{k},\hat{n}_{k})$ denotes the (randomly sampled) center of the $k$th island and $r^{(k)}_{mn}:=\sqrt{(m-\hat{m}_{k})^{2}+(n-\hat{n}_{k})^{2}}$ is the distance of entry $(m,n)$ from that center, then we define the contribution to $V^{0}_{mn}$ from the $k$th island as

[TABLE]

where ‘$\mathrm{mod}$’ denotes the modulo operation and $\lfloor x\rfloor=q$ is the greatest integer $q$ such that $q\leq x$ (“floor function”). Islands may overlap, in which case the contributions are simply added, i.e.,

[TABLE]

Hence in this ensemble, the sparsity is irregular (random), but the distribution of nonvanishing entries shows a regular pattern. For a visualization of a single realization of this ensemble, see Fig. 4(c).

In Fig. 7, we display the numerical estimate (64) of the local density of states for a perturbation with $M_{\mathrm{isl}}=2^{14}=16384$ islands of size $s_{\mathrm{isl}}=15$ and for two different perturbation strengths $\lambda=0.007$ and $\lambda=0.01$. The resulting average sparsity $\bar{r}=4.1$ is similar as in the first artificial example, i.e., approximately every fourth matrix element $V^{0}_{mn}$ is nonvanishing. The parameter $\sigma_{v}^{2}$ is computed as the sample variance of all perturbation matrix elements, so that as before there is no free parameter in the theory (62). Again, numerics and theory show good agreement.

We also explored several further examples of artificially “structured” or “correlated” perturbations $V$ along the same lines as for those depicted in Fig. 4; for instance, patterns generated by cellular automata or by concentric rings (“corner arcs”) around the upper left corner of the matrix, etc. We also investigated separately the impact of various “patterns” along the matrix diagonal. Usually, the results were similar to those in Figs. 5-7 and are therefore not reproduced here. Notably different results were only obtained for examples with extremely “little” randomness or extremely strong correlations, for instance, when choosing deterministically alternating signs along the off-diagonals in Fig. 4(b). (In this case, the only random ingredient is the Poissonian level statistics of the unperturbed spectrum.)

IV Derivation of Eqs. (m15)-(m18)

In this section, we deduce Eqs. (m15)-(m18) from the three assumptions above Eqs. (m15), (m17), and (m18).

In doing so, our principal tool is the Taylor expansion formula for functions of noncommuting operators derived in Ref. kum65 . It states that for two arbitrary linear (but not necessarily Hermitian) operators $A,B$ and any complex analytic function $\phi(z)$,

[TABLE]

where $\phi(A)$ is defined via the power series representation of $\phi(z)$, and where $\phi^{(n)}(z)$ denotes the $n$th derivative of $\phi(z)$. Moreover, the operators $C_{n}(A,B)$ are given by the recurrence relation

[TABLE]

with the square brackets denoting the commutator.

A second key point is the observation that $\gamma_{\tau\mu}(t)$ from (m8) can be rewritten in terms of the function $\phi(z):=e^{-izt/\hbar}$ as

[TABLE]

where $t$ is considered as arbitrary but fixed. Using the Taylor expansion formula (72) with $A:=H_{0}$ and $B:=V$, and noting that the above defined $\phi(z)$ is an analytical function with $\phi^{(n)}(z)=(-it/\hbar)^{n}\phi(z)$, we obtain

[TABLE]

Furthermore, (73) now takes the form

[TABLE]

More explicitly, the first few $C_{n}$’s are

[TABLE]

As exemplified by $C_{2}$ and $C_{3}$, the operators $C_{n}$ are thus not necessarily Hermitian. Finally, the definitions (m8) and (m16) imply that $g_{\mu}(t)=\gamma_{\mu\mu}(t)$, and with (75) we thus obtain

[TABLE]

IV.1 Derivation of Eqs. (m15) and (m16)

The goal of this subsection is to show that the assumption above Eq. (m15) implies

[TABLE]

whenever $\mu\not=\tau$. From this result and Eq. (m8), one readily recovers Eqs. (m15) and (m16).

In order to verify (79), it is according to (75) sufficient to show that

[TABLE]

for all $\tau\neq\mu$ and any given $n\in{\mathbb{N}}_{0}$. To this end, we take for granted the assumption above Eq. (m15). In particular, the two unperturbed basis vectors $|\tau\rangle_{\!0}$, $|\mu\rangle_{\!0}$ appearing in the average $\overline{\,{}_{0}\!\langle\tau|C_{n}|\mu\rangle_{\!0}}$ may thus be multiplied by arbitrary factors $\sigma_{\tau},\sigma_{\mu}\in\{\pm 1\}$ without changing the value of that average.

Due to the trivial fact that the operators $C_{n}$ in (76), (77) are independent of the basis vectors $|\tau\rangle_{\!0}$, $|\mu\rangle_{\!0}$ and thus of the factors $\sigma_{\tau},\sigma_{\mu}$, it follows that

[TABLE]

for arbitrary $\sigma_{\tau},\sigma_{\mu}\in\{\pm 1\}$. If $\tau=\mu$, then $\sigma_{\tau}\sigma_{\mu}=\sigma_{\tau}^{2}=1$, hence (81) is always fulfilled. If $\tau\not=\mu$, the factor $\sigma_{\tau}\sigma_{\mu}$ may assume both values $\pm 1$, hence (81) implies (80).

IV.2 Justification of Eq. (m17)

The objective of this subsection is to show that the function $\overline{g_{\mu}(t)}$ is independent of $\mu$ and hence $g(t)$ from Eq. (m17) is well-defined, provided the assumption above Eq. (m17) is fulfilled.

The assumption above Eq. (m17) consist of two parts: The first part requires that $E_{n+1}^{0}-E_{n}^{0}=D$ for all $n$ and thus

[TABLE]

for all $E_{m}^{0},E_{n}^{0}$. The second part requires that all statistical properties of the matrix elements $V^{0}_{mn}$ do not depend separately on $m$ and $n$, but only on the difference $m-n$. In particular, this implies that

[TABLE]

for arbitrary indices $\nu_{0},...,\nu_{k}$ and any $k\in{\mathbb{N}}$, where the right hand side indicates that the left hand side is given by some function $v_{k}$ which only depends on the differences $\nu_{1}-\nu_{0},\nu_{2}-\nu_{1},...,\nu_{k}-\nu_{k-1}$.

In view of Eq. (78), the objective from the beginning of this subsection will be achieved if we can infer from (82) and (83) that $\overline{\,{}_{0}\!\langle\mu|C_{n}|\mu\rangle_{\!0}}$ is independent of $\mu$ for all $n$. The latter will be shown in what follows.

From (76) or (77) it can be seen that the $C_{n}$ are composed of powers of $V$ and commutators of these powers with $H_{0}$. Working in the eigenbasis of $H_{0}$ and inserting complete sets of states between factors of $V$, we find that any diagonal matrix element $\,{}_{0}\!\langle\mu|C_{n}|\mu\rangle_{\!0}$ can be written as a sum of terms of the general form

[TABLE]

with $s_{i}\in\mathbb{N}_{0}$ and with the property that

[TABLE]

The explicit expression of $\,{}_{0}\!\langle\mu|C_{n}|\mu\rangle_{\!0}$ in terms of the functions (84) which satisfy (85) is given in Ref. kum65 , but does not matter for the following arguments.

From Eqs. (83) and (84) we can conclude that

[TABLE]

As usual, the value of the sum on the right hand side does not change when a summation index $\nu_{i}$ is shifted by an arbitrary integer. In particular, we thus may substitute each index $\nu_{i}$ by $\nu_{i}+\mu$ ($i=1,...,k-1$). Exploiting (82) it follows that

[TABLE]

where the right hand side is manifestly independent of $\mu$. Since all $\overline{\,{}_{0}\!\langle\mu|C_{n}|\mu\rangle_{\!0}}$ are sums of such terms (see above), it follows that they must be independent of $\mu$.

IV.3 Derivation of Eq. (m18)

The goal of this subsection is to show that the assumptions above Eqs. (m15), (m17), and (m18) imply

[TABLE]

Since also (m17) is known to apply under those assumptions, Eq. (m18) then readily follows.

In passing we note that (m16) implies $g_{\mu}(-t)=[g_{\mu}(t)]^{\ast}$, and with Eq. (m17) it follows that $g(-t)=[g(t)]^{\ast}$. Together with (88), we thus obtain

[TABLE]

which is a quite interesting result in itself.

In order to verify (88), it is according to (78) sufficient to show that the averaged diagonal matrix elements $\overline{\,{}_{0}\!\langle\mu|C_{n}|\mu\rangle_{\!0}}$ are purely real for even $n$ and purely imaginary for odd $n$. To this end, we recall that $\,{}_{0}\!\langle\mu|C_{n}|\mu\rangle_{\!0}$ is given by a sum of terms of the form (84). Hence, in order to verify (88), it is sufficient to show that

[TABLE]

where $n$ is determined by $k$ and $s_{1},\ldots,s_{k}$ according to (85).

Taking for granted the assumption above (m17), it follows, as in the previous subsection, that Eqs. (82), (83), and (87) apply.

Without loss of generality we can and will choose the energy scale so that $E_{0}^{0}=0$ and hence $E_{-m}^{0}=-E_{m}^{0}$ according to (82). Together with (85) we thus can conclude that

[TABLE]

Furthermore, the assumptions above (m15) and (m18) guarantee that the statistical properties of $-V^{0}_{mn}$ are identical to those of $V^{0}_{mn}$ for all $m,n$, implying that the functions $v_{k}$ from (83) satisfy

[TABLE]

Since $V$ is Hermitian, $V^{0}_{mn}=(V^{0}_{nm})^{*}$, it follows from (83) that

[TABLE]

and with (92) that

[TABLE]

Similarly as below (86) one sees that all summation indices $\nu_{i}$ on the right hand side of Eq. (87) may be replaced by $-\nu_{i}$, yielding

[TABLE]

Upon introducing (91) and exploiting (94) with $\nu_{0}=\nu_{k}=0$, it follows that

[TABLE]

Together with (87) we thus recover Eq. (90).

V Statistics of the random matrix elements $V^{0}_{mn}$

In this section, the statistical properties of the the random matrix elements $V^{0}_{mn}:=\,_{0}\!\langle m|V|n\rangle_{\!0}$ are discussed in more detail, especially those mentioned in the main text below Eq. (m23).

To begin with, we observe that the random distribution of any given matrix element $V_{mn}^{0}$ is captured by the probability density

[TABLE]

Here, it is understood that for random matrix ensembles with purely real elements $V_{mn}^{0}$, only real arguments $v$ are admitted in (97). On the other hand, for ensembles with complex $V_{mn}^{0}$ (which is only possible for $m\not=n$), the arguments $v$ in (97) are understood to be complex for $m\not=n$ and real for $m=n$. We also recall that

[TABLE]

for complex arguments $z:=V_{mn}^{0}-v$ on the right hand side of (97). It follows that

[TABLE]

both for real and complex arguments $z$.

Next, the following two conditions

[TABLE]

will be deduced from the three conditions above Eqs. (m15), (m17), and (m18). These results (100), (101) readily imply the statements below (m23), namely that the statistics of any given matrix element $V_{mn}^{0}$ only depends on $m-n$, and that $V_{mn}^{0}$ and $-V_{mn}^{0}$ are equally likely.

The assumption above Eq. (m15) (see also Sec. IV.1) together with the remarks below (97) readily imply

[TABLE]

for arbitrary $\sigma_{m},\sigma_{n}\in\{\pm 1\}$. If $m=n$, then $\sigma_{m}\sigma_{n}=\sigma_{m}^{2}=1$, hence (102) is always fulfilled. If $m\not=n$, the factor $\sigma_{m}\sigma_{n}$ in (102) may assume both values $\pm 1$. With (97) and (99) it thus follows that (102) is equivalent to

[TABLE]

for all $m>n$. (Unlike in (100), the case $m=n$ is still missing.)

Similarly, from the condition above Eq. (m17) (see also Sec. IV.2) one can readily deduce (101) for all $m\not=n$.

Finally, the condition above Eq. (m18) (see also Sec. IV.3) in combination with (97) and (99) implies

[TABLE]

for all $n$. Hence, (104) together with (103) yields (100).

In conclusion, the three conditions on $V^{0}_{mn}$ above Eqs. (m15), (m17), and (m18) imply the two conditions (100) and (101) on $p(v,m,n)$.

From now on we focus, as in the main text, on the simplest and most common case that all $V_{mn}^{0}$ with $m\geq n$ are statistically independent of each other. More precisely, since $V_{nm}^{0}=(V_{mn}^{0})^{\ast}$, we tacitly assume that $m\geq n$ in the above and the following statements.

As a consequence, all statistical properties of the random matrix are (by definition) fully captured by the distributions of the single elements $V^{0}_{mn}$ in (97). If those distributions satisfy (100) and (101), it is straightforward to invert the above line of reasoning with the final conclusion that the two conditions (100) and (101) imply the three conditions above Eqs. (m15), (m17), and (m18).

In other words, under the additional assumption that the matrix elements $V^{0}_{mn}$ are statistically independent for all $m\geq n$, the conditions (100) and (101) are actually equivalent to the three conditions above Eqs. (m15), (m17), and (m18).

VI Derivation of Eq. (m25)

Throughout this section, we adopt the same notation as in the main paper, except that all indices “[math]” are omitted for notational convenience.

Accordingly, we consider two arbitrary density operators $\rho(t)$ and $\tilde{\rho}(t)$ with identical initial conditions,

[TABLE]

but whose time evolution is governed by two different Hamiltonians, namely

[TABLE]

In other words, the eigenvectors $|n\rangle$ of $H$ and $\tilde{H}$ must be identical, while the eigenvalues $E_{n}$ and $\tilde{E}_{n}$ may be different.

In case the dynamics is governed by $H$, the expectation value of any given observable $A$ at time $t$ can be written, similarly as in Eq. (m1), in the form

[TABLE]

The right hand side of (110) is understood as usual:

[TABLE]

Likewise, if $\tilde{H}$ governs the dynamics, the expectation value of $A$ at time $t$ can be written as

[TABLE]

The last identity in (114) relies on the fact that $H$ from (106) and $\tilde{H}$ from (107) commute. Together with (108)-(111) it follows that $\tilde{\rho}(t)={\cal U}_{t}^{\prime}\rho(t)({\cal U}_{t}^{\prime})^{\dagger}$ and due to the cyclic invariance of the trace that

[TABLE]

where the dependence of $\delta$ on $t$ has been omitted for the sake of simplicity.

The main goal of this section is to show that

[TABLE]

for arbitrary $t$ and $A$, where $\Delta_{\!A}$ is the difference between the largest and smallest eigenvalues of $A$. Eq. (m25) then follows upon taking into account: (i) As mentioned at the beginning of this section, all missing indices “[math]” have to be restored. (ii) As said below Eq. (m16), only $n\in\{1,...,N\}$ actually count in (118).

In order to verify (118), we first evaluate the trace in (116) by means of the eigenbasis of $\rho(t)$, yielding

[TABLE]

where the maximization in (119) is over all normalized vectors $|\psi\rangle$.

For an arbitrary but fixed vector $|\psi\rangle$ of unit norm we can rewrite (120) with (117) as

[TABLE]

With the definition

[TABLE]

we can conclude that

[TABLE]

Eqs. (121), (124), and (126) imply

[TABLE]

From the definition (125) and the Cauchy-Schwarz inequality it follows that

[TABLE]

Since we assumed that $|\psi\rangle$ is normalized, also $|\psi^{\prime}\rangle$ in (122) will be normalized and the last factor in (129) can be upper bounded by $\|A^{2}\|=\|A\|^{2}$, where $\|A\|$ is the operator norm of $A$ (largest eigenvalue in modulus). Exactly the same upper bound can be obtained for $d_{2}$ in (127). With (128) we thus arrive at

[TABLE]

Obviously, $\delta$ in (116) remains unchanged when adding an arbitrary real constant $c$ to $A$. Hence, the inequality (130) with $\|A+c\|$ instead of $\|A\|$ on the right hand side remains valid for arbitrary $c$. The minimum over all $c$ is assumed when the largest and smallest eigenvalues of $A+c$ are of opposite sign and equal modulus, yielding

[TABLE]

where $\Delta_{\!A}$ is the difference between the largest and smallest eigenvalues of $A$.

Rewriting $|\psi\rangle$ as $\sum_{n}c_{n}\,|n\rangle$ with $c_{n}:=\langle n|\psi\rangle$, the normalization of $|\psi\rangle$ takes the form $\sum_{n}|c_{n}|^{2}=1$. Furthermore, we can infer from (106), (107), (115), and (122) that

[TABLE]

and from (123) that

[TABLE]

One readily verifies that $|1-e^{ia}|=2|\sin(a/2)|\leq|a|$ for arbitrary $a\in{\mathbb{R}}$, yielding

[TABLE]

By introducing (135) into (131) we can conclude

[TABLE]

Observing that this bound is independent of $|\psi\rangle$, and taking into account Eqs. (119) and (133), the announced final result (118) is recovered.