Energy spreading,equipartition,and chaos in lattices with non-central forces

Arnold Ngapasare, Georgios Theocharis, Olivier Richoux, Vassos Achilleos, and Charalampos Skokos,†

1Nonlinear Dynamics and Chaos Group,Department of Mathematics and Applied Mathematics,University of Cape Town,Rondebosch 7701,South Africa

2Laboratoire d’Acoustique de l’Universitdu Mans(LAUM),UMR 6613,Institut d’Acoustique-Graduate School(IA-GS),CNRS,Le Mans Universit ,France

We numerically study a one-dimensional,nonlinear lattice model which in the linear limit is relevant to the study of bending (flexural) waves.In contrast with the classic one-dimensional mass-spring system, the linear dispersion relation of the considered model has different characteristics in the low frequency limit.By introducing disorder in the masses of the lattice particles, we investigate how different nonlinearities in the potential (cubic, quadratic, and their combination)lead to energy delocalization,equipartition,and chaotic dynamics.We excite the lattice using single site initial momentum excitations corresponding to a strongly localized linear mode and increase the initial energy of excitation.Beyond a certain energy threshold, when the cubic nonlinearity is present, the system is found to reach energy equipartition and total delocalization.On the other hand,when only the quartic nonlinearity is activated,the system remains localized and away from equipartition at least for the energies and evolution times considered here.However,for large enough energies for all types of nonlinearities we observe chaos.This chaotic behavior is combined with energy delocalization when cubic nonlinearities are present,while the appearance of only quadratic nonlinearity leads to energy localization.Our results reveal a rich dynamical behavior and show differences with the relevant Fermi–Pasta–Ulam–Tsingou model.Our findings pave the way for the study of models relevant to bending(flexural)waves in the presence of nonlinearity and disorder,anticipating different energy transport behaviors.

Keywords: Anderson localization,energy spreading,energy equipartition,chaos

1.Introduction

The work of Anderson on wave dynamics in disordered systems has been hailed as one of the major milestones of 20th-century physics.[1]Since the publication of that pioneering work, several studies of Anderson localization (AL),which is the exponential localization of almost all eigenmodes of random realizations of mass and/or force constant distributions, have been performed.[2–10]The localization of wave functions(or eigenmodes)leads to the absence of energy transport in disordered systems.Over the years, AL has been extended to other models, apart from the original electronic systems, and related phenomena continue to reveal interesting results across many branches of physics.[11–16]Another scientific landmark of the previous century,which has culminated in a whole new branch of physics, is the work of Enrico Fermi and his collaborators.[17,18]Indeed, their work on the now famous Fermi–Pasta–Ulam–Tsingou(FPUT)problem gave birth to what we now know as computational physics(see Ref.[19]and references therein).

There are two main mechanisms through which waves can be localized in lattice chains, namely(i)disorder and(ii)nonlinearity,resulting in interesting physical phenomena such as AL,[1]solitons,[20]and discrete breathers.[21,22]Solitons are shape-preserving propagating localized formations,while discrete breathers are spatially localized time periodic structures.The interplay of nonlinearity and disorder has been an area of intense active research over the years (see Ref.[23] and references therein for a recent review on the topic of disorder and nonlinearity).In earlier studies, the actual role of nonlinearity on disorder remained not so clear.However,nowadays the theory of the interplay of disorder and weak nonlinearity in lattices such as the disordered Klein–Gordon(DKG) lattice of anharmonic oscillators and the disordered discrete nonlinear Schrdinger equation (DDNLS) has been well developed.[24,25]It is now known that weak nonlinearity leads to subdiffusive energy transport in these models.On the other hand,recent studies of discrete nonlinear models describing mechanical lattices have exploited the strongly nonlinear limit and have revealed superdiffusive energy transport in the presence of disorder.[16,26,27]

