Process of equilibration in many-body isolated systems: diagonal versus thermodynamic entropy

As recently manifested [], the quench dynamics of isolated quantum systems consisting of a finite number of particles, is characterized by an exponential spreading of wave packets in the many-body Hilbert space. This happens when the inter-particle interaction is strong enough, thus resulting in a chaotic structure of the many-body eigenstates considered in the non-interacting basis. The semi-analytical approach used here, allows one to estimate the rate of the exponential growth as well as the relaxation time, after which the equilibration (thermalization) emerges. The key ingredient parameter in the description of this process is the width Γ of the local density of states (LDoS) defined by the initially excited state, the number of particles and the interaction strength. In this paper we show that apart from the meaning of Γ as the decay rate of survival probability, the width of the LDoS is directly related to the diagonal entropy and the latter can be linked to the thermodynamic entropy of a system equilibrium state emerging after the complete relaxation. The analytical expression relating the two entropies is derived phenomenologically and numerically confirmed in a model of bosons with random two-body interaction, as well as in a deterministic model which becomes completely integrable in the continuous limit.


Introduction
The problem of thermalization in isolated quantum systems remains a hot topic in the field of modern statistical mechanics. It has been shown since long that thermalization can emerge without the presence of a heat bath, even if the number of interacting particles is small [2][3][4][5][6][7]. One of the main problems in this field is to establish the conditions under which a given system manifests strong statistical properties, such as the relaxation to equilibrium. Due to a remarkable progress in the study of this and related problems, much is already understood in theoretical and numerical approaches (see [8][9][10][11][12] and references therein) as well as in experimental studies of interacting particles in optical traps [13][14][15][16][17], even if many basic problems still need further intensive efforts.
As was shown in reference [18], the mechanism for the onset of statistical behavior for a quantum isolated system is the chaotic structure of the many-body eigenstates in a given basis defined in absence of inter-particle interaction. As for the properties of eigenvalues (such as the Wigner-Dyson type of the nearest-neighbor level spacings distribution), it was shown [19,20] that they have a little impact on the global statistical properties of the wave packet dynamics.
The concept of chaotic eigenstates originates from random matrix theory (RMT), which was suggested by Wigner in order to explain local properties of energy spectra of heavy nuclei, observed experimentally [21,22]. Being unable to describe global properties of energy spectra of complex physical systems, the random matrices turned out to be effective models for the description of the local statistical properties of number of many-body states, this process terminates when all states are excited, which create an energy shell in the Hilbert space (for details see [18]).
Numerical data obtained for the TBRI model with a finite number of bosons occupying a number of single-particle energy levels [1], as well as for models of a finite number of interacting spins-1/2 in a finite length chain [19,20], clearly demonstrated that, in presence of chaotic eigenstates, the global time behavior of the systems is very similar. One has to note that even if one of the spin models is integrable, with a Poisson-like level spacing distribution, this does not influence the quench dynamics. These results confirm the prediction that for many-body systems the type of energy level fluctuations is less important than the chaotic structure of the many-body states.
In this paper we continue the study of the quench dynamics, paying attention to the new question of the relevance between the diagonal entropy related to an initially excited state, and the thermodynamic entropy emerging in the process of thermalization. We show numerically and describe semi-analytically that there is a one-to-one correspondence between them, with some corrections due to the different size between interacting and non-interacting energy spectra. This remarkable result holds both for the TBRI model with finite number of bosons, and for the model originated from the celebrated Lieb-Liniger (LL) model [43][44][45]. The latter model, which has no random parameters, was proved to be integrable with the use of Bethe ansatz.
The LL model describes one-dimensional (1D) bosons on a circle interacting with a two-body point-like interaction. It belongs to a peculiar class of quantum integrable models solved by the Bethe ansatz [46,47]; in particular, it possible to show that it has an infinite number of conserved quantities. Apart from the theoretical interest, this model is important in view of various experiments with atomic gases [48][49][50]. For a weak inter-particle interaction the LL model can be described in the mean-field (MF) approximation. Contrarily, for a strong interaction, the 1D atomic gas enters the so-called Tonks-Girardeau (TG) regime in which the density of the interacting bosons becomes identical to that of non-interacting fermions (keeping, however, the bosonic symmetry for the wave function) [45]. The crossover from one regime to the other is governed by the ratio n/g between the boson density n and the interaction strength g. The latter constant is inversely proportional to the 1D inter-atomic scattering length and can be experimentally tuned with the use of the Feshbach resonance (see, for example, [51] and references therein). Specifically, the MF regime occurs for n/g 1 and the TG regime emerges for n/g 1 [45]. In our numerical study, we consider a finite many-body Hilbert space by fixing the total number of momentum states and the number N of interacting bosons. This truncated Lieb-Liniger model (t-LL) allows one to correctly obtain both the eigenvalues and many-body eigenstates of the original LL model, in a given range of energy spectrum [52]. This method (truncation of the infinite spectrum) can be considered as the complimentary one to the recently suggested way [53], according to which a finite number of eigenstates involving into the quench dynamics is used.

