Table Of ContentGlassy dynamics in disordered oscillator chains
Alen Senanian∗ and Onuttom Narayan†
Physics Department, University of California, Santa Cruz CA 95064
(Dated: February 1, 2017)
Theescapeofenergyinjectedintoonesiteinadisorderedchainofnonlinearoscillatorsisexamined
numerically. When the disorder has a ‘fractal’ pattern, the decay of the residual energy at the
injection site can be fit to a stretched exponential with an exponent that varies continuously with
the control parameter. At low temperature, we see evidence that energy can be trapped for an
7 infinitetime at theoriginal site, i.e. classical many body localization.
1
0
2 The Fermi Pasta Ulam (FPU) model of coupled non- whereclustersofnearbyoscillatorsareresonantresulting
n linear oscillators [1] has served as a testing ground for inchaoticspotsinthedynamics[10],thefrequencies ωi
{ }
a basis ideas in statistical mechanics for more than half haveto be chosenjudiciously. We choosethe frequencies
J a century. After the initial result that energy did not in a ‘fractal’ manner: a large gap between the frequen-
1 distribute itself efficiently between the oscillator modes, ciesofsitesneareachother,andasmallgapbetweenthe
3 therewasahugebodyofwork[2],includingthedevelop- frequencies of sites far away from each other, with the
mentofrelatedintegrablemodels[3,4]. Thefactthatthe gapsize decayingas a powerofthe distance betweenthe
]
n FPU model is approximately integrable is generally be- sites. This is made quantitative later in this paper.
n lievedtoresultinverylongequilibrationtimes,oftenbe- AsthenearestneighborcouplingconstantJ islowered
- yondwhatisobservablenumerically[5]. Athighenergies, at low temperature, we find that the decay of the excess
s
i equilibrationproceedsefficientlyandsuchmetastablebe- energy ∆EJ(t) slows down substantially. The curves for
d
haviorisnotseen[6]. TheFPUmodelhasalsobeenused ∆E(t)forvariousvaluesofJ canbesuperimposedontop
.
t to study heat conduction in low dimensional systems, of each other if lnt and ∆E are scaled for each J. Thus
a
m where a heat conductivity that diverges in the thermo- ∆E(t)=A(J)∆ˆE(tβ(J)). The numerical results are con-
dynamic limit is seen for various models [7]. sistentwithastretchedexponentialformforthefunction
d- Even when the FPU model does not equilibrate ef- ∆ˆE;stretchedexponentialdynamicsareoftenseeninex-
n ficiently, the normal modes in which the the energy is periments on glassy systems [11], and there are various
o concentrated are delocalized. It is possible to construct theoretical models with traps and a range of time scales
c disordered and linearized versions of the FPU model for that obtain similar behavior [12]. (A power law form
[
which the normal modes are localized, but the normal with a cutoff can also be fit to the data.)
1 modes are all decoupled from each other. Following the The smallest J value shows an essentially flat ∆E(t),
v great progress in the field of many body localization for and so we vary the temperature at this J and measure
9 quantum statistical systems [8], it is natural to ask if the decay of the excess energy. The behavior of ∆E(t)
2
localizationofenergycanbeachievedininteractingclas- is found to be a non-monotonic function of T: it decays
0
sicalsystemstoo. Thestudyofheatconductionindisor- slowly at both high temperature and low temperature,
9
dered nonlinear oscillator chains suggests otherwise [9]: with a more rapid decay at intermediate temperatures.
0
. evenasmallamountofnonlinearityisfoundnumerically In the low temperature regime, the numerical results in-
1
to result in normal heat conductivity for large chains. dicate thatthe systemfreezesata non-zerotemperature
0
Although this implies that a localized energy packetwill T , with∆E (t )=0.Toourknowledge,thisis
7 f T<Tf →∞ 6
1 spreadout,itisstillpossiblethatafractionoftheenergy the first evidence of classical many body localization.
: packet will remain at its original location. The Hamiltonian ofthe chainwith ‘fractal’disorderis
v
Inthispaper,weconsideraone-dimensionalringoflin-
i N N
X ear oscillators, each with a different frequency, in which m
r nearest neighbors are coupled together with a nonlinear H = 2 X[x˙2l +ωl2x2l]−JXcos(xl−xl−1) (1)
a potential. The system is initialized by equilibrating at l=1 l=1
a temperature T. Thereafter, a packet of energy is de- with periodic boundary conditions. The particles in the
positedatonesite,andthesystemisevolvedusingmolec- chain all have equal masses and are tethered to their
ulardynamics. Theexcessenergyatthesiteismeasured equilibriumpositionsbyharmonicsprings. Thetethering
asafunctionoftimet.Ifthesystemthermalizes,thedif- ensures that momentum conservation is destroyed and
ference of excess energy ∆E(t) should vanish as t . there is no anomalous transport [13, 14]. The frequency
→∞
In order to avoid accidental resonances between oscil- of each of the tethering harmonic oscillators is different.
lators that are far away from each other, or rare regions When J = 0, the energy is obviously localized at each
lattice site. For J =0, the oscillators are coupled, but if
6
J is small, one might try to use perturbation theory. If
we define
∗ [email protected]
† [email protected] a (t)=x˙ (t) iω x (t) (2)
l l l l
−
2
then the dynamical equations can be expressed as handsideofEq.(3)isweakcomparedtothe iωa term.
l l
−
Thisisnotaproofthatthesetermsareunimportant,and
J
a˙ = iωa + [sin(x x )+sin(x x )] (3) therefore we turn to numerical simulations.
l l l l+1 l l−1 l
− m − − Becauseenergytransportinthissystemisatbestvery
slow, one has to be careful to bring it to thermal equi-
with the supplementary equation
libriumbeforea packetofenergyis injected; simply run-
1 ningmoleculardynamicsforalongtime is notsufficient.
x = [a∗(t) a (t)]. (4)
l 2iω l − l WeinitializedthevelocitiesfromaGaussiandistribution,
l
and equilibrated the coordinates using Monte Carlo dy-
This is now a set of coupled first order differential equa- namics (with acceptance rate 0.5). Equilibrium was
tions. The solution to zeroeth order in J is trivial. considered to be achieved whe≈n the virial theorem was
To first order in J, al is forced by terms that are of satisfied to within 10%.
the form ∼ exp[i(mωl±1 + nωl)t], where m and n are Once in equilibrium, a heat packet of magnitude E+
integers. Thetermswithm=0,n= 1shiftthenatural was injected in the system in the form of equal and op-
−
frequencyoftheoscillatorωl byanamountthatisO(J). posite momentum between two neighboring sites. The
Other terms yield O(J) corrections to al, of the form system was then evolved dynamically with the Forrest-
Ruth algorithm, using a time step of h = 0.01. If we
J
exp[i(mω +nω )t]. (5) denote E(t;k) as the energy of sites l and l+1 at time
l±1 l l
∼ mω +(n 1)ω
l±1 − l t when the heat packet was injected at sites k and k+1
at t=0, the residual energy at l is defined as
TonextorderinJ,eachsiteisinfluencedbyitsnextnear-
est neighbors, and so on. If there is a near degeneracy
1
betweenthe frequenciesoftwositesl andl+r thatarer ∆El(t)= [El(t;l) El(t;l+N/2)] (8)
2 −
stepsapart,forsmallamplitudes(i.e. lowtemperatures),
the leading correction to al is of the order of sothat,ifthe systemweretoequilibrate,∆El( )would
∞
be zero. For the measurements reported here, the mass
Jr
. (6) ofeachsitewasm=1,thesizeoftheringwasN =64or
(ω ω )(ω ω )...(ω ω )
l l+r l+1 l+r l+r−1 l+r 128,and the extra energy injected was ∆E(t=0)=5.0.
− − −
Figure 1 shows the residual energy ∆E averaged
If the ωi ’s are chosen randomly, accidental near de- h lil
{ } over l, as a function of time, for the fractal oscillator
generacies can result in a small denominator in (and a
model with various values of the coupling constant J.
breakdownof)theperturbationexpansion. Toavoidthis,
The system consisted of N = 64 sites held at tempera-
we choose the frequencies in the following manner. Let
ture T = 0.25. The brackets ... denote averages over
N = 2n. First, the frequencies at all the odd and even h i
initial conditions as well as averagesoversites where the
sitesaresetto1and3respectively. Next,thefrequencies
energywasinjectedandmeasured. Thesimulationswere
of successive pairs of sites are increased or decreased by
averaged over 30 64 runs, with 30 runs for each site.
λ,sothatthefrequenciesare 1+λ,3+λ,1 λ,3 λ,... . ×
{ − − } As J was reduced, the dynamics become steadily slower
At the next step, the frequencies of successive groups of
and slower. The curves collapse onto one another if the
4 sites are increased or decreased by λ2. This is carried
horizontalaxisisscaledandtheverticalaxisisshiftedfor
outntimes,whenallthe ω ’sarenon-degenerate. The
i
{ } eachcurvebyadifferentamount,andtheresultfitsnicely
frequencygapbetweentwositesthatareanoddmultiple
to a stretched exponential. This implies that all the
of 2k apart satisfies
curves are of the form ∆E (t) =A(J)exp[ (t/t )β(J)],
l 0
δω 2λk[1 λ λ2...]=2λk1−2λ =2λk+1 (7) with t0 = 3 × 105. Hhowever,ithe curve fo−r J = 0.25
l,l+2k ≤ − − 1 λ is essentially flat, making it very difficult to determine
− whether it fits the same stretched exponential form or if
for N if λ = (5 √17)/2. Then it is easy to see thedecaytakesinfinitelylong: ∆E(t )=0.Inorder
thatif→J <∞<2λ2,the e−xpressioninEq.(6)issmallforall to elucidate this further, we hold the→co∞upl6ing constant
r. fixed at J =0.25 and vary the temperature.
Althoughthis is necessaryforthe perturbationexpan- Figure2showstheresidualenergyasafunctionoftime
sion to be well behaved, it is not sufficient. Regardless for the same system (but with N =128) at various tem-
of the relationship between two frequencies, it is always peratures with J = 0.25. Unlike in Figure 1, the energy
possible to find values of m and n in Eq.(5) for which wasonlyinjectedatthesitesl =(0,1).ForT >0.50,the
mω′ +(n 1)ω is as small as one wishes. This is the measurements were averaged over 1200 runs, while only
−
problemofsmalldenominatorsthatmakestheKAMthe- 900 runs were realized for the lower temperature curves.
orem[15]difficult. However,thesetermsarehigherorder In addition, the timestep for T > 5.0 were reduced by a
in the amplitudes of the oscillators. One might expect factor of 10 to account for the faster dynamics. Unex-
them to be important when the amplitude of an oscil- pectedly, the decay of the residual energy is slow at low
lator happens to be large, but in that case, because the and hightemperatures,butnotatintermediatetempera-
sine function is bounded, the coupling term on the right tures. Forthehightemperaturebehavior,wearguedear-
3
101
Temperature
0.4 0.75
2.6 0.45 0.9
0.5 1.25
2.4 0.6 1.5
il 0.7 1.75
El 100
∆ 2.2
h 0.25 J 2.0 Ei E0=2.1
∆
1.25 2.25 h 2.0
1.5 3.0
1.75 1.8 Temperature
10−1 0.4
10−1 100 101 102 103 104 105 106 0.6
t 1.75
1.6 5.0
101 32.0
10−1 100 101 102 103 104 105
1.4
102 103 104 105 106
t
il 0 1 2 3
∆El 100 0.4 2.0
h
0.3
)
βJ( 0.2 1.5
0.1 E)0
0.0 τ( 1.0
J
10−1
10−1 100 101 102 103 104 105 106
t 0.5
0.0
FIG. 1. (Top) Residual energy as a function of time for the 0.0 0.5 1.0 1.5 2.0
fractal oscillator modelwith cosineinteractions, with various T
values of J at T =0.25. The curves from bottom to top are
fordecreasingJ decreasing. (Bottom)Rescaledversionofthe
FIG. 2. (Top inset) Residual energy at sites l = 0,1 as
same. A stretched exponential curve is shown with a black
a function of time for the fractal oscillator model with co-
line. (Bottom inset) Time-scaling factor β(J) vsJ.
sine interactions, at J =0.25 and various values of T. (Top)
Residual energy for the low temperature curves from the in-
set. Filled circles indicate the points of intersection tI(T)
lierthatthecouplingbetweenoscillatorsisweak,andthe
of the curves with the dashed line. The vertical range of
slow decay of the residual energy is not surprising. and
thisplotanditsinsetareidentical. (Bottom) Theparameter
at low temperatures, the linear disordered model (which τ(T)=1/(lntI(T)−lnA) as a function of T, with A chosen
is localized) has small corrections, with the same result. to yield linear dependenceat low temperature.
Because ofthe possibility of metastability in oscillator
chains [5], one has to be careful whether the slow de-
cay really indicates energy being trapped for an infinite
exponential (except for large t and T). This implies
time. Therefore, the times t (T) at which the residual
I ∆E Aexp[ (t/t )β(T)], with the best fit values of
energy drops to E0 =2.1, approximately 80% of the en- hA =i2.∼8 and t −= 7 0 103. The data is consistent with
ergy originally injected, are calculated, and found to fit β(T) T0.2. T0hus β(×T 0)appearsto be zero,i.e. one
wthiethVAog=el0F.u47lche1r0[51,6B] fo=rm0.5t7I(aTn)d=TfA=ex0p.3[B4./T(Thi−s iTnfd)i-] canno∼t show that the dy→namics freeze at T 6=0.
catesthatafin×ite residualenergyremainsatthe original Finally, in Figure 4 we compare ∆E(t) when energyis
pair of sites when T <T , i.e. ∆E (t ) E . injectedatvariouspointsinthering,showingslowdecay
f l 0
To confirm our argument thhat the→bo∞undie≥d form of for some sites and fast decay for others. When energy
the interaction potential is (at least in part) responsi- is injected into the sites (0,1), the case discussed so far,
ble for energy localization, we carried out simulations the decay of ∆E(t) is one of the slowest.
on the fractal oscillator chain but with nearest neighbor In conclusion, we have studied energy trapping in dis-
potential V(u) = J(u2/2 + u4/4). With N = 64 and ordered classical oscillator chains, with disorder chosen
J = 0.25, the measurements were averaged over 1200 to avoid resonances. If the coupling between oscillators
runs for T 1.0 and 4800 runs for T > 1.0. The re- has a cosine form, we see evidence that, at low temper-
≤
sults for this system are in Figure 3. All the curves atures, energy can be trapped at a site for infinite time,
collapse on top of each other if the t-axis is scaled dif- indicating classical many body localization. This is not
ferently for each T, and the result fits to a stretched seen when the coupling is polynomial.
4
101 We thank Richard Montgomery, David Huse and Sid
Nagel for helpful discussions, and David Huse for sug-
gestingthisproblem. O.N.thankstheInternationalCen-
ter for Theoretical Sciences (ICTS) for their hospitality
during the program Non-equilibrium Statistical Physics
i (ICTS/Prog-NESP/2015/10), where some of this work
E 100
∆ was carried out.
h
Temperature
0.25 2.0 3.5
0.6 3.0
1.0 5.0 3.0
10−1
10−1 100 101 102 103 104 105
t 2.5
101
2.0
i SiteNo.
E
∆ 3 6
h 1.5
1 38
i 10−1 100 101100 1.0 0 54
∆E 100 10 45
h 0.5 42 29
)
J( 26 61
β
0.0
10−1 100 101 102 103 104 105 106
10−1 t
T
10−1
10−1 100 101 102 103 104 105
t FIG.4. Plotsof∆E(t)whentheenergyisinjectedatvarious
sites in the ring, with N = 64,T = 0.25 and J = 1.0. 30
runs were averaged for each plot. The curve labeled l has
FIG.3. (Top)∆E(t)atsitesl=0,1asafunctionoftimefor energy being injected at the sites (l,l+1). The l = 0 curve,
the fractal oscillator model with polynomial interactions, at correspondingtoenergyinjectionatsites(0,1)whichwehave
J =0.25 andvariousvaluesofT.Thecurvesfrom bottom to studied so far, is third from thetop, i.e. one of theflattest.
toparefordecreasingtemperature(Bottom)Rescaledversion
of the same data. A stretched exponential is shown with a
black line. (Bottom inset) Log-log plot of β(T) vsT.
[1] E. Fermi, J. Pasta, and S. Ulam, Studies of nonlinear [8] D.M. Basko, I.L. Aleiner, B.L. Altshuler, Ann. Phys.
problems (Los Alamos Document LA-1940, 1955). 321,1126(2006).ForarecentreviewseeR.Nandkishore
[2] Fora review, see J. Ford,Phys. Rep.213, 271 (1992). andD.A.Huse,Ann.Rev.Cond.Mat.Phys.6,15(2015).
[3] N.J. Zabusky and M.D. Kruskal, Phys. Rev. Lett. 15, [9] A.DharandJ.L.Lebowitz,Phys.Rev.Lett.100,134301
240 (1965). (2008).
[4] M. Toda, J. Phys.Soc. Jpn,22, 431 (1967). [10] V. Oganesyan, A. Pal, D. A. Huse, Phys. Rev. B 80,
[5] E.Fucito,F.Marchesoni,E.Marinari,G.Parisi,L.Peliti, 115104 (2009); D.M. Basko, Ann. Phys. 326, 1577
S. Ruffo and A. Vulpiani, J. Phys. 43, 707 (1982); L. (2011).
Berchialla,L.GalganiandA.Giorgilli,Discr.Cont.Dyn. [11] R. B¨ohmer, K.L. Ngai, C.A. Angell and D.J. Plazek, J.
Syst. A 11, 855 (2004); G. Bennetin, A. Carati, L. Gal- Chem. Phys.99, 4201 (1993).
ganiandA.Giorgilli, inTheFermiPasta Ulamproblem: [12] J.C. Phillips, Rep.Prog. Phys. 59, 1133 (1996).
a status report, G. Gallavotti ed., (Springer, Berlin Hei- [13] T.ProsenandD.K.Campbell,Phys.Rev.Lett.84,2857
delberg 2007). (2000).
[6] F.M.IzrailevandB.V.Chirikov,Sov.Phys.Dokl.11,30 [14] O. Narayan and S. Ramaswamy, Phys. Rev. Lett. 89,
(1966);F.Bocchieri,A.Scotti,B.BearziandA.Loinger, 200601 (2002).
Phys.Rev.A 2, 2013 (1970). [15] A.N. Kolmogorov, Dokl. Akad. Nauk. SSSR 98, 527
[7] For reviews see F. Bonetto, J.L. Lebowitz and L. Rey- (1954); J. Moser, Nachr. Akad. Wiss. G¨ottingen Math.-
Bellet, Fourier’s law: a challenge to theorists, in Mathe- Phys. Kl. II, 1 (1962); V.I. Arnold, Uspehi. Mat. Nauk.
matical Physics 2000,A.Fokas,A.Grigoryan, T.Kibble 18, 13 (1963).
and B. Zegarlinski, eds., (Imperial College Press, Lon- [16] H. Vogel, Phys. Z. 22, 645 (1921); G.S. Fulcher, J. Am.
don, 2000); S. Lepri, R. Livi and A. Politi, Phys. Rep. Ceram. Soc. 8, 339 (1925).
377, 1 (2003); A. Dhar, Adv.Phys. 57, 457 (2008).