Regarding mechanical lattices (seee.g., Refs.[28,29]),one of the most fundamental models describing the longitudinal displacement of point masses connected with springs is the FPUT model.[17,18]For the linear regime of this model it is known that energy propagation in the presence of disorder is dictated by the form of the dispersion relation in the low frequency limit,which starts from zero with a linear slope.[30]The nonlinearity is then found to, generally, enhance energy transport although a clear understanding of the mechanisms is not yet established.[31,32]Here,we wish to extend these studies, by considering a discrete model of a multidimensional conservative Hamiltonian system featuring a different dispersion relation in the low frequency regime.The proposed model is inspired by the recent interest on incorporating other degrees of freedom in mechanical lattices(e.g.,rotations[33])and also on the idea of AL of flexural waves.[34]We study a toy model whose linear limit describes lattice bending vibrations.Our main goal is to explore the interplay of nonlinearity and disorder and identify differences in energy transport with other models (e.g.the DKG and/or the FPUT models), which feature different linear spectra.[23,25,26]In our investigations we consider single site excitations on lattice sites which are not too close to the boundaries.This excitation choice typically results in the excitation of a single localized mode,due to the presence of strong mass-disorder.In particular,we endeavor to answer the question of the conditions under which the energy initially given to a single mode at a single initial excitation will spread to the rest of the lattice and whether the system will eventually reach energy equipartition or not, as well as study its chaotic dynamics.Here we focus on the chaotic behavior of conservative lattice systems,but the role of chaos in various setups of continuous and discrete time models is also a field of active research,seee.g.,Refs.[35–38].

The rest of this paper is organized in the following manner.In Section 2, we describe the proposed one-dimensional(1D) discrete model.In the next two sections, we discuss some representative results about the effects of a momentum single site initial excitation on energy localization, equipartition (Section 3) and chaos (Section 4).In Section 5, we present a statistical analysis of the localization, equipartition and chaotic behaviors of the model for a range of energy values.Finally,in Section 6,we summarize our work and present its conclusions.

2.The Hamiltonian model and its equations of motion

We consider a discrete nonlinear model described by the following Hamiltonian

where ζn=(un+1−un)−(un−un−1)andPn=mn˙un.In our notation[˙]denotes derivative with respect to time andmncorresponds to the mass of then-th site.In the linear limit,i.e.,whenK3=K4=0, equation(1)describes the transverse displacements of a mass spring chain as illustrated in Fig.1.[39]It is worth noting from the sketch in Fig.1 that for sufficiently small rotation angles θn,the relation

holds, withabeing the distance between neighboring sites.For simplicity in our study we seta=1.

Fig.1.A schematic representation showing the relation between the rotation angles θn,and the relative transverse displacement of the point mass at un of a lattice with distance a between neighboring masses(in our study we set a=1).In the small angle approximation the equations of motion for the lattice is the linear limit of the Hamiltonian(1).

For the purposes of our work,we consider three different cases regarding the nonlinear terms.One case has only cubic nonlinearity (K4=0) which we refer to asH3, an additional case with only quartic nonlinearity(K3=0)is,from now on,referred to asH4,and finally the general case where both nonlinear terms are present(K3/=0,K4/=0),which we refer to asH34.Note also that in this workK2>0 and we consider disorder only in the masses.The coefficientsKi,i=3,4 (when non-zero) are kept constant and the corresponding equations of motion for the generalH34model are

2.1.Dispersion relation and eigenmodes

Firstly,let us consider some aspects of the linear periodic case corresponding toK3=K4=0 andmn=m.We assume Bloch-like solutions of Eq.(2) of the formun=Aei(qn−ωt).This solution models traveling waves in a lattice chain,so that each particle oscillates about its equilibrium position with the same frequency ω,and amplitudeA,whileqdenotes the Bloch wave number.Using this ansatz we obtain the system’s dispersion relation,which is given by the following equation:

and is also plotted in Fig.2(a).The dispersion relation of this model is quite different from that of the harmonic chain with central forces,leading to some interesting dynamics.The uniqueness stems from the particularity of zero group velocity in the low frequency regime.Note that for low frequencies the dispersion relation of Eq.(3)resembles the one of a wave equation with fourth order spatial derivative describing bending waves in elastic beams.[39]

In what follows we consider finite chains with a size ofN=28particles.In addition, we introduce disorder in the particle masses such that the random masses are taken from a uniform distribution,i.e.,mn∈[0.1,1.8]while we fix the linear potential coefficient toK2=1.The nonlinear coefficientsK3,K4,take values of either 1 or 0 such that whenever any of the nonlinear coefficients is not zero,it is then set to 1,e.g.,for the systemH4we setK3=0 andK4=1.In Ref.[31]it was shown that the results do not significantly differ if either mass or stiffness disorder is implemented hence for this reason we will consider only the case of mass disorder.