Quench dynamics
Below, we study the quench dynamics in two different models described by the Hamiltonians written as Here the first term H 0 describes N non-interacting identical bosons occupying M single-particle levels specified by the energies s , while the second term V stands for the two-body interaction between the particles, Heren s =â † sâs is the number of particles in the corresponding s−level, withâ † s and â s as the creation/annihilation operators acting on the sth level, and the two-body matrix elements V s 1 s 2 s 3 s 4 specify the type and strength of the inter-particle interaction. The interaction conserves the number of bosons and connects many-body states that differ by the exchange of at most two particles. Due to the two-body nature of the interaction, the matrix elements H ij = i|H|j are non-zero only when the two non-interacting many-body basis states |i and |j have single-particle occupations which differ by no more than two units. This means that the Hamiltonian matrix is sparse which is a common property of realistic many-body systems.
Our aim is to compare the dynamical properties of two models: (a) the TBRI model with completely random and independent values of V s 1 s 2 s 3 s 4 (b) the t-LL model which is deterministic (without random parameters) and integrable in the limit M → ∞.
The non-interacting many-body eigenstates |k (basis states) of are obtained by creating N bosons in M single-particles energy levels, so that |k = a † s 1 . . . a † s N |0 , where 1 s 1 , ..., s N M. In this way, we determine the non-interacting many-particles basis, in which the Hamiltonian is a diagonal matrix. As a result, the many-body eigenstates |α of the total Hamiltonian H are represented in terms of the basis states |k as where C α k are obtained by exact numerical diagonalization. Below we use Latin letters when referring to the non-interacting basis of H 0 , while Greek ones stand for the basis of H.
After specifying the non-interacting basis |k , one can study the wave packet dynamics in this basis, after switching on the two-body interaction V. Our analytical and numerical results refer to the situation when initially the system is prepared in a particular eigenstate of H 0 , and evolved according to the full interacting Hamiltonian H. To do this, we focus on the time dependence of the effective number N pc (t) of principal components of the wave function, known in literature as the participation ratio.
In terms of eigenvalues and many-body eigenstates of the Hamiltonian H the participation ratio can be presented as where and are the diagonal and off-diagonal parts of As was recently shown for different models with chaotic behavior [1], the number N pc increases exponentially fast in time, provided the eigenstates involved in the dynamics, are strongly chaotic. After some relaxation time t s , the value of N pc fluctuates around the saturation value due to the complete filling of a portion of the total Hilbert space, called energy shell. Such a wave packet dynamics occurs when the many-body eigenstates of H are fully delocalized in the energy shell. This scenario explains the basic properties of the quench dynamics, before and after the saturation.

