Random-matrix model for thermalization

We show that for a system governed by a random-matrix Hamiltonian (a member of the time-reversal invariant Gaussian Orthogonal Ensemble (GOE) of random matrices of dimension N), all functions Tr(Aρ(t)) in the ensemble thermalize: For N→∞ every such function tends to the value Tr(Aρeq(∞))+Tr(Aρ(0))g2(t) . Here ρ(t) is the time-dependent density matrix of the system, A is a Hermitean operator standing for an observable, and ρeq(∞) is the equilibrium density matrix at infinite temperature. The oscillatory function g(t) is the Fourier transform of the average GOE level density and falls off as 1/|t|3/2 for large t. With g(t)=g(−t) , thermalization is symmetric in time. Analogous results, including the symmetry in time of thermalization, are derived for the time-reversal non-invariant Gaussian Unitary Ensemble of random matrices. Comparison with the ‘eigenstate thermalization hypothesis’ of (Srednicki 1999 J. Phys. A: Math. Gen. 32 1163) shows overall agreement but raises significant questions.


I. INTRODUCTION
Does an isolated quantum system thermalize?If so, under which conditions does that happen?Is thermalization linked to time-reversal invariance?These questions are important for the understanding of fundamental aspects of quantum statistical mechanics.Thermalization is seen as the quantum analogue of ergodicity in classical systems, see the reviews in Refs.[2][3][4].Therefore, thermalization continues to be discussed intensely until today [5,6].Thermalization cannot occur in quantum systems that are integrable or show manybody localization, and the discussion is focused on chaotic quantum systems.In an isolated system, thermalization cannot be investigated experimentally by definition.In ultracold atoms [5], for instance, the difficulty is circumvented by opening the system after some time and taking data.
In Refs.[2,4], thermalization is defined in terms of the time-dependent function Tr(Aρ(t)).Here ρ(t) is the timedependent density matrix of an isolated quantum system, and A is a Hermitean operator representing an observable used to test thermalization.The system is said to thermalize if for t → ∞, Tr(Aρ(t)) tends asymptotically to Tr(Aρ eq ) where ρ eq is the time-independent density matrix of statistical equilibrium.Thus, thermalization provides a test of the general hypothesis that in the long-time limit, expectation values of operators tend to their equilibrium values.
A central role in the discussion plays the "eigenstate thermalization hypothesis" introduced by Srednicki [1].In the basis of eigenstates α, β, . . . of the chaotic Hamiltonian with energy eigenvalues E α , E β , . .., the operator A is represented by the matrix A αβ .The hypothesis postulates that A αβ has the form Here A(E) and F (E, ω) are smooth functions of E = (1/2)(E α + E β ) and ω = E α − E β , S(E) is the thermodynamic entropy at energy E, and R αβ is a zero-centered * Electronic address: haw@mpi-hd.mpg.deGaussian-distributed random variable with unit variance.On the basis of Eq. ( 1) it is argued that in the semiclassical regime of high excitation energy, thermalization is generic for chaotic quantum systems [1].
That is an amazing assertion, for the following reason.The state of a quantum system is characterized by the Hermitean time-independent statistical operator Π that describes the distribution of the system over the states in Hilbert space.With H the Hamiltonian governing the system, the time-dependent density matrix ρ(t) is given by In the eigenstate representation of H, the diagonal matrix elements ρ αα of ρ(t) are independent of time, and the nondiagonal elements ρ αβ with α = β carry the time-dependent phase factors exp{−i(E α − E β )t/ }.These oscillate periodically in time.The operator A is independent of time.Therefore, the time-dependent function Tr(Aρ(t)) has a very simple structure, irrespective of whether the system is chaotic or not.How can that simplicity be reconciled with the statement that Tr(Aρ(t)) thermalizes, i.e., tends for time t → ∞ to the limit Tr(Aρ eq )?Why do chaotic systems thermalize while others do not?Eq. (1) indicates that the answer given in Ref. [1] is a probabilistic one.The random fluctuations of a chaotic quantum system encapsulated in the last term of Eq. ( 1) are used in Ref. [1] to show that with overwhelming probability, Tr(Aρ(t)) tends to Tr(Aρ eq ) in the long-time limit.
In the present paper we investigate thermalization from a different point of view.The random fluctuations of eigenvalues and eigenfunctions of chaotic quantum systems agree locally in energy with random-matrix predictions [7].That fact has motivated us to investigate the thermalization of a quantum system for which the Hamiltonian H is a random matrix.Throughout most of the paper, H is taken to be a member of the time-reversal invariant Gaussian Orthogonal Ensemble (GOE) of random matrices of dimension N ≫ 1 while Section V is devoted to the Gaussian Unitary Ensemble (GUE) that breaks time-reversal invariance.By definition, H is paradigmatic for the fluctuation properties that characterize chaotic quantum systems.Clearly, however, the structureless random-matrix Hamiltonian H which treats all states of the system on an equal footing, cannot represent a real chaotic quantum system with all its information content.We are going to show that, nevertheless, our study offers useful insights into the conditions for and characteristic properties of thermalization.We investigate the conditions under which Tr(Aρ(t)) thermalizes for a GOE Hamiltonian.We show that under the same conditions, Tr(Aρ(t)) thermalizes in very similar form also for a GUE Hamiltonian.In comparing our results with those of Ref. [1] we illuminate the difference between the two approaches.
Each realization of the random-matrix Hamiltonian H defined in Section II produces another time-dependent test function Tr(Aρ(t)).Jointly, these test functions form a "stochastic process", i.e., a random ensemble of time-dependent functions.We determine analytically mean value (Section III) and correlation function (Section IV) of that stochastic process in the limit N ≫ 1 of large matrix dimension.We calculate terms of order zero in N as well as corrections of order N −1 .We thereby explore the conditions under which every member Tr(Aρ(t)) of the stochastic process thermalizes.In Section V we investigate thermaliztion for the GUE in the same manner.In Section VI we compare our approach and our results with those obtained in the framwork of the eigenstate thermalization hypothesis.Section VII contains a discussion and the conclusions.
In discussing thermalization we keep in mind these different possibilities.
The operator Π is independent of time.The time evolution of the density matrix in Eq. ( 2) is entirely determined by the unitary operators exp{±iHt/ }. (We assume that Π = Π eq as otherwise ρ(t) = Π eq is independent of time, and the system remains in statistical equilibrium forever.)Except for Section V, the Hamiltonian H is a member of the GOE.The matrix elements H µν are real zero-centered Gaussian random variables with second moments Big angular brackets indicate the ensemble average, and λ defines the width of the spectrum.In terms of its eigenvalues E α and eigenfunctions O µα with α = 1, . . ., N , the Hamiltonian is given by For N ≫ 1 the matrix elements O µα are real zero-centered Gaussian random variables [8] with second moments The factor (1/N ) accounts for normalization.The eigenvalues E α obey Wigner-Dyson statistics.Eigenvectors and eigenvalues are statistically uncorrelated [9].The unitary timeevolution operator U = exp{−iHt/ } is given by The operator U carries the random variables O µα and E α and is, therefore, a matrix-valued stochastic process.
In analogy to the operator Π in Eq. ( 3), the operator A with eigenvalues A j and orthonormal eigenfunctions |j is written as and we have The operators A and Π are independent of time and not statistical.Both are Hermitean.Therefore, Tr(Aρ(t)) * = Tr(Aρ(t)) is real.The operator U (t) defined in Eq. ( 7) is both timedependent and stochastic and so is, therefore, also Tr(Aρ(t)).
We investigate thermaliztion by calculating the ensemble average and the correlation function of Tr(Aρ(t)).To that end we calculate in Appendix 1 the low moments of U (t) and in Appendix 2 the time-dependent functions that control thermalization.

