Analyzing Chaos in Higher Order Disordered Quartic-Sextic Klein-Gordon Lattices Using $q$-Statistics

TL;DR

This paper investigates the chaotic behavior of disordered Klein-Gordon lattices with quartic and sextic anharmonicities, showing that the system remains strongly chaotic over long times even with higher-order terms.

Contribution

It extends previous work by demonstrating that adding a sextic term does not diminish chaos, confirming persistent strong chaos in higher-order disordered Klein-Gordon lattices.

Findings

$q$-Gaussian distributions approach Gaussian in long time limit

System remains strongly chaotic for times up to 10^9

Chaos persists despite higher-order anharmonicity

Abstract

In the study of subdiffusive wave-packet spreading in disordered Klein-Gordon (KG) nonlinear lattices, a central open question is whether the motion continues to be chaotic despite decreasing densities, or tends to become quasi-periodic as nonlinear terms become negligible. In a recent study of such KG particle chains with quartic (4th order) anharmonicity in the on-site potential, it was shown that $q-$Gaussian probability distribution functions of sums of position observables with $q > 1$ always approach pure Gaussians ($q=1$) in the long time limit and hence the motion of the full system is ultimately "strongly chaotic". In the present paper, we show that these results continue to hold even when a sextic (6th order) term is gradually added to the potential and ultimately prevails over the 4th order anharmonicity, despite expectations that the dynamics is more "regular", at least in the regime of small oscillations. Analyzing this system in the subdiffusive energy domain using $q$-statistics, we demonstrate that groups of oscillators centered around the initially excited one (as well as the full chain) possess strongly chaotic dynamics and are thus far from any quasi-periodic torus, for times as long as $t=10^9$.

Equations

$$H=\sum _{l=1}^N\frac{p_l^2}{2}+\frac{{\widetilde{ϵ}}_l}{2}{x}_l^2+\frac{|x_l{|}^{\sigma +2}}{\sigma +2}+\frac{1}{2W}(x_{l+1}-x_l{)}^2,$$

$$z_l=\frac{E_l}{{\sum }_{i=1}^N{E}_i}$$

$$m_2=\sum _{l=1}^N(l-\bar{l}{)}^2{z}_l,\bar{l}=\sum _{l=1}^{N}l{z}_l$$

$$P=\frac{1}{{\sum }_{l=1}^N{z}_l^2},$$

$$\alpha =\frac{1}{\sigma +1}$$

$$S_M^{(j)}=\sum _{i=1}^M\eta (t_i{)}^{(j)},j=1,\dots ,N_{\text{mbox}ic},$$

$$s_M^{(j)}\equiv \frac{S_M^{(j)}-⟨S_M^{(j)}⟩}{{\sigma }_M},$$

$$\mathtt{P}(s_{M}^{(j)})=a\exp_{q}({-\beta s_{M}^{(j)2}})\equiv a\biggl{(}1-(1-q)\beta s_{M}^{(j)2}\biggr{)}^{\frac{1}{1-q}},$$

$$\beta=a\sqrt{\pi}\frac{\Gamma\Bigl{(}\frac{3-q}{2(q-1)}\Bigr{)}}{(q-1)^{\frac{1}{2}}\Gamma\Bigl{(}\frac{1}{q-1}\Bigr{)}},$$

$$H=\sum _{l=1}^N\frac{p_l^2}{2}+\frac{{\widetilde{ϵ}}_l}{2}{x}_l^2+A\frac{x_l^4}{4}+(1-A)\frac{x_l^6}{6}+\frac{1}{2W}(x_{l+1}-x_l{)}^2,$$

Full text

1 Introduction