Semi-analytical approach
The semi-analytical approach developed in [4] for TBRI matrices with an infinite number of interacting Fermi-particles and modified in [1] for finite Bose-systems, allows one to obtain simple estimates for two important characteristics of the quench dynamics. The first characteristic is the rate of the exponential growth for the N pc (t). The second characteristic is the time scale t s on which the exponential growth of N pc (t) occurs.
The two above characteristics can be estimated with the use of the semi-analytical approach originally developed in [54] for an infinite number of particles. To start with, we write down an infinite set of probability conservation equations, used for the description of the probability flow W(t) in the many-body Hilbert space of H 0 . To this end, let us consider an initial non-interacting state |k 0 , and the probability W 0 (t) to be in this state at the time t. Correspondingly, we define the set M 0 = {|k 0 } consisting of the initial state alone. In order to describe how the probability spreads in time we consider all states |k i directly coupled to the initial one via the two-body interaction V. We call this set M 1 = {|k i for which | k 0 |V|k i | = 0}. From this set we define W 1 (t) as the probability to be in M 1 at the time t. As time increases, the probability flows onto another set M 2 (with probability W 2 ) consisting of those states, coupled with the states of M 1 by non-zero Hamiltonian matrix elements and so on (see figure 1). Neglecting the backward flow the infinite set of rate equations reads, where we have introduced the parameter Γ describing the exponential decay rate of the initial probability W 0 (t) ≡ P k 0 (t). In reference [54] it was shown that the above equations can be solved exactly: with W 0 (t) = exp(−Γt). The solution (12) allows one to derive the expressions for various observables.
Defining the connectivity K of this network as the number of the elements in the first set M 1 , in [54] the following approximate expression has been obtained, It should be stressed that the analogy with an infinitely large 'tree' is correct in the thermodynamic limit only. For a finite Hilbert space, as in our case, the number of the effective subsets M k is finite. As a result, as t → ∞ each of the probabilities W k (t) converges to a non-zero value. Moreover, the set of equation (11) is practically restricted by a few sets (in our simulations, by W 1 (t) and W 2 (t) only, since the number of elements in M 2 is already of the same order as the dimension of the total Hilbert space). For the case of a small number of subsets, in reference [1] the equation (11) have been modified and used to explicitly obtain an extremely detailed quench dynamics. In this case one can show that, on an intermediate time scale, which coincides with equation (13) for large K values. From the above analysis one can see that the key parameter in the quench dynamics is the parameter Γ. This parameter is closely related to the width of the LDoS. Indeed, the LDoS is defined by, where |k 0 is an eigenstate of H 0 . Thus, it is obtained by the projection of the initial state |k 0 onto the eigenstates of H. The concept of LDoS is extremely important in the analysis of the dynamical properties of many-body systems. The Fourier transform of the LDoS determines the survival probability of an initially excited many-body state, and it is effectively used in the study of fidelity in many applications. In particular, the inverse width 1/Γ of LDoS gives the characteristic time scale, which is associated with the depletion of the initial state, thus representing an early stage of thermalization [42]. Initially introduced in atomic [55] and widely used in nuclear physics [56], it also serves as an important characteristic in other physical applications. As shown in many different papers (see for instance [57] and references therein), for systems with well defined classical limits, a classical analog of the LDoS can be defined and directly computed from the Hamilton equations of motion.
For isolated systems of interacting particles, described by a Hamiltonian H = H 0 + V, the form of LDoS changes on increasing the interaction strength V [18]. If for a weak, but not negligible, interaction the LDoS is typically a Lorentzian (apart from the tails which are due to the finite width of the energy spectrum), for a strong interaction (i.e. when the interaction energy becomes of the same order of the non-interacting one) its form becomes close to a Gaussian. Correspondingly, the width Γ of the LDoS can be estimated either using the Fermi golden rule, Γ ≈ 2πV 2 ρ f where ρ f is the density of the many-body states directly connected by V, or by the square root of the variance of LDoS, σ = H 2 ij for i = j. When studying the TBRI model, this crossover was found to serve as the condition for the onset of strong quantum chaos, defined in terms of a pseudo-random structure of many-body eigenstates. For this reason, instead of Γ in our case one can use σ since the latter is much easier to estimate than the former. Thus, when comparing our data in figures 4 and 8 with the predicted exponential dependence (14), we use the following expression: As will be shown below, the simple expression (14) nicely corresponds to numerical data demonstrating the exponential increase of N pc in time. Now let us discuss another important characteristic of the relaxation for finite systems, namely, the time scale over which the exponential growth of N pc lasts. To do that, we have to estimate the saturation value N ∞ pc of N pc after the relaxation of the system to equilibrium. This value, can be obtained by the time average performed after the saturation time t s , An analytical expression, assuming non-degenerate energy levels, in terms of the eigenstates can be written as, This expression determines the total number of non-interacting many-body states inside the energy shell, excited in the process of equilibration. Following [1], let us estimate the average value of N pc after the relaxation, for the situation when the eigenstates are strongly chaotic. First of all let us notice that the second term in the rhs of equation (18) is roughly 1/D times smaller than the first one, where D is the dimension of the many-body Hilbert space. This can be seen by taking uncorrelated components C α k (1/ √ D)e iξ α,k , where ξ α,k are random numbers. Thus, one gets, As a result (taking the first of the above terms only), we arrive at the estimate, Next, we assume a Gaussian shape for (i) the LDoS, (ii) the density of states for H 0 , and (iii) the density of states of H. These are realistic assumptions for chaotic many-body systems with random two-body interactions, such as the TBRI model close to the middle of the energy spectrum. Obviously, in the tails of the spectrum the Gaussian approximation is not valid, see figure 2. Since our interest is in the energy region far from the bottom of the spectrum, our assumptions are valid, at least when we are interested in global characteristics of the dynamics, which in the first line depends on the width of the above shapes, and not on their details.
(a) For the LDoS we then assume where Γ is the width of the LDoS and E k 0 is the energy of the non-interacting state. We also assume that Γ is independent of k 0 . Note that the LDoS is normalized, F k (E)dE = 1. (b) The Gaussian shape for the non-interacting density of states ρ 0 (E) of width σ 0 is written as (c) The Gaussian density of states, characterized by a width σ, is such that where for simplicity we set the middle of the spectrum at the energy E = 0. Both densities of states are normalized to the dimension of the Fock space, Numerical data confirm the Gaussian form of ρ(E 0 ) and ρ(E), for the TBRI model in the middle of the spectrum, see figure 2. As one can see, due to the interaction the energy spectrum increases its total width.
The above assumptions imply that in the continuum, one has which is defined only for 2σ 2 > Γ 2 . We can then approximate Taking into account that σ 2 0 = σ 2 − Γ 2 , equation (26) gives The time t s , determining the onset of the saturation, can be estimated from the relation, Considering the case for which M 2N, one gets an estimate for the maximal time [1] This is one of the important results, obtained in the frame of the discussed approach. As one can see, the characteristic time t s is N times larger than the time t Γ ≈ 1/Γ describing the early decrease of the return probability. The key point is that the time t Γ has to be associated only with an initial process towards the true thermalization. In contrast, the latter emerges when the flow of probability fills all the subsets W k that create the energy shell available in the thermalization process. This result can have important implications for addressing other issues such as the scrambling of information and the quantum butterfly effect.