III. ENSEMBLE AVERAGE
We average separately over matrix elements O µα and eigenvalues E α .We begin with the matrix elements.The average is indicated by an index O on the big angular brackets.We use Eq. ( 9) and separate uncorrelated and correlated factors U (t) to write We use the averages given in Eqs. ( 35) and (36) of Appendix 1 and Tr(Π) = 1 to obtain Here The upper index T in expression (11) denotes the transpose, the angular brackets denote the remaining average over the eigenvalues.The function f (t) is obviously of order zero in N and so are the first and the second term on the right-hand side of Eq. ( 11) whereas the third term is of order N −1 .That can be verified using the four forms of Π listed below Eq. ( 3).The second term on the right-hand side of Eq. ( 11) can be written as Tr(Aρ eq (∞)).Here (ρ eq (∞)) αβ = (1/N )δ αβ = (Π eq ) αβ is the time-independent density matrix that describes statistical equilibrium at infinite temperature.If |f (t)| 2 tends to zero for large time, Eq. ( 11) implies thermalization to leading order in N for the ensemble average of Tr(Aρ(t)).The rate of thermalization depends only on |f (t)| 2 , i.e., on the eigenvalues E α of H, and is independent of the test operator A and of the statistical operator Π.
We express |f (t)| 2 as the sum of the uncorrelated and the correlated parts, Both parts are calculated in Appendix 2. The resulting expressions involve the time scales Here 2λ defines the width of the average GOE spectrum, and d is the average GOE level spacing.At the center E = 0 of the spectrum, d has the value d = πλ/N .Since 2λ is the largest and d is the smallest natural GOE energy scale, τ λ is the shortest and τ d is the largest natural GOE time scale.The ratio is The uncorrelated part is proportional to the Fourier transform of the average GOE level density in Eq. ( 41).The oscillatory function g(t/τ λ ), given in detail in Appendix 2, obeys g(0) = 1, is symmetric in time so that g(t/τ λ ) = g(−t/τ λ ), and for large time |t| is bounded by is proportional to the Fourier transform of the GOE two-level correlation function Y 2 (y) defined in Eq. ( 46).The factor 1/N is due to the fact that Y 2 (y) is essentially confined to values y ≪ N of the argument.
In summary we have shown that The first term on the right-hand side represents the statistical equilibrium value Tr(Aρ eq (∞)) of Tr(Aρ(t)).The second term describes thermalization.Since both g(t/τ λ ) and |f (t)| 2 corr are symmetric under the operation t → −t, thermalization is symmetric in time.The function g 2 falls off with time as 4τ 2 λ /t 2 until it reaches the level of the statistical fluctuations given by |f (t)| 2 corr in Eq. ( 16).Since the fluctuations are of order 1/N ≈ τ λ /τ d , that happens when g 2 ≈ 1/N or, with |g| ≈ τ λ /t, at time t ≈ √ τ λ τ d intermediate between the shortest time τ λ and the longest time τ d .Beyond that time, the fluctuations take over.These are governed by τ d and, thus, are extremely slow.The last term in Eq. ( 17) is of order (1/N ).It provides an offset for the fluctuations which prevents for finite N the test function Tr(Aρ(t)) from ever thermalizing completely.The dependence on time both of g and of the correlated part in Eq. ( 16) is determined entirely by the statistical properties of the GOE and is independent of the test operator A and of the statistical operator Π.These two operators determine the values of the last two traces in Eq. ( 17) and, thus, initial and final values of thermalization.
The phase factors exp{±iE α t/ } in expression (10) represent solutions of the Schrödinger equation.That equation is invariant under a time shift t → t + t 0 .We are free to use that shift in expression (10) and, in consequence, in g(t/τ λ ).The resulting function g((t + t 0 )/τ λ ) is a valid alternative to g(t/τ λ ).The maximum is at t = −t 0 .Depending on the sign of t 0 , the maximum occurs at negative or positive values of t.In the latter case thermalization has the unexpected feature that, starting at t = 0, g 2 first grows (with oscillations) with increasing t, reaches the absolute maximum g = 1 at t = t 0 , and decreases thereafter.
As mentioned in the Introduction, thermalization is understood [2] as a test of the hypothesis that in the long-time limit, expectation values of operators tend to their equilibrium values.Thermalization does not imply that the density matrix ρ(t) itself thermalizes, i.e., tends to the time-independent statistical equilibrium distribution.That is demonstrated explicitly by the freedom to shift time t → t + t 0 which shows that the onset in time of thermalization can be chosen arbitrarily.Choosing t 0 > 0 shows that ρ(t) does not thermalize in the time interval 0 ≤ t ≤ t 0 .We first perform the average over the matrix elements O µα .Following Appendix 1 we distinguish contributions due to the first moment of U and to the totally correlated parts of the second, third, and fourth moment of U (t). Nonvanishing contributions to expression (18) are due to the following possibilities: (i) In each trace in the first term of expression (18), each of the two factors U is averaged individually.Then the only non-vanishing contribution to the correlation function ( 18) is due to the correlations between eigenvalues.(ii) In each trace one factor U is part of the correlated part of the second moment, the other factor is replaced by its average (four cases); (iii) in each trace each factor U is paired with a factor U in the other trace to form the correlated part of the second moment (two cases); (iv) two factors U in one trace combine with one factor U in the second trace to form the correlated part of the third moment, the remaining factor U is averaged separately (four cases); (v) the four factors U form the correlated part of the fourth moment (one case).
For case (i) Eq. ( 35) shows that the contribution to the correlation function is given by Here Tr(AΠ) is of order zero in N .The correlation function of |f (t)| 2 involves the GOE two-point correlation function in Eq. ( 16) as well as the GOE three-point and four-point correlation functions.As pointed out below Eq. ( 16) and in Appendix 2, the factor 1/N in Eq. ( 16) is due to the fact that the integration variable y is essentially confined to values of order unity.In the GOE three-point (four-point) functions, two (three) integration variables are similarly confined.That gives rise, respectively, to factors 1/N 2 and 1/N 3 multiplying these functions.In evaluating the correlation function in expression (19) we confine ourselves to contributions of order 1/N .These are due to the two-point function and given by Here f (t) is given in Eq. ( 15) and the two-point correlation function is given in Eq. ( 16).All terms are of order 1/N and are of second order in g and, thus, disappear with t −2 for large |t|.
For case (ii) the results carry the factor 1/N , factors f (2t)(f * (t)) 2 or |f (t)| 2 , and traces of the form Tr(A 2 Π 2 ), Tr(ΠAA T Π T ), or Tr(ΠA) 2 .To leading order in N , all timedependent factors are proportional to g k with k ≥ 2 and tend to zero for large |t|.The factor Tr(A 2 Π 2 ), representative for all terms, is of order 1/N if the rank of Π is of order N and is of order unity if the rank of Π is of order unity.That changes if we request that A be almost diagonal, with diagonal elements A j differing by terms of order 1/N , or that the inequality indicating a ratio of order  2 .All these terms are small of order 1/N or smaller, without any constraints on A. For case (iv) the results carry a factor (1/N ) 2 |f (t)| 2 and factors Tr(ΠAA T ), Tr(A)Tr(AΠΠ T ).These are small of order 1/N .For case (v) the results carry a factor 1/N 3 and factors Tr(A 2 ), Tr(AΠ) Tr(A), or Tr(Π 2 ) [Tr(A)] 2 .These are small of order 1/N .
We have, thus, shown that all terms in the correlation function (18) are of order 1/N or smaller.Without any constraints imposed on A that statement holds if the rank of Π is of order N .If, on the other hand, the rank of Π is of order unity, the statement holds provided that A obeys the constraint (21).
Summarizing the results of the previous and of the present Section, we have shown that for matrix dimension N → ∞, all members Tr(Aρ(t)) of the stochastic process defined by the GOE Hamiltonian H tend to a common limit, The first term on the right-hand side represents the statistical equilibrium value Tr(Aρ eq (∞)) of Tr(Aρ(t)) at infinite temperature.The oscillatory function g 2 has its maximum g 2 = 1 at t = 0, is bounded for large |t| from above by (2τ λ /t) 2 and, thus, tends to zero for large time |t|.Thus, each of the functions Tr(Aρ(t)) thermalizes.Since g(t/τ λ ) = g(−t/τ λ ), thermalization is symmetric in time.The result (22) holds in full generality if the rank of Π is of order N and with the constraint (21) on A if the rank of Π is of order unity.