We now describe the localization properties of the linear normal modes of this system.For this purpose, we consider solutions of the linear part of the system ofNmasses given in Eq.(2) to be of the form u(t) = U e−iωt, where U =[U1,U2,...,UN]Twith superscript T denoting the transpose of a matrix.Hence,we obtain the eigenvalue problem

where ωkare the eigenfrequencies and Ukthe corresponding eigenvectors.M is a diagonal matrix containing the mass elements,whilst K is a sparse diagonal matrix of the coefficients.We show an eigenfrequency spectrum for a representative disordered chain of sizeN=28in Fig.2(b).First we observe that the very low frequency part of the spectrum is similar to the ordered case(where all masses are equalmn=m=1),while the upper part has changed significantly and the cutoff frequency has been increased.

Fig.2.(a) Dispersion relation (3) of the lattice bending waves where a is the lattice spatial constant (which is set to a=1 in our study), q the wave number and ω the frequency.(b)The eigenfrequency spectrum of a representative disorder chain of size N=28 with the arrows indicating eigenfrequencies for an extended (blue arrow) and a localized (red arrow) normal mode.The eigenfrequencies have been sorted and indexed from lowest to highest frequency.

Next, we are interested to characterize the localization properties of the modes which can be done by computing the participation number

wherehn=Hn/His the normalized energy of each particle(seee.g., Ref.[25] for a further discussion of this quantity).First we note that, as is the case for the classic mass-spring model with central forces with disorder,[40]the lower part of the spectrum comprises of extended modes while the upper part of mostly localized ones.To illustrate that, in Fig.3(a)we plot the profiles for the two characteristic modes indicated in Fig.2(b).To have a more global picture and for the sake of comparison, in Fig.3(b) we also plot the averaged over 500 disorder realizations,mean participation number〈P〉,for linear chains of two types.In particular we contrast for the same disorder parameters(the random masses are taken from a uniform distribution,i.e.,mn∈[0.1,1.8])the harmonic chains with central forces,for which the equations of motion becomemn¨un=un+1+un−1−2un(blue curve)and non-central forces obtained by settingK3=K4=0 in Eqs.(1)and(2)(red curve).We note some of the most glaring differences which are: (i)modes with the lower indices(k20)have higher〈P〉for the chain with central forces than for the chain with non-central forces and(ii)the chain with central forces has more localized modes as compared to the chain with non-central forces at any lattice size[〈P〉≈2 fork160 in Fig.3(b)].

Fig.3.(a) Profiles of an extended (blue curve), as well as a localized (red curve) normal mode respectively indicated by the blue and red arrows in Fig.2(b), for a mass disordered harmonic chain with non-central forces.(b) Mean participation number 〈P〉 over 500 disorder realizations of linear eigenmodes for harmonic chains with non-central forces(red curve)and harmonic chains with central forces(blue curve)of size N =28.In panel(b),the modes have been sorted and indexed from the lowest frequency to the highest frequency modes.One standard deviation in the computation of〈P〉is indicated as the lightly shaded areas in panel(b).

3.Evolution of linear modes,energy delocalization and equipartition

In this section, we consider the fate of a single site momentum initial excitation under the influence of the different kinds of nonlinearities(see the three different kinds of Hamiltonians presented in Section 2).In Fig.3(a)we depict,using a red curve,an almost single site localized mode,namely modek=252,for which only siten=187 is significantly displaced from its equilibrium position.Note that the sizeNof the lattice determines the number of highly localized modes.In our study,we set the size of the chain to beN=28in order to have a significant number of highly localized modes and implement fixed boundary conditions(u0=un+1= ˙u0= ˙un+1=0).The dynamics is studied by numerically solving the different types of equations of motion given in Section 2 using the ABA864 symplectic integrator.[41–43]This integration scheme allows the accurate conservation of the total energyH,and keeps the relative energy error