Thermodynamic entropy versus diagonal entropy
Now let us discuss an important relation which has been discovered by analyzing the quench dynamics leading to thermalization. This relation links two entropies, S th and S diag . Here is the thermodynamic entropy characterizing the system after its relaxation to equilibrium and N ∞ pc is the average number of basis states (eigenstates of H 0 ) in the stationary distribution. This number can be associated with the occupied 'volume' V s (E) of the Hilbert space : where δ 0 = Δ H 0 /D is the non-interacting energy spacing, Δ H 0 is the effective width of the energy spectrum of H 0 and D is the dimension of the many-body Hilbert space. As for the diagonal entropy S diag , discussed in view of its relation to the Von Neumann entropy [58], it is given by, Note that the diagonal entropy is the Shannon entropy of the set of probabilities w k 0 (E α ) = |C α k 0 | 2 obtained by the projection of the non-interacting state |k 0 of H 0 onto the eigenstates of H. With the Shannon entropy we can built the entropic localization length giving the number of eigenstates of H excited by the initial basis state [54]. Thus, the volume occupied by the initial state is V i (E) H δ, where δ is the energy spacing estimated as δ Δ H /D. Since the two volumes are equal, N ∞ pc Δ H 0 = H Δ H , we arrive at the following relation: where Δ H and Δ H 0 are the widths of the energy spectra of H and H 0 , respectively. A similar correction due to the difference between Δ H and Δ H 0 also appeared in other context [59].

