Thermalization of closed chaotic many-body quantum systems

We investigate thermalization of a closed chaotic many-body quantum system by combining the Hartree–Fock approach with the Bohigas–Giannoni–Schmit conjecture. The conjecture implies that locally, the residual interaction causes the statistics of eigenvalues and eigenfunctions of the full Hamiltonian to agree with random-matrix predictions. The agreement is confined to an interval Δ (the correlation width). The results are used to calculate Tr(Aρ(t)) . Here ρ(t) is the time-dependent density matrix of the system, and A represents an observable. In the semiclassical regime, the average ⟨Tr(Aρ(t))⟩ decays on the time scale ℏ/Δ toward an asymptotic value. If the energy spread of the system is of order Δ, that value is given by Tr(Aρeq) where ρeq is the density matrix of statistical equilibrium. The correlation width Δ is the central parameter of our approach. We argue that Δ occurs generically in chaotic quantum systems and plays the same central role.


Introduction
Thermalization addresses the time evolution of the expectation value Tr(Aρ(t)) of a classical observable (represented in Hilbert space by a Hermitian operator A) in a closed chaotic quantum system S with density matrix ρ(t).Thermalization postulates that asymptotically (time t → ∞) the expectation value Tr(Aρ(t)) of A tends toward its equilibrium value, Tr (Aρ (t)) → Tr (Aρ eq ) . ( Original Content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence.Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.
Here ρ eq is the density matrix of S in statistical equilibrium.The asymptotic relation (1) is expected to hold in the semiclassical regime of very high excitation energy E and for systems with very small energy spread δE ≪ E. Thermalization is seen as the analog in quantum statistical mechanics of ergodicity in classical statistical mechanics.Reviews are given in [1][2][3].
Starting with [4,5], thermalization is commonly derived [2,3] in the eigenbasis of the Hamiltonian H of S with states labeled α, β where Tr (Aρ (t)) = ∑ αβ A αβ ρ βα (t) . ( It is assumed that the matrix elements of the operator A αβ are random variables while the matrix elements ρ βα (t) of the density matrix are taken to be fixed (i.e.non-statistical).If the statistical properties of A αβ reflect the statistical properties of the Hamiltonian H, that assumption is not convincing.Consistency requires that both A αβ and ρ βα (t) be taken as stochastic variables.That was demonstrated in [6].The distribution of A αβ and ρ βα (t) should follow from the spectral fluctuation properties of the Hamiltonian H of the system S (i.e. from the joint probability distribution of eigenvalues and eigenfunctions).
The question addressed in this paper is based on that premise.We ask: how to build a model of thermalization of a chaotic quantum system that is based on statistical properties of the Hamiltonian?We focus attention on a chaotic many-body quantum system S. Using a combination of a mean-field approach and the Bohigas-Giannoni-Schmit (BGS) conjecture [7], we determine the generic spectral fluctuation properties of the Hamiltonian of S. That makes it possible to calculate Tr(Aρ(t)) explicitly.We thereby determine the time scale of thermalization.We quantify the conditions 'high excitation energy' and δE ≪ E that are required for thermalization.We show what happens to the asymptotic relation (1) if the system is not in the semiclassical regime, or if the condition δE ≪ E is violated.
The present paper follows [6] wherein a random-matrix approach to thermalization has been developed.We closely follow the overall layout of that paper.Except for occasional references we do not, however, repeat here all of the central arguments nor some of the technical steps of [6].
The paper is structured as follows.In section 2 we describe our approach.In section 3 we formulate the Hartree-Fock (HF) approach and determine the generic statistical properties of the Hamiltonian of a chaotic many-body quantum system S.These are used to define in section 4 a random Hamiltonian ensemble.In terms of that ensemble, Tr(Aρ(t)) turns into a time-dependent stochastic process.In sections 5 and 6 we calculate ensemble average and correlation function, respectively, of Tr(Aρ(t)).In section 7 we extend the treatment to timereversal noninvariant systems.Section 8 contains a summary and the discussion.

Method
With H the Hamiltonian, the density matrix is given by ρ (t) = exp {−i Ht/ℏ} Π exp {i Ht/ℏ} . ( Here Π is the time-independent statistical operator of S. It completely specifies the system S by giving the distribution of S over the states in Hilbert space.The matrix representation of Π depends, of course, on the choice of basis in Hilbert space.Normalization implies All eigenvalues π κ , κ = 1, 2, . . . of Π are positive or zero.Together with equation (4), that implies 0 ⩽ π κ ⩽ 1 and ∑ κ π κ = 1.The system S being closed, the operator Π is independent of time.The entire time dependence of Tr(Aρ(t)) is due to the unitary time-evolution operator exp{−i Ht/ℏ} and its Hermitian adjoint in equation (3).Thermalization is expected to occur for most of the initial conditions on S, i.e. for almost any choice of the statistical operator Π provided only that the system is in the semiclassical regime and δE ≪ E applies.(In the context of the present paper, the expressions 'most' and 'almost' are given precise meaning in section 8).In [5] it is assumed, for instance, that Π possesses only a single nonzero eigenvalue and, thus, describes a system very far from statistical equilibrium.The operator A may be constrained by the condition that A possesses a viable classical counterpart.
As stressed in the literature [2,3], thermalization is not universal.Integrable quantum systems and systems that show many-body localization do not thermalize [2].Thermalization is expected to occur generically, however, in chaotic many-body systems [2].From equation (3) it is not evident how that difference in behavior comes about.Indeed, writing the observable Tr(Aρ(t)) explicitly in the eigenbasis of H with orthonormal eigenvectors denoted by the letters α, β, . . .and eigenvalues labeled E α , E β , . . .and using the Einstein summation convention we have That expression is universal and holds irrespective of the particular dynamical properties of the system.The matrix elements A αβ and Π βα and the diagonal elements ρ αα of the density matrix in equation ( 5) are independent of time, while the nondiagonal elements ρ βα (t) with β ̸ = α carry the phase factors exp{i (E α − E β )t/ℏ} and, thus, keep oscillating in time forever.These statements, too, hold alike for integrable, chaotic, and for systems that show many-body localization.
The fundamental difference between chaotic systems on the one hand and integrable systems or systems that show many-body localization on the other, lies in the spectral fluctuation properties of H.For chaotic quantum systems the BGS conjecture states that these coincide locally with those of the random-matrix ensemble in the same symmetry class (orthogonal, unitary, or symplectic) [7].No such correspondence exists for integrable systems or for systems that show many-body localization, and statistical assumptions cannot be used there.Therefore, we expect equation (1) to follow from the statistical properties of the Hamiltonian of a chaotic quantum system.That is in line with the standard approach to thermalization [2,3,5].As mentioned in section 1, in that approach it is shown that thermalization follows and can be worked out explicitly if it is assumed that the elements A αβ of the observable A in equation ( 4) are Gaussian random variables while the matrix elements Π βα of the statistical operator are taken to be nonstatistical.For reasons given in section 1 and, in more detail, in [6] we follow here a different route.
Using an arbitrary fixed basis of states labeled (µ, ν, . ..) in Hilbert space we write equation (5) in the form Here U(t) is the unitary time-evolution operator To work out thermalization we consider the matrices A µν and Π ρσ as fixed, i.e. independent of time and nonstatistical.The entire time dependence and the statistical properties of Tr(Aρ(t)) reside in U(t) and its Hermitian adjoint.We show that the statistical properties of the Hamiltonian of the chaotic many-body quantum system S and, thus, of U(t) can be determined via a combination of a mean-field approach and of the BGS conjecture [7].As mentioned above, the conjecture states that the spectral fluctuation properties of the Hamiltonian H of a chaotic quantum system coincide locally with those of the random-matrix ensemble in the same symmetry class.The word 'locally' refers to the existence of an energy interval ∆.The interval ∆ limits the range of energies wherein the spectral fluctuation properties of a chaotic quantum system coincide with those of random-matrix theory.The interval ∆ occurs generically in chaotic quantum systems.It is called the Thouless energy in diffusive quantum systems [8], the localization length in the theory of banded random matrices [9], it shows up in nuclear shell-model spectra [10], it has been identified for the quantum Sinai billiard [11] and for a quantized Cantori system [12], and it is given by the pion mass in chiral randommatrix theory of quantum chromodynamics [13].The existence of ∆ is an integral part of the BGS conjecture.Without it, the Hamiltonian would be a member of the Gaussian orthogonal ensemble (GOE) of random matrices.That is physically impossible since that ensemble carries no information.In what follows we refer to ∆ as the width of random-matrix correlations, in short the correlation width.In nuclear physics the expression 'spreading width' is often used instead [10].We do not follow that custom because it refers to wave-function mixing only, not to eigenvalue statistics.
The BGS conjecture forms a cornerstone of our approach.That is consistent with the central role played by the conjecture in the theory of chaotic quantum systems, see, for instance, the book by Haake [14].The conjecture has been tested successfully on a vast variety of numerical models of classically chaotic systems, and in the low-energy regime of several physical systems.We cite [15][16][17] as examples.While an analytical proof of the conjecture is still lacking, its universal validity for chaotic quantum systems has, for the case of the level-level correlation function, been demonstrated in a physically convincing manner with the help of Gutzwiller's periodic orbit theory, see [18,19].Periodic-orbit theory is semiclassical by construction.That is particularly fitting in the present context wherein the semiclassical regime plays an important role.We believe that the BGS conjecture is now universally accepted, with well-understood exceptions lending additional support to its generic validity, see, for instance, [20].

HF approach
We establish the generic spectral fluctuation properties of the Hamiltonian H of a chaotic many-body quantum system by combining the mean-field approach with the BGS conjecture.In the main body of the paper, we confine ourselves to time-reversal-invariant systems.Violation of time-reversal invariance is addressed in section 7. Without further mention we focus attention on a fixed set of conserved quantum numbers.We assume that the many-body Hamiltonian is governed by two-body interactions.The HF approach yields an approximation H HF to H.The total Hamiltonian H can then be written as where V is the residual interaction.The eigenstates of H HF are labeled m, n, . .., and we have The The residual interaction V is dominated by one-and two-body forces.That is what we focus on.The one-body matrix elements of V connect many-body HF states that differ only in the occupation of a single single-particle HF state.The two-body matrix elements of V connect only many-body HF states that differ in the occupation of not more than two single-particle HF states.The one-body matrix elements of V involve the overlap of two single-particle HF eigenfunctions.The bigger the difference in quantum numbers characterizing these two functions, the smaller the overlap.Therefore, there exists an effective cutoff that limits the number of many-body HF states connected to a given one by the one-body matrix elements of V.The same arguments apply to the two-body matrix elements of V.With the cutoff, the number of nonvanishing matrix elements of V in each row and column of the matrix H in equation ( 8) is finite.The matrix representing H is sparse and banded.The band width is defined by the matrix elements of V farthest from the main diagonal.The band width increases with increasing excitation energy because the average level density does.The average number of zeros separating neighboring nonvanishing elements V increases.In the sense of [21] the matrix H is local and is, therefore, a viable candidate for thermalization.
In spite of its being sparse in the representation of equation ( 8), the residual interaction V is capable of mixing the HF states sufficiently strongly to produce eigenvalues and eigenfunctions that locally follow random-matrix predictions.That has been demonstrated, for instance, by extensive numerical calculations of spectra and eigenfunctions for nuclei in the middle of the s − d shell in [10].The shell-model Hamiltonian used in [10] is the sum of singleparticle energies and a two-body residual interaction and, therefore, belongs to the class of Hamiltonians considered in equation (8).Confining Hilbert space to the states of the s − d shell, the authors diagonalized Hamiltonian matrices of dimensions up to several thousand, defined by parity, total spin and isospin.In the centers of the spectra, the results for eigenfunctions and eigenvalues showed close agreement with GOE predictions.
At the same time, the results of [10] indicate limits of that agreement.For the eigenvalues, the nearest-neighbor spacing distribution agrees with the GOE prediction.However, the spectral rigidity agrees with the logarithmic dependence predicted by the GOE for small energy intervals only.At some point the spectral rigidity shows a sudden upbend and increasing disagreement with the GOE result beyond that point (see, for instance, figure 27 of [10]).Likewise, the eigenfunctions of the shell-model Hamiltonian, although characterized by a Gaussian distribution of their components, are not a total mixture of all basis states.Correspondingly, the residual interaction does not spread the unperturbed many-particle states m of the shell model uniformly over the eigenfunctions of the shell-model Hamiltonian.Rather, their distribution is centered at the unperturbed energy E m and has a finite width ∆ in energy (the 'spreading width'), smaller than the range of the total spectrum of the shellmodel Hamiltonian.In [10], two forms of the distribution of the unperturbed states play a role.If the matrix elements of V, although sparse, cover uniformly the HF matrix (8), the matrix V is not banded, and the distribution has the shape of a Lorentzian with width ∆ given by the standard expression Here ⟨V 2 mn ⟩ denotes the mean square matrix element, the average with respect to m and n being taken over an interval centered at E that comprises a sufficiently large number of HF states labeled (m, n).Actually, equation ( 10) is a variant of Fermi's golden rule.Without averaging over m it describes the spreading of the HF state |m⟩ over the eigenstates of the system that is due to mixing with the states n ̸ = m.The average over m is used to attain independence of ∆ of the particular state |m⟩.Expressions of the form of equation ( 11) have been extensively derived, discussed, and used in the literature, see, for instance, [22,23].If, on the other hand, the matrix V is banded, the shape of the distribution is Gaussian.The results of [10] are in between these two limiting cases, see figure 46.If shells higher than the s − d shell were to be included in the calculation, the band structure of V would be pronounced more clearly, and we would expect the distribution to approach the Gaussian form more closely, see also [1].
In summary, the work of [10] has shown that the residual interaction V is capable of producing strong mixing of the basis HF states that lead to spectral fluctuations of the GOE type.These fluctuations are local.Their range in energy is bounded.That is true both for the statistics of eigenvalues and for the statistics of eigenfunctions.In the present paper we assume that the energy interval ∆ confining statistical agreement with GOE predictions is the same for eigenvalues and eigenfunctions.To the best of or knowledge, that assertion has not been tested yet in realistic cases.
We accordingly use the following formulation of the BGS conjecture.At each energy E, there exists an interval ∆ (called the correlation width and semiquantitatively given by equation ( 10)) centered on E within which the action of V in equation ( 8) causes spectral fluctuations of the GOE type.Such fluctuations exist also at an energy E ′ far removed from E (so that |E − E ′ | ≫ ∆) but are not correlated with the ones around E.
We are going to use our assumption in the semiclassical regime.It is, therefore, important to understand the dependence of ∆ on excitation energy E. A simple scaling argument shows that ∆ as given by equation ( 10) depends only weakly on E. For an energy E ′ much bigger than E we account for the exponential growth of the level density by the scaling parameter α ≫ 1 so that ρ(E) → ρ(E ′ ) = αρ(E).That causes the number of zeros separating two adjacent nonvanishing matrix elements of V to grow, too.The growth is proportional to α.In consequence, the average ⟨V 2 mn ⟩ is multiplied by the factor 1/α, and ∆ remains the same.Although the estimate for the spreading width of a banded matrix differs from equation ( 10) it seems reasonable to expect a qualitatively similar behavior in that case.We conclude that the correlation width ∆ does not possess the nearly exponential dependence on energy that characterizes ρ(E) (although ∆ may display a polynomial dependence on energy).Therefore, the number of states N ∆ = ∆ρ(E) in the interval ∆ increases strongly with energy.In the semiclassical regime we have N ∆ ≫ 1, and we may use an expansion in inverse powers of N ∆ , keeping only leading-order terms.

Statistical properties
The Hamiltonian H in equation ( 8) is diagonalized by an orthogonal transformation with elements O mα and has eigenvalues E α , so that As done for the eigenvalues of the HF Hamiltonian, we order the nondegenerate eigenvalues E α so that E α < E β for α < β.By assumption, the Hamiltonian H describes a chaotic many-body quantum system.Therefore, the BGS conjecture holds in the form formulated in section 3.
Obviously that requires V to be sufficiently strong so as to cause strong local mixing of the HF states |m⟩.
According to the BGS conjecture, there exists for each energy E an energy interval ∆ centered on E within which the eigenvalues E α and eigenvectors O mα approximately obey GOE statistics.The approximation improves with the number N ∆ = ∆ρ(E) of states located within the interval ∆ and becomes excellent in the semiclassical regime.For convenience we assume ∆dρ(E)/dE ≪ ρ(E) so that in the interval ∆ the average level density ρ(E) is effectively constant.
It is technically difficult to use the statistical properties of H directly to work out mean value and variance of Tr(Aρ(t)).It is preferable to construct an ensemble of Hamiltonians that all have the same statistical properties as H and to show that the resulting expression for the ensemble average of Tr(Aρ(t)) holds for every member of the ensemble and, therefore, also for H.That is the route we follow.We construct the ensemble by focusing attention on HF states with energies E m and eigenstates of H with eigenvalues E α that are located within a correlation length ∆ of each other.In the limit N ∆ → ∞, the eigenvalues E α and matrix elements O mα are uncorrelated random variables.The eigenvalues obey Wigner-Dyson statistics.The matrix elements O mα are zero-centered real Gaussian random variables with second moments The big angular brackets denote the ensemble average.The function F(x) in equation ( 12) limits the range of GOE correlations of the matrix elements O mα in terms of the correlation width ∆, and E α is the ensemble average of the E α with fixed α [24].Thus, F(x) must be positive for all x ⩾ 0 and must fall off steeply with increasing  12) implies that the residual interaction spreads the HF eigenstates over an energy interval of width ∆.Conversely, each eigenstate at energy E α is a linear combination of HF states located within a neighborhood of width ∆ centered at E α .We note that for the GOE itself, the correlations extend over the entire range of the spectrum.In equation ( 12) the function F(x) is, thus, replaced by 1/N where N is the GOE matrix dimension, signaling complete mixing of all states of the spectrum.With λ defining the width of the GOE spectrum, the role of F(x) is taken over by the GOE average level density ).In the present case, such complete mixing is restricted to states within the correlation width ∆.While our random-matrix ensemble may not reflect physical reality very closely near the ground state (where the inequality N ∆ ≫ 1 may not apply), it certainly does apply in the semiclassical regime with its very high average level density.
The orthogonality relations For N ∆ ≫ 1 we use the continuum limit, replacing summations by integrations.Then the orthogonality relations read The range of integration comprises the entire spectrum but actually is limited by the sharp falloff of F(x) with increasing x.We assume that the strong local mixing of HF states due to V does not affect the large-scale average level density.Then the average density ρ(E α ) of the average eigenvalues E α of H equals the average level density ρ(E m ) of the HF states Here E is an energy somewhere within the interval ∆.
The analytical form of F(x) has been addressed in section 3.For the standard model (which yields equation (10) for the correlation width ∆) the function F(x) has Lorentzian shape,

F
[ ( More realistic in the present case is the Gaussian form Equation ( 16) was also obtained recently for a system of interacting spins [25].We use equation ( 16) in the main part of the paper.In section 5 we show how the use of equation ( 15) affects the time dependence of thermalization.
Collecting results we obtain In the HF basis, the time-evolution operator U(t) has matrix elements The eigenvalues E α and the matrix elements O mα are random variables with distributions specified in and below equation (12).While A and Π are independent of time and nonstatistical, U(t) is time dependent and, via its dependence on the matrix elements O mα and eigenvalues E α , is stochastic.These properties of U(t) cause Tr(Aρ(t)) to be a time-dependent stochastic process.
It would not make sense to use the statistical assumptions on the eigenvalues E α and matrix elements O mα in equation ( 18) right from the beginning, i.e. for all t ⩾ 0. Starting at t = 0 from an arbitrary fixed basis of states (µ, ν), it takes a finite time t 0 for the stationary statistical correlations of eigenfunctions and eigenvalues of H defined around equation (12) to be dynamically observable.We are interested here in the long-time development of Tr(ρ(t)) and do not address the transients that occur in the time interval 0 ⩽ t ⩽ t 0 .We use the statistical assumptions (12) in equations ( 17) and (18) only after the transients have subsided, i.e. for times t > t 0 , without mentioning that restriction hereafter.
To characterize the semiclassical regime of large energy E and to display explicitly the condition δE ≪ E, we calculate mean value and second moment of H.We identify the energy E of the system with ⟨ Tr(HΠ) ⟩ .We use equations ( 10), ( 12) and ( 16) and find The energy E is defined in terms of the HF Hamiltonian only.On average the statistical fluctuations do not contribute to E. For E to be in the semiclassical regime, the diagonal matrix elements of Π in the HF basis must be concentrated on a corresponding range of large HF energies.
For the calculation of the second moment we assume that ∆ is constant (i.e.independent of energy) over the range of HF energies where the matrix elements Π mm are essentially different from zero and use Tr(Π) = 1.We find An estimate of the square of the energy spread δE is given by Equations ( 19)- (21) show that for the inequality δE ≪ E to hold, the range of HF states for which the diagonal elements Π mm essentially differ from zero, must be sufficiently narrow.
The big round bracket in equation ( 21) is positive or zero.That implies δE ⩾ ∆ for every choice of Π.The big round bracket vanishes when in the HF basis Π is restricted to a single nonzero diagonal element so that Π mm ′ = δ mm0 δ m ′ m0 .Even then, the statistical fluctuations of the random-matrix ensemble cause a minimum energy spread of size ∆.The condition δE ≪ E then reads ∆ ≪ E.

Ensemble average
In the present section we calculate the ensemble average ⟨Tr(Aρ(t))⟩ of Tr(Aρ(t)).As done in the derivation of equation ( 20), we assume that ∆ is independent of energy in the range of HF energies where the matrix elements of Π in the HF basis are essentially different from zero.Without that assumption, the increased complexity of the notation would only obscure the basic simplicity of our arguments and results, without furnishing additional insights.
Using expression (18) we calculate first the ensemble average over the elements O mα of the orthogonal matrix, keeping the eigenvalues E α fixed.The elements O mα are zero-centered Gaussian-distributed random variables with second moments given in equations ( 12) and (16).We find Here and in equation ( 23), angular brackets denote only the average over the matrix elements.The covariances are defined as A term carrying four Gaussian-distributed matrix elements O mα is averaged by 'contraction', i.e. by applying equations ( 12) and ( 16) to all pairs of matrix elements, a total of four possibilities.We find The upper (lower) sign holds, respectively, for v 1 (for v 2 ).Using equation ( 22), its complex conjugate, and the result (24) for v 1 in equation (17) gives The big angular brackets denote the remaining average over the eigenvalues.The average consists of two contributions resulting, respectively, from averaging the sums over α and β separately, and from the correlation function of E α and E β .The latter is given by the two-point GOE eigenvalue correlation function [26].That is a function of (E α − E β )ρ(E α ) which differs significantly from zero only over the range of a few mean level spacings.Therefore, the correlated part of the average in equation ( 25) is small of order 1/N ∆ .That is shown explicitly in [6].We keep the uncorrelated part only, i.e. we perform the sums over α and β independently.In the sum over α we write While (E m − E α ) extends over the range ∆, the fluctuations (E α − E α ) are confined to a few mean level spacings.Therefore, the term (E α − E α ) is small of order 1/N ∆ compared to (E m − E α ), and is neglected.We proceed analogously for the sum over β.All sums over α and β in equation ( 25) are then carried out in the continuum limit.We find We simplify the terms on the right-hand side of equation ( 26), neglecting terms of order 1/N ∆ .We use the eigenbasis of the operator Π with orthonormal eigenvectors ⟨n|κ⟩ and with eigenvalues π κ that obey ∑ κ π κ = 1 and 0 ⩽ π κ ⩽ 1 for all κ.Then Tr(AΠ) = ∑ κ A κκ π κ .The constraints on the eigenvalues π κ imply that Tr(AΠ) is of order unity.We show that the same is true for the term proportional to A mm Π nn while the term proportional to A mn Π nm will be shown to be of order 1/N ∆ .The term proportional to A mm Π nn carries a double sum over m and n.For our estimate we divide the summation over m into energy intervals of length ∆ each, labeled consecutively by k = 1, 2, . ... For fixed m ∈ k, the sum over n can be read as the partial trace of Π extended over the HF states in a neighborhood of width √ 2∆ of E m .For all states m ∈ k we approximate that partial trace of Π by an expression p k that is independent of m and given by Here 0 ⩽ p k ⩽ 1, and ∑ k p k = 1.The approximation ( 27) may be off by a numerical factor but not by a factor of order N ∆ .The term proportional to A mm Π nn takes the form Here ρ k is the average level density in the interval labeled k, and Tr k denotes the trace taken over the HF states in that interval.Since Tr k (A) is of order N ∆ , the term is of order unity.Using similar arguments we show in the appendix that the last term on the right-hand side of equation ( 26) is proportional to 1/N ∆ .That term is neglected.
In summary, we have shown that for 1/N ∆ ≪ 1, the ensemble average is given by ⟨ Tr (Aρ (t)) Equation ( 28) rests on the assumption (27).That assumption is not absolutely necessary.For an accurate treatment it is possible to retain the original form of the term proportional to A mm Π nn .However, the assumption (27) simplifies that term and makes the result physically more transparent.We keep the form (28) for the remainder of the paper.The central approximation used in deriving equation ( 28) is N ∆ ≫ 1.As shown in section 3, that condition confines the statistical operator to high energies.That is consistent with the condition 'energy E in the semiclassical regime' stated in the literature [2,3,5] as necessary for thermalization.
We discuss the result (28).The first term on the right-hand side describes the relaxation of the term ⟨Tr(Aρ(t))⟩ as function of time t.Starting at time t = 0 with Tr(AΠ), the term vanishes for t → ∞.The governing factor for relaxation is the Gaussian with characteristic time scale ∆.The Gaussian is the square of the Fourier transform of the correlation function F(x) in equation (16).That is consistent with the result of [6] where the time dependence of relaxation was found to be given by the square of the Fourier transform of the average GOE level density It is easily verified that use of the Lorentzian (15) for F(x) results in an exponential factor exp{−∆t} instead of the Gaussian.We conclude that the correlation width ∆ determines the time scale for relaxation irrespective of the form of F(x).The actual analytical form of the relaxation function depends upon the choice of the function F(x).
For time t → ∞, ⟨Tr(Aρ(t))⟩ approaches asymptotically the last term in equation ( 28).If the system thermalizes, that term must be equal to Tr(Aρ eq ).For an unrestricted statistical operator Π, that is, in general, not the case.The statistical operator may be spread over several or many intervals k.Then a corresponding number of the coefficients p k in equation ( 28) have values that differ significantly from zero.For an arbitrary distribution of the coefficients p k the asymptotic value of ⟨Tr(Aρ(t))⟩ differs from the equilibrium value Tr(Aρ eq ).Thermalization occurs only if the statistical operator Π is confined to a single interval k 0 or, less strictly, if Then the approximation (27) becomes nearly exact, we have p k ≈ δ kk0 , and the last term in equation ( 28) is given by (1/N k0 )Tr k0 (A).That is equal to the equilibrium value Tr(Aρ eq ) because, with Π confined to the interval k 0 , equilibrium corresponds to equal occupation probability of all states in k 0 .(For the narrow energy interval ∆, the Boltzmann factor is practically constant).We conclude that for thermalization to occur the condition δE ≈ ∆ must hold.That quantifies the usual condition δE ≪ E stated in the literature [2,3].Thermalization would also be attained if several or many of the coefficients p k differ from zero and if their distribution corresponds to statistical equilibrium.But that is the case only if Π itself is close to equilibrium, whereas thermalization as formulated in expression ( 1) is supposed to hold generically, i.e. for almost every form of the statistical operator subject to the conditions 'large E' and δE ≪ E.

Correlation function
The stochastic process Tr(Aρ(t)) is real because the operators A, H, Π are Hermitian.Therefore, the correlation function is In calculating expression (29) we order the terms by the number ν of pairs of matrix elements O mα that are 'cross-contracted', with one element of the pair located in one trace and the other element in the other.The elements O mα are Gaussian random variables.Nonvanishing contributions may, therefore, arise only for ν = 0, 2, 4. Each term in expression (29) contains two factors U(t) and two factors U † (t).Contraction of the matrix elements O mα causes the appearance of ⟨U(t)⟩, of ⟨U † (t)⟩, of v 1 and v 2 defined in equations ( 23), and of the thirdand fourth-order cumulants of U(t).Our results are written in terms of these cumulants.As in section 5 we use a division of the energy scale into intervals of equal length ∆ that are consecutively labeled k = 1, 2, . ... For ν = 0, all four matrix elements O mα in either trace are contracted with each other, and nonzero contributions to the correlation function (29) arise only from the correlations of the eigenvalues E α .That situation is realized if in either trace both factors U and U † are replaced by their mean values (22) or by their covariances v 1 in equations (24), resulting in a total of four possibilities.When all four factors U, U † are replaced by their mean values, the Gaussian factors and the Kronecker deltas in equation ( 22) confine all summation indices in either trace to the same intervals labeled k 1 and k 2 , respectively.We find The index 'corr' indicates the correlated part.The uncorrelated part would give the product of two factors, each of the form of the first term on the right-hand side of equation ( 26), taken, respectively, at times t 1 and t 2 .Both these factors are of order unity in N k1 , N k2 .Correlations between the eigenvalues E α , E β ∈ k 1 and E γ , E δ ∈ k 2 vanish unless k 1 = k 2 .Within the interval k 1 = k 2 contributions to the big angular bracket arise from the GOE two-point correlation function [26] of pairs of eigenvalues , from the product of two such correlation functions, and from the GOE four-point function.As mentioned in section 5, the GOE two-point correlation function is small of order 1/N k1 .The four-point function is even smaller [26].Therefore, the contribution (30) to the correlation function ( 29) is negligible.The same argument shows that the remaining three of the four possibilities mentioned above expression (30) are likewise negligible.For ν = 2 and for ν = 4, cross-contraction of pairs of matrix elements O mα does give nonvanishing contributions.In each of these, the leading-order term in 1/N ∆ arises when correlations between eigenvalues E α , E β , E γ , E δ are neglected.Therefore, the summations over these eigenvalues may be carried out in the expressions for ⟨U(t)⟩ and in the higher-order cumulants of U prior to using these in the correlation function (29).We do that using the approximation introduced below equation (25).For ⟨U(t)⟩ in equation ( 22) and for v 1 , v 2 in equation ( 24) that gives While the index m in the expression for ⟨U(t) mn ⟩ is unrestricted, the indices (m, n) in the expressions for v 1 and v 2 are confined to the lie in the same interval of width ∆.In addition to the terms in equation ( 31), cross-contraction of the matrix elements O mα gives rise also to the appearance of the third-and fourth-order cumulants of U(t).In the appendix we show that contributions of these higher-order cumulants to the correlation function (29) are negligible.
In the present section, we confine ourselves to the contributions (31).We simplify the resulting expressions by suppressing the phase factors exp{±i E m t/ℏ} as these are immaterial for the counting.The resulting equations are still formidable.We further simplify the presentation.We only give the right-hand sides of the resulting equations.Thus, each one of expressions (33), ( 35) and (36) gives the value of expression (29) for the case defined in the text.For ν = 2, there are four possibilities.One factor U or U † in either trace is replaced by its average, the remaining two factors U or U † being replaced by v 1 or v 2 as the case may be.As a representative example we consider the replacements Using equations (31) we obtain Here and in the remainder of the section, the symbol [AΠ] mn = ∑ l A ml Π ln stands for the matrix product AΠ with unrestricted summation over intermediate states l.Without additional assumptions, the term (33) cannot be shown to be small of order 1/N ∆ .The term (33) has the same time dependence as the first term on the right-hand side of equation (28).It gives the statistical fluctuations in the ensemble of matrices H around that term.Analogous results are obtained for the other three possibilities beyond the choice (32).
For ν = 4, all four factors U, U † are replaced either by two factors v 1 or by two factors v 2 , the factors U or U † contributing to v 1 or v 2 being taken from either trace.A representative example is given by the replacements Every matrix element of the operators A and Π appearing in expression (35) connects states in the two intervals (k 1 , k 2 ).The second and third term in round brackets can be written as a single trace in the interval k 1 .These terms are of order N k1 and, multiplied by 1/[(2π) 2 ρ k1 ρ k2 ∆ 2 ], are negligible.The last term is equal to ( ∑ mn A mn A nm )( The first factor is of order N k1 .The second factor is of order unity.To see that we use Π = Π † and write the factor as ∑ ps |Π ps | 2 .That expression is bounded from above by Tr k1 (Π 2 ) < 1. Together with the factor 1/[(2π) 2 ρ k1 ρ k2 ∆ 2 ] the last term is negligible.That leaves us with Without additional assumptions, the contribution (36) to the correlation function (29) cannot be shown to be negligible.The term (36) is independent of time and describes the statistical fluctuations of the ensemble of Hamiltonians H around the last term in equation (26).
We conclude that of all the terms that contribute to expression (29), only the terms (33) and (36) are not negligible.Both these terms vanish, however, too in the large N limit if either the operator A or the operator Π are essentially diagonal in HF space.That condition says that for all intervals |A mn | 2 , and correspondingly for Π.For A, that condition is similar to a condition formulated in [6] and to the 'eigenvalue thermalization hypothesis' of [5].The condition for Π coincides with the condition for thermalization formulated in section 5. We conclude that if the system thermalizes, the statistical fluctuations around the terms in equation ( 28) are vanishingly small.That statement holds without constraining assumptions on the operator A. Then all systems described by the stochastic process Tr(Aρ(t)) thermalize in the same manner, including the system described by the Hamiltonian defined in the first part of section 4.

Violation of time-reversal invariance
With proper modifications, the arguments used in section 4 apply also to the case of timereversal non-invariant systems, and the results in sections 5 and 6 hold as well.The BGS conjecture now postulates local agreement of the spectral fluctuation properties of the Hamiltonian H with those of the Gaussian unitary ensemble (GUE) of random matrices [26].Proceeding as in section 4, we replace equation (10) by The matrices U are unitary.Both the eigenvalues E α and the eigenfunctions U mα obey local GUE statistics.Eigenvalues and eigenfunctions are uncorrelated.In an interval of width ∆, the eigenvalues follow the level statistics of the GUE, and the eigenfunctions U mα are zerocentered complex Gaussian random variables with second moments The function F(x) is given by equation (16).The contraction rules used in calculating the ensemble average and the correlation function of Tr(Aρ(t)) are changed.In the orthogonal case, every matrix element O mα is contracted with every other such matrix element.Because of equations (38), only contractions of U with U † occur in the present case.That reduces the number of contraction terms.For instance, the term v 2 in equations (23) does not arise, and the term proportional to A nm Π nm in equation ( 26) is lacking.The result (28) remains unchanged, however.The absence of a number of terms does not affect the structure of the results in section 6.The arguments used and conclusions drawn there remain the same.We conclude that relaxation and thermalization are completely similar for systems that are and systems that are not, invariant under time reversal.

Results, discussion, summary
We have investigated the thermalization of a closed chaotic many-body quantum system by combining the HF approach with the BGS conjecture.The HF approach provides the scaffolding for the spectrum of the many-body Hamiltonian H.It defines the integrable HF Hamiltonian H HF , the average HF level density ρ(E) as a function of excitation energy E, and it allows expectation values of physical operators to take different values in different parts of the spectrum.We use, in particular, that within the HF approximation the level density ρ(E) increases almost exponentially with energy.We assume that the total Hamiltonian, given by H = H HF + V, describes a chaotic quantum system.We accordingly use the BGS conjecture.
We consider the conjecture to be universally valid and to be synonymous with quantum chaos.The conjecture states that locally, the spectral fluctuation properties of H coincide with those of the random-matrix ensemble in the same symmetry class.The word 'locally' means that for every energy E of the system, the coincidence is restricted to an energy interval of width ∆ (the correlation width) centered on E. In the framework of the HF approach, the BGS conjecture requires that within the interval ∆, the mixing of HF states due to the residual interaction V is so strong that the resulting eigenstates of H approach characteristic statistical properties of GOE eigenstates.Such strong local mixing of the HF states requires V to be sufficiently strong.At the same time we have assumed that the overall structure of the spectrum and, in particular, the average level density ρ(E), are not affected by V. Pictorially speaking, V puts a thin veneer on the scaffold of the spectrum provided by the HF approach.
Combining the HF approach and the BGS conjecture we, thus, postulate that locally the statistical properties of eigenfunctions and eigenvalues of H coincide with those of the GOE in an interval defined by a Gaussian of width ∆.Actually, the statistical properties of the GOE that we use hold in the limit of infinite matrix dimension.The required coincidence (and, therefore, the validity of the BGS conjecture) depends on the number N ∆ = ∆ρ(E) of states in the interval ∆ and grows with N ∆ .Using essentially Fermi's golden rule, we have estimated the correlation width ∆.We have shown that ∆ depends on energy E much less strongly than ρ(E).Since ρ(E) grows essentially exponentially with E, the condition N ∆ ≫ 1 is best fulfilled in the semiclassical regime.In that regime, we may neglect terms of order 1/N ∆ .
Using the statistical properties of H in the evaluation of Tr(Aρ(t)) poses a challenge.We solve the problem with the help of the following construction.Changing the eigenvalues and eigenfunctions of H into random variables with GOE properties, we define an ensemble of Hamiltonians, all with the same statistical properties as the original Hamiltonian H. Written in terms of that ensemble, the time-dependent function Tr(Aρ(t)) turns into a stochastic process.We calculate statistical average and correlation function of that process in the limit N ∆ ≫ 1.The average ⟨Tr(Aρ(t))⟩ is given by our central result equation (28).Vanishing of the correlation function (29) guarantees that equation ( 28) holds for all members of the ensemble and, thus, for our original Hamiltonian.The condition for that to happen is addressed below.
Equation (28) shows that ⟨Tr(Aρ(t))⟩ relaxes in time with a Gaussian factor of width ∆.That factor is due to the Gaussian form of the correlation function (16).Other forms lead to different time-dependent functions.For example, the Lorentzian in equation ( 15) yields an exponential.All these forms have in common that the time scale ℏ/∆ is set by the correlation width ∆.
For time t → ∞, the right-hand side of equation ( 28) relaxes to the time-independent term.The term depends upon the operator A and on the time-independent statistical operator Π which defines the system S.The form of Π determines whether or not ⟨Tr(Aρ(t))⟩ thermalizes.For a general form of Π, the coefficients p k in equation ( 28) may take any value, and ⟨Tr(Aρ(t))⟩ differs from Tr(Aρ eq ).Moreover, evaluation of the correlation function in section 6 shows that in that case there exist statistical fluctuations (both time-independent and time dependent and of Gaussian form) around the mean value (28).
Thermalization occurs only if constraining conditions are imposed on Π.The standard requirement [2,3] is that the average energy E of the system be in the semiclassical regime, and that the spread δE in energy around that average be small.In the present context, the first condition requires that the diagonal elements of Π in the HF basis essentially differ from zero only in the domain of large HF energies.That validates the inequality N ∆ ≫ 1.The spread δE in energy, expressed in terms of Π, is estimated as the root-mean-square deviation of the energy from its average value E. The result shows that δE ⩾ ∆ so that ∆ is the minimum energy spread in the ensemble of random Hamiltonians.Thermalization is seen to occur universally (irrespective of the detailed form of Π) only if Π is confined to HF states within an energy interval of width ∆.Then the last term of equation ( 28) is equal to Tr(Aρ eq ).If that condition is met, the statistical fluctuations around the mean value ⟨Tr(Aρ(t))⟩ vanish.All members of the stochastic process Tr(Aρ(t)) (including the one containing our original Hamiltonian) tend asymptotically in time towards the equilibrium value Tr(Aρ eq ).These conclusions hold both for systems that are and for systems that are not invariant under time reversal.
From a technical point of view, our results apply not for all but only for almost all members of the random-matrix ensemble defined in section 4. Exceptions form a set of measure zero.The measure is defined by the integration measure for the evaluation of ensemble averages.An exception is obtained, for instance, when one chooses Π = |ψ ⟩⟨ψ | with |ψ ⟩ an eigenfunction of one particular realization of the Hamiltonian H drawn from the random-matrix ensemble.In that case, the density matrix ρ(t) in equation ( 3) is independent of time, the system does not thermalize.However, the averages taken in sections 5 and 6 extend over all realizations of the random-matrix ensemble.The function |ψ ⟩ is not an eigenfunction of H for most other realizations of H, and the case just considered carries vanishing weight in the averages.
Our approach leaves open several questions.We have assumed that the average level density of the full Hamiltonian is the same as that of the HF Hamiltonian.That puts an upper bound on the strength of V which we have not identified.We have not investigated collective effects.In atomic nuclei, for instance, the presence of collective states strongly modifies the average level density at low energy, see, for instance, [27].Does such an effect persist in the semiclassical regime?How do scars in the spectrum [28] affect our approach?We have assumed that two energy intervals (the one that limits agreement with the ∆ 3 statistics of eigenvalues and the one that limits complete mixing of eigenfunctions) are equal.That assumption requires further scrutiny.We have assumed that the spectral correlation function F(x) is Gaussian.Other forms have been considered in the literature.Which form applies for a given system?The answer determines the time dependence of relaxation.
Our results support, strengthen, and quantify some of the conclusions on thermalization drawn in [1], see especially section 7.In that paper, thermalization is mainly discussed in terms of chaotic eigenstates and time averages.Semiquantitative arguments are supported by numerical examples.Thermalization is shown to require, in our notation, ∆ ≈ δE.That agrees with our findings.Our work goes beyond [1].The statistical approach allows us to derive exact conditions and exact results for thermalization.In particular, we obtain an explicit expression for the time dependence of thermalization.We determine the characteristic time scale.
Our approach differs from, and our results go beyond those of [2,3,5].We deduce thermalization from constraints on the operator Π and from the statistical properties of the Hamiltonian.The latter, in turn, follow stringently from the BGS conjecture.The conjecture holds within the correlation width ∆.That parameter plays the central role in our approach.It determines the characteristic time scale for relaxation.It also determines the value of the energy spread δE required for thermalization.The parameter ∆ plays no role in [2,3,5].Our approach makes a prediction that can, in principle, be tested experimentally: The correlation width ∆, defined by the spectral properties of the Hamiltonian, determines the time scale ℏ/∆ for relaxation of Tr(Aρ(E)).
We believe that our approach is not confined to chaotic many-body systems but applies to chaotic quantum systems in general.Certainly, the BGS conjecture holds, and the limiting interval ∆ occurs, generically, cf the examples in section 2. The possibility to construct an integrable Hamiltonian that defines the scaffolding of the spectrum is, likewise, probably generic.Consider, for instance, the quantum Sinai billiard, a square with an inscribed concentric circle for which the eigenfunctions of the Laplacian operator in two dimensions vanish on both surfaces.Replacing the circle by a polygon with n corners and the same area as the circle one obtains an integrable system the spectrum of which for n ≫ 1 provides the scaffolding for the spectrum of the Sinai billiard.The geometric difference between polygon and circle is the analogue of the residual interaction V in the present paper, mixing close-lying eigenstates of the polygon within the interval ∆.If generic, the possibility to find for every chaotic quantum system an integrable Hamiltonian that differs from the chaotic one only by a weak residual interaction, would qualitatively explain why the mixing of the eigenstates of the integrable Hamiltonian is local and confined to the correlation width ∆.
The correlation width ∆ is present universally in all chaotic quantum systems that obey the BGS conjecture.We are not aware of any other such common energy scale that would qualify for defining the time scale for relaxation.That raises the question whether ∆ possesses a classical counterpart that is present equally universally in chaotic classical systems.In [11], the upbend of the ∆ 3 statistics for the quantum Sinai billiard was related to the largest Ljapunov coefficient of the classical Sinai billiard.We believe that this is an interesting result.It raises the question whether the observed relation is generic.terms give negligible contributions implies that the terms that actually do occur are negligible as well.Each trace is at most of order k.Combined with the factor 1/(ρ k ∆) 3 that shows that terms in expression (4) which carry a single trace or two traces, are negligible.In the terms carrying three traces, one trace has either the form Tr k Π = p k , or the form Tr k (Π) 2 < p k .Both traces are of order unity.Combined with the factor 1/(ρ k ∆) 3 , the terms carrying three traces are negligible.In the term carrying four traces we also use Tr k Π = p k .The remaining two traces are of order (ρ k ∆) each.Combined with the factor 1/(ρ k ∆) 3 they are negligible.We conclude that contributions of the fourth-order cumulant of U to the correlation function (29) are negligible.

ORCID iD
Hans A Weidenmüller  https://orcid.org/0000-0001-8032-9524 eigenvalues E m of H HF with m = 1, 2, . . .are ordered so that E m ⩽ E n for m < n.Partial degeneracies are removed by diagonalization of V in the subspace spanned by the degenerate HF states.Without changing labels, we then have E m < E n for all m < n.