Anderson localization [1], i.e. the absence of wave diffusion in disordered media, is a phenomenon affecting many physical processes, such as for example the conductivity of materials, the behavior of granular matter, the dynamics of Bose-Einstein condensates, etc. The effect of nonlinearity on the localization properties of wave packets in disordered systems has attracted the attention of many researchers, both experimentally [31, 19, 30] and theoretically [28, 18, 12, 36, 39, 13, 33, 10, 21, 40, 24, 23, 27, 8, 7, 4, 9, 25, 20, 26, 35, 2, 22, 16, 11]. Recent studies of nonlinear disordered variants of two typical one-dimensional Hamiltonian lattice models, namely the Klein-Gordon (KG) oscillator chain and the discrete nonlinear Schrödinger equation, revealed the statistical characteristics of energy spreading and showed that nonlinearity destroys localization [12, 36, 33, 21, 8, 7]. In these papers, the existence of different dynamical behaviors in different energy density regimes was established, their particular dynamical characteristics were determined and their appearance was theoretically explained.

However, important questions regarding the asymptotic behavior of wave-packet spreading, and the persistence of chaos in such systems, still remain unanswered. Some researchers have conjectured [17, 3] that wave packets will eventually approach torus-like structures in phase space, exhibiting at the same time less chaotic behavior which eventually leads to the halt of energy spreading in the chain. Although most numerical investigations on nonlinear disordered lattices show that wave packets continue spreading chaotically, at least up to times accessible to computer simulations, some numerical indications of a possible slowing down of spreading for particular models have been reported in the literature [27, 26]. The diversity of the models studied so far includes systems with different numbers $N_{IoM}$ of integrals of motion or conserved quantities, and the degree of anharmonicity $\sigma$ which is related to a corresponding $n$-body interaction or equally to the number of interacting normal modes mediated by the anharmonicity. $N_{IoM}=1$ for the mentioned case of KG lattices, while $N_{IoM}=2$ for the DNLS case. In most studied cases $\sigma=2$ which corresponds to quartic anharmonicity and two-body interactions. Here we will also study $\sigma=4$ which corresponds to sextic anharmonicity and three-body interactions.

To investigate the potential asymptotic approach to regular (or irregular) dynamics, the properties of the motion in the subdiffusive regime of the KG model with $\sigma=2$ were recently studied [35, 2]. In [35], the computation of the maximum Lyapunov exponent (MLE) showed that although chaotic dynamics slows down as expected from a subdiffusive process, it does not cross-over to a regime of regular behavior. In [2] the dynamics of a disordered quartic (termed KG4) lattice was studied using $q$-statistics [38] by analyzing probability distribution functions (pdfs) of position observables describing the evolution of wave packets initiated by exciting the central lattice particle. In that work it was shown that the overall motion displays strongly chaotic behavior for long times, again with no signs of the dynamics relaxing on quasi-periodic tori.

In this paper, we apply the above methodology followed in [2, 35] to a KG model which includes a gradually increasing sextic anharmonicity which is assumed to show an asymptotic slowing down of the wave packet spreading - if existing - at earlier times, since the impact of this anharmonicity is weaker than the quartic one, as long as small densities are considered. In line with what we discovered in the KG4 model [2], both the full lattice, as well as groups of particles around the initially excited one at the center remain strongly chaotic and show no signs of approaching quasi-periodicity, even after very long integration times (at least up to $t=10^{9}$, with time scales dictated by the linear equations being of order one). On the other hand, individual particles far from the center, after interacting with the wave-packet behave at first weakly chaotically, but later also tend to strong chaos for times as long as $t=10^{9}$. We find that the closer the particles are to the center of the lattice, the more strongly chaotic their behavior is (with pdfs closer to the Gaussian $q=1$ case), and that the dynamics of the whole lattice is always strongly chaotic, with pdfs obeying Boltzmann-Gibbs thermostatistics. Thus, all our results indicate that the wave packet spreading is a truly chaotic process which does not show any tendency to become more regular. Our results are also supported by the computation of the time dependence of the largest Lyapunov exponent which again show no cross-over into a regime of regular behavior, similar to the KG4 case [35].

