Quantum Gibbs distribution from dynamical thermalization in classical nonlinear lattices

We study numerically time evolution in classical lattices with weak or moderate nonlinearity which leads to interactions between linear modes. Our results show that in a certain strength range a moderate nonlinearity generates a dynamical thermalization process which drives the system to the quantum Gibbs distribution of probabilities, or average oscillation amplitudes. The effective dynamical temperature of the lattice varies from large positive to large negative values depending on energy of initially excited modes. This quantum Gibbs distribution is drastically different from usually expected energy equipartition over linear modes corresponding to a regime of classical thermalization. Possible experimental observations of this dynamical thermalization are discussed for cold atoms in optical lattices, nonlinear photonic lattices and optical fiber arrays.


Introduction
The problem of thermal distribution for photons led to the discovery of the Planck constant and Planck law [1]. Further development of quantum mechanics generalized the Gibbs thermal distribution [2] to the quantum case leading to the quantum Gibbs distribution in a quantum system with discrete energy levels (see e.g. [3,4]). Thus, the problem of thermalization has always fascinated scientists, starting from the famous dispute between Boltzmann and Loschmidt on time reversibility and statistical description (see e.g. [4]).
The thermalization in a given system is based on the ergodicity of motion which can be produced by noise from a heat bath or by internal dynamical chaos. The mathematical and physical foundations of dynamical chaos are now well established and are described in [5][6][7][8].
The first numerical investigations of the onset of ergodicity and dynamical thermalization in a nonlinear lattice of coupled oscillators was performed for the Fermi-Pasta-Ulam problem [9][10][11][12] with an expectation to find energy equipartition over linear oscillator modes. Surprisingly, for a typical set of parameters the equipartition was absent, even if in certain cases signs of nonperiodic behavior were visible. The absence of ergodicity stimulated a great interest in the Fermi-Pasta-Ulam problem even if later it became clear that this model is rather close to the integrable Toda lattice and, hence, it does not belong to a class of generic models (see discussions in [8,11,12]).
Another approach to the investigation of the onset of ergodicity over linear oscillator modes in nonlinear lattices had been proposed in [13] by analyzing the effects of nonlinearity on the Anderson localization [14] in systems with disorder or systems of quantum chaos. It was found that below a certain critical nonlinearity spreading over modes is suppressed or is exponentially slow while at moderate nonlinearity subdiffusive spreading continues up to times a millions time larger than a typical timescale of oscillations. This result has been confirmed and significantly extended by further investigations [15][16][17][18][19][20][21][22][23][24][25][26]; however full understanding of the problem is still lacking. Thus, the results [24] indicate that at large times the spreading continues along certain chaotic but nonergodic layers. The mathematical studies [27][28][29] demonstrate all the complexity of this problem where the pure-point spectrum of the linear system generates intricate resonances induced by nonlinearity. Interest in the problem is also supported by experiments with disordered nonlinear photonic lattices [30,31] and Bose-Einstein condensates of cold atoms placed in a disordered optical lattice [32].
Recently, it was argued that in the discrete Anderson nonlinear Schrödinger equation (DANSE) a process of dynamical thermalization takes place leading to a statistical equilibrium in a finite disordered lattice at a moderate nonlinearity [33]. It was shown numerically that the Gibbs energy distribution takes place over linear eigenmodes. This work generated a certain interest in the process of dynamical thermalization in weakly nonlinear lattices [34]. It was also pointed out that such a thermalization is necessary for the emergence of Kolmogorov turbulence in finite size systems [35].
Here we extend the studies of dynamical thermalization in disordered lattices with weak or moderate nonlinearity. We especially stress the situation when the energies of linear modes grow linearly with the index of the linear modes corresponding to a static Stark field or finite density of levels in a unit energy (frequency) interval. Such a case is typical for the Kolmogorov (or weak wave) turbulence in finite systems [36,37]. As an example of such a system we can name the nonlinear Schrödinger equation in the Sinai billiard (or any other chaotic billiard) as discussed in [35]. It is also important to note that the DANSE with a static field is also characterized by a subdiffusive spreading [38].
In this work, we extend the research line of dynamical thermalization in nonlinear disordered lattices investigating a large number of models. Surprisingly, our results show that in lattices with weak or moderate nonlinearity there is emergence of a quantum Gibbs distribution over the energies of the linear eigenmodes. In some sense, the weak nonlinearity acts as a dynamical thermostat creating a quantum Gibbs distribution. We discuss the conditions under which such a quantum Gibbs replaces the usually expected energy equipartition over the linear modes predicted by the classical thermalization theory [3][4][5][6][7][8].
The paper is constructed as follows: in section 2 we describe all nonlinear lattice models investigated in this work, in section 3 we introduce the quantum Gibbs ansatz, results for 1D M1, M2 models and 2D M3, M4 models are presented in sections 4 and 5, the results for the Klein-Gordon (KG) lattice are given in section 6, the discussion of the results is presented in section 7.

