Thermalization in the Two-Body Random Ensemble
Using the ergodicity principle for the expectation values of several types of observables, we investigate the thermalization process in isolated fermionic systems. These are described by the two-body random ensemble, which is a paradigmatic model to study quantum chaos and specially the dynamical transition from integrability to chaos. By means of exact diagonalizations we analyze the relevance of the eigenstate thermalization hypothesis as well as the influence of other factors, like the energy and structure of the initial state, or the dimension of the Hilbert space. We also obtain analytical expressions linking the degree of thermalization for a given observable with the so-called number of principal components for transition strengths originated at a given energy, with the dimensions of the whole Hilbert space and microcanonical energy shell, and with the correlations generated by the observable. As the strength of the residual interaction is increased an order-to-chaos transition takes place, and we show that the onset of Wigner spectral fluctuations, which is the standard signature of chaos, is not sufficient to guarantee thermalization in finite systems. When all the signatures of chaos are fulfilled, including the quasi complete delocalization of eigenfunctions, the eigenstate thermalization hypothesis is the mechanism responsible for the thermalization of certain types of observables, such as (linear combinations of) occupancies and strength function operators. Our results also suggest that fully chaotic systems will thermalize relative to most observables in the thermodynamic limit.
Keywords: Connections between chaos and statistical physics, Matrix models, Quantum chaos
The success of thermodynamics is based on the fact that the state of a macroscopic system behaves as it were independent of the microscopic details of the system, and thus it can be determined from a few universal laws. In classical isolated systems this thermodynamic universality can be explained satisfactorily if one assumes that the dynamics is ergodic and mixing. These properties define classical deterministic chaos  and lead to the equiprobability of equal-volume regions of the available phase space, which is the starting point of the microcanonical ensemble formulation. This implies that, on average, the state of the system is independent of the initial conditions as well as independent of the measurement time. It is thus generally accepted that the system equilibrates into a thermal state, described by the microcanonical ensemble, when its dynamical regime is ergodic. On the other hand thermalization is expected to be inhibited for integrable or quasi integrable systems. In these systems (almost) every trajectory in phase space lies in one of the so-called invariant tori, which foliate the phase space. Therefore, being restricted the trajectories to these structures, the equiprobability equal-volume regions of the available phase space does not hold and the system will not thermalize, at least in the sense explained above.
In recent years the study of the equilibration and thermalization mechanisms in isolated quantum systems has attracted a great interest partly because the non-equilibrium dynamics, after an external perturbation is applied, has become experimentally accessible for ultra-cold quantum gases [2, 3] and electrons in solids . The technology makes it possible to induce sharp changes in the parameters controlling the system and then observe the subsequent time evolution, which is essentially unitary because on short and intermediate time scales the perturbed system is almost isolated from the environment. Thus, one can experimentally study if an isolated system equilibrates after a sharp perturbation and in this case whether it thermalizes or retains memory of the initial conditions. Experimental studies of the non-equilibrium dynamics in one-dimensional ultra-cold Bose gases have given, up to the moment, contradictory results [5, 6]
From a theoretical point of view, these questions have been addressed using different methods and points of view. The results from classical mechanics can not be translated directly into quantum mechanics. On one hand the concept of quantum integrability is not well defined, though it is generally considered as synonym for exact solvability [7, 8]. On the other hand the unitary time evolution of quantum states leaves no room for a dynamical regime similar to deterministic chaos. Actually, the name quantum chaos stands for the different type of signatures that certain quantum systems exhibit depending on the large time scale behavior of their classical analogues [9, 10]. These signatures appear on the statistical behavior of eigenenergies and eigenfunctions and is not clear at all how they can influence the evolution of the system on large-time scales.
The generally accepted assumption that integrable systems do not thermalize is corroborated in several models [11, 12, 13, 14, 15], where the non-equilibrium regime extends over large time scales or where a long-time steady state with non standard thermal properties is reached. The approach to thermal equilibrium of generic systems has been studied by several authors. It has been proven that almost any system in interaction with a large heat bath will equilibrate and thermalize [16, 17]. For isolated systems it has been shown that “typical” Hamiltonian and observables will be in thermal equilibrium for most times [18, 19], but it seems a very difficult task to decide whether a specific Hamiltonian or observable belongs to this class or not. The thermalization of specific 1D and 2D fermionic and bosonic systems has also been studied [20, 21, 22, 23, 24]. Surprisingly, thermalization was not obtained in all the cases. Several reasons have been reported to explain the lack of thermalization in these systems, but it seems that the so-called eigenstate thermalization hypothesis (ETH) plays a fundamental role [23, 24]. The ETH states that thermalization occurs at the level of individual eigenstates [25, 26], whenever they satisfy Berry’s conjecture on chaotic eigenfunctions , i.e., whenever they behave as (quasi) random superpositions of the basis states. For this and other reasons, the role played in the viability of thermalization by quantum chaos in general, and by the properties of chaotic wave functions in particular, has began to receive some attention [28, 29]. It seems that the number of principal components (NPC) or inverse participation ratio (IPR), which keeps track of the progressive eigenstate delocalization through the integrability to chaos transition, might be directly related to the deviations of the steady expectation values from the corresponding statistical values [30, 31]
In order to get a deeper understanding of the role played by quantum chaos in these processes we study the thermalization of isolated fermionic systems modeled by the simplest two-body random ensemble . Two-body random ensembles are paradigmatic models to study quantum chaos and specially the dynamical transition from integrability to chaos. They have also been used in the past, together with some related models, such as nuclear shell model and interacting spin models, to perform different studies on thermalization. The thermalization criteria were based on the equivalence between different definitions of entropy [33, 34, 35] and temperature ; representability of occupancies by Fermi-Dirac distribution [37, 38, 39] (Bose-Einstein distribution for bosons); and calculation of expectation values using the canonical distribution . A brief review of some of these studies can be found in Refs. [40, 41]. However, in the present work, as well as in most recent papers the focus is put on the Ergodicity principle [18, 19] which is the cornerstone for thermalization, and clearly more precise and general that the aforementioned criteria. We study the relevance of several factors in the thermalization process like ETH, the dimension of the Hilbert space, the structure and energy of the initial state and its proximity to the ground state. We also analyze the importance of the degree of chaos as measured by the different chaos markers.
The rest of the paper is organized as follows. Sec. 2 briefly introduces the embedded Gaussian ensembles of random matrices generated by two-body interactions, with emphasis in the embedded Gaussian orthogonal ensemble, and some of the well established main features of the order to chaos transition in this ensemble. However, the criteria for determining the transition points as the two-body interaction strength is increased, are somewhat different from those used in the past. Sec. 3 deals with the thermalization of embedded Gaussian orthogonal ensembles. Here we give all the new results of the paper, where the transition points of Sec 2.2 play a central role. After introducing some basic definitions in Sec. 3.1, Sec. 3.2 reports the main numerical results linking thermalization with the type of spectral fluctuations, the delocalization of the wave functions, the structure and the proximity of the initial state to the ground state, or the dimension of the Hilbert space. In Sec. 4 we gather together some analytical results that establish a connection between thermalization, relative to a given observable , and the value of the NPC for the transition strengths originated at a given eigenstate, or between thermalization and the dimension of the whole Hilbert space, the dimension of the microcanonical energy shell and the correlations generated by . Finally, Sec. 5 contains the conclusions.
2 EGOE(1+2) model: order-to-chaos transitions
2.1 EGOE(1+2) model
As previously mentioned, we try to analyze the relationship between the order-to-chaos transition and the thermalization of fermionic systems, which can be modeled by an appropriate two-body random matrix ensemble. Here we briefly introduce these ensembles and their relation with quantum chaos. Recent and comprehensive reviews can be found in references [40, 41].
There is a clear relationship between the energy level fluctuation properties of a quantum system and the large time scale behavior of its classical analogue. The spectral fluctuations of a quantum system whose classical analogue is fully integrable are well described by Poisson statistics, i.e., the spacings between successive energy levels are not correlated . According to Bohigas, Gianoni and Schmit , the fluctuation properties of generic quantum systems, which in the classical limit are ergodic, coincide with those of the Gaussian ensembles (GE) of random matrices . This statement, initially supported by many experimental data and numerical calculations, has been finally proven for quantum systems with few degrees of freedom, where the semiclassical approximation is valid . The large time scale behavior of the classical analogue also determines the properties of the wave functions. Extensive reviews of later developments can be found in [10, 41] and references there in.
Quantum many-body systems like complex atoms and atomic nuclei are usually considered to be chaotic if their spectral fluctuations are those of the Gaussian orthogonal ensemble (GOE), which is the appropriate GE for systems with time-reversal invariance and rotational symmetry. However, real quantum systems are usually well described by real or effective one- plus two-body interactions in the mean-field basis, whilst GE represent systems with multi-body interactions. The embedded Gaussian ensembles (EGE) of random matrices were introduced to tackle this problem, and to provide a more realistic picture of many-body quantum systems. Moreover, in the present context EGE are interesting because random interactions can illustrate the effects on thermalization caused by generic interactions, which lead the system from integrability to chaos.
The EGE(1+2) ensembles consider fermions or bosons distributed in single-particle states , interacting via the following Hamiltonian
where the single-particle energies and the two-body matrix elements (properly symmetrized or antisymmetrized) behave as independent Gaussian random variables. In this expression gives the strength of the two-body interaction, and and create and destroy a fermion (or a boson) in the th single-particle state. In fermionic systems only the strict inequalities and are valid. It has been numerically shown that for EGE(2) ensembles () the spectral fluctuations do agree with those of GE, provided that energy scale is redefined appropriately. Introducing the notation , where stands for the trace operation, and is the dimension of the Hilbert space, the centroid and the energy span must be and , respectively. 
2.2 The order-to-chaos transition in EGOE(1+2)
One of the most significant aspects of EGOE(1+2) is that as increases, starting from , the system undergoes a transition from a regular to a chaotic regime that deeply affects the state density, level fluctuations, and wave functions. This change of dynamical regime is characterized by three chaos markers [40, 41, 46, 47]. There is a first marker that signals the transition from Poisson to GOE spectral fluctuations. This transition occurs when the interaction strength is of the order of the spacing between the basis states that are directly coupled by the residual two-body interaction, a result that came out of nuclear structure calculations by Åberg [48, 49]. More developments in determining the marker are given for example in [50, 51]. An important outcome of Åberg criterion is that for EGOE(1+2), which is well verified in  and in the results presented below.
As increases further from , the structure of the eigenstates undergoes a deep transformation. First the strength functions change from Breit-Wigner to Gaussian at a transition point denoted by . Beyond we find a third chaos marker which defines the center of a region where different definitions of thermodynamic variables, such as entropy, temperature, specific heat, etc., give the same results, as it occurs for infinite systems. As far as the statistical entropy and the Shanon entropy are concerned we can understand the meaning of as follows. The former is proportional to the logarithm of the state density, which in our case has essentially Gaussian form, though its mean and variance depend on . At , where the eigenstates are fully localized in the mean-field basis, information entropy is , whilst for sufficiently large values of the interaction strength, the eigenstates are quite similar to those of GOE. In our case this occurs for and then , except perhaps near the spectrum edges. Changing from to the wave functions become more and more delocalized in the mean-field basis, being this process faster in the middle than in the spectrum edges, and signals the precise value of the residual interaction strength in which the energy dependence of and is quite similar. Thus, if we define an appropriate “distance” between the two entropies, it should have a minimum at . Another property is that signals the duality point between and basis, i.e., the point where the eigenstates become equally delocalized in the two basis [47, 50]. Beyond the eigenstates become quickly similar to those of GOE systems, for which they are essentially Gaussian random superposition of the basis states. This property resembles Berry’s conjecture about the ergodic structure of chaotic wave functions in phase space . Nevertheless, a word of caution should be included here. Since is basis dependent, for it would be more appropriate to define in the basis.
For our purposes the markers and are the most relevant. We detect the change in the spectral fluctuations, and thus the position of , using the nearest-neighbor spacing distribution, denoted . The spacings for generic integrable systems obey the Poisson distribution, i.e, , while the Wigner surmise, , provides a very good approximation for GOE-like systems . The Brody distribution ,
where is usually called the Brody parameter and is a normalization constant, it is used to asses how close the fluctuations are to the Poisson limit, corresponding to , or to the Wigner surmise with . Although Eq. (2) is only a heuristic formula, and the Brody parameter has no definite meaning for Hamiltonian systems, it has been employed in a variety of studies since its introduction. Very recently a physical foundation of has been found by Sakhr and Nieminen in the context of self-similar fractals . As the interaction strength in Eq. (1) increases from to a sufficiently large value, the Brody parameter changes from to . In what follows, the position of the first chaos marker is fixed by the condition . Similar conditions can be found in the literature [51, 54].
To locate the value of the third marker we consider the values of three different entropies:
Thermodynamic entropy, , where is the density of states.
Information entropy in the mean-field basis , where the coefficients are the eigenstates components in the mean-field basis.
Single-particle entropy , with the occupancy of the th single-particle state at energy .
where , the value of corresponds to the minimum of because this ensures that the values of the different entropies will be very close to each other.
Fig. 1 displays the ensemble averages and along , for a member EGOE(1+2) with and . In order to enlarge the region where the order-to-chaos transition takes place we use a logarithmic horizontal axis. The two vertical dashed lines indicate the respective positions of and . It can be seen that and , values that are consistent with theoretical estimates (see [40, 41, 46]). For we obtain , which is very close to the actual GOE result. However, the structure of these states is still very different to those of a GOE system because . Only when one finds that , signaling that the dynamics has become fully chaotic.
3 Thermalization definitions and numerical results
3.1 Basic definitions
We study the thermalization properties of finite fermionic systems with time-reversal and rotational invariance. To be precise we consider fermions distributed in independent particle states, interacting via the Hamiltonian (1), where and are independent real Gaussian random variables. These systems are usually called EGOE(1+2) ensembles. In this work we use , , , , and the energy scale is such that and , regardless of the value of the interaction strength . The results presented below have been obtained by fully diagonalizing member EGOE(1+2) systems with , and . The corresponding dimensions are given in Table 1.
Let be the initial state of the system, that we decompose as
in the eigenstate basis of the Hamiltonian (1). Then, given a certain observable we define:
The instantaneous value , where .
The time average . When , , where . The steady state average can also be written as , with .
The statistical average , where is the density operator corresponding to an appropriate statistical ensemble.
We say that the system thermalizes if for any relevant observable and almost any state , it is satisfied that
In practice, to asses whether Eq. (5) do approximately holds for a specific observable, one can use the relative error
which is well suited to compare the degree of thermalization of different systems, or the thermalization of a single system relative to distinct observables.
Let us consider that the system is prepared in a non-equilibrium state for which the energy is essentially constant. Hamiltonian eigenvectors are excluded as they are stationary states, but we can assume that , and , with sufficiently small compared to the energy spectrum span, but large enough to contain many energy eigenstates. In such a case the microcanonical ensemble is the preferred statistical ensemble. Denoting by the corresponding energy shell, i.e., , the microcanonical density operator is
where is the dimension of the subspace and the symbol means that the sum is restricted to eigenstates belonging to . Thus, the corresponding microcanonical average is given by
To examine whether the energy plays a role in the thermalization of the system we shall use three different energy intervals with and and . We denote the corresponding energy shells by , and respectively. Moreover, to see how the initial conditions affect the process we let the system evolve from three different types of initial states, defined as:
, where is the projector onto , and is the mean-field state with energy .
, with the coefficients Gaussian random variables with mean zero and variance equal to one, i.e., .
, where the expansion coefficients are .
The states and are random superpositions of the eigenstates belonging to , but due to the Gaussian factor the distribution of the coefficients is wider for . As we shall see below the distribution of the state amplitudes inside is one of the factors that may affect the thermalization process.
Before we turn to the main results, let us introduce the four types of generic observables for which we study the validity of Eq. (5):
diagonal one-body operators ,
one-body operators ,
two-body operators ,
strength function operators ,
where the parameters , and are taken as random variables.
3.2 Numerical results
Let us consider the ensemble averages , , , and for a member EGOE with . Fig. 2 displays their evolution with the strength when the system is prepared at in a state with energy ( ). Because of the large differences between the relative errors of these operators and to properly visualize their evolution in the short interval , a log-log scale is used. As for Fig. 1 the two vertical lines give the positions of and . In all the cases becomes smaller as the interaction strength increases up to . It is very important to realize that the transition from Poisson to GOE spectral fluctuations, which is considered the most relevant signature of quantum chaos, occurs at and does not modify this trend. On the contrary, for the relative errors either remain essentially constant or the decreasing rate is much smaller. Recall that defines a region where the three entropies , and take essentially the same values, and signals the point at which the wave functions start to become very delocalized in the mean-field basis. Beyond , becomes clearly smaller than one percent only for two operators, namely and . Their errors are % and %, respectively. Thus, as long as the the system is prepared in an initial state and , Eq. (5) approximately holds for the observables and and then we can assert that the system thermalizes relative to these observables. This is not the case of the observables and . It is worth noting that perhaps the main difference between , in one hand, and , in the other, is that the latter have meaningful smoothed form for large , as given by spectral distribution methods [41, 56, 57].
The results corresponding to other choices of are shown in Fig. 3. For simplicity we have represented only the results for and . Curves in black, red and green correspond to , and , respectively. In all cases the initial state belongs to the energy shell . We see that the choice of the initial conditions does not affect the main trend: , and diminish progressively as the strength is increased, and once their values remain essentially constant. However, the precise values are quite different. When , the initial states and give rise to very similar results while the error corresponding to is clearly larger.
To shed some light on the mechanisms leading to these results lets us note that the difference
becomes very small when: the distribution of the coefficients is nearly flat, i.e, , or the matrix elements almost do not fluctuate inside the energy shell , and therefore . This condition is known in the literature as the eigenstate thermalization hypothesis (ETH), which has been conjectured to hold in chaotic quantum systems [22, 26]. The standard deviation of the initial state components in the energy eigenbasis, , is plotted in Fig. 4. When is very small takes very different values in the three cases, being quite larger for . However, in this case undergoes a sharp decreasing up to , where its values are very similar to those of and half those of .
Defining the standard deviation of the expectation values inside an energy shell as
the ratio provides us a measure of the fluctuations of these matrix elements. Fig. 5 plots for a randomly selected member of the ensemble inside . Throughout the interval to the ratio is reduced sharply by a factor ranging between 5 and 40, but even then we get for operators like and . Remarkably, for or , meaning that the fluctuations of their expectation values are small compared to in the energy shell .
The extraordinary resemblance of Fig. 5 with Figs. 2 and 3 shows that the fluctuations of the eigenstate expectation values determine for which observables thermalization occurs. In other words, the thermalization of the system is correlated with the degree of compliance of the ETH. Moreover, given a certain observable the exact degree of thermalization depends on the width of the initial state: the wider is the initial state the worse thermalizes the system, a fact that is in agreement with some previous claims .
Thermalization may be affected by other factors, such as the proximity of to the edges of the spectrum, and in particular to the ground-state energy. In order to explore this question we have compared the behavior of in three energy shells , and , of width and and , respectively. Although the energies of the largest eigenstates oscillate between and depending on the value of , it is very difficult to get closer to the spectrum end. The reason is that outside the central interval the state density is so scarce that it is impossible to obtain a narrow energy shell which contains a large number of states , but small compared to the dimension of the whole space. We plot in Fig. 6 the variation with of the two averages and , calculated for an EGOE system with , which is initially prepared in a state of type . One sees immediately that behaves very differently for these operators. Once the thermalization region () has been reached, the relative error grows as the energy approaches the end of the energy spectrum. By contrast, the evolution of is more irregular, but still this error is clearly larger for and . Therefore, the proximity of the initial state energy to the spectrum edges inhibits the thermalization of the system.
Before closing this section we shall briefly consider whether the dimension of the Hilbert space hinders or enhances thermalization. For illustration we consider systems with and , with dimensions (see Table 1) and , respectively. Fig. 7 shows the evolution of for and as the dynamical regime changes from regularity at to full chaos for . The growth rate of is very small, but it is enough to induce a gentle decrease in , while maintaining the same qualitative trend. Similar results are obtained for the other two observables and . If we also take into account the analytical results of the next section we can conclude that the system will always thermalize in the thermodynamic limit (the difference between microcanonical and diagonal predictions decreases with system size), but there may be relevant differences for finite dimensions, which is important since many physical systems of interest are mesoscopic.
4 Analytical results for thermalization
Our goal now is to relate the relative error with the fluctuation properties of the eigenstates and with the correlations generated by the observable . To this end let us focus on the quantity
with . It can be considered as a random variable because we model the physical system by means of an appropriate ensemble and thus the external products appearing in change from one member of the Hamiltonian ensemble to another. Moreover, being interested in “typical” properties we can introduce a fictitious ensemble of initial states, , with similar energies. Assuming that they are uniformly distributed in the unit sphere in , the fluctuation properties of the coefficients obey the Porter-Thomas (P-T) distribution  whenever is large enough. Thus
where means averaging over the fictitious ensemble of initial states. Therefore implying that . Then, is the variance of (9) and we may agree on defining the “typical” value of as
Introducing the shifted operator , with , it is easy to see that . Moreover, a straightforward calculation gives
and taking into account the fluctuation properties of the coefficients we obtain
Introducing this result in (11), the value of is given by the expression
which can be further simplified because the averaged forms are smooth functions of the energy that do not change substantially inside the energy window . Thus, we arrive to
Since is a free parameter we can choose , which gives rise to
For strength function operators, like , the value of can be related to the NPC (also called IPR) in transition strengths originating from the central eigenstate . Before we proceed it seems suitable to define the usual NPC for eigenstates and its extension for transition strengths. The former is a measure of the eigenstate complexity in the mean-field basis of Slater determinants with and . Expanding the eigenstates in this basis, i.e.,
where the amplitudes satisfy that , the number of principal components is defined as
Note that NPC, as a function of the amplitudes , attains its absolute minimum whenever the wave function is localized in a single basis state, and its absolute maximum when the wave function is completely delocalized and all the amplitudes are equal to . Generally speaking NPC gives the effective number of basis states that build up Hamiltonian eigenstates. Thus, it is small for localized states, while for chaotic states, which are quite delocalized in the mean-field basis, it takes values comparable but somewhat smaller than because system symmetries and orthogonality prevent reaching its maximum. For instance, the average value for GOE is .
Since the analysis of eigenvector amplitudes is largely equivalent to dealing with transition strengths we can extend the previous discussion. Using the standard notation of spectral distribution methods, let us introduce the so called locally normalized strength , and the normalized stregth generated by the action of the operator on a certain eigenstate as
where is the total strength sum. In our notation . Note that the initial and final spaces connected by do not need to be the same, but for simplicity we shall consider here a single space. If we normalize the state vector and expand it in the Hamiltonian eigenbasis
then the NPC for the strength distribution generated by on the eigenstate is defined as
Contrary to usual NPC in wavefunctions it only depends on the operator and the initial state, but it is basis independent. The NPC for strength distributions gives the effective number of eigenstates over which the strength generated by the action of the operator on a given eigenstate is spread. For operators generating collective states (then and are highly correlated) NPC should be small, while for operators generating chaotic states it should take large values.
Before turning back to the calculation of we give an expression for the smoothed NPC that will be used below. If we assume that the fluctuations in the locally renormalized strengths follow the P-T distribution then , , and
An EGOE(1+2) formula for NPC in transitions strengths and not too far from is given in Ref. . The typical error can be written as
The step from the first to the second line in the previous expression follows from the independence of the strengths due to the P-T assumption, and the third line is obtained by adding and subtracting terms with . Using again that , this result can be simplified as
This expression was reported for the first time in  (see also ), and its predictions have been compared with shell model results in [40, 60] for the width of the strength sums fluctuations. Finally, inserting (26) into (23), one gets
which establish a connection between the thermalization of the system, relative to an observable , and the value of the NPC for the transition strengths generated by acting on the eigenstate with energy . An important outcome is that for chaotic systems the NPC is expected to be large and hence these system will thermalize, while for regular systems NPC has to be small and thus thermalization will be hindered.
It is also possible to derive an expression for generic operators when . In this case the eigenstates become more and more delocalized in the mean-field basis. Finally we reach a regime were they behave as (quasi) random combinations of the basis states, the squared amplitudes abide the P-T distribution and (). In this regime the expression (16), valid for a generic observable, and (27) for , take on a very simple form. Using Eq. (17) to expand in the mean-field basis, one gets
where the usual P-T formulas for the moments, and , have been applied.
Gathering these expressions together with (16) one finally arrives to
that links the value of with the dimension of the whole Hilbert space, the dimension of the microcanonical energy shell and the correlations generated by the operator , which are measured here by the ratio . Fig. 8 displays the values of calculated for the four observables and in EGOE(1+2) systems with and . These are prepared in three different initial states