V. GAUSSIAN UNITARY ENSEMBLE
The symmetry in time of the thermaliztion process in Eq. ( 22) is evident already in Eq. (11).Since To test whether that symmetry is caused by the time-reversal invariance of the GOE, we study thermalization for the Gaussian Unitary Ensemble (GUE) of random matrices.That ensemble is not invariant under time reversal.The GUE consists of Hermitean matrices H U of dimension N ≫ 1.The elements H U;µν are complex zero-centered Gaussian random variables with second moments Here λ U defines the width of the GUE spectrum.In terms of its eigenvalues E U;α and eigenfunctions U µα with α = 1, . . ., N , the Hamiltonian matrix H U is given by The elements U µα of the unitary matrix U are complex zerocentered random variables that carry random phases.Therefore, U µα U µ ′ α ′ = 0 while The eigenvalues E U;α obey GUE Wigner-Dyson statistics.Eigenvectors and eigenvalues are statistically uncorrelated [9].The unitary time-evolution operator U U (t) = exp{−iH U t/ } is given by With the replacement U (t) → U U (t), the stochastic process Tr(Aρ(t)) is given by Eq. ( 9).We average Tr(Aρ(t)) first over the eigenvectors U µα .That is done by contracting each of the two factors U µα with one of the two factors U * µα and use of Eq. ( 25).There are two contraction patterns.The result is Here f U (t) is defined as in Eq. ( 12) except for the replacement  16) by replacing the GOE two-level correlation function by its GUE counterpart.The steps in Section IV can be followed to show that the GUE analogue to the correlation function in Eq. ( 18) vanishes for N → ∞.Thus, thermalization is very similar for the GUE and for the GOE.Since The time dependence of thermalization for the GUE is symmetric about t = 0. We conclude that the symmetry of thermalization with respect to the replacement t → −t is not a consequence of time-reversal invariance.It follows from the fact that ρ(t) is Hermitean.That causes ρ(t) to have the forms U (t)ΠU † (t) for the GOE and U U (t)ΠU † U (t) for the GUE.These forms, in turn, cause the appearance of the time-symmetric factors |f (t)| 2 in Eq. ( 22) and |f U (t)| 2 in Eq. ( 27).