Description of nonlinear lattice models
To investigate the phenomenon of the emergence of a quantum Gibbs distribution we study several models of linear lattices with disorder and additional weak or moderate nonlinear terms. These models represent one-dimensional (1D) and two-dimensional lattices (2D) which in the absence of nonlinearity can be reduced to the Anderson model of noninteracting electrons (see e.g. [39]) on a disordered lattice in 1D and 2D, respectively. The main DANSE model [16] is described by the equation In the following, we use dimensionless units withh = V = 1, the Boltzmann constant is taken to be unity so that we have all the dimensionless variables. In total, we consider the lattice with N sites and periodic boundary conditions. For E n = 0 and long wave limit, the system is reduced to the nonlinear Schrödinger equation which is also known in the field of cold atoms as the Gross-Pitaevskii equation [40]. At β = 0 and random values of E n distributed in the interval −W/2 E n W/2, the system (1) represents the 1D Anderson model with the localization length ≈ 96/W 2 [39]. For this distribution of E n and nonzero β equation (1), named the DANSE model, was discussed and investigated in [13,[15][16][17][18] and other papers.
The Hamiltonian of the DANSE has the form with ψ n and ψ * n being the conjugated variables. The energy and the probability norm n |ψ n | 2 = 1 the exact integrals of motion. The Hamiltonian (2) can be rewritten in the basis of the linear eigenmodes ϕ nm related to ψ n = C m ϕ nm . In the eigenmode representation the Hamiltonian is with m |C m | 2 = 1, and V mm m 1 m 1 ∼ −3/2 being the transition matrix elements [13] (the dependence on is given by assuming the random matrix estimate for the eigenstates overlap). From this representation it is especially clear that spreading takes place only due to the nonlinear β coupling.
In 1D we consider the extensions of the DANSE model given by the following replacements in equation (1): Here E n have the same random distribution as in DANSE, n 0 = (N + 1)/2 marks the center of the lattice and the periodic conditions link sites N and 1. This is the M1 model with the static Stark field f which models the constant density of the states in energy as is the case in the quantum Sinai billiard [35]. We also study the M2 model which is obtained from M1 by the following replacement of the nonlinear term: In the M2 model, the nonlinear term grows with the level number that often happens for nonlinear wave interactions in wave turbulence (see e.g. [36,37,41]). We also analyze the 2D DANSE lattice studied in [19]: i ∂ψ n x n y ∂t = E n x n y ψ n x n y + β|ψ n x n y | 2 ψ n x n y + (ψ n x +1n y + ψ n x −1n y + ψ n x n y +1 + ψ n x n y −1 ).
Periodic boundary conditions are used for the N × N square lattice with −N /2 n x , n y N /2. However, here we use the extended version of this model assuming that This is the M3 model with random values of energies δ E n x n y in a given interval. In addition we study the M4 model obtained from the model M3 by the replacement This is the 2D analogue of the M2 model. Since the term (n 2 x + n 2 y ) grows linearly with index k = |n x | + |n y | we can consider the M3 model as the model for the nonlinear Schrödinger equation in the Sinai billiard (see equation (6) at F = 0 in [35]). Indeed, in a Sinai billiard the energy levels are randomly and homogeneously distributed over the energy axis, as is the case in the M3 model at f > 0, and also the nonlinear term has a similar form coupling the linear modes. The advantage of the M3 model is that it is significantly easier for numerical simulations compared to the case of Sinai billiard. The M4 model has a stronger nonlinear interactions at high wave vectors that is typical for the weak wave turbulence [36,37].
We note that the 2D M3, M4 models also can be written in the form (3) with more complex matrix elements induced by the nonlinear coupling on 2D lattice.
The above M1, M2, M3 and M4 models have two integrals of motion being energy and the wavefunction norm. The latter is generally absent in nonlinear lattices. For this reason we consider the KG lattice (KG model) described by the Hamiltonian where˜ l are taken as random in the interval [1/2, 3/2] (see e.g. [18]). This KG model was studied in [18] and it was shown that it has the same type of subdiffusive spreading as DANSE.
We keep the same notations as in [18] (see equation (6) there) but we introduce the nonlinear coefficient β (it is taken at β = 1 in [18]) and we add a static field f replacing˜ l →˜ l + f |l − l 0 | keeping the random distribution in the same interval (in [18] f = 0). We use l 0 = (N + 1)/2 and periodic conditions linking sites l = 1 and N . As shown in [18], the linear part of the Hamiltonian at f = 0, β = 1 can be reduced to the 1D Anderson model. The time evolution of the M1, M2, M3 and M4 models was integrated numerically using the symplectic integration scheme as described in [19]; the KG model was integrated by SABA 2 C method described in [18]. The time average is performed over the time interval δt in the vicinity of time t. The integration time step was fixed at δt = 0.05 for all of the models but we checked that its decrease by a factor 10 did not affect the results of numerical simulations.

Quantum Gibbs ansatz
For the DANSE and M1, M2, M3 and M4 models we make a quantum Gibbs conjecture that the nonlinear terms act like some kind of dynamical thermostat which creates the quantum Gibbs distribution over quantum states with linear mode eigenenergies m . Then according to the standard relations of statistical mechanics [3,4] we find the probabilities ρ m = |C m | 2 and the statistical sum Z of the system Here, T is a certain temperature of our isolated system which depends on the initial energy given to the system. As usual for any quantum system with energy levels m we have the total probability m ρ m = 1 and total energy E = ρ m m (here we neglect a small nonlinear term correction to energy). The norm conservation can also be taken into account using the standard approach of statistical mechanics with the chemical potential and conservation of number of particles (or norm) [3,4] that is equivalent to the normalization used in (10). We note that possibilities of thermalization has been discussed in nonlinear chains starting from the FPU problem [9][10][11][12] and continuing even for nonlinear breathers [42][43][44]. However, here we consider the case of weak or moderate nonlinearity when the nonlinear terms are relatively small comparing to linear quadratic terms. In this case the classical system is expected to reach energy equipartition over linear modes [3,4,45,46].
The entropy of the system can be expressed via the average probability ρ m on level m via the usual formula where the overline means time averaging. The entropy S, energy E and temperature T are related to each other via the standard thermodynamics expressions [3] This value of the entropy yields the maximal possible equipartition for a given initial energy.
In an implicit way, the value of energy E determines the temperature T of the system and its entropy, or by varying temperature T in the range (−∞, +∞) we obtain the variation E(T ), S(T ) and implicitly the curve S(E). The advantage of the variables E, S is based on the fact that they both are extensive variables [3,4] and thus they are self-averaging and hence in numerical simulations they have significantly smaller fluctuations compared, e.g., to temperature T . It is important to note that the above quantum Gibbs relations can also be obtained from the condition that the entropy S takes the maximal value at variation of probabilities ρ m . In fact, the quantum Gibbs ansatz was introduced in [33] for the DANSE and it was shown that it works at moderate nonlinearity β and not very strong disorder W (see also discussions in [44]). However, in [33] the striking paradox of the quantum Gibbs ansatz was not pointed out directly. Indeed, the nonlinear classical lattice is expected to have energy equipartition over linear modes that is in drastic contrast to the quantum Gibbs distribution described above.
The examples of dependence S(E) and T (E) produced by the quantum Gibbs ansatz for the M1 and KG models are shown in figure 1. We use one disorder realization with eigenvalues m for M1 and m = ω 2 m /2 for KG (more details on the KG model are given in section 6). To compare the numerical data obtained from time evolution with the Gibbs ansatz we use the exact eigenenergies m obtained from exact matrix diagonalization of the linear problem at a given disorder realization. Examples of dependence of m on index m are shown in figure 2 for the M1 model. We also can use the average dependence m ≈f m/2 (1 m N ) withf = 2( max − min )/N , which in an approximate manner takes into account the disorder fluctuations with the maximal max and minimal values min of linear eigenenergies. This approach of an effective average densityf gives a good description of numerical data (see figure 1). It gives a slight shift of the maximum of S(E) curve which is more sensitive to disorder and is not of principal importance. We return to the discussion of the KG model in section 6. In contrast to the quantum Gibbs distribution the classical thermodynamics implies the energy equipartition over all the modes [3,4] that gives where E min is a certain minimal energy of the system and C 0 is a numerical constant. The results of figure 1 show the drastic difference between the predictions of quantum and classical thermodynamics.
The dependence S(E) has one maximum and according to the standard thermodynamics relations (12) the system has a negative temperature T < 0 at the right branch of S(E) curve. It is known that such situations can appear in quantum systems with energies located in a finite bandwidth [3,4]. We note that recent experiments with cold atoms in optical lattices [47] allowed us to realize finite quantum systems at negative temperatures.
We should stress that the quantum Gibbs distribution we find has close similarities to the thermal quantum distribution in real quantum systems; however, it appears as a result of dynamical thermalization in weakly nonlinear classical coupled oscillators without any second quantization. This Gibbs distribution results from dynamical thermalization and entropy maximization over linear modes without real quantum Planck constant entering the game. In this respect our physical interpretation is very different from the one developed in [48] where the authors discussed the appearance of the real quantum Planck constant in the thermal equilibrium of classical nonlinear lattices. In our consideration we have an effective Planck constant which may be effectively introduced in a system of weakly coupled nonlinear oscillators (e.g. as a typical frequency difference between frequencies for DANSE or KG models). Below we present the numerical results on the detailed verification of the quantum Gibbs ansatz for various lattice models.

Results for one-dimensional lattice models
The dependences S(E) for the 1D lattice M1, M2 models are shown in figures 3 and 4.
For M1 we see that at β = 2 the quantum Gibbs works well at f = 0.5 and W = 2. At fixed W , an increase of f leads to the appearance of a significant number of nonthermalized modes at f = 2. Indeed, at large f the average distance between the linear modes is growing ω ≈ f and a nonlinear frequency broadening δω becomes too small so that the nonlinear coupling between linear modes starts to be perturbative and the integrability sets in for larger and larger number of initially excited modes. A similar situation for three oscillators with a nonlinear coupling had been discussed in [49]. An increase of disorder from W = 2 up to 4 reduces the localization length and the number of coupling terms between linear mode drops. This leads to a larger number of nonthermalized modes.
For the M2 model in figure 4 we take a relatively small value of nonlinearity β = 0.2. Thus, at m < m c ≈ N /3 we have a local effective β eff ≈ β|m − n 0 | < 1, thus the dynamicals remains mainly integrable and the dynamical thermalization is absent for low-energy modes. However, at m > m c ≈ N /3 we have the onset of dynamical thermalization and the Gibbs law works for high-energy modes. With the increase of time we see the increase of number of thermalized modes at m > m c .

Results for two-dimensional lattice models
The results for the 2D lattice M3 and M4 models are presented in figures 5 and 6. We note that the M3 model can also be viewed as a model for a nonlinear interaction of laser modes in optical fibers in which z-propagation along the fiber is analogous to the time evolution in our model. At present the nonlinear dynamics of modes in laser fiber arrays attracts significant interest from the optics community (see e.g. [50][51][52]). We take here a relative large value of a static field f = 1 having in mind to model the evolution of the nonlinear Schrödinger equation in a Sinai billiard. Of course, the M3 model is only an approximation of this physical system. The obtained results resemble those found for 1D models. At weak nonlinearity we have a large fraction of nonthermalized modes while for β 2 (in M3) and β 1 (in M4) we find that practically all initial conditions with linear eigenmodes follow the S(E) curve given by the quantum Gibbs ansatz.
A more detailed comparison between the numerically obtained probabilities ρ m and the probabilities given by the Gibbs ansatz is shown in figure 6. For a given disorder realization we start from linear eigenmode m and numerically determine the time-averaged probability ρ m in each of N linear modes. In addition ρ m are averaged over ten disorder realizations. The numerical results at β = 1, 4, W = 2 and N = 32 are compared with the theoretical probabilities of the Gibbs ansatz obtained for the same disorder realizations (figure 6). We see that for β = 1 (left panel) there is a significant probability to find nonthermalized modes well visible as a high density near the diagonal. However, for β = 4 we have good agreement with the probability distribution of the Gibbs ansatz.

Results for the one-dimensional Klein-Gordon lattice model
The above lattice models have two exact integrals of motion being the energy of the system and the total probability. These models are obtained from a nonlinear Schrödinger equation and hence the appearance of the quantum Gibbs distribution can be viewed as somewhat natural result with the dynamical thermalization over quantum linear modes produced by a moderate nonlinearity. Due to that it is interesting to study the case of the KG model which has only one energy integral.
To understand the properties of the KG model we note that the eigenmodes of its linear part are described by the same linear equations as the 1D Anderson model (see the correspondence description in [18]). To explore this correspondence in a deeper way we determine the eigenmodes of displacements u lm with eigenfrequencies ω 2 m . The time evolution of the nonlinear KG equation (9) is solved numerically up to times t = 10 8 for different disorder realizations. At the initial time t = 0 we start with u l (t = 0) = u lm and p l (t = 0) = 0 ( l u 2 lm = 1). During the time evolution we compute the expansion coefficients C m (t) = l u l (t)u lm . From them we determine the time averaged norm κ = m |C m | 2 where the averaging is performed over a time interval δt around time t. The dependence of κ on time for various initial eigenmodes m is  shown in figure 7. We see that even at very large times κ remains approximately constant with variations remaining on a level of 1-2%. On average we have κ ≈ 1/2 since half of the energy is concentrated in the kinetic part which is taken as zero at t = 0. Since κ remains an approximate integral of motion we define the probabilities ρ m = |C m | 2 /κ so that their sum is normalized to unity at a given moment of time m ρ m = 1. With such a definition of ρ m we compute the entropy S of the KG model via the usual relation (11). Of course, this normalization does not The dependence on effective disorder strength W is shown in figure 9 for the standard parameters of KG model β = 1 and f = 0. We see that at disorder W = 1 there is a significant fraction of nonthermalized states. We attribute this to the fact that, in the 1D Anderson model at such W , the localization length ≈ 96 becomes much larger than the system size and the linear modes cross the system practically in a ballistic way leading to a different onset of chaos. For W = 2, 4 we find good agreement with the Gibbs ansatz.
The dependence of S(E) on f for a fixed β = 1, W = 2 is shown in figure 10. As for the DANSE-type models discussed above we find that at large f the fraction of nonthermalized modes becomes significant. The physical reasons are the same: the average spacing between linear modes becomes larger than the nonlinear coupling and the system starts to approach an integrable regime.
The dependence of S(E) on nonlinear parameter β at fixed f = 0.125, W = 2 is shown in figure 11. At β = 0.5 we still have nonthermalized modes in the Kolmogorov-Arnold-Moser integrability regime. The numerical data are in good agreement with the Gibbs ansatz at β = 1 while at β = 2 the deviations become slightly visible. The deviations become larger for β = 4 (data not shown). This happens since at large β the nonlinear part of Hamiltonian is not weak or moderate and, hence, it leads to the appearance of significantly nonlinear effects including breathers and other phenomena. It is possible that the classical energy equipartition over linear  modes will appear at such larger nonlinearities. Thus, we find that the quantum Gibbs ansatz is valid within a certain finite range of nonlinearity β min < β < β max . At moderate nonlinearities β = 1 we find that not only the curve S(E) is well described by the Gibbs ansatz but also the probabilities ρ m (m ). This fact is illustrated in figure 12 where the probability distributions ρ m (m ), shown in color, are in good agreement with the quantum distribution of probabilities given by the theoretical Gibbs distribution. At small β = 0.5 we have nonthermalized states with a higher density at the diagonal similar to figure 6.

Discussion
In this work, we studied numerically the time evolution in various types of classical lattices with moderate nonlinearities. We show that at moderate values of nonlinear parameter β min < β < β max and at large timescales the nonlinear interactions between linear lattice modes create a steady-state quantum probability distribution over energies of linear modes. This steady-state probability distribution is well described by the quantum Gibbs ansatz (10) being drastically different from the classical steady state energy equipartition over the linear modes expected from the classical thermalization picture. In a certain sense, the nonlinear term generates a dynamical thermalization in the whole system with the emergence of the quantum Gibbs distribution.
The appearance of such quantum statistics takes place not only in the lattices with a discrete Schrödinger equation (DANSE type), where energy and norm are both conserved, but also in other types of lattices which have only one exact integral of energy. We argue that in the latter case there is approximate conservation of the norm that again makes such nonlinear lattices similar to the DANSE type case. The emergence of the quantum Gibbs ansatz in nonlinear lattices with only one energy integral of motion allows us to make a conjecture that the quantum Gibbs ansatz is a generic phenomenon typical for many-mode lattices with weak or moderate nonlinearities. Indeed, a system of linear oscillators is effectively equivalent to a certain effective Schrödinger equation and thus a nonlinear interaction of modes can drive a generic lattice to the quantum Gibbs distribution via dynamical thermalization. We think that the further analysis of dynamical thermalization in nonlinear classical lattices is of fundamental importance for deeper understanding of the onset of ergodicity and thermalization in such systems.
We hope that the phenomenon of dynamical thermalization described here can be tested in experiments with cold atoms in optical lattices (e.g. like in [32,47]), nonlinear photonic lattices (e.g. like in [30,31]) or optical fiber arrays [50][51][52] which seems to us to be especially promising.