A Ngapasare

G Theocharis

O Richoux

Ch Skokos

V Achilleos

Chaos and Anderson Localization in Disordered Classical Chains: Hertzian vs FPUT models

We numerically investigate the dynamics of strongly disordered 1D lattices under single particle displacements, using both the Hertzian model, describing a granular chain, and the α + β Fermi-Pasta-Ulam-Tsingou model (FPUT). The most profound difference between the two systems is the discontinuous nonlinearity of the granular chain appearing whenever neighboring particles are detached. We therefore sought to unravel the role of these discontinuities in the destruction of Anderson localization and their influence on the system's chaotic dynamics. Our results show that the dynamics of both models can be characterized by: (i) localization with no chaos; (ii) localization and chaos; (iii) spreading of energy, chaos and equipartition. The discontinuous nonlinearity of the Hertzian model is found to trigger energy spreading at lower energies. More importantly, a transition from Anderson localization to energy equipartition is found for the Hertzian chain and is associated with the"propagation" of the discontinuous nonlinearity in the chain. On the contrary, the FPUT chain exhibits an alternate behavior between localized and delocalized chaotic behavior which is strongly dependent on the initial energy excitation.

I. INTRODUCTION

The role of nonlinearity in disordered systems which exhibit Anderson localization [1,[START_REF] Abrahams | 50 Years of Anderson Localization[END_REF] is a topic that triggered a vast amount of theoretical, numerical [START_REF] Kopidakis | [END_REF][4][5][6][START_REF] Flach | Nonlinear Optical and Atomic Systems[END_REF] and experimental studies [8][9][10][11]. The two principal questions under consideration are (a) does the energy carried by localized wave-packets eventually spread or not and (b) what is the route to equipartition? Among different nonlinear models, large theoretical work and progress has been made especially for the Klein-Gordon (KG) system and the Discrete Nonlinear Schrödinger (DNLS) equation with disorder. For these systems, it has been found that the combined influence of disorder and nonlinearity leads to sub-diffusive energy transport [5]. It is also now understood that whether nonlinear Anderson localization persists or is destroyed has probabilistic features and is directly associated with chaos [12]. Additionally a variety of different physical settings have been exploited in order to study this interplay between nonlinearity and disorder, especially in optical and atomic systems [9,13].

Recently a classical lattice i.e. the granular chain described by the Hertzian contact force [START_REF] Johnson | Contact Mechanics[END_REF], has also attracted much attention in the same context [START_REF] Pikovsky | [END_REF][16][17][18][19][20]. The considerable interest in the Hertzian chain can be attributed to the strong nonlinearity of the system which is however easily tuned (usually by the pre-compression of the chain). The Hertzian contact forces also allow access to wave propagation in an almost linear system up to the case of a lattice where only nonlinear waves propagate ("sonic vacuum") [START_REF] Theocharis | Daraio in Acoustic Metamaterials and Phononic Crystals[END_REF][START_REF] Ponson | [END_REF]. An additional interesting dynamical feature of the granular chain is that the power law nonlinearity, due to the Hertzian force, coexists with a non-smooth nonlinearity describing de-tached particles [23][START_REF] Mitchell | Fundamentals of Soil Behavior[END_REF][START_REF] Andreotti | Granular media: between fluid and solid[END_REF][START_REF] Yunker | [END_REF]. Recent works on both uncorrelated and correlated disorder granular chains showed that the system traverses from a sub-diffusive regime for sufficiently weak nonlinearities to a super-diffusive regime for increasing nonlinearity [19]. In strongly disordered granular chains it was found that localization coexists with chaos and equipartition is reached for finite times [20].

Furthermore, the granular chain in the weakly nonlinear regime provides an experimental setting to study the Fermi-Pasta-Ulam-Tsingou (FPUT) model with both α and β type terms [27,[START_REF]The Fermi-Pasta-Ulam Problem: A Status Report[END_REF]. In contrast to the DNLS and KG models, which have been studied for disordered systems, the FPUT system has been mostly studied in the homogeneous case, although some studies regarding disorder also exist e.g. in Refs. [START_REF] Dhar | [END_REF][30][31].

In fact, the phenomenon of equipartition for the homogeneous case is a long standing problem, which originates from the pioneer work of Fermi, Pasta, Ulam and Tsingou [27,[START_REF]The Fermi-Pasta-Ulam Problem: A Status Report[END_REF], although substantial progress has been made on the subject [32,33]. Very recent studies both in α-FPUT (but also in the KG model) periodic lattices showed that the thermalization is reached through high order resonant interactions leading to large timescales for equipartition [START_REF] Onorato | Proc. Natl. Acad. Sci. USA[END_REF][START_REF] Pistone | [END_REF]. It was also found that the fluctuations of the entropy after the system reaches equipartition are characterized by sticky dynamics close to qbreathers for the FPUT model and discrete breathers for the KG model [36].

In this work we aim to expose the role of different nonlinearities in the destruction of Anderson localization, the chaoticity of the system but also the timescales to reach equipartition. To do so we perform a detailed comparison between the granular chain model and the FPUT system. It is beyond the scope of our work to study the differences between the two models on general grounds. Our interest is to study the fate of strongly localized modes and focus on the role of discontinuous nonlinearities in their dynamics. A statistical analysis of the linear limit, which is common for both models, shows that for sufficiently strong disorder, the system acquires a significant number of strongly localized, almost single particle, modes. From this ensemble we choose a representative realization to illustrate its nonlinear behavior. Our goal is to identify the mechanisms that lead to energy spreading of an initially excited localized mode. Additionally, we use chaos indicators [37,38] to quantify the total systems' chaotic behavior. We provide information about chaos propagation in the lattice enabling us to differentiate localized and extended chaos. By tracking the mode distribution during the dynamics' evolution we provide insights regarding equipartition.

The paper is organized as follows. In Section II we introduce the two models and establish the disorder strength for the lattice which ensures strongly localized modes in the linear limit. Selecting a single configuration and focusing on a highly localized mode near the center of the lattice, in Section III we study the mode's evolution in both models for increasing excitation energy. A thorough analysis of the lattice dynamics is performed focusing on the spreading of the initially localized mode, on the chaoticity of the system and on monitoring the appearance of particle detachments. Finally we show results illustrating how and for which energy the two systems reach energy equipartition. In section IV we summarize our findings and discuss their significance.

II. HERTZIAN AND FPUT MODELS WITH DISORDER

Both models studied here, namely the granular chain with Hertzian interactions and the FPUT system, are considered to be energy preserving (i.e. without losses). Their total energy for a chain with N spherical homogeneous beads of radius R n and mass m n (n = 1, 2, 3, . . . , N ) is given by the following Hamiltonian

H = N ∑ n=1 H n = N ∑ n=1 p 2 n 2m n + V (Hz,F ) n . ( 1
)
Here, p n = m n un and u n denote respectively the momentum and displacement from equilibrium for each particle, (˙) denotes the first order time derivative, while the random radii R n are uniformly chosen in the interval

[min(R n ), max(R n )].
The Hertzian potential V Hz n for each bead due to the nearest neighbor coupling is defined as

V Hz n = [V Hz (u n ) + V Hz (u n+1 )]/2 where V Hz (u n ) = 2 5 A n [δ n + u n-1 -u n ] 5/2 + - 2 5 A n δ 5/2 n -A n δ 3/2 n (u n-1 -u n ). (2)
The static overlap δ n between two neighboring beads n-1 and n is given by δ n = (F 0 /A n ) 2/3 where F 0 is the precompression force. The coefficient A n for spherical beads is given as

A n = (2/3)ε √ R n-1 R n /(R n-1 + R n )/(1 -ν 2 )
where ε and ν are the elastic modulus and the Poisson ratio respectively [START_REF] Johnson | Contact Mechanics[END_REF]. The plus sign in [•] + describes the fact that this term is present as long as δ n +u n-1 -u n > 0 and is absent otherwise, since then the particles are no longer in contact. This is the non-smooth nonlinearity which substantially differentiates the two models.

The FPUT model is described by Eq. ( 1) with a potential

V F (u n ) = 4 ∑ k=2 K (k) n (u n -u n-1 ) k . ( 3
)
Accordingly, the potential of the n th particle is written as

V F n = [V F (u n ) + V F (u n+1 )]/2.
For the rest of this work we consider a chain of N = 40 particles. In our simulations we choose units corresponding to a mean radius of R = 0.01m, and a static precompression force F = 1N. The mean radius is used as a reference to the uniform system with particles of radius R = (α + 1) R/2. The disorder strength, is quantified by the parameter α = max(R n )/ min(R n ). This choice of disorder naturally leads to a random distribution of both the masses and stiffness coefficients [16]. In all calculations we use fixed boundary conditions with dummy beads on both ends such that u 0 = u N +1 = 0 and p 0 = p N +1 = 0. The corresponding equations of motion for the Hertzian model [START_REF] Abrahams | 50 Years of Anderson Localization[END_REF] are

m n ün = A n [δ n + u n-1 -u n ] 3 2 + -A n+1 [δ n+1 + u n -u n+1 ] 3 2 + (4)
while for the FPUT model we obtain

m n ün = 4 ∑ k=2 [ K (k) n+1 (u n+1 -u n ) k-1 -K (k) n (u n -u n-1 ) k-1 ] . (5)
A direct connection between the two models is made by taking the Taylor series expansion of Eq. ( 4) up to fourth order (assuming small dipslacements) u n /δ n,n+1 ≪ 1. Doing so we recover Eq. ( 5) with coefficients K (2) [START_REF] Theocharis | Daraio in Acoustic Metamaterials and Phononic Crystals[END_REF]. Below we normalize our units such that for the linear homogeneous chain with α = 1 the frequency cut-off is ω max = √ 4K m = 1 with K = K (2) n and m = m n = 1.

n = (3/2)A n δ 1/2 n , K (3) n = -(3/8)A n δ -1/2 n and K (4) n = (3/48)A n δ -3/2 n

A. Linear mode analysis of the disordered chain

In this work we are interested on the fate of strongly localized modes. Thus, we first identify the sufficient disorder strength able to sustain a significant amount of localized modes. To do so, we perform a statistical analysis of the linearized equation of motion m n ün = K (2) n+1 (u n+1 -u n ) -K (2) n (u n -u n-1 ) (6) which is common for both models. Assuming harmonic solutions of the form U(t) = U 0 e iωt , where U 0 is a column matrix with elements U n , n = 1, 2, 3, . . . , N . We then solve the corresponding eigenvalue problem

-ω 2 MU 0 = KU 0 . ( 7
)
The matrix M is a diagonal matrix with elements m n and K is a sparse diagonal matrix containing the stiffness coefficients K

(2)

n . To quantify the localization properties of the disorder system, we calculate the participation number [START_REF] Flach | Nonlinear Optical and Atomic Systems[END_REF] of the wave-packet P = 1/ ∑ h 2 n where h n = H n /H. This quantity is defined in a way so that its maximum value equals the total number of particles (extended mode) and its minimum value equals to 1 when only one particle is participating in a mode (strongly localized mode).

In Fig. 1(a), we show the mean value ⟨P ⟩ of the participation number of the eigenmodes for different disorder strengths, using an ensemble of 1000 disorder realizations. The modes are sorted with descending values of P for each realization. For relatively weak disorder (e.g. for α = 2) ⟨P ⟩ largely deviates for the homogeneous case (α = 1) and some localized modes appear in the system. On the other hand for values of α ≥ 4 the averaged participation number reaches a limiting curve with about 10 strongly localized modes with ⟨P ⟩ ≈ 2. The above analysis, provides clear evidence that a single disorder realization with α = 5 is sufficient for the chain to possess several strongly localized modes.

To monitor the spreading we calculate the time evolution of the energy density h n and the participation number P . At the same time we identify and quantify chaos in the system using the maximum Lyapunov characteristic exponent (mLCE) [37,38], which is obtained by numerically integrating the corresponding variational equations [39]. The two sets of equations where integrated using the so called "Tangent Map" method with a fourth order optimal integration scheme with a marching step of 5 × 10 -4 in all our simulations [39,40]. The variational equations govern, at first order of approximation, the time evolution of a deviation vector ⃗ v(t) = [δu 1 , δu 2 , . . . , δu N , δp 1 , δp 2 , . . . , δp N ] where δu n , δp n , n = 1, 2, . . . , N are respectively small perturbations in positions and momenta (see e.g. [38]).

The mLCE is given by λ = lim t→∞ Λ(t), where

Λ(t) = 1 t ln ||v(t)|| ||v(0)|| , ( 8
)
is the so-called finite time mLCE [38]. Note that in Eq. ( 8), || • || denotes the usual Euclidean vector norm.

For chaotic orbits, Λ(t) eventually converges to a positive value, while for regular orbits it tends to zero following the power law Λ(t) ∝ t -1 [38].

To gain more insight about the spatial properties of chaos, we calculate the deviation vector density (DVD) given by

w n = δu 2 n + δp 2 n ∑ n (δu 2 n + δp 2 n ) . ( 9
)
The deviation vectors are known to align with the most unstable region in phase space. They have been employed in disorder nonlinear lattices in order to visualize the spatial evolution of the most chaotic regions [41][42][43]. Here, we make use of the DVDs in order to spatially characterize the chaoticity of the system, either as localized or extended chaos. The initial condition used for the deviation vectors ⃗ v(0) is a random uniform distribution of momentum perturbations δp i as for this choice, the time evolution of the finite time mLCE was found to converge faster to the Λ(t) ∝ t -1 law for regular orbits.

The numerical results shown in the rest of this work (unless stated otherwise), are performed using a representative single realization of the statistical ensemble for α = 5. The corresponding eigenfrequencies of this realization are shown in Fig. 1(b). Generally, low frequency modes extend over many particles, whilst high frequency modes are localized. We identify the 34 th as a strongly localized mode (P ≈ 2.5) located in the middle of the chain at site n = 21 as shown in the inset of Fig. 1

(b).

To study the effect of nonlinearity, we initially excite the 21 st site which results in the excitation of almost only the 34 th mode, and we monitor the evolution of both models as we increase the initial excitation energy.

III. DYNAMICAL EVOLUTION OF AN INITIALLY LOCALIZED MODE

A. Near linear limit

For sufficiently small energies H, we have numerically confirmed that the two models behave both qualitatively and quantitatively the same. An example is given in Fig. 2 which corresponds to H = 0.25. As shown in panels (a) and (b) the energy density for both models is completely localized around the initially excited site n = 21 as shown by the black solid line which indicates the mean position of the energy distribution. Localization is quantified by the almost constant value of P ≈ 1.8 for both models, shown in Fig. 2(c). The curves of the Hertzian (red) and the FPUT (blue) models almost overlap. The time evolution of Λ(t) is depicted in Fig. 2(d) and confirms that the dynamics is regular as Λ(t) follows the power law Λ(t) ∝ t -1 .

The spatiotemporal evolution of the corresponding DVD, plotted in Figs. 3(a deviation vector distribution in contrast to the localized, pointy shape that DVDs exhibit for chaotic orbits [41][42][43]. Accordingly, particular profiles of the DVDs taken at different times shown in Figs. 3(c) and (d) are found to be extended covering the whole excited part of the lattice in a relatively smooth way.

However, a difference between the two models is found by closely inspecting the corresponding participation number P D of the DVDs shown in Figs. 3(e) and (f). This quantity is calculated in a similar way as the energy density and it gives the number of sites that are significantly participating in the dynamics of the DVD. In Figs. 3(e) and (f) we observe that although up to t ≈ 10 3 both DVDs exhibit approximately P D ≈ 20, for the case of the Hertzian chain [panel (e)] it starts to drop to a smaller value. As discussed earlier, the tendency of the DVD to start to localize is a precursor of a chaotic spot that may appear in the dynamics over a longer timescale.

It is interesting to note that, although the two models behave almost identically for H = 0.25, this energy corresponds for the Hertzian model to a initial displacement of u 21 (0) = 1.01 with the neighboring static overlaps being δ 21,22 ≈ 1.06. These values are far from the small amplitude approximation (u n /δ n,n+1 ≪ 1). The two models however show no differences (at least for the studied time scales), mainly due to the fact that practically only a sin- gle mode is participating in the dynamics.

B. Chaos and destruction of localization

Energy density evolution and chaos

In Fig. 4(a) we show the energy density evolution for the Hertzian model with energy H = 0.5. The energy distribution for both models [results for the FPUT are similar to Fig. 2(b)], is still localized for H = 0.5. However, there is a difference during the last decade, better captured by the evolution of P as illustrated in Fig. 4(c), since the Hertzian chain exhibits a tendency to increase the number of highly excited particles.

The most intriguing feature for this particular case is found in the system's chaoticity as quantified by the time evolution of Λ(t) shown in Fig. 4(d). The red solid line, which corresponds to the Hertzian chain with H = 0.5, deviates from the Λ(t) ∝ t -1 curve, at the last decade, and attains an almost constant value. This signals that the system is chaotic. In contrast, for the same energy the FPUT model's orbit remains regular. The Hertzian model therefore exhibits localized chaos whilst the FPUT model is localized and regular.

In Fig. 4(b) we show that initially localized wavepacket for the Hertzian model at H = 1.8, gradually spreads throughout the lattice signaling the destruction of Anderson localization. In particular, up to t ≈ 2 × 10 2 the wave-packet remains localized [see Figs. 4 (b) and (c)] with a participation number P < 3, it then rapidly spreads until t ≈ 4 × 10 3 . At the last decade the participation number saturates to a value P ≈ 26. This is the maximum observed value of P in all our simulations. According to the corresponding Λ(t) shown in Fig. 4(d for H = 1.8, the system also becomes chaotic as early as t ≈ 2 × 10 2 acquiring an almost constant positive value of Λ(t) ≈ 10 -3 . Results for the FPUT are not shown for this energy since excitations were still found to be localized and regular.

It is important to note here that for the particular single site excitation, all energies H > 1.8 lead to a final chaotic and delocalized state. This suggests the appearance of an energy threshold beyond which the final state of the Hertzian model is delocalized and chaotic. Below we show that this is true for different modes of this realization but also for different realizations.

Spatiotemporal evolution of chaos

In order to better understand the onset of chaos in the aforementioned cases, we study more closely the behavior of the DVDs. In Figs. 5(a) and(c) we plot the DVDs for the Hertzian model for H = 0.5 and H = 1.8 respectively. Focusing on the case of H = 0.5 we see that initially, when the system behaves regularly, the DVD exhibits an extended smoothed profile. This is more clearly seen by the red (dotted) curve in Figs. 5(b). Thereafter, during a period up to t ≈ 4 × 10 3 the DVD gradually converges around site n = 21. A profile of the DVD in this era is shown with the magenta curve in Fig. 5(b). Finally for the rest of the simulation the profile of the DVD is strongly localized around site n = 21 as also confirmed by two different profiles during the last decade shown in Fig. 5(b). Other recent studies (i.e. Refs. [41][42][43]) also used the DVD to spatially characterize chaos. In these works, it was found that the profile of the DVD exhibits a peak that oscillates within a chaotic region, while in our case it remains attached to a single site indicating strongly localized chaos.

At a higher energy of H = 1.8 the DVD initially exhibits an extended and smooth profile and an example is plotted in Fig. 5(d) with the red (dotted) line. Further on, it concentrates around a region close to the center of the chain and beyond this point the system is chaotic. The evolution of the DVD during this chaotic era, is characterized by one dominant peak along with other smaller peaks usually two orders of magnitude smaller (at most) as illustrated in Fig. 5(d). We particularly choose three cases where the dominant peak is either at the center (magenta), closer to the right edge (blue) or at the left edge (black). This is to emphasize the fact that for this energy, the chaoticity of the system is extended featuring strongly chaotic spots throughout the whole lattice.

Here we would like to stress the importance of the DVD which enables us to differentiate between localized and extended chaos.

Chaos and delocalization for the FPUT model

In contrast to the Hertzian model, the dynamics for the FPUT model appears to remain localized and regular up to an energy excitation of H = 1.8 [see Figs. 4 (c) and(d)]. The first energy at which the FPUT model's wave-packet is delocalized, exhibiting also a chaotic behavior, is around H = 2.9 (first row of Fig. 6). After an initial transient time for which the wave-packet remains localized, in the last decade of the simulation, it eventually spreads as shown by the energy density and P in Figs. 6(a) and (b) respectively. The time evolution of Λ(t) shown in panel (d) significantly deviates from the Λ(t) ∝ t -1 line, indicating chaotic dynamics for t ≳ 4 × 10 3 . Furthermore the DVD shown in Fig. 6(c) exhibits peaks at different places within the lattice when the system is chaotic, similarly to Fig. 5(c) for the Hertzian model, which is associated with extended chaos.

The most striking difference between the two models is found by examining higher energy excitations. To our surprise, we found that increasing the energy for the FPUT model does not necessarily lead to delocalization. In other words there isn't an energy threshold beyond which the final state of the FPUT lattice is delocalized. For example, as shown in the second row of Fig. 6; for H = 4 the excited wave-packet remains well localized and the participation number hardly changes [compare Figs. 6(a) and (b) with Figs. 6(e) and(f)]. This is somewhat a surprising result and it highlights the complexity of the phase-space of a disordered FPUT lattice. We could qualitatively describe the results for the FPUT model as alternating between spreading an localization as the energy increases. To better visualize this alternate behavior, an example for H = 8.7381 in the bottom row of Fig. 6 is shown, which exhibits a delocalized and chaotic final energy profile. For this energy, the system behaves qualitatively the same as in the first row with H = 2.9.

Regarding chaoticity, the dynamics of the DVDs are shown in the third column of Fig. 6. For all cases, the initially localized DVD around n = 21 finally departs from this site and it oscillates within the lattice. Accordingly, Λ(t), shown in the last column of Fig. 6, initially follows the regular orbit slope but eventually signals chaotic dynamics by diverging from this line and acauiring a nonzero value. We have found that for all H > 2.8, the final state of the lattice is always chaotic irrespectively of the localized or delocalized nature of the wave-packet.

C. Role of the non-smooth nonlinearity and energy equipartition

In order to further track down the mechanisms responsible for the different behaviors between the two models we monitor the appearance of the non-smooth nonlinearity [i.e. whenever (u n -u n-1 ) > δ n ] for the Hertzian model, or in other words the appearance of gaps. Fig. 7(a), shows the position of gaps for the case H = 0.5, which corresponds to the panels of the first row of Fig. 4. We clearly see that on the left and right side of site n = 21 a gap often opens during the system's evolution triggering the appearance of the non-smooth nonlinearity. At this energy no more than one gap is open at any instant as observed in Fig. 7(c) where the total number of gaps as a function of time is plotted. Importantly, since for this energy the dynamics of both the Hertzian and FPUT models is equivalent, but the Hertzian model appears to be chaotic, we identify the non-smooth nonlinearity around n = 21 as the ingredient which induces chaos for the Hertzian model.

For the energy H = 1.8 shown in Fig. 7(b), we find that more gaps start to open"moving" away from site n = 21, covering eventually the whole lattice. In fact, for the energy region 0.5 ≲ H ≲ 1.8 the wave-packet starts to delocalize (as quantified by P ) at the same time that additional gaps start to move away from site n = 21. For H = 1.8, as shown in Fig. 7(b), this happens around t ≈ 3 × 10 2 which is the same time that P [see Fig. 4(b)] starts to increase and the wave-packet starts to delocalize. These results, indicate a direct connection between the spreading of gaps within the lattice and the energy threshold beyond which the Hertzian model always traverses to delocalized and extended chaos.

To complete the comparison between the two models we also calculate the so called "spectral entropy" [44] by monitoring the corresponding normal modes. We write the weighted harmonic energy of the kth mode as v k = E k / ∑ N k=1 E k where E k is the kth mode's energy. We thus obtain the spectral entropy at time t as:

S(t) = - N ∑ k=1 v k (t) ln v k (t). ( 10
)
with 0 < S ≤ S max = ln N . It is however more convenient to use the normalized spectral entropy η(t) which can be written as,

η(t) = S(t) -S max S(0) -S max . ( 11
)
The value of η is normalized such that 0 ≤ η ≤ 1. With this normalization, when η remains close to one the dynamics does not substantially deviate from the initially excited modes. On the other hand as more modes are excited, η decreases towards zero. For a system at equipartition, a theoretical prediction for the mean entropy ⟨η⟩ exists, which assumes that the modes at equipartition follow a Gibbs distribution when the nonlinearity is weak.

The analytical form of the mean entropy ⟨η⟩ is given by [36,45] ⟨η⟩ = 1 -C ln N -S(0) (12) with C ≈ 0.5772 being the Euler constant.

In Fig. 8 we plot the time evolution of η and of the normal modes for different values of the energy H. As shown in Fig. 8(a), for H = 0.25 where the dynamics for both models is localized, the normalized entropy initially has a value of η = 1 and only slightly decreases from that value. This indicates that the dynamics is dominated by the single mode initially excited along with some weakly excited low frequency modes. This is also very clear in Fig. 8(e) where the time evolution of the weighted modes is shown. Initially only mode 34 is visible, and after some brief transient phase, a set of extended (low frequency modes) are slightly excited. In fact after t ≈ 10 2 the amplitude of each mode remains approximately constant and so does the time evolution of the normalized entropy η(t).

Similar behavior for the Hertzian model is observed at H = 0. the mean entropy at equipartition as given by Eq. ( 12). The asymptotic value of η(t) approaches the theoretically predicted value of ⟨η⟩ with H ≥ 1.8 as indicated in panels (c) and (d). The fact that the final stages of these simulations are close to an equipartition state is also supported by the mode energy distribution which clearly shows that at the last decade all modes appear to participate in the dynamics.

The existence of an energy threshold beyond which equipartition is reached for the Hertzian model depends neither on the particular mode nor the chosen realization shown in Fig. 1(b). To illustrate this, we first identified the 9 most localized modes of the distribution shown in Fig. 1(b). We then excite these modes by using a single site excitation around the point of localization of each mode with an energy H = 3. The choice of energy is to ensure that it is above the threshold for each mode. The results are shown in Fig. 9(a) and it is clear in all cases that the system finally reaches equipartition. We also performed simulations using different disorder realizations with α = 5, and exciting them at the central site n = 20 with energy H = 3. Since we always excite the same site but for different realizations, we may or may not excite a single localized mode. In any case, as it is shown in Fig. 9(b), the system reaches to equipartition in all cases.

The absence of an energy threshold leading to equipartition for the FPUT model is also shown in Fig. 10. F H = 0.25 [panels (a) and (e)] the behavior is the same as for the Hertzian model: η(t) saturates to a finite value close to 1 and a dominant mode along with some low frequency modes are present. For a much higher energy excitation of H = 2.9 shown in Figs. 10(b) and (f), from the early stages of the evolution more modes are excited and the entropy exhibits a plateau at η ≈ 0.7. Note that such a plateau is well known and studied in homogeneous FPUT chains and is associated with a metastable phase [36]. Beyond this point the entropy abruptly falls at t ≈ 5 × 10 3 and finally reaches a minimum value which is found to be close to the analytical result for equipartition given by Eq. ( 12). As shown in Fig. 10 (f) this is associated with the excitation of almost all linear modes.

For a larger initial energy H = 4, i.e. the case presented in the second row of Fig. 6, the dynamics of η is quite surprising. As shown in Fig. 10(c) the entropy saturates for most of the evolution around a relatively large value η ≈ 0.82. For the last two decades it starts to decrease, but with a very small slope. This is unexpected (also in accordance to the homogeneous FPUT studies e.g. Ref. [36]) since for higher energy excitations we anticipate to have a shorter plateau (than the one for H = 0.25) and the system to be driven faster towards equipartition. However, here the dynamics suggests that the contribution of modes other than mode 34 remains weak. This is also seen in Fig. 10(g) where not all modes have been excited at the end of the simulation, and in particular the highest frequency ones are still "mute". However, it is expected, that for larger timescales the system will reach equipartition, and η will eventually drop.

To highlight the alternate behavior found for the disorder FPUT model, in Fig. 10(d) we show the entropy for an even higher energy excitation of H = 8.7381 which corresponds to the results presented in the third row of Fig. 6. Similarly to the case of H = 2.9 the entropy saturates for a long time interval at a value η ≈ 0.8. Then at t ≈ 10 4 , η starts to drop and at the end of the simulation reaches a minimal value well captured by the analytical prediction of Eq. (12). Accordingly in Fig. 10 (h) we observe that as time increases more modes participate in the dynamics, and at the final stages of the simulation all modes are present.

IV. SUMMARY AND CONCLUSIONS

In this work, we numerically studied energy localization/delocalization and the chaoticity of two onedimensional disorder models: the Hertzian model fea-turing a non-smooth nonlinearity, and the FPUT model. The two models share the same linear limit and thus the same linear eigenmodes. Statistics on 1000 disorder realizations demonstrate that beyond a sufficient disorder strength, the linear chain acquires a significant amount of strongly localized modes. Focusing on a single realization from the aforementioned ensemble, we show that the evolution of such a mode can be characterized by three different scenarios: (i) localization with no chaos; (ii) localization and chaos; (iii) spreading of energy, chaos and equipartition.

In particular, for sufficiently small energies the two models behave quantitatively similar, with excitations remaining localized and non chaotic, at least for the time scales of our simulations. For larger energy values, a transient energy region is found for which the Hertzian model exhibits localized but chaotic behavior. After an energy threshold, associated to the spreading of gaps in the lattice, the Hertzian model evolves into an equipartition, chaotic state independent of the particular value of the initial energy. The appearance of such a threshold is confirmed for other modes of the particular disorder realization but also for different disorder realizations.

On the other hand the dynamics of the FPUT model is substantially different from that of the Herztian model. Firstly, delocalization and chaos emerge for higher energies for the FPUT model. We find strong numerical evidences that this difference is attributed to the nonsmooth nonlinearity which is present only in the Hertzian model. Furthermore, for higher energy values, the FPUT system shows an alternating behavior between chaotic localized and chaotic extended dynamics lacking a par-ticular threshold beyond which equipartition is always reached. We can therefore conclude that, in contrast to the Hertzian model, the final state of a strongly disordered FPUT lattice under single site excitation, strongly depends on both the disorder realization and the initial excitation energy.

Our results provide further insights into the chaotic dynamics of strongly disordered chains. Using additional chaos indicating tools such as the deviation vector densities, we are able to clearly separate localized (in space) from extended chaotic behavior. In addition we show that non-smooth nonlinearities do not only induce the destruction of Anderson localization but also provide a mechanism to drive the system into equipartition. An interesting direction stemming from our results is to pursue a thorough statistical analysis in order to probe the interplay between disorder and nonlinearity and the resultant effect on the corresponding time scales for equipartition.

FIG. 1 .

 1 FIG. 1. (a) Mean (over 1000 disorder realizations) participation number ⟨P ⟩ of the eigenmodes for varying disorder strengths α, sorted in descending order k for each realization. The standard deviation at each point is shown by the error bars. (b) The eigenfrequencies of a particular disordered chain of 40 sites for α = 5 sorted by increasing frequency. The insert shows the profile of the 34 th mode.

FIG. 2 .

 2 FIG. 2. (a) and (b) The spatiotemporal evolution of the energy distribution for the Hertzian and FPUT chains respectively for H = 0.25. The black curves indicate the running mean position of the energy distributions. The color bars on the right sides of (a) and (b) are in logarithmic scale. (c) The locally weighted smoothed values of P as a function of time for the Hertzian chain (red curve) and the FPUT chain (blue curve). (d) The time evolution of Λ(t) for the Hertzian chain (red curve) and the FPUT chain (blue curve). Both lines practically overlap and the dashed line indicates the law Λ(t) ∝ t -1 .

  ) and (b), exhibits an extended

FIG. 3 .

 3 FIG. 3. (a) [(b)]: The spatiotemporal evolution of the deviation vector density (DVD) for the Hertzian [FPUT] disordered chain. The color bars on the right sides of (a) and (b) are in logarithmic scale. (c) [(d)]: Deviation vector profiles for three time instances of t ≈ 10 1 indicated by the blue (b) curve, t ≈ 10 3 indicated by the red (r) curve and t ≈ 10 5 indicated by the black (bl) curve. These times correspond respectively to the blue, red and black horizontal lines in panel (a) [b]. (e) [(f)]: The time evolution of the participation number PD of the DVD for the Hertzian [ FPUT] model. All results are obtained for H = 0.25.

FIG. 4 .

 4 FIG. 4. Panels (a) and (b) show the spatiotemporal evolution of the energy distribution for the Hertzian model with H = 0.5 and H = 1.8 respectively. Black curves indicate the running mean position of the energy distributions. The color [41-43]bars on the right sides of (a), (b) are in logarithmic scale. Panels (c) and (d) are the same as Figs. 2(c) and (d); for the Hertzian model H = 0.5 and H = 1.8 and for the FPUT model with H = 1.8.

)FIG. 5 .

 5 FIG. 5. Panel (a) shows the spatiotemporal evolution of the DVD for the Hertzian model at H = 0.5 whilst panel (b) shows the profiles of the DVDs at t ≈ 1.7 × 10 1 red (r) curve, t ≈ 1.7 × 10 4 magenta (m) curve and t ≈ 8.2 × 10 4 blue (b) curve. (c) Same as (a) but for H = 1.8. (d) Same as (b) but for H = 1.8 at t ≈ 1.7 × 10 1 red (r) curve, t ≈ 4.9 × 10 3 magenta (m) curve, t ≈ 3.5 × 10 4 blue (b) curve and t ≈ 4.8 × 10 4 black (bl) curve. The color bars in (a) and (c) are in logarithmic scale. .

FIG. 6 .

 6 FIG. 6. (a), (b), (c) and (d) depict the energy density, P , DVD, PD and Λ(t) respectively for the FPUT with H = 2.9. The second, and third rows correspond to energies H = 4 and H = 8.7381 respectively. The color bars on the right sides of panels (a), (c), (e), (g), (i) and (k) are in logarithmic scale.

FIG. 7 .

 7 FIG. 7. The spatiotemporal evolution of the gaps in the Hertzian model for energies H = 0.5 (a) and H = 1.8 (b). The yellow (lighter) color corresponds to the lattice points where (un(t) -un-1(t)) > δn. The instantaneous total number of gaps for the Hertzian model for energies H = 0.5 (c) and H = 1.8 (d).

FIG. 8 .

 8 FIG. 8. Top row: The time evolution of the normalized spectral entropy η(t) for the Hertzian model. The dashed horizontal line in panels (c) and (d) show the mean value ⟨η⟩ given by Eq. (12). Bottom row: The evolution of the weighted harmonic energy of eigenmodes as a function of time. The modes are sorted by increasing frequency [c.f. Fig. 1(b)]. The values of the energy are H = 0.25 (a)-(e), H = 0.5 (b)-(f), H = 1.8 (c)-(g), and H = 3 (d)-(h). The color bars on the right sides of panels (e)-(h) are in logarithmic scale.

FIG. 9 .

 9 FIG. 9. (a) Temporal evolution of η(t) for the 9 most localized modes of the distribution corresponding to Fig. 1(b). (b) Temporal evolution of η(t) obtained by exciting site n = 20 of 9 different disordered realizations with α = 5. In both panels, the results correspond to the Hertzian model with an energy H = 3. The dashed horizontal lines show the average (on the different initial conditions), mean entropy at equipartition ⟨η⟩.

  FIG. 10. Same as in Fig. 8 but for the FPUT model. The dashed horizontal line in panels (b) and (d) show the mean value ⟨η⟩ given by Eq. (12). The values of the energy in this case are H = 0.25 (a)-(e), H = 2.9 (b)-(f), H = 4 (c)-(g), and H = 8.7381 (d)-(h).

ACKNOWLEDGMENTS

Ch. S. was supported by the National Research Foundation of South Africa (Incentive Funding for Rated Researchers, IFFR and Competitive Programme for Rated Researchers, CPRR) and thanks LAUM for its hospitality during his visits when part of this work was carried out. We also thank the Center for High Performance Computing (https://www.chpc.ac.za) for providing computational resources for performing part of this paper's computations.