VI. THE EIGENSTATE THERMALIZATION HYPOTHESIS
We compare our procedure and results with those of Ref. [1].In the eigenbasis of the Hamiltonian H with eigenvectors |α and eigenvalues E α we write Eq. ( 8) as Here A αβ = α|A|β and Π βα = β|Π|α .Eq. ( 28) is completely general, i.e., holds for every Hamiltonian, and is not restricted to the model defined in Section II.In Ref. [1] it is assumed that the system is in a single state |ψ(t) = α c α exp{−iE α t/ }|α .That is equivalent to assuming that Π has the special form i.e., that Π has rank one.Moreover, it is assumed that A αβ has the form of Eq. ( 1), with R αβ a zero-centered Gaussian random variable.We interpret Eq. ( 1) and especially the occurrence of the random variable R αβ as due to the fact that in a chaotic quantum system, eigenfunctions and eigenvalues fluctuate randomly [7].In the eigenstate representation of the Hamiltonian, the matrix elements A αβ of the operator A are then random variables.These are parametrized in the form of Eq. (1).Consistency would seem to require that the elements Π βα of the operator Π in Eq. ( 29) and the eigenvalues E α in Eq. ( 28) are treated as random variables, too.That is not done in Ref. [1].The only random element in the approach of Ref. [1] is the variable R αβ .Changing that and taking, in addition to Eq, (1), the energies E α and the coefficients c α as random variables, would alter the results of Ref. [1].We show that explicitly below Eq. (33) for the function C eth (t).
Aside from using the eigenstate thermalization hypothesis Eq. ( 1), the approach of Ref. [1] investigates thermalization of generic chaotic quantum systems with the help of long-time averages of powers of Tr(Aρ(t)) defined for m = 1, 2, . . .as For m = 1 the long-time average gives the mean value of Tr(Aρ(t)), and powers of higher order are used to determine the fluctuations of Tr(Aρ(t)).We investigate that approach by applying it to Tr(Aρ(t)) defined in Eq. ( 9) in terms of the manifestly chaotic Hamiltonian (5), and by comparing the results with those of the random-matrix approach.We do so first for the mean values and then for the time dependence of thermalization.In Ref. [1], the long-time average of Tr(Aρ(t)) is shown to be equal to the thermal average of A at the appropriate temperature provided that the energy spread of the system is sufficiently small compared to the mean energy.That constrains the statistical operator Π.For the random Hamiltonian (5), calculation of the long-time average of Tr(Aρ(t)) in Eq. ( 9) Averaging over the random matrix elements O µα yields three terms, two of which are of order 1/N .With Tr(Π) = 1 the leading-order term is (1/N )Tr(A).That agrees with the first term on the right-hand side of Eq. ( 22) and, thus, with Tr(Aρ eq (∞), in agreement with the result of Ref. [1].(Actually, in the random-matrix approach there is no need for any constraint on Π since Eq. ( 11) is an immediate consequence of orthogonal invariance.)We note that in both approaches, the average terms are actually present for all times t > 0. In our case that is seen directly from Eq. ( 11).For Ref. [1] it follows from the fact that fluctuations around the long-time average are calculated as fluctuations of (Tr(Aρ(t)) − Tr(Aρ eq )).
We turn to thermalization.In Ref. [1] the time dependence of thermalization is determined in two steps.In the first step, higher moments of (Tr(Aρ(t))−Tr(Aρ eq )) are estimated with the help of Eq. ( 1).It is shown that in the semiclassical regime and provided that the spread in energy of the system is sufficiently small, the time evolution of (Tr(Aρ(t)) − Tr(Aρ eq )) is with overwhelming probability given by C(t)Tr(Aρ(0)).
Here C(t) is independent of the initial value Tr(Aρ(0)) and is given by where the long-time average is over t ′ .In the second step, the numerator on the right-hand side of Eq. ( 32) is written in the eigenstate representation of the Hamiltonian as Replacing A αβ by the last term in Eq. ( 1) and averaging over the ensemble results in with F (E, ω) introduced in Eq. ( 1).The index "eth" stands for the eigenstate thermalization hypothesis.As mentioned earlier, in deriving Eq. (33) it is assumed that the eigenvalues E α are ordinary nonstochastic variables.If the E α were taken as random variables, one would have to average } over the joint probability distribution of E α and E β .Averaging would have to take into account that the random variables E α and E β occur not only in the exponential but also as arguments of the function F (E, ω) in Eq. ( 1).The result would differ from Eq. (33).
We compare the first step taken in Ref. [1] with the randommatrix model.There it is obvious almost from the outset (see Eq. ( 11)) that the time dependence of thermalization is independent of the initial value of A and given by |f (t)| 2 .To compare with that result, we calculate C(t) from Eq. (32) and use Eq. ( 9).We focus attention on the numerator on the righthand side of Eq. (32).In averaging over the matrix elements O µα we take into account only contraction patterns that retain the dependence on t.That is consistent with Ref. [1].Using C(0) = 1 we find where the index "rm" stands for the random-matrix model.Eq. ( 34) is in agreement with Eq. ( 22), confirming the definition Eq. (32) of C(t) given in Ref. [1] and the claim that the time evolution of thermalization is given by C(t)Tr(Aρ(0)).
The second step taken in Ref. [1] leads to the explicit expression for C eth (t) in Eq. (33).That expression differs from C rm (t) in Eq. ( 34).The function C eth (t) in Eq. ( 33) is the Fourier transform of the non-negative function |F (E, ω)| 2 .Therefore, C eth (t) oscillates in time about the value zero.If the band width in ω of |F (E, ω)| 2 is finite, C eth (t) is asymptotically inversely proportional to t.In contrast, the function C rm (t) in Eq. ( 34) is non-negative for all times t.To leading order in 1/N , |f (t)| 2 is equal to the square of the Fourier transform g(t/τ λ ) in Eq. ( 15) of the average level density and, for large |t|, falls off as 1/t 2 .Thus, the time dependence of thermalization is fully determined analytically within the random-matrix model.Except for the assumption on smoothness, the dependence on E and ω of the function F (E, ω) in Eq. ( 1) is not defined nor is it obvious how information on that function might be obtained.Actually, in Ref. [1] the band width of the function F (E, ω) with respect to ω is used as a parameter and is estimated with the help of physical arguments.
The fundamental difference between the random-matrix model and the eigenstate thermalization hypothesis lies in the identification of the source of randomness.In the randommatrix model, it is the operator U (t) defined in Eq. ( 7) and appearing in Eq. ( 9) that carries the random variables and that causes the function Tr(Aρ(t)) to be a stochastic process.The operator A, in contradistinction, is written in the fixed basis of states (µ, ν) as in Eq. ( 8) and, therefore, is not random.The central result, Eq. ( 11), is valid for any operator A. The approach of Ref. [1] takes the converse point of view.The elements ρ αβ (t) of the density matrix of the system are treated as ordinary (non-statistical) variables.The stochasticity of the chaotic system is contained entirely in the eigenstate representation A αβ of the operator A.
We see that Eq. ( 1) is taylored to yield thermalization.From the point of view of formal logic, Eq. ( 1) formulates a sufficient condition for thermalization.It is the great merit of Ref. [1] to have in that way given a first analytical approach to and understanding of thermalization.A number of questions remains.Why is randomness attached only to the elements A αβ of the operator A? Does Eq. ( 1) impose a constraint on the operator A? If so, which operators obey that constraint?Do all operators that obey the constraint possess the same function F (E, ω), or does that function depend on the operator A under consideration?How does F (E, ω) (beyond being smooth) actually depend on the variables E and ω?
The answers should link Eq. ( 1) with the Hamiltonian of the chaotic quantum system.Without that link, Eq. ( 1) remains a hypothesis.
Quantum chaos and stochasticity of the Hamiltonian are closely linked [7].Therefore, it would seem more natural to analyze thermalization in terms of a stochastic model for the Hamiltonian rather than with the help of Eq. ( 1).The randommatrix model defined in Section II is too restricted for that purpose.The approach of Ref. [10] would seem a suitable candidate.Unfortunately and in contrast to Ref. [1], it does not offer the chance of an analytical treatment.That is the strength of the eigenstate thermalization hypothesis.