Two-body random interaction model
The TBRI model is characterized by the Hamiltonian of N interacting bosons occupying M single-particle levels, in which the single-particle levels are specified by random energies s with mean spacing s − s−1 = 1 setting the energy scale. This model belongs to the class of random models since the two-body matrix elements V s 1 s 2 s 3 s 4 are random Gaussian entries with zero mean and variance v 2 . Originally suggested to describe systems with Fermi particles, recently the TBRI model has been applied to randomly interacting bosons [26][27][28][29]. As mentioned above, the Hamiltonian matrix H has a band-like form and has many zero off-diagonal elements, due to the two-body nature of the inter-particle interaction V. It is also worthwhile to mention that being constructed from random entries V s 1 s 2 s 3 s 4 , the many-body matrix elements are slightly correlated since some of them are constructed from the same two-body matrix elements. How important are these tiny correlations, have been studied analytically and numerically in reference [60]. The typical structure of the Hamiltonian matrix H is shown in figure 3. The wave packet dynamics in the non-interacting basis |k after switching on the interaction V was numerically analyzed after finding all many-body eigenstates of H and their energy eigenvalues. As one can see, the numerical study is strongly restricted by the number of particles N and the numbers M of single-particle states since the total size of the Hamiltonian matrix increases exponentially with both N and M. In our study we have mainly considered the case when M ≈ 2N, which is akin to the study of one-dimensional lattices with M sites and N Fermi-particles. As is known, in many aspects some of the observables have the same properties when both N and M increase, keeping the ratio N/M fixed.
Our interest is in the time dependence of the number of principal components N pc (t) in the wave function considered in the basis of H 0 . As we discussed in the previous sections, the system is initially prepared in a particular non-interacting state |k 0 , which determines the quench dynamics. Having all eigenstates and eigenvalues, we obtain the number N pc (t) as is shown in figure 4. For sake of completeness we also plot in the same figure the curve 2Γt corresponding to the analytical prediction obtained beyond an initial time scale on which the perturbation theory is valid. As one can see, this prediction is fully confirmed by numerical data. Now, in order to check the relation between the diagonal and thermodynamic entropies, one needs to have an accurate estimate of the diagonal entropy in the same way as we get equation (27). Let  us start with the expression given in equation (32), associated with the initial state |k 0 using the same approximations as above, Taking only the dominant contribution one gets with C 2 = exp[(1 − Γ 2 /σ 2 )/2]. Comparing equation (27) with equation (37) one finally gets, This proves that the global functional dependence of both N ∞ pc and diag on the initial energy E k 0 is the same and they turn out to be proportional one to each other.
In order to check both the consistency of our physical argument and of the relation (38) we considered for the whole set of initial energies E k 0 the two entropies, see figure 5 panel (a). As one can see they have the same dependence even if they are slightly shifted due to the variation in the size of the energy spectra. Taking into account the correction due to the different spectral widths of H and H 0 , we can see that the two entropies are pretty much the same, apart from small deviations at the edge of the spectrum where eigenstates are expected to be non-chaotic, see panel (b) of the same figure.  . On x-axis we put the dimensionless time Γt where Γ is the width of the LDoS averaged over the degenerate initial states. The average (blue curve) has been done over all initial degenerate states having the same energy E 0 close to the band center (with j 0 ≈ 4050). Here we considered N = 9 particles in = 13 momentum levels, n/g = 10, for a fixed total momentum M = 1 (the matrix size is 8122). The normalized diagonal entropy S diag + ln(Δ H /Δ H 0 ) = 2.75 ± 1 does not correspond to the equilibrium value since for this interaction strength the eigenstates involved in the dynamics are not chaotic.