less than 10−4when the integration time step is τ =0.1.For illustrative purposes,we show in Fig.4 a representative example of the time evolution of ΔH(t)for a case of theH4system.From the presented results it is clearly seen that the used integration scheme keeps the error bounded for the whole duration of the simulation and that its upper bound does not increase in time.

Fig.4.Temporal evolution of the relative energy error ΔH(t)(6)of a single disordered realization with energy H4=0.2 and time step τ =0.1.

3.1.Spatiotemporal evolution of modes

To better understand the fate of single site initial excitations,let us discuss the evolution of normal modes.To do so,we transform the linear part of the Hamiltonian (1), and cast it in terms of the system’s normal modes.At any timet, we define the velocity vector V(t)=[˙u1(t),...,˙uN(t)]T.The projection of this vector onto the system’s normal modes is given by T =A−1V(t)with matrix A having the lattice eigenvectors as columns.In the same manner,the displacement vector U(t)=[u1(t),...,uN(t)]Tis projected onto normal modes to yield R=A−1U(t).TkandRkrespectively are elements of the projection vector corresponding to normal mode momenta and positions of thek-th mode (k=1,...,N).Thus the total energy of the linear modes is given as

We first consider the spatiotemporal evolution of the linear modes after a single site momentum excitation of the system has been implemented for the various nonlinearities.The evolution of normal modes forH3is shown in Fig.5 for both the near linear regime[panel(a)]as well as the nonlinear regimes[panels(b)and(c)].We shift between these regimes by changing the system’s energy.

More specifically, in Fig.5(a), we depict results for the evolution ofH3in normal mode space after a single site momentum initial excitation for the near linear regime where mainly the localized modek= 252, is significantly excited(≈90% of the total energy), although some other normal modes are weakly excited as well(≈10% of the total energy).However,no coupling of the normal modes takes place during the dynamical evolution.On the other hand,an increase in energy entails that the system becomes more nonlinear.We depict results for this nonlinear regime in Fig.5(b).There we see that the nonlinearity due to the cubic potential is able to facilitate coupling of the linear normal modes until the modes are well mixed especially aftert104.We shall investigate further this mixing of linear eigenmodes in Subsection 3.3 where we will discuss the route to energy equipartition.Surprisingly enough,the quartic nonlinearity(H4)does not appear to facilitate normal mode mixing as depicted by the results in Fig.5(c)which show that,although other normal modes besides modek=252 are initially excited,they do not mix at all.The results forH34(not shown here)are similar to those forH3.

Fig.5.Spatiotemporal evolution of normal mode energy content after a momentum single site initial excitation for the chain with non-central forces(a)H3 with energy H3 =10−4, (b) H3 with energy H3 =0.2, and (c) H4 with energy H4=0.2.The color bar is in logarithmic scale.

3.2.Energy delocalization

We now turn our attention onto another issue and investigate how the initially localized energy spreads in time to the rest of the lattice sites when the system is initialized by a single site momentum excitation.In particular,we follow the lattice dynamics for initially localized energy excitations for systems modeled by the HamiltoniansH3,H4, andH34for the same value of total energy.More specifically, we setH=0.2 and we show results for momentum single site initial excitations in Fig.6 for the three types of nonlinearities.In Fig.6(a),we show the spatiotemporal evolution of the initially localized wave-packet when the system is described byH3.The results in Fig.6(a)show the delocalization of the wave-packet which eventually leads to energy spreading to the whole lattice.This happens especially for the time intervalt104.Results for the dynamics when the system is governed byH4are depicted in Fig.6(b) showing that the wave-packet remains localized at least for the duration of our simulation.For a combination of these two types of nonlinearities,i.e.,the cubic and quartic nonlinearity, the dynamical behavior again reverts to showing energy spreading for the initially localized wave-packet around siten=187 [Fig.6(c)].To ascertain the energy localization we integrate the systems for the three cases up tot=106time units,which is well beyond the times we initially observe energy spreading.With this in mind we set the final integration time att=106for all our simulations.

Fig.6.Spatiotemporal evolution of the energy density after momentum single site initial excitation.Results are shown for systems(a)H3,(b)H4,and(c) H34.The total system energy for all three cases is H =0.2.The black curves in all plots depict the mean position of the energy density.The used color bar is in logarithmic scale such that yellower (lighter) regions have more energy than greener(darker)regions.