The paper is structured as follows: In Sec. 2 we present the KG model and outline the used statistical methods in the spirit of the Central Limit Theorem (CLT). In Sec. 3 we examine the mixed case of both quartic and sextic anharmonicities in the potential and consider the statistical properties of the dynamics as the sextic terms become increasingly more important. In Sec. 4, we focus on the purely sextic anharmonicity model and describe the results obtained when we excite only the central particle of a $500$ particle chain and monitor the chaotic evolution of individual particles, groups of particles about the central one, as well as the whole system. Our conclusions and a discussion follow in Sec. 5.

2 Model and numerical methods

The Hamiltonian of the disordered one-dimensional KG lattice studied in [33] is

[TABLE]

where $l$ is the lattice site index, $x_{l}$ and $p_{l}$ are respectively the generalized canonically conjugated coordinates and momenta (with $x_{N+1}=0$), $\sigma$ measures the degree of anharmonicity and $W=4$ controls the nearest neighbour interaction strength and thus the effective strength of disorder. Disorder enters through the onsite harmonic squared frequencies $\tilde{\epsilon}_{l}$ which are random uncorrelated numbers chosen uniformly from the interval $\left[\frac{1}{2},\frac{3}{2}\right]$. The total energy $E\equiv H\geq 0$ of the system serves as a control parameter of the nonlinearity for fixed disorder strength $W$. The case $\sigma=2$ corresponds to the typical quartic disordered KG4 model. Wave packet spreading in the KG4 case was studied in several papers [12, 36, 21, 8, 7, 35, 2]. The equations of motion follow as $\dot{x}_{l}=\partial H/\partial p_{l}$ and $\dot{p}_{l}=-\partial H/\partial x_{l}$.

The dynamics of wave packet spreading in the Hamiltonian (1) was analyzed in detail in [33] for several anharmonicity values $\sigma$, following the evolution of the normalized energy density

[TABLE]

of the site-energies $E_{l}=\frac{p_{l}^{2}}{2}+\frac{\tilde{\epsilon}_{l}}{2}x_{l}^{2}+\frac{|x_{l}|^{\sigma+2}}{\sigma+2}+\frac{1}{4W}\left[(x_{l+1}-x_{l})^{2}+(x_{l-1}-x_{l})^{2}\right]$. The main goal was to monitor the evolution of the second moment $m_{2}$

[TABLE]

which quantifies the wave packet degree of spreading, and the participation number

[TABLE]

which measures the number of the most strongly excited sites in the system. In [33] it was found that the wave packet spreads incoherently and subdiffusively with $m_{2}\propto t^{\alpha}$ and $P\propto t^{\alpha/2}$ with the exponent

[TABLE]

up to the largest computed times, indicating no loss of incoherent (chaotic) dynamics and no cross-over to coherent (regular) dynamics.

In the present paper we consider a hybrid model (termed KG46 here) which interpolates between the quartic $\sigma=2$ and sextic $\sigma=4$ cases by including sextic terms in (1). We expect that sextic terms are generating weaker nonlinear terms in the equations of motion, and could amplify a crossover to regular dynamics, if present. In all our models we follow the evolution of single-site excitations of the central particle in the subdiffusive regime defined in [33] (see Fig. 1 in [33]), by considering $N=500$ sites and setting the total energy of our lattice to some constant energy $E$. In particular, for a given set of $\tilde{\epsilon}$ values, we choose $x_{l}=0$ and $p_{l}=\sqrt{2E}\delta_{l,N/2}$, where $\delta_{i,j}$ is the Kronecker delta, thus exciting precisely one oscillator in the center of the system which has $N$ sites.

To study the resulting trajectories, we integrate numerically the equations of motion of Hamiltonian (1) by the 4th order Yoshida’s symplectic integrator [41]. In our simulations, we set the integration time step to $\tau=0.05$, which typically keeps the relative energy error at about $10^{-6}$. Furthermore, to obtain reliable statistical results that are independent of the particular realizations, we consider an ensemble of 64 disorder realizations, i.e. 64 random sequences of $\tilde{\epsilon_{l}}$ values in (1). Apart from the computation of $m_{2}$ and $P$, we also evaluate the MLE $\lambda_{1}$. For this purpose, we use the same symplectic integration scheme for the integration of the variational equations of system (1) according to the tangent map method [34, 15, 14]. The variational equations govern the evolution of small deviation vectors from the studied trajectory and are used for the evaluation of the MLE according to the so-called standard method [5, 6, 32].

