Fragile Many Body Ergodicity
Thudiyangal Mithun, Carlo Danieli, M. V. Fistul, B. L. Altshuler, Sergej Flach
FFragile Many Body Ergodicity
Thudiyangal Mithun,
1, 2
Carlo Danieli,
3, 2
M. V. Fistul,
2, 4, 5
B. L. Altshuler,
6, 2 and Sergej Flach
2, 7 Department of Mathematics and Statistics, University of Massachusetts, Amherst Massachusetts 01003-4515, USA Center for Theoretical Physics of Complex Systems,Institute for Basic Science, Daejeon 34051, Korea Max Planck Institute for the Physics of Complex Systems, N¨othnitzer Str. 38, 01187 Dresden, Germany Theoretische Physik III, Ruhr-Universit¨at Bochum, Bochum 44801, Germany Russian Quantum Center, National University of Science and Technology ”MISIS”, 119049 Moscow, Russia Physics Department, Columbia University, New York 10027, USA Basic Science Program(IBS School), Korea University of Science and Technology(UST), Daejeon 34113, Korea
Weakly nonintegrable many-body systems can restore ergodicity in distinctive ways depending onthe range of the interaction network in action space. Action resonances seed chaotic dynamics intothe networks. Long range networks provide well connected resonances with ergodization controlledby the individual resonance chaos time scales. Short range networks instead yield a dramatic slowingdown of ergodization in action space, and lead to rare resonance diffusion. We use Josephson junctionchains as a paradigmatic study case. We exploit finite time average distributions to characterize thethermalizing dynamics of actions. We identify a novel action resonance diffusion regime responsiblefor the slowing down. We extract the diffusion coefficient of that slow process and measure itsdependence on the proximity to the integrable limit. Independent measures of correlation functionsconfirm our findings. The observed fragile diffusion is relying on weakly chaotic dynamics in spatiallyisolated action resonances. It can be suppressed, and ergodization delayed, by adding weak actionnoise, as a proof of concept.
The conventional perception of evolving dynamical sys-tems with a macroscopic number of degrees of freedom(DoF) is them being in a state of thermal equilibrium,i.e. ergodic. This assumes all allowed microstates hav-ing the same probability. It goes along with trajectoriesvisiting the vicinity of all points of the available phasespace (i.e. the phase space subject to constraints due tointegrals of motion such as e.g. the energy), and infi-nite time averages equaling available phase space aver-ages [1]. Statistical physics approaches were paved byGibbs and Boltzmann and provide a straight connectionbetween microcanonical dynamics and the emergence ofcanonical distributions [1]. The more interesting it is tostudy cases when this connection is not evident, erod-ing, or even missing. This can happen (i) for dynamicsin the proximity of an integrable limit, (ii) for dynam-ics in the proximity of nonergodic sets of measure zero(such as periodic orbits), and (iii) for dynamics drivenout of ergodicity due to additional constraints (e.g. con-densation). Fermi-Pasta-Ulam-Tsingou problems [2–6]can be associated with (ii), and non-Gibbs states for in-teracting Bose lattice gases [7–11] with (iii). As for (i),the celebrated Kolmogorov-Arnold-Moser (KAM) theo-rem is available [12], but applies to systems with finitenumbers of DoF and dictates weakly non-integrable dy-namics to be non-ergodic on a finite measure set of in-variant tori, while being ergodic on the complementaryone (Arnold diffusion [13]). The KAM borders are as-sumed to quickly diminish with increasing DoF numbers[12]. What lies beyond those borders for macroscopic sys-tems? The expectation that Gibbs and Boltzmann takeover, was shattered by recent results on Many-Body Lo-calization (MBL) [14, 15] which show that certain quan- tum many-body systems can resist thermalization at fi-nite distance from integrable limits. With most analyti-cal results being non-rigorous, and computations notori-ously heavy due to exploding Hilbert space dimensions,the weakly touched field of ergodization and thermaliza-tion of corresponding classical many body systems is inthe focus of this work.Networks of weakly coupled superconducting grainsare one of the few paradigmatic examples of systemswhere the above scenaria have been considered [16–18].Also, related networks of interacting anharmonic oscil-lators were used to argue for and show the existence oftwo different classes of nonintegrable perturbations of anintegrable Hamiltonian H ( { J k } ), with a countable setof actions J k ( k being an integer) [19]. Nonintegrableperturbations H ( { J k , Θ k } ) typically span long range orshort range networks (LRN or SRN) in the action-anglespace (here Θ k are the canonically conjugated angles). Areference action J k in that network is coupled to R k × L k -tuples of other action-angle pairs. L k is typically a singledigit integer. R k however can be intensive (SRN) or ex-tensive (LRN) [19].Consider a typical LRNs, and translationally invari-ant two-body interactions L k = 3, R k ∼ N with N being the volume (system size) [19]. Chaotic dynamicscan develop in a given L -tuple on a time scale T Λ = 1 / Λwhere Λ is the typical (largest) Lyapunov exponent inthe system. Chaotic dynamics develops due to nonlin-ear resonances which take place when ratios of networkmatrix elements to certain frequency differences are large[20]. In the proximity to an integrable limit the networkmatrix elements scale with 1 /N [19]. Therefore the prob-ability for an L -tuple to be resonant will be π r /N (cid:28) a r X i v : . [ n li n . C D ] J un where π r (cid:28) R k ∼ N L -tuples in which one given action is involved. There-fore the probability Π k that the reference action is inresonance with at least one of its L -tuples (i.e. satisfy-ing the resonance condition) turns exponentially close tounity: Π k ≈ − (1 − π r /N ) R k ≈ − e − R k π r /N for sucha macroscopic LRN. It follows that LRNs thermalize ho-mogeneously in action space, see e.g. [19].On the contrary, SRNs show anomalously slow er-godization and thermalization dynamics in proximity toan integrable limit [18, 19]. Since R k is now intensiveand N -independent, the resonance probability Π k ∼ π r is small, resonances are rare, and thermalization is de-layed until resonances were able to migrate through thewhole system. Thermalization is expected to be a highlyinhomogeneous process in action space.In this work we quantitatively describe the dynamicsof thermalization by making use of finite time average(FTA) distributions. We (a) observe a novel regime ofaction diffusion, and extract the diffusion coefficient asa function of the proximity to the integrable limit, (b)show the connection between the dynamics of FTA dis-tributions and auto-correlation functions and predict andobserve algebraic decay of correlations in time, and (c)finally predict the diffusion delay through action noise de-stroying resonances and provide computational evidenceof the delay.We consider the Hamiltonian H ( q, p ) = N (cid:88) n =1 (cid:20) p n E J (1 − cos( q n +1 − q n )) (cid:21) , (1) describing the dynamics of a chain of N superconductingislands with nearest neighbor Josephson coupling in itsclassical limit. This model is equivalent to a 1D XY chainor simply a coupled rotor chain with rotor momenta p n and angles q n . E J controls the strength of Josephsoncoupling and will be compared to the energy density h = H/N . The equations of motion of Eq. (1) read˙ q n = p n , ˙ p n = E J (cid:2) sin( q n +1 − q n )+sin( q n − − q n ) (cid:3) . (2)We apply periodic boundary conditions p = p N +1 and q = q N +1 . The system has two conserved quantities:the total energy H and the total angular momentum P = (cid:80) Nn =1 p n . We will choose P = 0 without loss ofgenerality. A SRN limit is obtained for E J /h →
0, with H = (cid:80) Nn =1 p n and H = (cid:80) Nn =1 E J (1 − cos( q n +1 − q n )).The actions J n ≡ p n , and the angles Θ l ≡ q n . Notethat the opposite limit E J /h → ∞ (not further studiedin this work) yields a LRN network as discussed in theintroduction.The microcanonical dynamics of (1) explores the avail-able phase space Γ. The phase space average of an ob-servable f ( (cid:126)X ) = (cid:104) f (cid:105) ≡ Z (cid:82) f ( (cid:126)X ) d Γ , Z = (cid:82) d Γ. Here (cid:126)X is a point in Γ. The ergodicity property is tested quanti-tatively by showing that the infinite time average of anyobservable f ( (cid:126)X ) will be equal to its phase space average (cid:104) f (cid:105) . Lacking infinite times, we rather compute finite log T -6 -5 -4 -3 -2 -1 µ h=1 Figure 1. (Color online) The time dependence of thesecond moment µ ( T ) for various values of E J : E J =0 . , . , . , . , . , . , . h = 1, N = 1024 and R = 192. time averages, which depend on both the averaging time T , and the initial condition (cid:126)X : f T ( (cid:126)X ) = 1 T (cid:90) T f ( (cid:126)X ( t )) dt , (cid:126)X ( t = 0) = (cid:126)X . (3)For an ergodic system it followslim T →∞ f T ( (cid:126)X ) = (cid:104) f (cid:105) , (4)for any choice of (cid:126)X except for a subset of measure zero.Dense scanning of all initial points (cid:126)X over Γ yields the finite time average distribution ρ ( f ; T ) of the finite timeaverages f T ( (cid:126)X ). It is a function of f , parametricallydepends on T , and is characterized by its moments µ m ( T ) = (cid:90) f m ρ ( f ; T ) df , µ = 1 . (5)It follows that the first moment µ ≡ (cid:104) f (cid:105) is invariantunder variation of the averging time T . All higher mo-ments will in general depend on T . For an ergodic systemit follows lim T →∞ µ m ( T ) → , m ≥ . (6)In our studies ρ is close to a Gaussian distribution, whichallows us to focus on µ . log T -5 -4 -3 -2 -1 µ h=1, E J =0.7 T -0.5 T -1 Figure 2. (Color online) The time dependence of the sec-ond moment µ ( T ) for different system sizes N : N =2 , , , , , , , from bottom to top. Here, h = 1, E J = 0 .
7, and R = 192 . We use the actions (momenta) p n as the relevant slowobservables in the SRN proximity to the integrable limit E J → f ≡ p n . With P = 0 it follows µ = 0. Thesecond moment µ ( T ) is then simply the variance of ρ ,and further related to the momentum-momentum auto-correlation function R ( t ) = lim τ →∞ τ (cid:82) τ p l ( τ ) p l ( t + τ ) dτ as µ ( T ) = T (cid:82) T R ( t ) dt . Under the usual assumptionthat the correlation function will have an exponential de-cay at large enough times (with ways to even weaken therequirement) the conclusion is that µ ( T → ∞ ) ∼ /T .The details of the integration methods are outlinedin Ref. [18]. We numerically integrate R trajecto-ries using symplectic integrators [21], where each initialpoint was chosen by setting q n = 0, drawing p n from aMaxwell distribution, constraining P = 0, rescaling allmomenta such that the desired energy density h is ob-tained, and giving each trajectory a prethermalizationrun of t prethermal = 10 . Since all actions p n are statis-tically equivalent, we measure them all and add all datainto one pool which is used to compute ρ . Fig. 1 shows µ ( T ) for h = 1, N = 1024 and 0 . ≤ E J ≤
1. Thevariance µ ( T ) resists decay up to some characteristicergodization time scale T E , after which it turns decayingas expected, signalling restoration of ergodicity. Notethat this time scale T E was assessed in Ref. [18] and isan intensive time scale.A close inspection of the size dependence of µ ( T ) isshown in Fig. 2 for h = 1, E J = 0 . µ ( T E ≤ T ≤ T D ) loosely E J = - - - - T δ Figure 3. (Color online) The time dependence of the parame-ter δ = d ( log µ ) /d ( log T ) for the curves in Fig.2. System sizeincreases at δ = − . δ = − . δ = − follows a 1 / √ T diffusive decay, which is followed by theanticipated 1 /T decay for T D ≤ T . The new time scale T D ( N ) is evidently system-size dependent. To supportthat finding, we plot δ = d ( log µ ) /d ( log T ) versus T in Fig. 3. The curves clearly show intermediate satura-tion on a plateau with δ ≈ − .
5, and a subsequent decayat T D ( N ) down to δ = −
1. Since T D is increasing withsystem size, we conjecture that T D ( N → ∞ ) → ∞ , ex-tending the 1 / √ T decay in µ ( T ) to infinite times forinfinite size. In turn this implies that the correlationfunction R ( t ) ∼ / √ t without any exponential cutoff inthe same limit.Let us discuss possible mechanisms leading to the ob-served behavior of µ ( T ) for different system sizes. Weconsider the occurence of a chaotic resonance to takeplace when first order perturbation theory for the evo-lution of a given rotor at site n breaks down. A sim-ple calculation yields ∆ + n < E J and ∆ − n < E J with∆ ± n = | p n ( p n − p n ± ) | . The presence of such reso-nantly coupled triplets of grains along the network gen-erate chaotic dynamics and results in a Lyapunov expo-nent whose inverse yields a time scale T Λ ≤
10 on thestudied interval 0 . ≤ E J ≤ T Λ (cid:28) T E , T D . The resonance probability can be easilycomputed as π r ∼ ( E J /h ) . At variance to the LRNcases, the SRN resonances are rare and inhomogeneouslydistributed over the system at any time. The typical dis-tance between consecutive chaotic triplets l r ∼ ( h/E J ) grows with reducing E J turning the resonances moresparse and rare [22]. The assumption of partial ther-malization of actions involved in the rare resonances willstill not lead to any substantial observation of the on-set of ergodization, simply because resonances are rareand separated by non-chaotic (regular) regions. To onsetergodicity instead, chaotic resonances have to diffusivelymigrate throughout the entire system [18]. This presum-ably happens due to an incoherent detuning of the mo-menta of rotors in a neighbourhood of a given resonance.Once such a neighbouring rotor is sufficiently detuned,it could become resonant with its own neighbourhoodforming a new resonance.To confirm that we observe action diffusion, we rescale µ → µ N and T → T /N - as shown in Fig.4. For agiven value of E J we observe very good collapse of allcurves onto one master curve for T ≥ T E . The master (T/N )10 -4 -2 µ N h=1 ( T / N ) - E J = . E J = . E J = . E J = . E J = . E J = . E J = . ( T / N ) - . Figure 4. (Color online) The rescaled time dependence of µ N versus bT/N for various E J : E J = 1 . , . , . , . , . , . , . b =0 . , . , , , , , h = 1. curves in Fig.4 show the turnover from diffusive 1 / √ T to asymptotic 1 /T decay at the time T D . The diffusionprocess assumes that a diffusion coefficient D ∼ N /T D can be read off a fit of T D . We found the inverse D − from the intersection of the fit of the diffusive 1 / √ T andthe asymptotic 1 /T trends (marked with green dots inFig 4). The measured values of the diffusion coefficient D are then reported as a function of E J in Fig.5, whichappears to be reasonably close to a power-law over theanalyzed interval. To check whether the asymptotic be-havior of D may turn into an exponential behavior ratherthan power law [23–27] we plot the data in Linear-Logscale in the inset of Fig.5. From the available data weconclude that a power law is more close to the obtaineddata.Our analysis shows that the dynamics of resonancesstarts with a diffusion process between chaotic triplets,i.e. on a length scale l r ∼ √ DT E . After that the diffu-sion continues until all fluctuations stored in N/l r non-resonant patches each of the size l r , were exchanged and J -3 -2 -1 D ( E J ) J -3 -2 -1 D ( E J ) h=1 E J-5.5 e (-2.3/E J ) Figure 5. (Color online) Diffusion coefficient D vs E J for fixed h = 1 in Log-Log scale (main) and Linear-Log scale (inset). Thered dashed lines guide the eye and represent a power-law trend D ∼ E − . J (main) and an exponential trend D ∼ exp( − . /E J )(inset). reached a given location in the system. This happens for T D ∼ N /D [28].If the above scenario is correct, we can expect to de-lay the diffusion, relaxation, and ergodization process, ifwe manage to efficiently destroy resonances, before theyhad time to diffuse. To check this prediction, we take asystem with N = 512 and E j = 0 . h = 1. Every timeinterval T Λ ≈
10 we randomly pick a site n , and incre-ment or decrement its momentum p n by a given value ∆ p with equal probability [29]. On average a given site n isreached on a time T Λ N ≈ p ≈ E J /h we expectto efficiently detune and destroy a nonlinear resonance.The effect should become visible for T ≈ log T -4 -3 -3 -1 µ h=1, E J =0.7, N=512 Figure 6. (Color online) The time dependence of µ ( T ) in thepresence of a random kick process with ∆ p = 0 (magenta), 0 . . N = 512, h = 1 and E J = 0 . tributions to characterize the thermalizing dynamics ofactions. The slowing down of the thermalization dynam-ics upon approaching the integrable limit results in a de-creasing of an effective diffusion constant which is re-lated to heat conductivity. This slowing down appearsto follow a power law in the distance from the integrablelimit, rather than an exponential one. We identify anovel action resonance diffusion regime responsible forthe slowing down. The observed fragile diffusion is re-lying on weakly chaotic dynamics in spatially isolatedaction resonances. We were able successfully delay andsuppress it by adding weak action noise, as a proof of con-cept. Among a number of intriguing open questions, wemention the search for further distinct classes of nonin-tegrable action networks (neither short nor long ranged),and the impact of quantization on the fragile short rangenetwork dynamics in the vicintiy of an integrable limit.The authors are thankful to Ara Go and Hyeong JunLee for assistance on computational aspects of the work.The authors acknowledge financial support from IBS(Project Code No. IBS-R024-D1). M.V.F. acknowledgesthe partial financial support in the framework of IncreaseCompetitiveness Program of NUST ”MISiS” K2-2020-001. [1] K. Huang, Statistical mechanics (Wiley, 1987).[2] Joseph Ford, “The fermi-pasta-ulam problem: Paradox turns discovery,” Physics Reports , 271 – 310 (1992).[3] Thomas P Weissert, The Genesis of Simulation inDynamics: Pursuing the Fermi-Pasta-Ulam Problem(Springer-Verlag, New York, NY, 1997).[4] S. Flach, M. V. Ivanchenko, and O. I. Kanakov, “ q -breathers and the fermi-pasta-ulam problem,” Phys. Rev.Lett. , 064102 (2005).[5] Giovanni Gallavotti, The Fermi-Pasta-Ulam problem: astatus report, Vol. 728 (Springer, 2007).[6] C. Danieli, D. K. Campbell, and S. Flach, “Intermittentmany-body dynamics at equilibrium,” Phys. Rev. E ,060202 (2017).[7] K. Ø. Rasmussen, T. Cretegny, P. G. Kevrekidis, andNiels Grønbech-Jensen, “Statistical mechanics of a dis-crete nonlinear system,” Phys. Rev. Lett. , 3740–3743(2000).[8] S Iubini, R Franzosi, R Livi, G-L Oppo, and A Politi,“Discrete breathers and negative-temperature states,”New J. Phys. , 023032 (2013).[9] Thudiyangal Mithun, Yagmur Kati, Carlo Danieli, andSergej Flach, “Weakly nonergodic dynamics in the gross-pitaevskii lattice,” Phys. Rev. Lett. , 184101 (2018).[10] Alexander Yu. Cherny, Thomas Engl, and Sergej Flach,“Non-gibbs states on a bose-hubbard lattice,” Phys. Rev.A , 023603 (2019).[11] Stefano Iubini, Liviu Chirondojan, Gian-Luca Oppo, An-tonio Politi, and Paolo Politi, “Dynamical freezing ofrelaxation to equilibrium,” Phys. Rev. Lett. , 084102(2019).[12] Jurgen Moser, Stable and random motions in dynamicalsystems: With special emphasis on celestial mechanics,Vol. 1 (Princeton university press, 2001).[13] Florin Diacu and Philip Holmes, Celestial encounters:the origins of chaos and stability, Vol. 22 (Princeton Uni-versity Press, 1999).[14] D.M. Basko, I.L. Aleiner, and B.L. Altshuler, “Metalinsulator transition in a weakly interacting many-electronsystem with localized single-particle states,” Annals ofPhysics , 1126 – 1205 (2006).[15] Rahul Nandkishore and David A Huse, “Many-body lo-calization and thermalization in quantum statistical me-chanics,” Annu. Rev. Condens. Matter Phys. , 15–38(2015).[16] Dominique Escande, Holger Kantz, Roberto Livi, andStefano Ruffo, “Self-consistent check of the validity ofgibbs calculus using dynamical variables,” Journal of Sta-tistical Physics , 605–626 (1994).[17] Manuel Pino, Lev B. Ioffe, and Boris L. Altshuler, “Non-ergodic metallic and insulating phases of josephson junc-tion chains,” PNAS , 536–541 (2016).[18] Thudiyangal Mithun, Carlo Danieli, Yagmur Kati, andSergej Flach, “Dynamical glass and ergodization timesin classical josephson junction chains,” Phys. Rev. Lett. , 054102 (2019).[19] Carlo Danieli, Thudiyangal Mithun, Yagmur Kati,David K. Campbell, and Sergej Flach, “Dynamical glassin weakly nonintegrable klein-gordon chains,” Phys. Rev.E , 032217 (2019).[20] Boris V Chirikov, “A universal instability of many-dimensional oscillator systems,” Physics reports , 263–379 (1979).[21] C. Danieli, B. Many Manda, T. Mithun, and Ch. Skokos,“Computational efficiency of numerical integration meth-ods for the tangent dynamics of many-body hamiltonian systems in one and two spatial dimensions,” Mathematicsin Engineering , 447 (2019).[22] D. M. Basko, “Local nature and scaling of chaos in weaklynonlinear disordered chains,” Phys. Rev. E , 036202(2012).[23] C. Giardin`a, R. Livi, A. Politi, and M. Vassalli, “Finitethermal conductivity in 1d lattices,” Phys. Rev. Lett. ,2144–2147 (2000).[24] Vadim Oganesyan, Arijeet Pal, and David A. Huse, “En-ergy transport in disordered classical spin chains,” Phys.Rev. B , 115104 (2009).[25] Sergej Flach, Mikhail Ivanchenko, Nianbei Li, andBaowen Li, “Thermal conductivity of nonlinear waves indisordered chains,” The European Physical Journal B ,1007 (2011).[26] Yunyun Li, Nianbei Li, and Baowen Li, “Temperature dependence of thermal conductivities of coupled rotatorlattice and the momentum diffusion in standard map,”The European Physical Journal B , 182 (2015).[27] S Iubini, S Lepri, R Livi, and A Politi, “Coupled trans-port in rotor models,” New Journal of Physics , 083023(2016).[28] J T Edwards and D J Thouless, “Numerical studies oflocalization in disordered systems,” Journal of PhysicsC: Solid State Physics , 807–820 (1972).[29] Each time one momentum is changed by a given value∆ p , the other momenta are corrected such that the to-tal momentum P is again zero. After that all momentaare slightly rescaled to reach the original energy density hh