In order to quantify the energy localization, we consider the temporal evolution of the participation numberPgiven by Eq.(5)for the three cases.We find that the nonlinear dynamics in cases where the cubic nonlinearity is present (i.e., for the cases ofH3andH34),exhibit saturation of the participation number toP≈125 fort104.This is depicted by the blue and green curves in Fig.7(a)and the saturation ofPis clearly observed by the flattening of these two curves,especially during the last two decades of the evolution.However,forH4,where only the quartic nonlinearity is present, the dynamics reveal localized behavior with no energy spreading as depicted by a practically flat curve ofPin Fig.7(a) (red curve) [P(t)≈2].We believe that this absence of wave-packet spreading is possibly due to the non-excitation of other modes[see Fig.5(c)]which does not permit the spreading of energy to the rest of the lattice.It is interesting to note that this behavior is in contradiction to what is observed in the two versions of the homogeneous FPUT model,where it is known that the α-FPUT system(cubic nonlinearity)thermalizes much later than the β-FPUT one(quartic nonlinearity).[57]Even though our system is disordered,H3andH34reach equipartition faster thanH4.A combination of the two kinds of potentials(H34)shows that the dynamics is strongly dominated by the cubic potential as its behavior is similar to the one observed for theH3model.In fact, there are not many differences between the dynamics ofH3andH34.In fact by comparison of Figs.6(a) and 6(c)as well as the blue and green curves in Fig.7(a)it is clear that there are not many differences between the dynamics ofH3andH34.

Fig.7.Time evolution of(a)P (Eq.(5)),and(b)η (Eq.(9));for systems H3 (blue curves),H4 (red curves)and H34 (green curves).The total system energy for all three cases is H=0.2.The dashed horizontal line in panel(b)indicates the analytical value of〈η〉when energy equipartition is reached.

3.3.Energy equipartition

It is also interesting to determine whether the system eventually reaches energy equipartition for the three types of nonlinearities.In order to do so,we numerically compute the so called “spectral entropy” by monitoring the corresponding time evolution of the normal modes.[44]To this end we write the weighted harmonic energy of thek-th mode as

whereEkin Eq.(7) is thek-th linear mode energy.Then the spectral entropy at timetis computed as

where 0

Thus,through such normalization,0 ≤η ≤1.When η is close to one,the dynamics do not substantially deviate from the initially excited mode(s).However,as more modes are excited,η decreases and approaches zero.For a system at equipartition,a theoretical prediction for the mean entropy η exists, which assumes that the modes follow a Gibbs distribution when the nonlinearity is weak.The analytical form of the mean entropy〈η〉anin that case is given as

withC≈0.5772 being the Euler constant.[45,46]In Fig.7(b),we show the time evolution of the spectral entropy η in order to determine how the system approaches energy equipartition.When η is equal to the analytically predicted value of Eq.(10),we conjecture that energy equipartition is achieved.Again a similar behavior, in terms of how energy equipartition is approached, is observed betweenH3andH34where the cubic nonlinearity is active.This is shown in Fig.7(b)where there is close resemblance between the results forH3andH34, indicated by the blue and green curves.These curves also show that energy equipartition (indicated by η ≈〈η〉an=0.0930)is achieved at around the same time as energy delocalization.For the case where the system only has a quartic nonlinearity(H4),we find no energy equipartition,as shown by the corresponding time evolution of η,which remains fairly constant at almost η =1 in Fig.7(b).

4.Chaotic behavior

Chaoticity is another aspect that is interesting to study about the dynamics of the system.To study chaos, in this work we utilize a global indicator of chaos called the maximum Lyapunov exponent(mLE),which we compute using the so-called standard method.[47–49]The finite time mLE(ftmLE)is defined as

where w(t) is a vector of small perturbations from the phase space trajectory at timet, also called deviation vector, which we denote as