We use the solutions of the equations of motion of Hamiltonian (1) to construct pdfs of suitably rescaled sums of $M$ values of a generic observable $\eta(t_{i}),\;i=1,\ldots,M$, which depends linearly on the position coordinates of the solution. Viewing these as independent and identically distributed (iid) random variables (in the limit of $M\rightarrow\infty$), we evaluate their sum

[TABLE]

for $N_{\mbox{ic}}$ initial conditions and study the statistics of quantities (6) centered about their mean value $\langle S_{M}^{(j)}\rangle$ and rescaled by their standard deviation

[TABLE]

where $\sigma_{M}^{2}=\langle S_{M}^{(j)2}\rangle-\langle S_{M}^{(j)}\rangle^{2}$. Plotting the normalized histogram of the probabilities $\mathtt{P}(s_{M}^{(j)})$ as a function of $s_{M}^{(j)}$, we then compare the resulting numerically computed pdfs with a $q$-Gaussian of the form

[TABLE]

where $q$ is the so-called entropic index, $\beta$ is an arbitrary parameter and $a$ a normalization constant. Equation (8) is a generalization of the well-known Gaussian pdf, since in the limit $q\rightarrow 1$ we have $\exp_{q}(-\beta x^{2})\rightarrow\exp(-\beta x^{2})$. Moreover, it can be shown that the $q$-Gaussian distribution is normalized when

[TABLE]

where $\Gamma$ is the Euler $\Gamma$ function. Clearly, equation (9) shows that the allowed values of $q$ are $1<q<3$ for the normalization to be possible.

The index $q$ appearing in (8) is connected with the Tsallis entropy [38]. Systems characterized by the Tsallis entropy are said to lie at the “edge of chaos” and are significantly different from Boltzmann-Gibbs systems, in the sense that their entropy is non-additive and generally non-extensive [38].

Let us now consider the values of one (or a linear combination) of coordinates $x$ of the solution of system (1) at discrete times $t_{i},\;i=1,\ldots,M$, as realizations of $N_{\mbox{ic}}$ random variables $X^{(j)}(t_{i})$, $j=1,\ldots,N_{\mbox{ic}}$ and study in detail their statistics. Typically, the results presented in this work are obtained by setting $N_{\mbox{ic}}=10^{6}$ and $M=10^{3}$. According to the CLT [29], if these variables are random, their sum-distributions yield a Gaussian pdf, whose mean and variance are those of the $X^{(j)}$’s. This happens for many dynamical systems in strongly chaotic regimes, where correlations decay exponentially and the system obeys Boltzmann-Gibbs statistics. However, in weakly chaotic regimes, these pdfs do not rapidly converge to a Gaussian, but instead are approximated for long times, by a $q$-Gaussian distribution (8) with $1\leq q<3$.

3 Results for the hybrid quartic-sextic (KG46) case

Let us consider a Hamiltonian that includes quartic as well as sextic order terms in the on-site potential, i.e.

[TABLE]

where $0\leq A\leq 1$. Note that the introduction of the parameter $A$ in (10) allows us to investigate a family of Hamiltonians that ranges from the KG4 ($A=1$) to the KG6 ($A=0$) system.

Following [2], we have used energy values $E=0.4$ and $E=0.6$ in (10), which corresponds to the subdiffusive regime of the KG4 system. Although we studied the behavior of system (10) for various values of parameter $A$, we present in Fig. 1 illustrative results for $A=0.5$ only. Figure 1 shows the time evolution of the entropic index averaged over 64 realizations for four indicative cases: $\langle q_{21}\rangle$ corresponding to 21 particles about the center, $\langle q_{500}\rangle$ for the full chain, and $\langle q(200)\rangle$ and $\langle q(300)\rangle$ for the 200th and 300th particles. Here $\langle...\rangle$ indicates averaging of the quantity over all 64 disorder realizations. We first plot these indices in Fig. 1(a)-(d) for the hybrid quartic-sextic KG46 case $A=0.5$ and then repeat the calculation for the same indices in the purely sextic KG6 case $A=0$ in Fig. 1(e)-(h). The only difference is that particle groups have indices that rapidly converge to $q=1$, while individual particles behave weakly chaotically with indices significantly higher than unity, but also show a tendency to approach $q=1$ for very long times. From these results we conclude that the dynamics of the chain is strongly chaotic as the $q$-values are close to 1, especially when large parts of the lattice are considered, just as was observed in the purely quartic KG4 case studied in [2].

