Estimating ergodization time of a chaotic many-particle system from a time reversal of equilibrium noise

We propose a method of estimating ergodization time of a chaotic many-particle system by monitoring equilibrium noise before and after time reversal of dynamics (Loschmidt echo). The ergodization time is defined as the characteristic time required to extract the largest Lyapunov exponent from a system's dynamics. We validate the method by numerical simulation of an array of coupled Bose-Einstein condensates in the regime describable by the discrete Gross-Pitaevskii equation. The quantity of interest for the method is a counterpart of out-of-time-order correlators (OTOCs) in the quantum regime.

We propose a method of estimating ergodization time of a chaotic many-particle system by monitoring equilibrium noise before and after time reversal of dynamics (Loschmidt echo). The ergodization time is defined as the characteristic time required to extract the largest Lyapunov exponent from a system's dynamics. We validate the method by numerical simulation of an array of coupled Bose-Einstein condensates in the regime describable by the discrete Gross-Pitaevskii equation. The quantity of interest for the method is a counterpart of out-of-time-order correlators (OTOCs) in the quantum regime.
Quantitative characterization of ergodicity in manyparticle systems is a long-standing challenge for the foundations of statistical physics, which dates back to the Poincaré recurrence theorem [1] and Zermelo's paradox [2]. It was already pointed out by Boltzmann [3,4] and since then became fairly obvious for the practitioners in the field [5,6] that the ergodization time of many-particle systems, defined as the Poincaré recurrence time, is impractically long to be observable on experimental timescales. Instead, it is common to call many-particle systems "ergodic", when they establish the Boltzmann-Gibbs equilibrium on an experimentally observable timescale. But even with such a concept in mind, it still remains a challenge to define the corresponding ergodization time and to measure this time experimentally.
In this paper, we define the ergodization time of a chaotic system as the characteristic time one needs to monitor the system in order to extract its primary chaotic parameter, namely, the largest Lyapunov exponent, which characterizes the sensitivity of a system to infinitesimal perturbations, the so-called "butterfly effect". The advantage of this definition is that it is unbiased in the sense of not being coupled to any particular system's coordinate. Our goal is to theoretically propose and numerically test a method, which can be used to experimentally determine whether the system is ergodic or not, and if it is, then to extract the system's ergodization time. The method is based on monitoring the equilibrium noise of the system. It involves the time reversal of system's dynamics -the so-called "Loschmidt echo".
In Fig. 1, we outline the method. It consists of the following steps.
(i) Measuring equilibrium noise of observable X before Equilibrium noise of an observable X before and after an imperfect time reversal of a system's dynamics at time τ is denoted, respectively, as X(τ − ∆t) (green line) and X(τ +∆t) (red line), where ∆t = |t − τ |. In order to facilitate visual comparison, "Time" on the horizontal axis represents t before the time reversal and 2τ − t after the time reversal. The difference between the direct and the reversed noise ∆X(∆t) = X(τ +∆t)−X(τ −∆t) (thick black line) fluctuates around 0, while its amplitude grows, on average, exponentially as a function ∆t with a rate equal to the largest Lyapunov exponent λmax. The exponentially growing envelope of ∆X(∆t) is represented by dashed lines. and after slightly imperfect time reversal. The noise is to be denoted as X(τ − ∆t) and X(τ + ∆t), where τ is the time of the dynamics' reversal, and ∆t = |t − τ |.
(iii) Repeating the procedure for an ensemble of randomly chosen initial conditions on an energy shell.
(iv) Calculating two kinds of ensemble averages ln |∆X(∆t)| and ln |∆X(∆t)| . For ∆t → ∞, the former average approaches λ max ∆t, while the latter one approaches Λ∆t, where λ max is the largest Lyapunov expo- nent, and Λ is a parameter to be discussed later.
(v) Extracting the ergodization time τ erg , which, as we show below, is proportional to the difference between Λ and λ max .
The method is rather general. Here, we illustrate it for one-, two-and three-dimensional arrays of coupled Bose-Einstein condensates (BECs) in the regime describable by the discrete Gross-Pitaevskii equation (DGPE): where ψ j is the complex order-parameter, describing the condensate at site j = 1 . . . N , J is the hopping parameter, and β is the nonlinear on-site interactions parameter, respectively. The summation over k = 1 . . . N nn extends over the nearest-neighbors NN(j) of site j. The DGPE generates conservative dynamics corresponding to the Hamiltonian As a measurable quantity of interest we have chosen a set of on-site occupations , where D(t) = |δR(t)| is the distance between the two phase-space trajectories: the reference trajectory R 1 (t) and the slightly perturbed one R 2 (t) = R 1 (t) + δR(t) [28].
The ratio log D(t) D(0) fluctuates in time as the reference trajectory R 1 (t) explores the energy shell. We define instantaneous local stretching rates as λ(t) = d dt log D(t) D(0) . The largest Lyapunov exponent is the average of local stretching rates over a sufficiently long time: λ max = λ(t). We denote fluctuations of the local stretching rates by δλ (t) ≡ λ (t) − λ max , and their autocorrelator by We propose to use the convergence of λ(t) as an indicator of ergodization, and define the ergodization time as In numerical simulations, λ max and ϕ(t) can be obtained from the direct calculations of R 1 (t) and R 2 (t). However, such an approach is impractical experimentally, because it requires tracking all phase-space coordinates of the system. An alternative, more practical approach was proposed in Refs. [14,29]. That approach is based on monitoring the effect of Loschmidt echo on equilibrium noise of almost any observable (see Supplementary Material).
In the present setting, the Loschmidt echo is implemented by reversing the sign of Hamiltonian (2) at time τ , and simultaneously perturbing the state vector, where δψ i is a very small random perturbation. We track the equilibrium noise of the on-site occupations {n i (t)} before and after the time reversal, and introduce the deviation between the reversed and direct dynamics of the on-site occupations, ∆n i (∆t) ≡ n i (τ +∆t)−n i (τ −∆t). As sketched in Fig. 1, the deviations ∆n i (∆t) fluctuate with an exponentiallygrowing envelope.
We introduce two ensemble averages of ∆n i (∆t): (5) and where Λ ≡ 1 t ln exp t 0 λ(t )dt . The limit (5) was verified recently in Ref. [29]. Now, we concentrate on relation (6). The reason for the difference between parameter Λ (sometimes referred to as the generalized maximum Lyapunov exponent [30][31][32][33]) and λ max is the different order of operations of taking logarithm and ensemble averaging. This difference is controlled by the In order to demonstrate this, we first note that The average on the right-hand side can be calculated analytically on the basis of the assumption that variable t 0 δλ(t )dt is Gaussian, which gives (see Supplementary Material): Using this relation together with Eq. (4), we obtain Λ − λ max = ∞ 0 ϕ(t )dt ≡ δλ 2 τ erg . Therefore, the ergodization time can be expressed as The experimental use of Eq. (9) requires determining λ max and Λ from Eqs. (5) and (6) and, in addition, the knowledge of δλ 2 . While there might be ways of extracting δλ 2 from experimental time-series, here we resort to an empirical estimate where N 2 nn is the number of nearest neighbors for a lattice site. In Fig. 3, we substantiate the estimate (10) on the basis of our direct numerical simulations. Why this approximation works so well for the DGPE on large lattices and whether it works for a more general class of systems needs further investigation. A possible explanation of Eq. (10) is that, in our simulations, the Lyapunov eigenvector corresponding to λ max is usually localized at only a handful of sites, which is consistent with other observations of wandering localization of Lyapunov eigenvectors [34][35][36][37][38][39][40][41].
The estimate (10) leads to the following approximation for the ergodization time When the ergodicity of a system is about to break down, one obvious indicator of this is an anomalously large value of the ergodization time given by Eq. (9). One may wonder, however, whether the Loschmidt echo response contains other signatures of broken ergodicity. In an ergodic regime, the distribution of ln |∆X(∆t)| should be Gaussian (see Supplementary Material), and its variance σ 2 G (∆t) ≡ ln 2 |∆X(∆t)| − G 2 (∆t) is supposed to grow linearly in time: In the opposite case of a non-ergodic regime, the averages in G(∆t) and W (∆t) converge poorly, which in turn leads to a non-Gaussian distribution for individual realizations of ln |∆X(∆t)| [42], accompanied by a deviation from the linear growth of σ 2 G (∆t) given by Eq. (12). Thus, relation (12) can be used for an experimentally feasible test of ergodization.
For illustration, we chose three model systems: a onedimensional chain with N = 100 sites, a two-dimensional square lattice with N = 10 × 10 sites and a threedimensional cubic lattice with N = 4 × 4 × 4. We used J = 1, β = 0.01. The initial conditions corresponded to the total energy E total = 100N and the number of particles N p ≡ i |ψ i | 2 = 100N , so that the particles were distributed equally among all lattice sites n i (0) ≡ |ψ i (0)| 2 = 100 with random phases. The perturbation to the state vector at the moment of time reversal was ψ i (τ + 0) = ψ i (τ − 0) + δψ i , where δψ i is a random vector in the phase space subject to the constraint In order to test the relation (9), we calculated the two averages of Loschmidt echoes G(∆t) and W (∆t) for one-, two-and three-dimensional DGPE lattices. The results are presented in Fig. 2. The values of the characteristic exponents λ max and Λ extracted in each case are listed in Table I. We also collected long enough time-series of local stretching rates λ(t), then calculated the autocorrelation function ϕ(t) and extracted δλ 2 and τ erg . Table I compares three values of the ergodization time: the one calculated on the basis of the definition (4), the one given by Eq.(9) and the one given by the approximation (11). In Eq.(9), we used the directly calculated value of δλ 2 .  (9) and (11): λmax and Λ are extracted from Fig. 2; δλ 2 is extracted either directly from a time-series of local stretching rates according to Eq. (4) or from emprical estimate (10); the three values of τerg are obtained on the basis of the definition (4), from the Loschmidt echo relation (9), and from the approximate relation (11). For two-and three-dimensional lattices, the values of the ergodization time obtained from Eqs. (4) and (9) agree very well. And at the same time, we observe clear discrepancy between Eqs. (4) and (9) for the onedimensional lattice, which indicates that the system has not ergodized on the timescale covered by the Loschmidt echo. Non-ergodized fast growing samples in W (∆t) from Eq. (6) reach a plateau significantly earlier than others: an indication of this in Fig. 2 is an early departure of W (∆t) from the linear growth regime. Overall, the ergodization time decreases with the increasing lattice dimension, being the longest in the one-dimensional case. Slow ergodization of one-dimensional chains (for Fermi-Pasta-Ulam, Klein-Gordon chains and also DGPE) has been also noticed and investigated in Refs. [43,44].
In all three cases, we further observe that the values obtained from Eq. (11) give a satisfactory approximation to Eq. (4).
We also performed the ergodicity test associated with relation (12). The results are presented in Fig. 4, where the ratio is plotted as a function of the echo time ∆t. For the quickly ergodizing two-and three-dimensional systems, the above ratio levels off rather quickly around the expected value Λ − λ max , whereas for the slow-ergodizing one-dimensional case it never reaches the expected plateau.
Relation to quantum systems: (i) Quantum-mechanical description of Loschmidt echoes is known to involve out-of-time-order correlators [14,26,45,46]. It was recently shown in Ref. [17], that, by analyzing the growth rate of OTOCs, one can impose a temperature-dependent constraint on λ max for quantum systems. In fact, however, the growth rate of OTOCs is the counterpart of Λ defined in the present work, which is, in general, larger than λ max . This general observation was made in Ref. [20], where it was illustrated by the example of kicked rotator. Kicked rotator, however, is a system with a two-dimensional phase space. One could, therefore, hope, that, as the number of degrees of freedom in a system increases, the difference between λ max and Λ approaches zero. Our findings, however, indicate that, for a lattice of a given dimension (1D, 2D and 3D), Λ − λ max remains finite for rather large systems. Yet, this difference decreases with the increase of the lattice dimension from 1D to 2D to 3D.
(ii) The difference Λ − λ max originates from the fluctuations of Loschmidt echo amplitude, which is, as shown in the present work, sensitive to ergodicity breakdown in classical systems. The counterpart of this breakdown in quantum systems is the transition from an ergodic to a many-body localized phase. In a related study [26], the authors argued that the fluctuations of a Loschmidt echo in quantum systems are sensitive to the many-body localization transition.
To summarize, we proposed a method of estimating ergodization time of a chaotic many-particle system by monitoring equilibrium noise before and after time reversal of dynamics, and validated it numerically by simulations of the discrete Gross-Pitaevskii equation. We showed that the difference between the largest Lyapunov exponent and the growth rate of the classical counterpart of OTOCs is proportional to the ergodization time of a system. We also introduced a related test for the breakdown of ergodicity.
We acknowledge discussions with S. Flach and D. Campbell. This work was supported by a grant of the Russian Science Foundation (Project No. 17-12-01587). If an experiment can track all phase-space coordinates of a system, then it can obtain the largest Lyapunov exponent by identifying the phase-space direction δR along which the growth of a perturbation is the quickest, i.e. the eigenvector corresponding to the largest local Lyapunov exponent. However, a realistic experiment is limited to an observable X. In such a case the eigenvector is unlikely to belong to the subspace of the whole phase space that contains X, but it is overwhelmingly likely to have a non-zero projection onto that subspace. This means that where α(∆t) is the angle between the eigenvector and the direction corresponding to ∆X(∆t) in the manydimensional phase space.
Here we consider the growth of the initial difference ∆X(0) introduced by an imperfect time reversal, and justify the limits ∆t → ∞ for G(∆t) in Eq. (5) (cf. Ref. [14]) and for W (∆t) in Eq. (6) We use Eq. (A1) to express G(∆t) as G(∆t) = ln |∆X(0)| + ln |cos α(∆t)| + ln e ∆t 0 λ(t )dt , (A4) where the first term is constant, the second term remains limited from above after ensemble averaging over initial conditions, and the third term is the only one growing linearly with ∆t. The second term ln |cos α(∆t)| may appear problematic for ∆t corresponding to |cos α(∆t)| = 0. However, this singularity is integrable: it vanishes after ensemble averaging. Given the definition of λ max from the main text of the article, Eq. (A4) implies Eq. (A2).
To prove the limit (A3) for W (∆t), we assume that |cos α(∆t)| is uncorrelated with e ∆t 0 λ(t )dt and hence factorize the average This assumption is, presumably, appropriate for almost any non-local observable. It is supported by the extensive numerical experience, e.g. Refs. [12][13][14][15], showing that the eigenvectors corresponding to λ max exhibit rather erratic behavior. The above factorization leads to W (∆t) = ln |∆X(0) cos α(∆t)| + ln e Here we derive Eq. (8) by a stochastic-noise method analogous to the one developed by Anderson and Weiss [47] in a different context, namely, for the calculation of exchange-narrowed magnetic resonance linewidths. We represent the left-hand side of Eq. (B1) as where and P t (Y ) is the probability distribution of Y (t). We assume that the system fluctuates near equilibrium, and, therefore, the process δλ(t) is stationary, i.e. its probability distribution p(δλ(t i )) is independent of t i . If δλ(t) is a Gaussian random variable, then Y is also a Gaussian random variable for all times, i.e. P t (Y ) is Gaussian. If p(δλ) is not Gaussian, but the variable δλ(t) has a finite memory time τ erg , then P t (Y ) still becomes Gaussian for t τ erg (consequence of the central limit theorem).