with δun(t)and δpn(t)respectively indicating small perturbations in position and momentum for thenth lattice site.The mLE is defined as Λ =limt→∞λ(t).In Eq.(11),||·||denotes the usual Euclidean vector norm.For chaotic trajectories, λ attains a finite positive value,otherwise,λ ∝t−1for regular orbits.An efficient and accurate method to follow the evolution of w(t) is to numerically integrate the so-called variational equations,[50]which govern the vector’s dynamics, together with the Hamilton equations of motion using the tangent map method outlined in Refs.[51–53].Let us now investigate the chaotic behavior of the system for the three types of nonlinearity under study.The time evolution of λ reveals in all cases a chaotic response,as shown by the ftmLE which deviates from λ(t)∝t−1(see Fig.8).In fact, forH3andH34we attain the same final values of λ towards the end of the simulations as shown by the blue and green curves in Fig.8.In our study,we also utilize the so-called the deviation vector distribution(DVD) to analyze chaos further.The DVD is able to identify spots in real or phase space which are more sensitive to chaos.[54–56]The DVD is written as

We are thus able to describe the spatial distribution of chaos in a dynamical system by making use of the DVD.We now identify the most active chaotic spots in the lattice by following the time evolution of the DVDs.We use as a representative case the results ofH3indicated by the blue curve in Fig.8.The spatiotemporal evolution of the DVD is plotted in Fig.9(a)and the figure shows that,on average,the mean position of the DVD remains close to the excitation region during the early stages of the dynamicst103. Beyondt103,the mean position of the DVD starts to move throughout the chain.This is shown by the erratic motion of the mean position(solid black curve in the figure),which implies that there is a chaotic spot randomly moving throughout the lattice. In order to clearly show this behavior, we take some snapshots of the DVD especially towards the end of the simulation (t105).These snapshots show peaks at different positions implying that indeed the chaotic spots are moving throughout the chain as shown in Fig.9(b).

Fig.8.The time evolution of λ (11) for systems (i) H3 (blue curve), (ii)H4 (red curve), and (iii) H34 (green curve) for the same momentum initial excitation discussed in Section 3.The dashed diagonal line guides the eyes for slope −1.

Fig.9.(a)Spatiotemporal evolution of the DVD[Eq.(13)]and(b)instantaneous snapshots of the normalized DVDs at t =T1 ≈2.47×103 (black curve), t =T2 ≈2.08×105 (red curve), and t =T3 ≈2.477×105 (blue curve) respectively corresponding to black, red and blue horizontal dashed lines in panel (a).The results are for the momentum initial excitation discussed in Section 3 for H3=0.2.Panels(c)–(d)are similar to Panels(a)–(b)but for H4 =0.2.The snapshots were taken at t =T1 ≈103,t =T2 ≈105,and t=T3 ≈5×105.The black curves in panels(a)and(c)depict the mean position of the DVD.The color bar is in logarithmic scale such that yellower(lighter)regions have more energy than greener(darker)regions.

In the case of the combination of cubic and quartic potentials(H34)the results are qualitatively similar to those presented in Figs.9(a)–9(b)hence we do not show them here.The case which has only a quartic potential(H4)is also found to be chaotic,as shown by its λ(t)(red curve in Fig.8).Upon closer analysis of the DVDs,the chaotic behavior exhibited byH4is revealed to be different from the behavior already described for the other two types of nonlinearities which feature cubic potentials(H3andH34).In particular,the DVD indicates that chaos remains localized around the excitation region forH4as is shown by the DVD and the corresponding snapshots in Figs.9(c)and 9(d).Note that w(0)is centered slightly to the right of the excitation site in Fig.9.

5.Statistical analysis

We further study how the system behaves statistically by considering different mass disorder realizations for several energy levels fromH=0.05 up toH=0.3. In particular, we present results obtained over 30 realizations at each energy level. The strength of the disorder is such that there is a significant number of strongly localized modes and many of those are in the center (they are located all along the chain) similar to the one already presented in Section 2. The results forH3indicate that the mean participation number 〈P〉 saturates to 〈P〉≈125 (or ≈N/2) as the energy is increased forH3[Fig.10(a)]. We also find a threshold energy(H0.2)above which, on average, any realization for which a highly localized mode if taken as initial condition leads to this maximum value ofP≈125. For the quartic potential (H4), we find no energy delocalization and〈P〉remains very low as the energy is increased as shown by the red curve in Fig. 10(a) where〈P〉≈2. We studied the system for even higher energy excitations up toH=0.9 and we still found the system’s density distribution to remain localized.