Now, for any non-zero values of $A$ it is expected that the asymptotic behavior of the dynamics will be controlled by the quartic terms of the Hamiltonian (or equivalently by the cubic terms in the corresponding equations of motion), since the contribution of the sixth order part will eventually become negligible as the wave packet spreads and the values of $x_{l}$ decrease in magnitude. This is indeed what we observe for $A=0.5$ in Fig. 1(a),(b) since the overall dynamics is similar to the behavior of the KG4 model studied in [2]. This turns out to be more generally true as the results for several other non-zero values of $A$ (not presented here) are qualitatively similar to the ones shown in Fig. 1(a)-(d). For this reason, we focus our attention in the next section on the $A=0$ case.

4 Results for the purely sextic KG6 case

We set $A=0$ in (10), and present in Figs. 2 and 3 the subdiffusive indicators of single-site excitations of the central particle averaged over 64 disorder realizations. First, in Fig. 2(a),(b) we show the time dependence of the second moment $\langle m_{2}\rangle$ and participation number $\langle P\rangle$ respectively, shading by grey one standard deviation away from the mean. Here $\langle...\rangle$ indicates averaging of the logarithm of the quantity over all 64 disorder realizations. In Fig. 3(a) we show the MLE $\langle\lambda_{1}\rangle$, while in Fig. 3(b) we depict the normalized energy $z_{l}$ of all particles $l=1,\ldots,500$ for one of the 64 realizations, at $t=10^{7}$. As we see in Fig. 3(a), $\lambda_{1}$ first decays as $1/t$ during an initial transient phase, and after that continues to decay with a much weaker exponent as $1/t^{0.3}$, which indicates that the motion remains chaotic in time (since regular motion would correspond to $1/t$ [32]). This behavior is similar to the one observed in [35] for the KG4 model, with a slightly different decay exponent $1/t^{0.25}$. The results in Fig. 2(b) are in good agreement with the prediction of Eq. (5). In Fig.  3(b), on the other hand, it is evident that the particles at both edges of the chain are not excited at this time. We have also been careful to check that the normalized energies $z_{l}$ for particles close to the edges are always smaller than $10^{-10}$, at least up to $t=10^{9}$, to guarantee that the wave packet did not reach the edges of the chain during the integration times.

Let us now turn to the statistical properties of the KG6 lattice as the initial excitation of the central particle diffuses to its neighboring sites. In Fig. 1(e),(f) we present for $A=0$ the time evolution of the same two entropic indices as in Fig. 1(a),(b) for $A=0.5$, implying again that any collection of particles around the central one, including the full lattice, show strongly chaotic behavior as their averaged entropic indices are practically equal to one. On the other hand, Fig. 1(g),(h) show that single particle statistics is analogous to what one finds for $A=0.5$ in Fig. 1(c),(d). These results suggest that the wave packet spreading dynamics in the subdiffusive regime of the KG6 system also does not approach a KAM regime of invariant tori. This is in agreement with [2, 33, 35], where similar behaviors were found for KG4 systems, and with Sec. 3 where the hybrid KG46 models were considered, and thus supports the conclusions that the overall motion in these systems never approaches a KAM regime of quasi-periodic motion.