Truncated LL model (t-LL)
The Hamiltonian of the Lieb-Liniger model with N bosons occupying a ring of length L, in dimensionless units can be written as, Here g is the parameter determining the strength of interaction between particles, and the single-particle energy levels (for non-interacting bosons) are given by The δ−function in equation (39)   Here we considered N = 9 particles in = 13 momentum levels, n/g = 0.5, for a fixed total momentum M = 1 (the matrix size is 8122). Horizontal green thick line represents the normalized diagonal entropy S diag + ln(Δ H /Δ H 0 ), the thickness corresponds to one standard deviation due to different initial states. In computing Δ H and Δ H 0 we excluded a number of energies close to the band edges.
single-particle momentum states. Note that the total number of single-particle energies s is M + 1 since the states with momentum ±s are degenerate. We choose N and 2M to be approximately the same, in analogy with the case considered above for the TBRI model. A rough estimate, n = N/L ∼ g, for the crossover from the MF to the TG regime in connection with quantum chaos, is discussed in [61]. The structure of the Hamiltonian matrix at some fixed total momentum value M is shown in figure 6 and as one can see it is very similar to that presented in figure 3 for the TRBI model. Our numerical study of the quench dynamics of the t-LL model demonstrates that N pc (t) oscillates in time in the MF regime (n/g 1) as shown in figure 7. In contrast, it grows exponentially fast in the TG regime (n/g 1, see figure 8), after a short time where, due to the standard perturbation theory, the time-dependence is quadratic in time. This exponential growth lasts up to some time t s after which a clear saturation of N pc emerges, together with irregular fluctuations around its mean value. In order to reduce these fluctuations which are due to different initial conditions (various values of j 0 ), we have performed the average N pc over all those initial non-interacting basis states with the same energy.
First, the numerical data clearly manifest a very good correspondence to the analytical estimates of the exponential increase of N pc (t) occurring for t t s . Second, the relaxation time t s roughly corresponds to the estimate t s ≈ N/Γ (specifically, to Γt s ≈ N with N = 9). Finally, the relation (34) between the diagonal and thermodynamic entropies holds with a very good accuracy. All these results are highly non-trivial since they are obtained for the t-LL model which is non-random. Specifically, one can show that the non-diagonal matrix elements of the total Hamiltonian are strongly correlated. For the case shown in figure 8 diagonal and off-diagonal Hamiltonian matrix elements take the values from a set with very few elements, as one can see from the plots shown in figure 9. It is important to note that the validity of our analytical expression holds also in this case where the non-interacting density of states is a set of δ−functions (and certainly not a Gaussian).

Summary
We have studied two different models, the TBRI model with random two-body interaction and the truncated LL model (t-LL) with a finite number of particles in a finite number of momentum states which is originated from the completely integrable Lieb-Liniger model. By studying the quench dynamics in the region where the many-body eigenstates can be considered as strongly chaotic, we have found that the time evolution of both models is quite similar to that recently found for the random TBRI, as well as for the deterministic XXZ model of interacting 1D spins-1/2 [19,20]. Specifically, the number of components in the wave packets evolving in the Hilbert space, increases exponentially in time with a rate given by twice the width of the LDoS related to the initial state of the non-interacting many-body Hamiltonian H 0 . This growth lasts approximately until the saturation time t s which is much larger than the inverse width of the LDoS characterizing the initial decay of the survival probability. Both the rate of the exponential increase and the characteristic time t s for the onset of equilibration, nicely correspond to our semi-analytically estimates.
By studying the process of relaxation of the system to the equilibrium, we have discovered a remarkable relation between the thermodynamic entropy S th (defined in terms of the number of principal components N pc in the wave function) of the system after the equilibration, and the diagonal entropy S diag related to an initial many-body state. This relation (34) establishes a direct link between statistical and thermodynamical properties, and seems to be generic, valid for both deterministic and random many-body systems. Recently, another relation between the diagonal entropy and the GGE entropy was found in a fully integrable 1D Ising model in a transverse magnetic field [62,63]. Thus, a further study of the relevance of the diagonal entropy to the thermodynamic properties of many-body systems seems to be very important, especially, in view of possible experimental studies of the onset of thermalization in isolated systems.