Transport, correlations, and chaos in a classical disordered anharmonic chain
Manoj Kumar, Anupam Kundu, Manas Kulkarni, David A. Huse, Abhishek Dhar
TTransport, correlations, and chaos in a classical disordered anharmonic chain
Manoj Kumar , , Anupam Kundu , Manas Kulkarni , David A. Huse , and Abhishek Dhar International Centre for Theoretical Sciences, Tata Institute of Fundamental Research, Bengaluru – 560089, India. Centre for Fluid and Complex Systems, Coventry University, CV1 5FB, United Kingdom. Physics Department, Princeton University, Princeton, NJ 08544, USA, Institute for Advanced Study, Princeton, NJ 08540, USA (Dated: August 13, 2020)We explore transport properties in a disordered nonlinear chain of classical harmonic oscillators,and thereby identify a regime exhibiting behavior analogous to that seen in quantum many-body-localized systems. Through extensive numerical simulations of this system connected at its ends toheat baths at different temperatures, we computed the heat current and the temperature profile inthe nonequilibrium steady state as a function of system size N , disorder strength ∆, and temperature T . The conductivity κ N , obtained for finite length ( N ) systems, saturates to a value κ ∞ > N limit, for all values of disorder strength ∆ and temperature T >
0. We show evidence thatfor any ∆ > T in the ( T / ∆) → κ ∞ ∼ e − B | ln( C ∆ /T ) | fits our data. This form has earlier been suggestedby a theory based on the dynamics of multi-oscillator chaotic islands. The finite-size effect can be κ N < κ ∞ due to boundary resistance when the bulk conductivity is high (the weak disorder case), or κ N > κ ∞ due to direct bath-to-bath coupling through bulk localized modes when the bulk is weaklyconducting (the strong disorder case). We also present results on equilibrium dynamical correlationfunctions and on the role of chaos on transport properties. Finally, we explore the differences inthe growth and propagation of chaos in the weak and strong chaos regimes by studying the classicalversion of the Out-of-Time-Ordered-Commutator. I. INTRODUCTION
In the last two decades, there has been a consider-able amount of interest in understanding the transportproperties of systems in the presence of both disorderand interactions. It is well-known that disordered sys-tems described by quadratic Hamiltonians (e.g noninter-acting electrons in a disordered potential or disorderedharmonic crystals) exhibit the phenomena of Andersonlocalization [1], whereby the single-particle eigenstatesor normal modes (NMs) of the system form spatially lo-calized states. This has a profound effect on transport— in particular for one dimensional systems all statesare localized and one finds that the system is a thermalinsulator.A question of great interest is to ask what happenswhen one introduces interactions in such a system: Doesone need a nonzero critical strength of interactions tosee an insulator-to-conductor transition? For quantumsystems, this question has been extensively studied inthe context of many-body localization (MBL) [2, 3]. Itis now generally accepted that for one-dimensional quan-tum systems with a sufficiently strong disorder, the local-ized insulating state persists up to a critical interactionstrength. One can ask the same question in the contextof classical systems and this has been addressed some re-cent works [4–8]. The work in [4–7] leads one to believethat there is no classical analogue of an MBL phase, while[8] provides evidence that such a phase might exist in anonlinear oscillator chain, for a specially-designed real-ization of spring constants. Theoretical arguments in [7]indicate that the thermal conductivity of a disorderednonlinear system goes to zero with decreasing tempera- ture T faster than any power of T . The numerical studyin [5] is consistent with this finding, however Flach et al.in [6] found evidence for a power-law dependence. A re-cent study proved sub-diffusive transport in a disorderedchain with sparse interacting regions [9].Several other numerical as well as theoretical studieshave also investigated the phenomena of Anderson lo-calization, wave-packet diffusion, and transport in non-linear disordered media [4–7, 10–20]. Numerical studieshave shown that nonlinearity gives rise to the subdiffu-sive spreading of a wave packet in an otherwise emptylattice (thus zero temperature), implying the destructionof Anderson localization [11, 13, 20]. A theoretical expla-nation of the subdiffusive spreading is based on the factthat the nonlinearity results in non-integrability of thesystem because of which the wave packet evolves chaoti-cally, and this leads to an incoherent spreading [20–27]. Apossible mechanism of chaos generation and thermaliza-tion at nonzero temperature was discussed in [26], in thecontext of the disordered discrete nonlinear Schrodingerequation. Based on this picture it was estimated that theprobability of chaos generation scales at a low temper-ature as e − B | ln( C ∆ /T ) | , where B, C are constants, ∆ isdisorder, and T denotes the temperature. It is, therefore,argued that the conductivity follows the same scaling.The main aim of this paper is to extract, throughextensive numerics, the temperature dependence of thethermal conductivity of the disordered anharmonic chainin the T → c . For our study, a r X i v : . [ c ond - m a t . s t a t - m ec h ] A ug we consider a one-dimensional chain of harmonic oscil-lators with random frequencies and purely anharmonicnearest-neighbor coupling. This model lies in the classintroduced by Fr¨ohlich, Spencer, and Wayne [21], andtherefore we referred it as the FSW model. It is thestrong disorder limit of the so-called Klein-Gordon (KG)model [6, 25, 28, 29]. At zero temperature, the modeleffectively consists of disconnected oscillators at randomfrequencies and hence a small localization length at lowtemperatures. Thus the FSW model is suitable to studythe low-temperature behavior of conductivity since weexpect that relatively small system sizes can be used toobtain the asymptotic (infinite size) conductivity.We have performed extensive nonequilibrium simula-tions for a range of temperatures T and the disorderstrength ∆, and for different system sizes N . We showthat finite-size effects in both the low and high disorderregimes can be understood as arising from boundary ef-fects, and use finite-size scaling to extract the thermalconductivity in the infinite size limit. As one of ourmain results we find that our data at the lowest tem-perature fits well to the form κ ∞ = Ae − B | ln( C ∆ /T ) | ,which has earlier been argued on the basis of the dy-namics of multi-oscillator chaotic resonances [26]. In ad-dition to the nonequilibrium simulations, we have alsoexamined the form of equilibrium dynamical correlationfunctions for the weak and strong disorder, and our cru-cial observation is that the behavior of dynamical corre-lations is truly Gaussian for a weak disorder, while forthe strong disorder, the behavior becomes non-Gaussianbut has a diffusive scaling. Finally, we have lookedat the spatio-temporal propagation of chaos in the sys-tem by computing a classical analogue of Out-of-Time-Ordered-Commutator (OTOC) for our nonlinear disor-dered model.The contents of the paper are as follows. In Sec. II, wedescribe the Hamiltonian and the reservoirs (which aremodelled as Langevin baths). We also introduce impor-tant dimensionless units which transparently shows thattemperature is equivalent to nonlinearity strength. InSec. III, we present simulation results for the nonequilib-rium steady state heat current. We analyze various as-pects such as system size scaling, disorder, and tempera-ture dependence of the thermal conductivity. In this sec-tion, we also present results for the temperature profilesin the nonequilibrium steady state. Sec. IV is devotedto the energy correlations in space and time and, in par-ticular, their dependence on the strength of the disorder.In Sec. V, a classical analogue of Out-of-Time-Ordered-Commutator (OTOC) is investigated for our nonlineardisordered system. In particular, the behavior of theheat map, butterfly velocity, and Lyapunov exponentsas one changes disorder strength are analyzed. Finally,in Sec. VI, we conclude this paper with a brief discussionon our main findings. FIG. 1. Schematic diagram of the model. U is the nonlinearinteractions between oscillators given by U ( x l − x l − ) = ( x l − x l − ) / II. DEFINITION OF THE FSW MODEL ANDNONEQUILIBRIUM DYNAMICS
We start by taking a chain of N oscillators with masses m and random spring constants k i = mω ε i , with each ε i chosen uniformly in the interval [1 − ∆ , U of strength ν (see Fig. 1). The Hamiltonianof the system is given by H = N (cid:88) i =1 (cid:20) p i m + k i x i (cid:21) + N (cid:88) i =0 ν ( x i − x i +1 ) , (1)where { x i , p i } are respectively the positions and mo-menta of the oscillators in the chain and we set x = x N +1 = 0. The limit ∆ = 0 represents the pure case and∆ = 1 is the maximum disorder strength for this disorderdistribution.The chain of oscillators is attached to two thermalreservoirs at unequal temperatures T L and T R at the leftand right ends, respectively, so that a temperature gra-dient is generated, and there is a heat current along thechain [30, 31]. The two thermal reservoirs are modeled byLangevin equations. In dimensionless units, t → ωt and x → (cid:112) ν/ ( mω ) x , the equations of motion for 1 ≤ i ≤ N are given by¨ x i = − ε i x i − [( x i − x i − ) + ( x i − x i +1 ) ] − γ i ˙ x i + η i , (2)with η i = η L δ i, + η R δ i,N and γ i = γ ( δ i, + δ i,N ).The Gaussian white noise, η L,R , satisfies the fluctuation-dissipation relation (cid:104) η L,R ( t ) η L,R ( t (cid:48) ) (cid:105) = 2 γT L,R δ ( t − t (cid:48) )with (cid:104) η L,R (cid:105) = 0. Here the dissipation constant γ ismeasured in units of mω and temperature in units of m ω / ( νk B ), where k B is the Boltzmann constant. Theonly relevant dimensionless parameters in the problemthat remain with this scaling are the disorder strength∆, the temperature T (which is equivalent to the nonlin-earity strength ν ), dissipation constant γ , and the systemsize N . III. SIMULATION RESULTS FORNONEQUILIBRIUM STEADY STATES
We compute the heat current and the temperatureprofile in the nonequilibrium steady state (NESS), when T L > T R . The (scaled) heat current along the chain from -4 -3 -2 -1 κ ∆ = 0∆ = 0.1∆ = 0.2∆ = 0.3∆ = 0.4 κ ∆ = 0.1∆ = 0.2∆ = 0.3∆ = 0.4∆ = 0.5∆ = 0.6
10 100 1000 10000 N κ ∆ = 0.3∆ = 0.4∆ = 0.5∆ = 0.6∆ = 0.7∆ = 0.8 T = 0.08T = 0.04T = 0.01 (a)(b)(c)
FIG. 2. Plots of κ N vs N for different ∆ for a fixed value of γ = 1, and at different temperatures values, as specified inthe frames (a), (b), and (c). Point symbols are the measuredvalues of κ N ( T ), whereas the solid lines are plotted using Eq.(3) for ∆ < ∆ c and Eq. (4) for ∆ > ∆ c . The values ofthe parameters are summarized in Table. I. Data have beenshown only for those values of ( T, ∆ , N ), where a steadystate profile of the local heat current (cid:104) J l (cid:105) is obtained. left to right is given by J = (cid:104) J N (cid:105) = (cid:80) Nl =2 (cid:104) f l,l − ˙ x l (cid:105) / ( N −
1) where f l,l − = ( x l − − x l ) is the force exerted bythe ( l − th particle on the l th particle. We define T = ( T L + T R ) /
2. Then for a finite system we define athermal conductivity κ N (∆ , T ) = JN/ ( T L − T R ). For adiffusive system one expects a finite value for κ ∞ (∆ , T ) =lim N →∞ κ N (∆ , T ). In all our numerical studies we set( T L − T R ) /T = 0 . T L = 1 . T and T R = 0 . T ) and explore the system properties as wevary ∆ , T and N .We perform numerical simulations by using the veloc-ity Verlet algorithm, adapted for Langevin dynamics [32].To speed up relaxation to the steady state, the initialconditions are chosen from a product form distributioncorresponding to each disconnected oscillator being inequilibrium at temperature T i that varies linearly acrossthe chain. The system is evolved up to times ranging T ∆ r ξ A κ ∞ ±
10 - - 0.200(5)0.1 793 ±
36 - - 0.0270(4)0.2 30 ±
31 - - 0.00340(7)0.3 - 19 ± . ± . ± . ± . . ± . ± ±
13 0.008(5) 0.0204(2)0.6 - 27 ±
15 0.006(3) 0.0111(2)0.08 0.3 33 . ± . . ± . . ± . ± ±
38 0.003(1) 0.053(1)0.8 - 55 ±
42 0.005(3) 0.0383(4)TABLE I. A summary of exponents r , ξ , A , and κ ∞ deter-mined from the best non-linear fits of finite-size conductivitiesas shown in Fig. 2 with the form of Eqs. (3) and (4). Thenumbers in parentheses are the error estimates on the lastsignificant figures. All these data are for γ = 1. from 2 × − × time steps of step size dt = 0 . N , ∆, and with decreasing T . We also averageover 50 disorder samples, and our error bars representsample-to-sample variations. For T (cid:46) .
01, the conduc-tivity becomes very small ( (cid:46) − ) and reaching a steadystate becomes computationally challenging because thefluctuations become more pronounced. Therefore, for lowtemperatures, to reduce the impact of such fluctuationswe perform 10 times steps for a NESS and compute κ N by taking an average over the NESS measurementsfor another 10 time steps, which has been possible for N = 32 and 64.In Fig. 2 (a)-(c), we plot κ N against N for differentvalues of ∆ at a fixed value of γ = 1, and for tempera-tures T = 0 . , . , .
08. We observe that in all cases,the conductivity seems to converge with increasing sys-tem size to a nonzero κ ∞ (∆ , T ). However, the approachto κ ∞ with increasing N is different for the small andlarge disorder, demarcated by a characteristic disorderstrength ∆ c ( T, γ ) that depends on the temperature T and the coupling γ to the reservoirs at the ends of thechain. For ∆ < ∆ c ( T, γ ), we find that κ N is an increasingfunction of N , while for ∆ > ∆ c ( T, γ ), it is a decreasingfunction. At ∆ = ∆ c ( T, γ ) we find that the conductivityis almost independent of system size. This is illustratedin Fig. 3 which shows a plot of κ N vs ∆ for different N at T = 0 .
01 and γ = 1, and where the curves for different FIG. 3. A plot of κ N vs ∆ on a log-linear scale for different N at a fixed T = 0 .
01 and γ = 1. An inset shows the behaviorof ∆ c ( T ) vs T . ∆ -1 κ N = 16N = 32N = 64N = 128 ∆ -2 -1 κ N = 16N = 32N = 64N = 128N = 256N = 512N = 1024 ∆ -2 -1 κ N = 16N = 32N = 64N = 128 ∆ -2 -1 κ N = 16N = 32N = 64N = 128
T = 0.04, γ = 0.1
T = 0.04, γ = 1
T = 0.04, γ = 2
T = 0.04, γ = 3
FIG. 4. Plots of κ N vs ∆ on a log-linear scale with varying N at a fix T = 0 .
04, but for different values of γ . The curvesfor different N intersect at ∆ c ( T = 0 . , γ ), where ∆ c ( T =0 . , γ ) (cid:39) γ = 0.1, 1, 2, 3, respectively. N intersect at ∆ c (cid:39) .
2. The variation of ∆ c ( T, γ = 1)with temperature T is shown as an inset of Fig. 3. Fig. 4shows plots of κ N vs ∆ for different N at a fix T = 0 . γ -values. Clearly, the ∆ c ( T, γ ) for afixed T increase with an increase in γ . In the following,we present all the numerical results for a fix γ = 1.We now discuss the difference in the system-size de-pendence of κ N between weak and strong disorder. Forweak disorder ∆ < ∆ c ( T, γ ), the bulk of the chain has arelatively low thermal resistivity 1 /κ , and the boundaryresistance to the reservoirs at each end is high enoughto produce an increase in the apparent resistivity 1 /κ N of finite length chains. If we assume a total boundarythermal resistance r due to the two ends, then the bulkand boundary resistances add to give total thermal re-sistance R = ( N/κ ∞ ) + r . Hence the effective finite-size
128 256 512 1024 2048 4096 8192 N κ T = 0.005T = 0.0060.0070.0080.010.020.040.080.160.320.641.282.565.1210.24 ∆ = 0
FIG. 5. A plot of κ vs N on a log-log scale for ∆ = 0. Pointsymbols are the simulated values of κ whereas the solid linesare the best nonlinear fits of Eq. (3). conductivity is given by κ N (∆ , T ) = N/R = κ ∞ κ ∞ r/N ) , for ∆ < ∆ c . (3)For system sizes N (cid:28) rκ ∞ the boundary resistance dom-inates and one has κ N ∼ N , while for larger N (cid:29) rκ ∞ the heat transport is diffusive with κ N ∼ N .For the zero disorder (∆ = 0) case also, the finite-sizeeffects are well described by Eq. (3) and Fig. 5 shows aplot of finite size conductivities as a function of systemsize N with varying temperatures. For low system sizes,the heat transport is ballistic, i.e., κ ∼ N , while withan increase of system size, the anharmonic oscillator po-tential part in the FSW model leads to a saturation ofthe conductivity κ . The point symbols in Fig. 5 are thesimulated values of κ whereas the solid lines are the bestnonlinear fits of Eq. (3). Clearly, Eq. (3) fits accuratelyto the simulated data at all temperatures. From thesefits, the saturated values of conductivity, κ ∞ , can be ex-tracted and plotted as a function of temperature T . Wefind the dependence κ ∞ ( T ) ∼ T at low temperatures [seeFig. 6(a)].Coming to the disordered case, for strong disorder∆ > ∆ c ( T, γ ) and low enough temperature, the short-distance, short-time behavior of the chain is insulating(Anderson localized), with the thermal conduction due tochaos being relatively weak. This situation can be viewedas two parallel channels of conduction: One channel islinear conduction through Anderson-localized modes ofthe linearized system. These modes couple to both reser-voirs for the finite system. Since such states decay withdistance d as e − d/ξ (∆ ,T ) , where ξ is a localization length,their contribution to the current ∼ e − N/ξ (∆ ,T ) . The sec-ond channel is the conduction of energy between locallychaotic multi-oscillator nonlinear resonances via the pro-cess of Arnold diffusion [26, 33, 34]. This leads to a smallconductivity (system-size independent), which essentiallygives κ ∞ . Hence the contribution from these two parallelprocesses suggests the following net conductivity for the T -6 -4 -2 κ Ν N = 32 ( ∆ = 0.5)
N = 64 ( ∆ = 0.5)
N = 256 ( ∆ = 0.5)
N = 1024 ( ∆ = 0.5)
N = 256 ( ∆ = 0)
N = 1024 ( ∆ = 0)
N = 4096 ( ∆ = 0)
AAAA |ln(2.72 ∆/Τ)| -6 -4 -2 κ ∞ ( ∆ ) ∆ = 0.1 A A ∆ = 0.2∆ = 0.3∆ = 0.4∆ = 0.5∆ = 0.6∆ = 0.7∆ = 0.8∆ = 0.9
N = 32, ∆ = 0.5
N = 64, ∆ = 0.5 κ ∼ Τ κ ∼ Τ κ ∞ ∼ Τ (a)(b) κ ∞ = 0.65 exp[-0.082*|ln(2.72 ∆/Τ)| ] FIG. 6. (a) The plot of κ N vs. T for ∆ = 0 . . T , and at a low T with increasing N . We find κ ∼ T in our lowest attained temperature range for N =64. For the ordered case, a dependence κ ∞ ∼ T is seen atlow T . (b) The plot of κ ∞ as a function of ∆ /T shows agood collapse. The solid line is the fit to the form κ ∞ = A exp( − B | ln( C ∆ /T ) | ). Some finite-size κ N data at ∆ = 0 . finite system: κ N (∆ , T ) = AN e − N/ξ + κ ∞ (∆ , T ) , for ∆ > ∆ c . (4)As shown in Figs. (2), the forms in Eqs. (3,4) pro-vide excellent fits (shown by solid lines) to the finite-sizesimulation results (plotted as point symbols) in the twodifferent regimes (also see an appendix B). One of thefitting parameters gives the true thermal conductivity κ ∞ (∆ , T ). The parameters r, ξ, A, and κ ∞ , obtainedfrom our best nonlinear fits for the data of Figs. (2), aretabulated in table I. In this way, we fit κ N (∆ , T ) to thescaling forms [Eqs. (3,4], and obtain κ ∞ (∆ , T ) for manydifferent sets of (∆ , T ). Temperature dependence of κ ∞ : We next studythe temperature dependence of κ ∞ , particularly at low T . In Fig. 6(a) we plot κ N ( T ) vs T for ∆ = 0 . T = 0. As mentioned earlier, for theordered case, κ ∞ ∼ T at low T , while for the disorderedcase κ ∞ ( T ) appears to decrease at low T faster than anypower of T . If we fit the behavior to a power law, κ ∞ ∼ T a , over a narrow range of T , then around T ∼ = 0 . a ∼ = 4, as was also reported in t λ ( t ) t λ ( t ) i ε ( i ) i ε ( i ) (a) (b) T = 0.01T = 0.02T = 0.06T = 0.04 T = 0.01T = 0.02T = 0.04T = 0.005 T = 0.06T = 0.005
FIG. 7. Finite-time Lyapunov exponents as a function of timefor a single disorder sample { ε i } (shown on the top panel) forthe case (a) [the left panel] where there are no resonancesand the case (b) [the right panel] where a three-oscillatorresonance is inserted in the middle of the chain by setting ε N/ ± = ε N/ . In both cases, data correspond to N = 16with ∆ = 0 . ∆ λ T λ ∆ = 0.1∆ = 0.5∆ = 1 (a) (b) T = 0.06T = 0.1T = 0.01T = 0.02T = 0.04
FIG. 8. (a) Lyapunov exponent as a function of disor-der strength for various temperatures. (b) The Lyapunovexponent as a function of temperature for different disorderstrengths. For a fixed disorder strength, our numerics indi-cates that λ ∝ √ T . Here N = 16. Ref. [6]. Going down to T = 0 .
005 and N = 64 we find arapid increase of this effective exponent to a ∼ = 8, whichindicates that a might be even larger at this T for larger N . In Ref. [26] it has been argued that the behaviorat small T / ∆ should be of κ ∞ ∼ e − B | ln( C ∆ /T ) | , withconstants B and C ; we show in Fig. 6(b) that the datafit rather well to this form. Transport mechanism : There are several possibili-ties for the precise mechanism by which transport occursat low temperatures in the FSW model. One argument isbased on the formation of localized chaotic islands (CI),which could provide an effective channel for energy trans-port. It has been argued earlier [27] that the formationof such CIs requires three consecutive oscillators withresonant frequencies ( | ε i +1 − ε i | ∼ T ) and thus occurswith probability p ∼ T / ∆ . From our numerical stud-ies with three oscillators, we found, however, that if anyneighboring pair out of the three oscillators is in reso-nance, this seems sufficient to generate chaos [35]. Thiswould imply p ∼ T / ∆. Since the CIs form with proba-bility ∼ T / ∆, they are separated on average by distance d ∼ ∆ /T . They would then act as effective thermalreservoirs between which energy is transmitted via inter-mediate localized states. However, a detailed calculationalong the lines in [9] shows that one ends up with regionsof large resistance and eventually sub-diffusive transport.An alternate mechanism suggested in [26] for the dis-ordered nonlinear Schrodinger equation is that a moreefficient mechanism of chaos generation does not requirenearby pairs to be close in frequencies. Instead, it is pos-sible for n ∼ ln(∆ /T ) oscillators to satisfy a resonancecondition and be driven to chaos by nearby sets of oscilla-tors. Based on this picture it is estimated in [26] that theprobability of chaos generation scales as e − B | ln( C ∆ /T ) | and then one can argue that the conductivity follows thesame scaling. In fact, in Fig. 7, we show two scenariosfor N = 16. One case (left panel) has all frequencies off-resonant, and the other case (right panel) has three os-cillators in resonance. However, irrespective of whetherthere is resonance or not, we notice that the system ischaotic with almost the same value of the Lyapunov ex-ponent. Therefore, as argued by [26], we also find thatchaos generation does not require nearby oscillator pairsto be in resonance. In fact, the Lyapunov exponent turnsout to be independent of details of how the random fre-quencies are chosen. In Fig. 8 we show the dependenceof Lyapunov exponent on disorder strengths and tem-peratures. Interestingly, our numerics indicates that fora fixed disorder strength, λ ∝ √ T , which is similar towhat is seen in several other very different classical sys-tems [36, 37]. In Sec. (V) we present further results onchaos propagation in this system. Temperature profiles : The signatures of boundaryresistance, the strong temperature dependence of κ ∞ ( T ),and disorder can also be seen in the NESS temperatureprofiles. Note that using Fourier’s law j = − κ ( T ) dT /dx along with knowledge of the form κ ( T ) and the bound-ary conditions T = T L , T N = T R uniquely fixes the tem-perature profile in the steady state. In Fig. 9 we plotthe temperature profile T i = (cid:104) p i (cid:105) for different temper-atures and disorder strengths. In Fig. 9(a), which isin the low-disorder regime, the boundary resistance isclearly seen for small N , and the profile slowly converges(with increasing N ) to an asymptotic form that is con-sistent with the form κ ∞ ∼ T . For somewhat strongerdisorder and not too low temperature [Fig. 9(b)] we arenear ∆ = ∆ c ( T ), so the profile converges quickly to theasymptotic form which is now consistent with κ ∞ ∼ T . in this range of T . These two asymptotic forms are shownby black dashed lines in Figs. 9(a) and 9(b). At evensmaller temperatures and high disorder, a sufficientlysmall size system is effectively in the localized regime,and we expect a step temperature profile [38, 39]. Thereis some indication from our numerics that this is indeedthe case [see Fig. (9c)]. This is a signature for the clas-sical analogue of an MBL-like regime, which, however, T i N =16N = 32N = 64N = 128N = 256N = 512N = 1024a = 1 T i N =16N = 32N = 64N = 128N = 256N = 512N = 1024a = 3.3 i/N T i N =16N = 32
T = 0.04, ∆ = 0.1
T = 0.04, ∆ = 0.5
T = 0.01, ∆ = 0.5 (c)(a)(b)
FIG. 9. (a) Temperature profiles in the steady state with T L = 0 . T R = 0 .
03, and ∆ = 0 . T L , T R as in (a) but with ∆ =0 .
5. The analytical fits (black dashed lines) to the asymptoticprofiles were obtained by solving − κ ( T ) ∂ i T ( i ) = J . With theform κ ∼ T a , it can be solved exactly for T ( i ) and plotted asa function of i/N using a = 1 in (a), and a = 3 . T = 0 . .
5, which shows signatures of a step profile, as seenin MBL systems [38, 39]. does not survive in the thermodynamic limit.
IV. SIMULATION RESULTS FOREQUILIBRIUM DYNAMICAL CORRELATIONFUNCTIONS
Equilibrium dynamical correlation functions (unequalspace and unequal time) serve as another probe of trans-port properties and we now present results on the form ofthese correlations in different parameter regimes. Let usin particular focus on the spread of energy fluctuationscharacterized by the correlation function C ( i, t ) = N − N (cid:88) l =1 (cid:104) [ (cid:15) i + l ( t ) − (cid:104) (cid:15) i + l (cid:105) ] [ (cid:15) l (0) − (cid:104) (cid:15) l (cid:105) ] (cid:105) , (5)where (cid:15) i = p i / k i x i / ν ( x i +1 − x i ) / (cid:104) ... (cid:105) de-notes an equilibrium average. Here, to generate the equi-librium ensemble of initial conditions, we took the systemwith periodic boundary conditions and attach Langevinheat baths at temperature T to every oscillator, thus en-suring a fast equilibration. Using initial states preparedin this way, the heat baths are then removed and the sys-tem is evolved with the Hamiltonian dynamics to com-pute the time evolution of C ( i, t ) as defined in Eq. 5.In Fig. (10) we show the time evolution of C ( i, t )at T = 0 .
04 for four different disorder strengths. Wefind diffusive scaling of the correlations at all disorderstrengths, but with non-Gaussian scaling functions ex-cept for the ordered case ∆ = 0 (at least for the space-time ( i, t ) scales we have reached in our numerics). Apossible explanation would be that the system has a dis-tribution of local diffusivities, which can lead to such non-Gaussian forms, yet diffusive scaling (see, for example[40, 41]). However, we expect that, in the very long-timelimit (inaccessible in our numerics), the scaling form willeventually become a Gaussian for the disordered case, assuggested by the observation of the absence of MBL inthe previous section.We note that the diffusion constant can be indepen-dently obtained using D = κ N /c v where c v is the specificheat. We find the values of D = 0 . D obtained in this way is con-sistent with the diffusion constant obtained by fitting aGaussian in Fig. 10(a). Due to the fact that the diffusionconstant turns out to be very small for the disorderedcase, therefore one needs to go to extremely long-timesto see Gaussian behavior. From our studies of the equi-librium correlations we find that there is no qualitativedifference in their forms between the weak and strongdisorder regimes. This seems consistent with the picturethat the differences that we see in the nonequilibriumstudies basically arise from boundary effects. V. NUMERICAL RESULTS ON CHAOSPROPAGATION ANDOUT-OF-TIME-ORDERED-COMMUTATOR(OTOC)
Finally, we investigate the differences in chaos prop-agation in this system in the two regimes of weak andstrong chaos. In quantum systems, this has been studiedthrough the so-called Out-of-Time-Ordered-Commutator(OTOC) and it is seen that chaos propagates linearly intime with a finite velocity in the conducting phase whilein the MBL phase, the growth is logarithmic [42–44].As the classical analogue of the OTOC, one replaces thecommutator by the Poisson bracket ( {· · · } PB ) and thisleads to an observable [45] which essentially measureshow an initial perturbation at the site i = N/ -2 0 2 400.00020.00040.0006 √ t C ( i , t ) t = 500t = 1000t = 2000t = 3000t = 5000 -120 -60 0 60 120 i C ( i , t ) -2 -1 0 1 200.00050.0010.00150.002 √ t C ( i , t ) t = 500t = 1000t = 2000t = 3000t = 5000 -40 -20 0 20 40 i C ( i , t ) -1 0 100.0020.004 √ t C ( i , t ) t = 500t = 1000t = 2000t = 3000t = 5000 -10 0 10 i C ( i , t ) -1 0 1 i/ √ t √ t C ( i , t ) t = 500t = 1000t = 2000t = 3000t = 5000 -8 -4 0 4 8 i C ( i , t ) ∆ = 0∆ = 0.2∆ = 0.4∆ = 0.5 (a)(b)(c)(d) FIG. 10. Diffusive scaling of energy spreading for disorderstrengths ∆ = 0 , . , . , . T = 0 .
04. For this temper-ature ∆ ≈ . ae − b | x | c with b ≈ . , . , . , . c ≈ , . , . , .
924 for∆ = 0 , . , . , . We start with the Hamiltonian equations of motion ofthe system˙ x i = p i , ˙ p i = f i , i = 1 , , . . . , N , where (6) f i = − ∂ H ∂x i = − k i x i − ν (cid:2) ( x i − x i − ) + ( x i − x i +1 ) (cid:3) . Let us consider an infinitesimal perturbation { δx i (0) =0 , δp i (0) = δ i,N/ } at site i = N/ t = 0 to anyspecific initial condition of positions and the momenta ofthe oscillators ( X (0) = { x i (0) } , P (0) = { p i (0) } ). Ouraim is to study how this initial localized perturbationspreads and grows through the system both in space aswell as in time. In order to do so, we look at the OTOC FIG. 11. Heat maps showing the spatio-temporal growth ofthe OTOC for a system of size N = 256 in the (a) weak dis-order regime (∆ = 0 .
1) [left panel] and (b) strong disordercase (∆ = 0 .
6) [right panel]. The map shows the strength of (cid:104) ln D ( r, t ) (cid:105) where the average is over 10000 initial configura-tions drawn from an equilibrium distribution at temperature T = 0 . D ( r, t ) defined as D ( r, t ) = { p ( N/ r ) ( t ) , x N/ (0) } = (cid:18) ∂p ( N/ r ) ( t ) ∂p N/ (0) (cid:19) . (7)From the linearized form of the equations of motion inEq. (6), the evolution of the perturbation is given by˙ δx i = δp i ˙ δp i = − k i δx i − ν (cid:2) ( x i − x i − ) ( δx i − δx i − )+ ( x i − x i +1 ) ( δx i − δx i +1 ) (cid:3) , (8)for i = 1 , , . . . , N .The quantity of interest for measuring the spreading ofa localized perturbation given in Eq. 7 can be rewrittenas D ( r, t ) = [ δp ( N/ r ) ( t )] . (9)To obtain δp ( N/ r ) ( t ), we need to solve the system ofequations in (6) and (8). We solve these ODEs using afourth-order Runge-Kutta (RK4) numerical integrationscheme with a time step dt = 0 .
005 and with periodicboundary condition ( x = x N , x N +1 = x ). The ini-tial condition of X (0) , P (0) is chosen from the equi-librium Gibbs distribution, ρ ( X (0) , P (0)) = e − β H /Z ,where Z = (cid:82) dXdP e − β H is the partition function. Forour nonlinear model, this initial condition can easily begenerated by connecting all sites to the Langevin heatbaths at the same temperature T . In this way, the sys-tem equilibrates very fast and the distribution of { x i , p i } follow the equilibrium Gibbs distribution.For a chaotic system, it is expected that the signalshould arrive at the site r at a time t r = r/c where c gives the speed of chaos propagation (we define the ar-rival time through D ( r, t r ) = 1). At long times the signalwould eventually grow exponentially with time with Lya-punov exponent λ = (cid:104) ln D ( r, t ) (cid:105) / t . It is to be noted that (cid:104) ... (cid:105) denotes the average over equilibrated initial condi-tions ( X (0) , P (0)) and disorder realizations. For a given t -0.05-0.02500.0250.05 < l n D (r , t ) > / t r = 0r = 32r = 64r = 96r = 128 r t r ∆ = 0.1 FIG. 12. A plot of (cid:104) ln D ( r, t ) (cid:105) / t with time t for ∆ = 0 . r . The inset shows a linear behavior of t r with r , for which D ( r, t r ) = 1, and a solid red line is thebest linear fit. Here T = 0 . t -0.04-0.0200.020.04 < l n D (r , t ) > / t r = 0r = 32r = 64r = 96r = 128 r t r ∆ = 0.6 FIG. 13. Similar to Fig. 12, but for ∆ = 0 . disorder realization, the quantities ( c, λ ) which character-ize chaos propagation depend on initial conditions, andwe thus study the averaged quantity (cid:104) ln D ( r, t ) (cid:105) .In Figs. (11a,b) we display heat maps showing thespace-time evolution of (cid:104) ln D ( r, t ) (cid:105) in the weak and strongdisorder regimes respectively for a chain of size N = 256.Unlike the quantum case, here we do not see (even atearly times) any signature of logarithmic growth in thestrong disorder case. We see ballistic propagation in bothcases with a notable difference in the magnitude of thespeed and the Lyapunov exponent.Fig. 12 shows a plot of (cid:104) ln D ( r, t ) (cid:105) / t with time t at different values of r for ∆ = 0 .
1. The quantity (cid:104) ln D ( r, t ) (cid:105) / t gives the Lyapunov exponent λ in the limit t → ∞ . The numerical data are averaged over 10 suchequilibrated initial conditions at a temperature T = 0 . t , the curves for different r saturates ata value λ = 0 . t r as a time atwhich the Lyapunov exponent vanishes, i.e., D ( r, t r ) = 1.The inset in Fig. 12 shows a plot of t r vs r , and asolid line is the best linear fit. From the slope of thisfit, the speed of growth of the perturbation is given by c = 1 / . (cid:39) . .
6. Here we found λ = 0 . c = 1 / . (cid:39) . c, λ ) ≈ (0 . , . . c, λ ) ≈ (0 . , . . T = 0 .
04 and N = 256. We see that as one increases dis-order strength, both the butterfly velocity and Lyapunovexponent decrease. Note that the Lyapunov exponentsare somewhat larger than the ones reported in Fig. (8)in Sec. (III). This is because of the smaller system sizestudied there ( N = 16). VI. CONCLUSIONS AND OUTLOOK
We studied the transport properties of a nonlinearchain with the weak and strong disorder, and lookedfor signatures of classical many-body localization. Fromour numerical studies of the nonequilibrium steady state,we find an interesting cross-over behavior whereby thesystem-size scaling of conductivity ( κ N ) is qualitativelydifferent above and below a characteristic disorder (∆ c ),that depends not only on temperature but also on thecoupling strength to the baths. We find that the finite-size effects in the thermal conductivity are consistentwith boundary effects. On the one hand, there is aregime of weak enough disorder where the system can beviewed as thermal resistors in series, one for the lengthof the oscillator chain and the others for the couplingsbetween the ends of the chain and the heat reservoirs. Inthis regime the boundary resistances suppress the mea-sured thermal conductance. On the other hand, at lowenough temperature the nonlinearity and thus the chaosare weak, and for strong enough disorder the system canbe approximately realized as linearized, resulting in An-derson localized modes. In this regime short chains canbe viewed as having two parallel channels for thermalconduction, one directly from reservoir to reservoir viathe localized modes of the chain and the other throughthe weakly chaotic diffusive transport within the bulk ofthe chain. In this regime the extra conductance via thelocalized modes enhances the measured thermal conduc-tance of short chains.Our finite-size scaling analysis leads to estimates of thethermodynamic limit conductivity κ ∞ and we find evi-dence that for strong disorder κ ∞ is a function of thescaled variable ∆ /T . Our data are described well bythe form κ ∞ ∼ e − B | ln( C ∆ /T ) | which is consistent withRef. [26] for the discrete nonlinear Schrodinger equation.As argued in this reference, our numerical studies alsosuggest that chaos results from many-particle resonancesrather than a few particle ones. The form of κ is quitenon-trivial and a similar form for time-scales associatedwith the spreading of perturbations in disordered nonlin-ear media was suggested earlier in [46]. We also inves-tigated the temperature profile in NESS, and for strongdisorder and low temperatures, we found hints of step- like profiles, a feature that is expected in systems withlocalization [38, 39].We do not see signatures of the weak-strong chaoscross-over in the form of equilibrium correlation func-tions which exhibit diffusive scaling in both the weakand strong disorder regimes, as expected since the cross-over is dominated by boundary effects. We find strongnon-Gaussianity which we expect would go away in thelong time limit. A study of the OTOC in the tworegimes shows that chaos propagation is always ballisticthough the butterfly speed and the Lyapunov exponentare smaller in the strong disorder regime. We observed a T / dependence of the Lyapunov exponent on tempera-ture for both strong and weak disorder.Our study naturally leads to asking similar questions(such as spread of the perturbations) in models in whichthe oscillators in space are coupled even in the absence ofnonlinearity. This could provide more insight into many-body localization transition in classical systems. Futurework also includes understanding transport and OTOCin a model where disorder has a fractal pattern [8] whichhas been proposed as the closest classical analogue tomany-body localization. Needless to mention, a rigorousunderstanding of the transport mechanism at ultra-lowtemperatures in a nonlinear disordered many-body clas-sical system remains an interesting open problem. ACKNOWLEDGMENTS
We thank F. Huveneers and W. De Roeck for manyuseful discussions and for pointing out errors in an ear-lier analysis. We also thank C. Dasgupta for useful dis-cussions. Manoj Kumar would like to acknowledge ICTSpostdoctoral fellowship and the Royal Society - SERBNewton International fellowship (NIF \ R1 \ Mario
HPC at ICTS-TIFR and a
Zeus
HPC ofCoventry University.0 i -6 × -4 -4 × -4 -2 × -4 × -4 × -4 < j i > i -6 × -5 -4 × -5 -2 × -5 × -5 × -5 × -5 < j i > i -2 × -5 -1 × -5 × -5 × -5 < j i > i × -6 × -6 × -6 × -6 × -5 < j i > t = 10 t = 10 t = 10 (a) (b)(c) t = 10 (d) FIG. 14. Energy current profiles, (cid:104) j i (cid:105) versus i , of the system ofsize N = 128 at ∆ = 0 .
5, and T = 0 .
04. We plot these profilesfor different amounts of time t as specified in the panels (a)to (d). i -4 × -5 -2 × -5 × -5 < j i > i -2 × -6 -1 × -6 × -6 × -6 × -6 < j i > i × -6 × -6 < j i > i × -7 × -7 × -7 < j i > t = 10 t = 10 t = 10 (a) (b)(c) t = 10 (d) FIG. 15. Analogous to Fig. 14, but for T = 0 . Appendix A: Nonequilibrium steady state of theheat current
For our one-dimensional system of N oscillators con-nected to two heat baths at its end, the time derivativeof energy (cid:15) l associated with l th particle or oscillator, interms of current j l,l − from l − l , is given as [31]˙ (cid:15) = − j , + j ,L , ˙ (cid:15) l = − j l +1 ,l + j l,l − for l = 2 , , · · · , N − , ˙ (cid:15) N = j N,R + j N,N − , (A1)where j l,l − = 1 /
2( ˙ x l − + ˙ x l ) f l,l − with f l,l − = − ∂U ( x l − − x l ) /∂x l = ( x l − − x l ) . j ,L and j N,R arethe instantaneous energy current from the left and rightreservoirs into the system, respectively. These are givenas j ,L = f L ˙ x = ( − γ ˙ x + η L ) ˙ x and j N,R = f R ˙ x N =( − γ ˙ x N + η R ) ˙ x N .In the steady state, if we denote the time average as (cid:104)· · · (cid:105) , the Eq. (A1) then demands the equality of current -2 -1 / κ Ν − / κ ∞ ∆ = 0.1, Τ = 0.01∆ = 0.2, Τ = 0.04 N -6 -5 -4 -3 -2 -1 κ Ν − κ ∞ ∆ = 0.3, Τ = 0.01 ∆ = 0.5, Τ = 0.04 ∆ > ∆ c (T) ∆ < ∆ c (T) (a)(b) FIG. 16. (a) The plot of 1 /κ N − /κ ∞ versus 1 /N (on a log-log scale) for ∆ < ∆ c ( T ), where ∆ c ( T = 0 . (cid:39) .
2, and∆ c ( T = 0 . (cid:39) .
4. (b) The plot of κ N − κ ∞ versus N (ona log-linear scale) for ∆ > ∆ c ( T ). The point symbols are thesimulated values for κ N (∆ , T ), which have been shown onlyfor those N , to which the steady state is obtained. The solidlines are the fits of two different forms, presented in Eqs. (3)and (4) for ∆ < ∆ c ( T ) and ∆ > ∆ c ( T ), respectively. flowing between any neighboring pair of particles, i.e, (cid:104) j ,L (cid:105) = (cid:104) j , (cid:105) = (cid:104) j , (cid:105) = · · · (cid:104) j N,N − (cid:105) = −(cid:104) j N,R (cid:105) , (A2)with (cid:104) j l,l − (cid:105) = (cid:104) ˙ x l f l,l − (cid:105) upon using (cid:104) ˙ x l − f l,l − (cid:105) = (cid:104) ˙ x l f l,l − (cid:105) . Thus, in order to reach the steady state ofthe system, we determine the energy current (cid:104) j l,l − (cid:105) be-tween all neighboring pair of particles, and examine thebehavior of (cid:104) j l,l − (cid:105) versus l for different amounts of time.A nonequilibrium steady state (NESS) is reached whenEq. A2 holds, i.e., the current profile of the system, for (cid:104) j l,l − (cid:105) or in simple notation (cid:104) j l (cid:105) versus l , is showing es-sentially a flat behavior.In Fig. 14, we show the current profiles of the sys-tem of size N = 128, for a particular disorder sampleat ∆ = 0 .
5, and at T = 0 .
04. We compute these energycurrents independently for various values of time, as men-tioned in the panels (a) to (d). Notice from Fig. 14 thescale of fluctuations in energy current, flowing betweeneach neighboring particle, which decays as the time t israised, and eventually a steady state is reached in time t of the order of 10 as seen in the panel (d). In order1to check the effect of changing the disorder sample onrelaxation, we repeated this same analysis for differentdisorder realizations drawn from the same ∆-value, andfound that a steady state is reached in typically the sameorder of relaxation time t eq . Therefore we emphasize thatthe relaxation time does not depend upon a disorder sam-ple. Instead, it depends on parameters like N , ∆, and T .With lowering T , t eq increases rapidly as shown in aFig. 15, where we plotted the energy currents for thesame values of N = 128 and ∆ = 0 .
5, but at T = 0 . t of O (10 ). Comparing Figs. 14 and 15,the relaxation time t eq increases about 10 times in low-ering the temperature T = 0 .
04 to T = 0 .
02. Thus, itis the reason that at much lower temperatures T (cid:46) . t eq (cid:39) . With this procedure of reaching asteady state, we then started our measurement to com-pute NESS averaged heat current and also averaged itover several disorder samples. Appendix B: Finite-size scaling corrections of thethermal conductivity
To look for any dominant finite-size corrections in thescaling of conductivity given in Eqs. (3,4), we replot someof the data of Fig 2 in different manners, as shown inFig. 16. In particular, we plotted the residual-like quan-tities, 1 /κ N − /κ ∞ against 1 /N for the weak disorder∆ < ∆ c ( T ), and κ N − κ ∞ against N (on a semi-log scale)for the strong disorder ∆ > ∆ c ( T ). The point symbolsdenote the simulated values, whereas solid lines are thefitting lines. In panel (a), such lines are linear fits, repre-senting 1 /κ N − /κ ∞ ∼ /N , while the lines in panel (b)are exponential fits of the form κ N − κ ∞ ∼ exp( − N/ξ ).Clearly, the fits in both panels agree very well to thesimulated points, implying that our data do not showthe presence of any finite-size scaling corrections. Hence,Eqs. (3) and (4), respectively, for ∆ < ∆ c and ∆ > ∆ c ,precisely describe the system size scaling of conductivity κ N . [1] P. W. Anderson, Physical review , 1492 (1958).[2] R. Nandkishore and D. A. Huse, Annu. Rev. Condens.Matter Phys. , 15 (2015).[3] F. Alet and N. Laflorencie, Comptes Rendus Physique , 498 (2018).[4] A. Dhar and J. L. Lebowitz, Phys. Rev. Lett. , 134301(2008).[5] V. Oganesyan, A. Pal, and D. A. Huse, Phys. Rev. B , 115104 (2009).[6] S. Flach, M. Ivanchenko, and N. Li, Pramana , 1007(2011).[7] F. Huveneers, Nonlinearity , 837 (2013).[8] A. Senanian and O. Narayan, Phys. Rev. E , 062110(2018).[9] W. De Roeck, F. Huveneers, and S. Olla, Journal ofStatistical Physics , 1 (2020).[10] R. Bourbonnais and R. Maynard, Phys. Rev. Lett. ,1397 (1990).[11] A. S. Pikovsky and D. L. Shepelyansky, Phys. Rev. Lett. , 094101 (2008).[12] G. S. Zavt, M. Wagner, and A. L¨utze, Phys. Rev. E ,4108 (1993).[13] C. Skokos, D. O. Krimer, S. Komineas, and S. Flach,Phys. Rev. E , 056211 (2009).[14] Y. S. Kivshar, S. A. Gredeskul, A. S´anchez, andL. V´azquez, Phys. Rev. Lett. , 1693 (1990).[15] P. Devillard and B. Souillard, Journal of statisticalphysics , 423 (1986).[16] B. Doucot and R. Rammal, EPL (Europhysics Letters) , 969 (1987).[17] R. Knapp, G. Papanicolaou, and B. White, Journal ofstatistical physics , 567 (1991).[18] M. J. McKenna, R. L. Stanley, and J. D. Maynard, Phys.Rev. Lett. , 1807 (1992).[19] T. Laptyeva, J. Bodyfelt, and S. Flach, EPL (Euro-physics Letters) , 60002 (2012). [20] S. Flach, D. O. Krimer, and C. Skokos, Physical ReviewLetters , 024101 (2009).[21] J. Fr¨ohlich, T. Spencer, and C. E. Wayne, Journal ofstatistical physics , 247 (1986).[22] G. Kopidakis, S. Komineas, S. Flach, and S. Aubry,Phys. Rev. Lett. , 084103 (2008).[23] S. Flach, Chemical Physics , 548 (2010).[24] S. Tietsche and A. Pikovsky, EPL (Europhysics Letters) , 10006 (2008).[25] T. Laptyeva, J. Bodyfelt, D. Krimer, C. Skokos, andS. Flach, EPL (Europhysics Letters) , 30001 (2010).[26] D. Basko, Annals of Physics , 1577 (2011).[27] D. M. Basko, Phys. Rev. E , 036202 (2012).[28] T. Laptyeva, M. Ivanchenko, and S. Flach, Journal ofPhysics A: Mathematical and Theoretical , 493001(2014).[29] J. D. Bodyfelt, T. V. Laptyeva, C. Skokos, D. O. Krimer,and S. Flach, Physical Review E , 016205 (2011).[30] S. Lepri, R. Livi, and A. Politi, Physics reports , 1(2003).[31] A. Dhar, Advances in Physics , 457 (2008).[32] M. P. Allen and D. J. Tildesley, Computer simulation ofliquids (Oxford university press, 2017).[33] B. V. Chirikov et al. , Physics reports , 263 (1979).[34] B. Chirikov and V. Vecheslavov, Journal of statisticalphysics , 243 (1993).[35] A. Dhar, Unpublished notes (2019).[36] T. Bilitewski, S. Bhattacharjee, and R. Moessner, Phys.Rev. Lett. , 250602 (2018).[37] D. Kumar, S. Bhattacharjee, and S. S. Ray, arXivpreprint arXiv:1906.00016 (2019).[38] C. Monthus, Journal of Statistical Mechanics: Theoryand Experiment , 043303 (2017).[39] W. De Roeck, A. Dhar, F. Huveneers, and M. Sch¨utz,Journal of Statistical Physics , 1143 (2017).[40] A. V. Chechkin, F. Seno, R. Metzler, and I. M. Sokolov,Physical Review X , 021002 (2017). [41] M. V. Chubynsky and G. W. Slater, Physical review let-ters , 098302 (2014).[42] I. H. Kim, A. Chandran, and D. A. Abanin, arXivpreprint arXiv:1412.3073 (2014).[43] Y. Huang, Y.-L. Zhang, and X. Chen, Annalen derPhysik , 1600318 (2017).[44] K. Slagle, Z. Bi, Y.-Z. You, and C. Xu, Physical ReviewB , 165136 (2017). [45] A. Das, S. Chakrabarty, A. Dhar, A. Kundu, D. A. Huse,R. Moessner, S. S. Ray, and S. Bhattacharjee, Physicalreview letters , 024101 (2018).[46] G. Benettin, J. Fr¨ohlich, and A. Giorgilli, Communica-tions in mathematical physics119