As an illustrative example of how the corresponding pdfs look like, we present in Fig. 4(a), for $A=0$, a pdf with $q$ slightly smaller than 1. This distribution refers to a group of 121 particles symmetrically located around the particle $i=250$ for one disorder realization of the KG6 model with $E=0.6$. The pdf is constructed from values of the observable $\eta_{121}=\sum_{i=190}^{310}x_{i}$ for $t\mbox{ in }[9\cdot 10^{6},10^{7}]$. The analysis yields $q_{121}=0.95\pm 0.017$, which is very close to the $q=1$ case of Boltzmann-Gibbs strong chaos. Entirely analogous results are obtained for the hybrid KG chain with $A=0.5$ in Fig. 4(b) for 21 particles around the central one.

The significance of $q<1$ values is as yet unclear. Thus, we have performed an analysis with different summands for the calculation of the pdfs as a function of time and have found they exhibit smaller than 1 entropic indices for long time intervals. This suggests that $q<1$ may not be due to numerical errors as it is persistent for long enough times, but does not entirely exclude this possibility (e.g. finite time and size effects). Conversely, it may be due to a deeper reason related to a similar result reported recently in [37], where large systems of coupled logistic maps near their chaotic threshold were also found to exhibit pdfs with $q<1$ when subjected to a different type of randomness, namely that of additive noise in the equations of motion. Consequently, we believe that more work on this fascinating behaviour should be devoted to identify the underlying reason.

5 Conclusions

In this paper we have analyzed the chaotic properties of wave-packet spreading in KG disordered lattices with higher than 4th order nonlinearity in their on-site potential, using the methodology of $q$-statistics and non-extensive Statistical Mechanics [38]. Our aim was to provide further evidence to strengthen the conclusions of earlier investigations [2, 35] that wave packet spreading remains strongly chaotic for longer and longer integration times. Studying the entropic $q$ index of suitably chosen probability density functions, we showed numerically that the lattice in parts, but also as a whole, behaves strongly chaotically, with $q\rightarrow 1$ as time increases to $t=10^{9}$. This suggests that the overall motion in these systems does not approach the quasi-periodic regime of invariant tori as conjectured by some authors, due to progressively smaller nonlinear effects at each lattice site as time goes to infinity.

To further establish this conclusion, we examined a class of on-site disordered potentials with quartic and/or sextic anharmonicities, and analyzed wave packet dynamics generated by exciting the central particle and spreading subdiffusively over a 500-site one-dimensional chain. We found that all these KG systems remain strongly chaotic, in parts and as a whole, showing no sign of approaching quasi-periodic behavior, even after times as long as $t=10^{9}$. We also verified that individual particles close to the center of the lattice, after interacting with the wave, exhibit first weakly chaotic motion with $q$ values significantly larger than unity, but eventually also show a tendency towards strong chaos, just as was reported in [2] for a KG4 model with purely quartic anharmonicity. Focusing on the strictly sextic KG6 model, which is expected to be closer to linear dynamics for small energy densities, we noted that here also strongly chaotic behavior prevails as time becomes very large.

In conclusion, the indisputable fact is that when we consider the full lattice or groups of particles around its center, the long-time dynamics is always strongly chaotic with entropic index $q=1$. Interestingly, in the purely sextic KG6 case, as well in the hybrid KG46 lattice, there are particle groups that attain $q$ values slightly less than one (see Fig. 4)! This is reminiscent of a similar result in a recently published $q$-statistical analysis of coupled weakly chaotic logistic maps in the presence of noise [37]. Of course, noise in [37] enters additively in the evolution equations, while in our system disorder is spatial and present only in the parameters. Thus, the occurrence of $q<1$ in such models requires further study to shed more light on the different types of chaotic behavior occurring in these multi-dimensional Hamiltonian systems.

6 Acknowledgements

Ch. S. acknowledges support by the National Research Foundation of South Africa (Incentive Funding for Rated Researchers, IFRR and Competitive Programme for Rated Researchers, CPRR). This work was supported by the Institute for Basic Science of South Korea, project IBS-R024-D1. We are grateful to Professor L. B. Drossos for many useful discussions and for helping us use the computational facilities of his HPCS-DL Lab at the Technological Educational Institute of Western Greece.