VII. DISCUSSION AND CONCLUSIONS
We have shown that in the random-matrix model, every test function Tr(Aρ(t)) thermalizes, i.e., attains the limit (22) for N → ∞.Thermalization is due to the orthogonal invariance of the GOE: Orthogonal invariance causes the matrix elements O µα in Eq. ( 5) to be zero-centered Gaussian random variables that obey Eq. ( 6).That is the only fact used in deriving Eq. (11).In particular, orthogonal invariance determines the form of the equilibrium term (1/N )Tr(A).The last term in Eq. ( 22) gives the time dependence of thermalization.One might have expected that the time dependence of thermalization is determined by destructive interference of the phase factors exp{±iE α t/ } in Eq. ( 2) for the density matrix.That is not the case.The leading-order contribution to thermalization is given by the function f (t) , the Fourier transform g of the average GOE level density.That function does not contain interference terms.The GOE two-point level correlation function only yields a term of order 1/N .Both functions are determined by orthogonal invariance.Corresponding statements hold for the GUE where unitary invariance takes the role of orthogonal invariance.
Thermalization as defined in Section I is a property of the test function(s) Tr(Aρ(t)).It characterizes the system but is not an obvious property of the density matrix ρ(t) in Eq. ( 2) itself.Thermalization is caused by the statistical properties of the time-evolution operator U (t) in Eq. (9).Thermalization is symmetric in time.That symmetry holds for both, the time-reveral invariant GOE and the time-reversal non-invariant GUE.
Comparison with the eigenstate thermalization hypothesis Eq. (1) shows that long-time averages agree with randommatrix results.The predicted time dependence of thermalization is qualitatively correct.The statistical assumptions used in Ref. [1] need justification.It would be desirable to link form of and parameters in Eq. (1) to the Hamiltonian of the chaotic system.
Thermalization is very different from equilibration in open quantum systems [11].The time dependence of equilibration is described by a master equation that violates time-reversal invariance.The statistical operator of the system changes with time and tends toward the statistical equilibrium distribution.The process is exponential in time.In contrast, thermalization occurs in an isolated system.The statistical operator of the system does not change with time, the elements of the density matrix of the system keep oscillating in time forever.Only the test function Tr(Aρ(t)) thermalizes with a characteristic dependence on time t given by 1/t 2 .
In Ref. [5], the eigenstate thermalization hypothesis is formulated without reference to Eq. ( 1) and in more general terms by stating that "every eigenstate in an energy shell represents thermal equilibrium".We believe that our randommatrix models are paradigmatic for the eigenstate thermalization hypothesis in that more general form.The energy shell is the entire Hilbert space.All eigenstates |α = µ O µα |µ are random superpositions of the basis states |µ .The statistical weight of every component vector |µ is the same, independent of the eigenvalue E α to which the state belongs.
Thermalization tests the hypothesis that in an isolated system and in the long-time limit, expectation values of operators tend to their equilibrium values.That hypothesis is seen as the quantum analogue of the ergodic hypothesis in classical statistical mechanics [2].It is not surprising that in the random-matrix models, thermalization follows from orthogonal or unitary invariance.Such invariance implies that all states in N -dimensional matrix space are equivalent.That is similar in spirit to classical ergodicity.
IV. CORRELATION FUNCTIONAs a measure of the size of the fluctuations of individual members Tr(Aρ(t)) around the ensemble average in Eq. (17) we calculate the correlation function Tr(Aρ(t))Tr(Aρ(t)) − Tr(Aρ(t)) Tr(Aρ(t)) .(18) 1/N .The postulate (21) resembles the eigenstate thermalization hypothesis of Eq. (1).It does not, however, invoke any statistical assumption on A. If condition (21) is met, the term (1/N )Tr(A 2 Π 2 ) is of order 1/N if the rank of Π is of order unity and is considerably smaller than that if the rank of Π is of order N .That same statement holds for the other contributions to case (ii).For case (iii) the results carry a factor 1/N 2 , a factor unity or |f (2t)| 2 , and factors µν A 2 µν Π 2 µν , Tr(A 2 )Tr(Π 2 ), or [Tr(AΠ)] Proceeding as in Appendix 2 we see that f U (t) is proportional to the Fourier transform of the average GUE level density, and that the correlated part |f U (t)| 2