The statistical results of the mean spectral entropy, 〈η〉show contrasting behaviors forH3andH4.H3is found to reach energy equipartition as the predicted analytical value is reached (〈η〉 ≈〈η〉an) for what appears to be an energy threshold of aroundH0.2, as shown by the blue curve in Fig. 10(b). On the other hand,H4does not reach energy equipartition even at higher energy excitations than those presented in Fig. 10. The results forH4are shown by the red curve in Fig.10(b)for which〈η〉≈1 for all initial energies indicating no energy equipartition. In fact the smaller error bars on the blue curves for energiesH0.2 in both Figs. 10(a)and 10(b) validates these conclusions. This result is in some sense the inverse of what was observed for the homogeneous α-FPUT model, which attains equipartition later than the β-FPUT model.[32]In the latter model, one discovered route to equipartition is through modulational instability which has an energy threshold.[32,57]

As for the system’s chaoticity,the final value of the mean ftmLE 〈λ〉 does not practically depend on the system’s total energy;for as depicted in Fig.10(c)where〈λ〉obtains values around 1.5×10−3for bothH3andH4,at least for the energy range considered in this study, although a slight increase of〈λ〉 is seen as energy grows.H4is found to be in general less chaotic thanH3as shown by the lower values of〈λ〉[red curves in Fig.10(c)] forH4.Note that the statistical results forH34are qualitatively similar to those forH3presented in Fig.10.

Fig.10.Mean(a)participation number〈P〉,(b)entropy〈η〉,and(c)ftmLE,〈λ〉for H3(blue curves)and H4(red curves)against system’s energy H,over 30 disorder realizations for which highly localized modes are taken as initial excitations.The error bars represent one standard deviation.The dashed line in panel (a) guides the eye to the level =125, whilst a similar line in panel (b) guides the eye to the analytical mean entropy at equipartition[Eq.(10)].

6.Conclusions

In this study, we numerically investigated the dynamics of a disordered 1D nonlinear lattice model which is relevant to the study of bending waves,focusing on the energy delocalization and equipartition,as well as the system’s chaotic behavior.For this study, we chose a localized wave-packet which is almost similar to a highly localized mode as the initial energy excitation.We considered three kinds of nonlinear configurations of the Hamiltonian: (i)with only a cubic nonlinearity named asH3, (ii) with only a quartic nonlinearity calledH4,and finally (iii) a combination of the cubic and quartic terms namelyH34.

The dynamics shows that the energy delocalization strongly depends on the system’s type of nonlinearity.A threshold exists beyond which any initial excitation will lead to energy delocalization for the cases ofH3andH34in which the cubic potential dominates the dynamics.The dynamics forH4remains localized throughout the whole simulation even for the large values of the energy.Similarly,energy equipartition is observed only when the cubic terms are present,even from a statistical viewpoint when the system is studied at different energy levels.

Additionally, for the same initial conditions, we studied the system’s chaotic response using the finite time mLE and the related DVDs.Here all three types of nonlinearities induce chaotic responses.Statistically,the chaotic behavior appears to be slightly stronger when cubic terms are included in the model (systemsH3andH34), with the strength of the chaoticity moderately increasing as the system energy grows.The chaotic behavior of the systems as analyzed using the DVDs shows both localized and delocalized chaos.Localized chaos is exclusive toH4even for energy excitations as high asH=0.9.

We find that our numerical results have not only revealed interesting nonlinear dynamical behaviors but also open the road for further theoretical studies of chains with non-central forces.This is due to the fact that the proposed model suggests fundamentally different dynamics than the ones of wellstudied discrete systems like the DKG and FPUT models,for which we have a better understanding.In particular,our results regarding the very different behavior between the cubic and quartic nonlinearities,poses the question of what are the mechanisms responsible for these observations especially in view of the comparison with the α and β-FPUT models.Along the same lines, an interesting question is how the energy spreading is obtained in infinite chains and whether the increase of nonlinearity leads(or not)to subdiffusive behaviors.

Our work furthermore contributes to the debate about the potential localization or spreading of initially localized excitations in disordered nonlinear systems when such processes are juxtaposed to the system’s chaotic strengths.

Acknowledgements