© 2018 al-Farabi Kazakh National University Printed in Kazakhstan IRSTI 27.29.23

A. Bountis

Nazarbayev University, Department of Mathematics, School of Science and Technology,

Astana, Kazakhstan, e-mail: anastasios.bountis@nu.edu.kz

**Complex Dynamics and Statistics in Hamiltonian **
**1-Dimensional Lattices **

**Abstract. **In this paper I review in a brief and introductory way some important developments in the
analysis of the dynamics and statistics of N – dimensional Hamiltonian systems, in which my research
team and I have played an important role over the last two decades. The results I describe here have
helped us understand the surprising importance of simple periodic orbits and their local stability
properties in revealing crucial dynamical and statistical properties of the systems as a whole. This has led
us to introduce the concepts of “strong” and “weak” chaos that are expected to play a significant role in
better understanding the complexity of these multi-dimensional systems, which have important
applications in solid state, field theory, superconductivity and more recently nonlinear optics.

**Key words:** N – dimensional Hamiltonian systems, superconductivity, nonlinear optics.

**Introduction **

Let us consider the 2- degree of freedom Hamiltonian system:

� = ��+ ���=

=^{�}_{�}��_{�}^{�}+ �_{�}^{�}) +^{�}_{�}��^{�}+ �^{�}) + ��^{�}�^{�} (1)

Its solutions for ε = 0 (the uncoupled case), plotted as intersection points on a Poincare surface of section of the 4 – dimensional space, are shown as a family of smooth closed curves in the graph on Figure 1.

**Figure 1** – A surface of section plot on the plane x(t k), px(tk), (every time y(tk)=0).

All orbits are obtained for the same constant value of the Hamiltonian, H=E)

International Journal of Mathematics and Physics 9, №2, 21 (2018) It is important to note that in this integrable

(separable) case, all solutions correspond to periodic and quasiperiodic orbits with two frequencies ω1

and ω2 belonging to the two uncoupled oscillators in (1). Thus, if ω1 and ω2 are rationally related, i.e. ω1 / ω2 = m/n (m,n positive integers) the corresponding closed curves in Figure 1 are filled with periodic orbits, while if ω1 / ω2 = irrational, the curves are filled by a single initial condition by an orbit that never closes and is called a quasiperiodic solution.

Let us now make ε > 0, e.g. ε=0.02. What we discover is remarkable: While many smooth closed curves corresponding to ω1 / ω2 = irrational survive, those that correspond to ω1 / ω2 = rational have disappeared giving their place to chains of islands having a stable periodic orbit at their center and a chain of “saddle – like” regions where two pointed regions meet at an unstable periodic orbit of the same period, see Figure 2 below. Now if we magnify the region close to one of those unstable

periodic points shown in Figure 2(a) by an arrow, an amazingly complex network of islands, as well as small scale chaotic region that we will later identify with “weak chaos”, are revealed in Figure 2(b) below [1].

If we now further increase the value of ε to ε = 0.2, we observe in Figure 3, that the small scale chaotic regimes observed near unstable periodic orbits (saddle points) of Figure 2 have now grown considerably into domains that we will refer to as domains of “strong chaos” later in this paper.

Although we don’t show it here, as one can imagine, increasing further the value of ε, the islands of stable periodic motion will diminish in size, while the strongly chaotic regimes will further increase, showing a tendency to occupy most of the available phase space. The same effect will occur, if we fix the value of ε > 0 and start increasing the total energy E at which the surfaces of section are computed.

(a) (b)

**Figure 2 (**a) – The surface of section of Figure 1 showing orbit intersections in the plane x(tk), px(tk),
for ε=0.02. (b) A magnification of the region shown by the arrow in (a), where,

besides the chains of small islands having stable periodic orbits at their center, one observes a region of randomly scattered points which constitute a “weakly chaotic’ domain.

**Figure 3** – The surface of section of the orbits for ε=0.2 and the same energy value
E as in Figure 2. Note how the weakly chaotic domains of Figure 2 have grown

to a much greter size forming regions that we will later

International Journal of Mathematics and Physics 9, №2, 21 (2018)
**Simple periodic orbits, weak and strong **

**chaos **

We study Hamiltonian dynamical systems of* N*
degrees of freedom (dof), in an 2N–dimensional
phase space of position and momentum coordinates,
whose equations of motion are written in the form

���

�� =_{��}^{��}

�,^{��}_{��}^{�} =_{��}^{��}

�, k = 1,2, . . . N, (2)
where *H* is the Hamiltonian function. For more
details on the results that follow in the present paper
the reader is invited to consult [2].

If *H* does not explicitly depend on the time *t*, it
represents a first integral, whose value gives the
total energy of the system *E*. We will assume that
the Hamiltonian can be expanded in *power series as *
*a sum of homogeneous polynomials* of degree *m* ≥2,
so that the origin is a stable equilibrium point of the
system:

� = �_{�}(�_{�}, � , �_{�}, �_{�}, � , �_{�}) +

+��(��, . . . , ��)+. . . = �. (3)
We now assume that *H**m* = 0 for all m > 2 and
that the linear equations resulting from (2) and (3),
yield a matrix, whose eigenvalues all occur in
conjugate imaginary pairs, ±iωq, and thus provide
the frequencies of the so-called normal mode
oscillations of the linearized system.

�_{�}= � �_{�}

�

���

= �, �_{�} =1

2��_{�}^{�}+ �_{�}^{�}�_{�}^{�} �,

� = 1,2, . . . , �, (4) where Pq, Qq are the normal mode coordinates.

Then, according to a famous theorem by Lyapunov, if none of the ratios of these eigenvalues, ωj/ωk is rational, for any j, k = 1,2,…,N, j ≠ k, all linear normal modes continue to exist as periodic solutions of the nonlinear system.

If the frequencies for *H**m*≠0 are close to those of
the linear modes, the continuation of the linear
modes are examples of simple periodic orbits
(SPOs) of the nonlinear system, where all variables

oscillate with the same frequency. We shall mention below the spectrum of Lyapunov exponents, and will discuss how its properties are connected to the emergence of strongly (large scale) chaotic behavior in the solutions.

We will also describe the method of the
Generalized Alignment Indices GALIk,
k=1,2,…,2N, which efficiently identify domains of
chaos and order in *N* dof Hamiltonian systems and
2*N*-dimensional (2N-D) symplectic maps.

**Indicators of regular and chaotic dynamics **
One of the most important questions in
Hamiltonian dynamics concerns the connection
between the* * local (linear) stability properties of
simple periodic solutions of Hamiltonian systems,
with the more “global” dynamics. We will examine
this question using the one-dimensional lattice (or
chain) of coupled oscillators called the Fermi Pasta
Ulam β-model described by the *N* dof Hamiltonian

� =1
2� �_{�}^{�}

�

���

+ �1

2 (�^{���}� �_{�})^{�}

�

���

+

+^{�}_{�}∑ (�^{�}_{���} ���� �_{�})^{�}= �, (5)
where *x*j are the displacements of the particles from
their equilibrium positions, and *p**j** = dx**j** /dt* are the
momenta, β is a positive real constant and *E* is the
total energy. Let us now consider some examples of
simple periodic solutions (SPOs), which have well-
defined symmetries and are known in closed
form:

(I) the out of phase (pi-mode)

���(�) = ������(�) = ��(�), � = 1,2, . . . , � (6) where N is even, under periodic boundary conditions;

(II) the SPO1 mode, where every 2 particles one is stationary and those on its either side move out of phase;

(III) the SPO2 mode, where every 3 particles one is stationary and the two on either side move out of phase both under fixed boundary conditions (fbc).

International Journal of Mathematics and Physics 9, №2, 21 (2018)
**Figure 4 –** Examples of SPOs that we have called the Out of Phase (or pi-)

Mode (above), the SPO1 orbit (middle) and the SPO2 orbit (below)

Applying Lyapunov's Theorem to the FPU system we can prove the existence of SPOs as continuations of the linear normal modes of the system, whose energies and frequencies are

�_{�} =1

2 ��^{�}^{�}� �_{�}^{�}�_{�}^{�}���

��������_{�}= 2��� �_{�(���)}^{��} � ����� = 1�2� � � � � ������ (7)
It is thus easy to verify that SPO1 and SPO2
orbits, as NNMs, are identified by the indices q =
(N+1)/2 and q = 2(N+1)/3 respectively.

In fact, it is possible to formulate a semi – analytical criterion for “weak” chaos:

We have verified numerically that the above NNMs first destabilize at energy densities of the form

* E**c**/N~ 1/N*^{ α}, α=1,2, as *N*→∞. (8)
In agreement with an analytical criterion by
Flach and co-workers [3], Ec/N ~ π^{2}/6βN^{2}, we find
that for α = 2 orbits (like SPO2) instability implies

“weak” chaos and the breakup of FPU recurrences.

On the other hand, if α = 1, for which the SPO1 mode destabilizes we find what we shall later call

“strong” chaos. Indeed, we believe that (8) may be true for other NNM solutions as well, but so far no proof of this statement is available.

**Lyapunov exponents and “strong” chaos **
Chaotic behavior is usually studied by
evaluating the spectrum of Lyapunov exponents,* L**i**, *
*i=1,….2N, *(LEs) *L**1**=L**max**>L**2**>….>L**2N* , defined as
follows:

�_{�}^{�}=1

� ��

�|��(�)|�

�|��(0)|������������� ���

�������_{�} = ���_{���}�_{�}^{�}��������� = 1�2� � � � �2� (9)
where �_{�}^{�} represent rates of separation from the
reference orbit of small deviations �_{�}(�) along the
2N directions in phase space. If the maximum of
these exponents *L**max* > 0, the orbit is chaotic, while
if *L**max* = 0 the orbit is stable. In the thermodynamic
limit, where *E*→∞ and *N *→∞ (with *E/N* fixed), the
Lyapunov spectrum near unstable NNMs tends to a
smooth curve, see Figure 5(a) above [2].

For our two orbits SPO1 and SPO2, at low energies when they are unstable, we find that their Lyapunov spectra are distinct see Figure 5(b).

Raising the energy, however, we observe in Figure
5(c) that the Lyapunov spectra converge to the same
exponentially decreasing function *L**i **(N) ~ exp(-*
*αi/N)*, ,i=1,2,…,N, thus providing evidence that the
orbits explore the same chaotic region.

International Journal of Mathematics and Physics 9, №2, 21 (2018)
**Figure 5** – (a) The spectrum of Lyapunov exponents near an out of phase orbit

of the β – FPU model as E and N grow (E/N=3/4). In (b) and (c)

the Lyapunov spectra of solutions starting near unstable SPO1and SPO2 orbits converge, as the energy grows from E= 2.1 in (b) to E=2.6 for (c),

indicating that the chaotic regions about these orbits have merged

**Beyond Lyapunov exponents: The **
**Generalized Alignment Indices (GALI) **

As mentioned above, the Lyapunov exponents
are always computed with respect to a *single *
*deviation vector* from the reference orbit. More than
a decade ago, however, several researchers [4, 5]

introduced an alternative approach by defining as
*the GALI indicators* quantities which take into
account simultaneously 2, 3 or more deviations from
the reference orbit, obtaining thus more
comprehensive results, enabling us to: (a) detect the
chaotic nature of the orbits more rapidly than other
methods and, (b) identify quasiperiodic motion
providing also the dimension of the torus.

The GALIk, k = 2, 3,…., N indicators are
defined, through the evolution of *k* initially linearly
independent deviation vectors wi(0), as the volume

of a k-parallelepiped given by the wedge product

����_{�} � �|��_{�}(�) � ��_{�}(�) � � ��_{�}(�)|��

�� � ���� � � � � � (10)
whose k edges are the unitary deviations ��_{�}(�) �

�_{�}(�)��|�_{�}(�)|�� Thus, it is evident that if at least
two of the deviation vectors become linearly
dependent, the volume of the k – parallelepiped
represented by the wedge product in (10) becomes
zero, and the GALIk vanishes. Thus, as expected, for
chaotic orbits, deviation vectors tend to become
linearly dependent in the direction defined by the
maximal Lyapunov exponent. As an example of this
effect, we show schematically in Figure 6 below
how this happens in the case of the GALI2 indicator

International Journal of Mathematics and Physics 9, №2, 21 (2018)
**Figure 6 –** Behavior of GALI2 for chaotic motion

In fact, one can show analytically that in the case of chaotic orbits, all GALIk tend to zero exponentially for large t, following a detailed asymptotic argument analyzing the quantities (10 ) in determinants and keeping the most dominant terms as t → ∞ [2]. To see the main idea of how this is done, we show in the next subsection that GALI2(t) ~ exp[-(σ1 – σ2 )t] →0 , σ1 > σ2 , being approximations of the two largest Lyapunov exponents L1 > L2 .

**Asymptotic analysis of the GALI****2**** for chaotic **
**motion **

The evolution of one deviation vector from a chaotic orbit can be approximated by the expression:

�_{�}(�) = � �_{�}^{(�)}�^{�}^{�}^{�}

��

���

�̂_{�} ≈

≈ �_{�}^{(�)}�^{�}^{�}^{�}�̂_{�}+ �_{�}^{(�)}�^{�}^{�}^{�}�̂_{�}+. .. (9)
where σ1 > σ2 > … are approximate values of the
Lyapunov exponents up to the time t of integration.

Thus, dividing this deviation by its magnitude we derive a leading order estimate for w1(t):

��= �_{�}(�)

�|�_{�}(�)|�=�_{�}^{(�)}�^{�}^{�}^{�}�̂_{�}+ �_{�}^{(�)}�^{�}^{�}^{�}�̂_{�}

��_{�}^{(�)}� �^{�}^{�}^{�} =

= ±�̂�+_{��}^{�}^{�}^{(�)}

�(�)��^{�(�}^{�}^{��}^{�}^{)�}��̂�, (10)
and an analogous expression for w2(t):

�� = �_{�}(�)

�|�_{�}(�)|�=�_{�}^{(�)}�^{�}^{�}^{�}�̂_{�}+ �_{�}^{(�)}�^{�}^{�}^{�}�̂_{�}

��_{�}^{(�)}� �^{�}^{�}^{�} =

= ±�̂�+_{��}^{�}^{�}^{(�)}

�(�)��^{�(�}^{�}^{��}^{�}^{)�}��̂�, (11)
Taking their cross product gives the following
result:

�����= �|��(�) � ��(�)|� ≈

≈ ��^{�}^{�}^{(�)}

�_{�}^{(�)}±^{�}^{�}^{(�)}

�_{�}^{(�)}�� �^{�(�}^{�}^{��}^{�}^{)�}, (12)
which clearly demonstrates what we referred to
above as exponential decrease of the GALIs to zero
as t goes to infinity.

**Behavior of GALI****2**** for regular motion **

It is also of great importance to analyze how
GALIs will behave in time if they represent
deviations of “ordered” or “regular” orbits, which
have zero Lyapunov exponents and lie on tori of N –
dimensional quasiperiodic motion in N dof
Hamiltonian systems. Remarkably enough, all
deviation vectors in that case become tangent to the
torus, and, for a k -dimensional torus, the GALI of
the associated k linearly independent vectors *will *
*not go to zero* since the volume of the corresponding
parallelepiped will *not vanish *(see a pictorial
representation of this in Figure 7 below in the case
of a 2 – dimensional torus).

We now make the following *very important *
*observation:* As we just explained, in the case of
regular orbits lying on s-dimensional tori, all

International Journal of Mathematics and Physics 9, №2, 21 (2018) deviation vectors tend to fall on the tangent space of

the torus. As a result, if we start with k ≤ s, the deviation vectors will remain linearly independent on the tangent space of the torus and the GALIk will be approximately constant, different from zero.

*Hence, for quasiperiodic motion, we find GALI**2**(t) ≈ *
*const. for all t > 0*.

Now, what is interesting is that if we start with k

> s deviation vectors, since only s of them will in
the end be linearly independent, the GALI will
again go to zero, but this *time following a power *
*law! *Clearly, this will be of great help in **identifying **
**the actual dimension** of the torus, for which no
other such criteria are available.

**Figure 7** – Behavior of GALI2 for regular motion occurring
on a 2-dimensional torus.

Summarizing, we have shown by asymptotic
analysis that: (a) for a chaotic orbit, all deviation
vectors tend to align in the direction of *L**1* , and all
GALIk tend to zero exponentially following the law

�����∞ exp (−(��− ��� ��− ��� � � ��− ��)�)(13)
where *L**i* are the k largest LEs. (b) On the other
hand, for k > s, all GALIk approach 0 as *t *→∞

following power laws,

����_{�}∞ 1

�^{���}� � � � � �� − ��

����� �

�^{�(���)}� �� − � � � � ��� (14)
since some deviation vectors will eventually become
linearly dependent. In Figure 8 below we display
some applications of the above theory to the study
of tori in a 2 dof and a 3 dof Hamiltonian system. In
the former case the tori are 2 – dimensional and for
this reason GALI2 tends to a constant while higher
GALIs decay by power laws given by (14), while in
the latter the tori are 3 – dimensional and, therefore,
not only GALI2 , but also GALI3 go to zero
following the power laws (14) presented above.

**Figure 8** – The GALI indices for a Hamiltonian system of (a) 2 degrees of freedom
and (b) 3 degrees of freedom. In case (a), since only GALI2 is constant the motion lies

on a 2-dimensional torus, while in (b), where both GALI2 and GALI3 are constant, the torus is 3-dimensional

International Journal of Mathematics and Physics 9, №2, 21 (2018)
**Localization in 1-dimensional lattices **

**Localization in Fourier space **

In 1955, by E. Fermi, J. Pasta and S. Ulam (FPU) used the computers available at the Los Alamos National Laboratory to integrate a chain of 31 nonlinear oscillators, coupled to their nearest neighbors, and investigate how energy was shared by all normal modes of the system. Starting with

initial conditions placed on the q=1 linear normal mode, they discovered, for small energies, a near- recurrence to their initial state after relatively short times exciting very few other modes, see Figure 9 (a) (for more explanations, the reader is again invited to consult the relevant chapters of [2]).

This remarkable observation ran contrary to the expectation of energy sharing among all modes predicted by equilibrium statistical mechanics and was termed the “paradox” of FPU recurrences.

**Figure 9 –** (a) Localization in modal space in the form of FPU recurrences,
discovered by Fermi Pasta and Ulam, for a lattice of 31 particles

"Energy localization” here implies localization
in Fourier q-modal space, as the FPU recurrences
were observed when all the energy was placed in the
*q*=1 mode.

Flach and his co-workers, in 2005 [3] introduced the concept of q-breathers, as exact periodic solutions of the problem. They showed that if we excite a single low q-breather mode the total energy remains localized only within a few of these low frequency modes, also called metastable states or natural packets.

A more complete interpretation of the FPU
paradox was provided by our group [6], where
we introduced the concept of* *q-tori, reconciling

q-breathers with the metastable packets of low- frequency modes. Now we shall use the GALI indices to study the stability of these q–tori and the breakdown of the associated FPU recurrences.

More specifically, in Figure 10 below we show
that it is indeed possible in an N = 8 dof case to
study *a torus of dimension 2* by selecting initial
conditions exciting a continuation of 2 linear modes
(see Figure 10(a)). Then, because this torus is stable,
its dimensionality for long times is verified by only
the GALI2 being constant while the higher order
GALIs decay by algebraic power laws given by (14)
(see Figure 10(b)).

International Journal of Mathematics and Physics 9, №2, 21 (2018)
**Figure 10 –** FPU with 8 particles: (a) Only the E1 and E3 modes are excited.

Observe that the associated q-torus is 2-dimensional, since (b) only GALI2 =const.

and all other GALIk decay by power laws

On the other hand, if the torus is unstable, and initial conditions nearby are going to wander in the weakly chaotic region that surrounds it, the GALIs are going to show it by falling exponentially to zero.

In Figure 11(a) below we exemplify this by showing the behavior of GALI2 , while all higher order

GALIs (not shown here) also fall to zero exponentially according to the laws (13). Note that, if one wanted to study this phenomenon by tracking the energies of the two excited modes, he would discover it much later in time through the excitation of other modes, s shown in Figure 11 (b)

**Figure 11 –** FPU with 8 particles and initial conditions near a q=2 torus:

(a) The evolution of GALI2 shows already at t ≈ 1000 that the orbit diffuses away from the torus weakly chaotically.

(b) This becomes visible, as the FPU recurrences break down at t ≈ 14000 through the excitation of different modes

**Localization in configuration space **

It is also very interesting to apply the above approach to the localization of spatial coordinates in nonlinear lattices, through the occurrence of a fascinating type of exponentially localized periodic oscillations, called discrete breathers [7, 8]. These solutions have been verified analytically and

numerically on a variety of lattices, like the Klein- Gordon (KG) chain

���� � ��^{�}(��) + �(����� ���+ ����),

�(�) �^{�}_{�}��^{�}+^{�}_{�}�^{�} (15)

International Journal of Mathematics and Physics 9, №2, 21 (2018)
- ∞ < *n* < ∞ , where *V*(x) is an on-site potential and

α > 0 is a coupling parameter. One example of such a breather solution is shown here in Figure 12 below. Expanding in Fourier series, one finds that

discrete breathers are directly related to homoclinic orbits of invertible maps, through which one can prescribe a numerical procedure for constructing them to arbitrarily high accuracy [7, 8].

**Figure 12** – Localization in configuration space in the form of a discrete breather
of a harmonic nearest neighbor chain with on site nonlinear potential of the Klein Gordon type

Indeed, keeping only the leading term
*x*n(t)=*A*ncos(ωbt) in such an expansion we obtain the
map ����� ������ ���� �^{��}���,

� � �� � �� � �_{�}^{�}���, (16)
which provides a very good approximation for the
amplitudes An, as homoclinic orbits lying at the
intersections of the invariant manifolds of the saddle
point at the origin of (16), at |C|>2.

Discrete breathers constitute one more example of what we call Simple Periodic Orbits, with all particles oscillating with frequency ωb outside the phonon band of NNMs.

An interesting question here would be to identify whether discrete breathers are surrounded by low-dimensional tori when they are stable. If so,

it would be pertinent to study the dimensionality of these tori and their stability using our GALI indices to determine if these localized solutions will eventually break down as time evolves.

As we observe in Figure 13 below, in the case of
a 31 particle lattice where a stable breather is
followed for very long times, the oscillations of the
central particle (larger band in Figure 13(a)), its two
adjacent particles (middle size band) and the very
small oascillations of all the others (small size band)
remain practically constant up to t =1.4x10^{6 }time
units. Using the GALIs, however, one does not need
to integrate over such extended time intervals. As is
already evident in Figure 13(b), for times as short as
t = 10^{4} time units, the stability and dimensionality of
the breather are detected by the constancy of GALI2,
and the fact that all other higher order GALIs decay
by power laws.

International Journal of Mathematics and Physics 9, №2, 21 (2018)
**Figure 13 –** (Stable torus): (a) The oscillations of the central three particles
of a KG chain of N=31 particles do not break down, forming a quasiperiodic breather.

(b) The torus is 2-dimensional, since only GALI2 remains constant, while all other GALIk decrease by power laws

**Figure 14** – (Unstable torus): (a) The oscillations of the central 3 particles,

starting further from the breather, appear quasiperiodic for very long times. (b) The solution, however, is chaotic and the “torus” eventually breaks down since the GALIs decay exponentially

On the other hand, when we increase the energy
of the system somewhat, the breather oscillations
become irregular (see Figure 14 (a)) and the
breather collapses after t=1.4x10^{6} time units, while
the GALIs fall exponentially much sooner (see
Figure 14 (b)), declaring after only a few thousand
time units that this solution is dynamically unstable
and will eventually break down!

**Complex statistics of chaotic dynamics **

To study and understand the statistical properties of chaotic behavior in Hamiltonian systems, it is important to recall first some basic facts of equilibrium thermodynamics also reviewed in [9]. As is well-known, in Boltzmann – Gibbs statistics, if a system can be at any one of i=1,2,...,W

states with probability *p**i* , its entropy is given by the
famous formula

��� = �� ∑� ������

��� ,

under the constraint

∑^{�}_{���}�_{�} = 1, (17)
where* k* is the Boltzmann's constant. The BG
entropy satisfies the property of additivity, i.e. if *A*
and *B* are two independent systems, their union
entropy is *S*BG(*A+B*) = *S*BG(*A*)+*S*BG(*B*). At thermal
equilibrium, and for a continuum set of states
depending on one variable, *x*, the probability density
that optimizes the BG entropy subject to the
constraints (17), zero mean and variance *V* is, of
course, the well-known Gaussian

International Journal of Mathematics and Physics 9, №2, 21 (2018)

�(�) = �^{��}^{�}^{/��}/√��. (18)
Another important property of the BG entropy is
that it is* extensive*, i.e. that *S*BG/*N* is finite in the limit
*N*→∞. But, many physically important systems
governed by long range interactions are neither
additive nor extensive, like self-gravitating systems
of finitely many mass points and ferromagnetic spin
models. For such systems the so-called Tsallis
entropy has been proposed [9]

�� = �^{��∑}_{���}^{�}^{���}^{�}^{�}^{�},
under the constraint

∑^{�}_{���}�_{�} = 1, (19)
depending on an index *q*. For a continuum set of
states* x*, the Tsallis entropy is optimized by the *q*-
Gaussian pdf

��(�) = ��_{�}^{���}^{�} =

= �(1 � (1 � �)��^{�})^{�/(���)}, (20)

where β=1/kT is a free parameter and *a* > 0 a
normalization constant. Expression (20) tends to a
Gaussian, as *q *→ 1 *e**q* →*e. *The Tsallis entropy is
not additive, and, in general, non-extensive. It offers
us the possibility of studying problems whose
correlations decay *not exponentially* but by power
laws, thus implying that the interactions within such
systems are of the long range type.

**The case of multi-degree-of freedom **
**Hamiltonian systems **

In the realm of multi – dimensional Hamiltonian systems analyzed in this paper, there are many situations where the dynamics is weakly chaotic and may, therefore, possess Tsallis statistics of the type described above.

For example, we have found that, in the β-FPU
model, near an unstable SPO1 orbit of a 5-particle
chain, orbits that remain “trapped” for very long
times in a thin chaotic region (see Figures 15, 16)
and are described by pdfs of the* *q-Gaussian type
with q ≈ 2.8 [11].

**Figure 15 –** Three different orbits with initial conditions very close to an unstable SPO1 orbit
of the 5 particle FPU-β chain: The black “figure eight” in the middle starts from a distance of 10^{-7} ,

the green one starts within 0.0001 and the red one extending over a much larger region, starts within 0.01 from the saddle point

These are what we call quasi-stationary states
(QSS) of the dynamics. Following these states for
long times, one typically finds that their pdfs pass
through* *QSS described by q-Gaussians of q > 1, see

Figure 16 (a,b,d,c), until they finally converge to Gaussians with q = 1, when the orbits escape to a much larger domain of strong chaos, see Figure 16 (c,e).

International Journal of Mathematics and Physics 9, №2, 21 (2018)
**Figure 16 –** Orbits starting at a distance of 1.0*10^{-7} from the unstable SPO1 orbit,

integrated for: (a) t=10^{5 }, (b) t=10^{7 }, until they eventually escape in the large chaotic sea (c) t=10^{8 }
(d – f): Plots of pdfs of position variables for a 5- particle FPU chain

and initial conditions close to an unstable SPO1 orbit. The QSS observed here are well described by q-Gaussians with (d) q = 2.78, then (e) q = 2.48, until the orbit drifts away to a wide chaotic sea

and the pdfs converge to (f) Gaussians with q=1.05

International Journal of Mathematics and Physics 9, №2, 21 (2018) Another example of a weakly chaotic orbit

located near an SPO2 orbit of a a multi – dimensional β – FPU Hamiltonian system has been

found to “stick” on a low – dimensional torus and remain for extremely long times near a type of quasiperiodic motion is shown in Figure 17 below

**Figure 17 –** (a) The dynamics near an SPO2 orbit ``sticks'' to a quasiperiodic torus,
at least up to t=10^{8} . The weakly chaotic nature of the motion is shown in (b),

where we have plotted the 4 largest Lyapunov exponents up to t=10^9.

Although they all decrease towards zero, at about t >10^{9} , the largest exponent shows a tendency
to converge to a positive value, indicating that the orbit is chaotic

As Figure 17 demonstrates this orbit is
dynamically very “stable”, as it remains for very
long times near a regular quasiperiodic torus. Its
chaotic nature, however, is clearly exemplified by
the fact that after nearly t= 10^{8 }time units its largest
Lyapunov exponents stops decreasing towards zero
and starts to converge to a positive value! This is
clear evidence that this orbit is not regular and can

therefore be characterized as weakly chaotic, according to the terminology used in this paper.

The interesting question that arises, therefore, is whether a statistical analysis of this orbit also shows that this orbit can also be characterized as weakly chaotic, by plotting the probability distributions of averaged sums of its coordinates, as was done above for the orbit shown in Figure 16.

**Figure 18 –** Left: The distribution of the normalized sum pdf of the orbit starting near SPO2,
for a total integration time t=10^{6}. Right: Final integration time t =10^{10}.

The pdf has converged to an almost analytical shape that is close to a q--Gaussian with q≈2.769 near the center, and seems to have analytical form

Indeed, as is clearly seen in Figure 18, the probability density functions (pdfs) associated with this orbit are not Gaussian but are well described by q- Gaussian distributions of the form (20). The

remarkable observation we make here is that, after
10^{10} time units, the pdfs appear to converge to a
*smooth* distribution, that departs from a q- Gaussian
type at large distances away from its mean value.

International Journal of Mathematics and Physics 9, №2, 21 (2018) This type of complexity provides one further

justification for the title of this paper.

**Weak chaos in 2–dimensional area – **
**preserving maps **

It is instructive to illustrate some of the phenomena we have described in this paper on a 2- dimensional area – preserving map, which represents a simple model of a Hamiltonian systems of 2 – degrees of freedom. The question we wish to address here is whether what we have called strong and weak chaos can be observed in systems of low dimensions as well. For this reason, we shall examine the behavior of a model called the McMillan map expressed by the following equations mapping the xn, yn plane to itself [11]:

����� ��,
�_{���}� ��_{�}+ ���_{�}

�_{�}^{�}+ 1+ ��_{�},

� � �,1,�,�, � �� (21)
For ε = 0 this system is integrable as it possesses
a constant of the motion of the form �^{�}+ �^{�}+

�^{�}�^{� }� ���� � ������ It is also easy to see that it
also possesses a saddle point at the origin for μ >1.

When ε > 0 and small (21) becomes non – integrable and chaotic orbits are expected to appear near its origin forming “figure eight” domains of the type shown in Figure 15. What we would like to investigate here is whether these orbits also display strongly and weakly chaotic properties similar to what we found earlier when we were discussing multi – dimensional systems.

The results we obtained in [11] we indeed quite interesting: We first noted that for several μ, ε values the chaotic orbits wandering around the saddle point at (0,0) formed indeed a “figure eight”

similar to the one of Figure 15 in a generally

strongly chaotic fashion [11]. In other words, when
followed for as many as 2^{20 }iterations (time units)
the pdfs produced by averaged sums of the
coordinates appeared to converge to Gaussian
distributions. In other words, the orbits wandered
around a “figure eight” domain chaotic domain and
their pdfs passed through a sequence of q-Gaussian
states, with q > 1, until they become true Gaussians,
with q = 1.

However, for certain choices of parameter values, the orbits exhibited a remarkable “diffusive’

behavior, as they began after a certain time to escape from the “figure eight” and wander about in the plane along a chaotic boundary surrounding chains of islands that encircled the central “figure eight” region! This produced a complex pattern of chaotic domains in which the orbits wandered about,

“sticking” often as it were in the vicinity of

“thinner” chaotic regions surrounding higher order
saddle points of the map extending over island
chains of stable periodic orbits of (21) that also
encircle the origin. It was precisely in these cases
where we discovered that the pdfs of our chaotic
orbits began to converge to a true q (>1)-Gaussian,
for n →∞ , as shown in Figure 19 and Figure 20
below, thus demonstrating their weakly chaotic
nature.** **

**The role of Long Range Interactions **

As was mentioned earlier, all the results
described above appear to suggest that the statistical
properties of chaotic motions in Hamiltonian
systems must depend on the *type of interactions*
characterizing these motions. More specifically, we
postulated that if these interactions are “short” (e.g.

exponential) this might account for the strong chaos characterizing states that are described by Gaussian pdfs. On the other hand, if the interactions are long ranged (e.g. decaying by power laws) this would entail that the pdfs obtained after sufficiently long times would be of the q – Gaussian type.

International Journal of Mathematics and Physics 9, №2, 21 (2018)

**Figure 19 –** Upper: Diffusive motion of orbits in a thin chaotic layer of the 2-d area – preserving map (21),
starting near the unstable fixed point at the origin, and evolving to N=2^{20} iterations

**Figure 20 –** The pdfs representing the normalized sum of averages of the xn coordinate of the map,

for the chaotic orbit of Figure 18 are seen to converge after 2^{20} iterations
to the q-Gaussian shown here, with q=1.6

International Journal of Mathematics and Physics 9, №2, 21 (2018) Thus, to test the validity of the above the above

conjecture we decided to extend our studies and consider a class of N – dimensional Hamiltonians that involve Long Range Interactions (LRI) of the kind exemplified by the following class of FPU models:

� =1
2 � �^{�}^{�}

�

���

+1

2 �(����� �_{�} )^{�}

�

���

+

+ _{�(�,�)}^{�} ∑ ∑ ^{(�}^{�}^{��}^{�}^{)}^{�}

�|���|�^{�}

������

���� = �(�), (22) where b>0 and α ≥ 0 is an important parameter introduced to “measure” the length of the interactions [12.13]. Note that to keep all energy terms in the Hamiltonian extensive, i.e. proportional to N, we have introduced before the quartic part of the potential the factor

�(�, �) = 1

� � �
1
(� � �)^{�}

�

�����

�

���

=

= _{�}^{�}∑^{�}_{���}^{�����}_{(���)}_{�} , (23)

What we found was indeed very interesting:

Noting first that α = ∞ represents the shortest type of
(nearest neighbor) interactions considered already in
sections 3 and 4 of this paper, we studied the full
range of α values all the way down to 0 ≤ α ≤ 1
representing the regime of the longest type of LRI
possible. As Figure 21 and Figure 22 below
demonstrate, as α became smaller than α = 1, a
surprising phenomenon of “regularization” of the
dynamics was observed (see Figure 21): The
maximum Lyapunov exponent λ > 0 was seen to
decrease to values that seem to tend to zero! In
other words, a kind of *weakly chaotic* behavior was
discovered, showing that LRI has a “stabilizing”

effect on the dynamics.

To see whether this “regularization” effect extends also to the statistical properties of the motion we studied the pdfs of sums of averaged momenta of our chaotic orbits and discovered that these also tended to become closely approximated by q – Gaussian distributions in the LRI regime 0 ≤ α ≤ 1 (see Figure 22). On the other hand, in the regime α > 1 where the interactions may be characterized of “shorter” type, the same pdfs quickly converged to q – Gaussians, demonstrating that the nature of the dynamics is characterized by a stronger form of chaos (see Figure 22)!

**Figure 21 –** LRI restores order out of chaos: For 0 < α < 1
the maximal Lyapunov exponent λ starts to decay to zero as

N increases and α “weaker” form of chaos is expected

International Journal of Mathematics and Physics 9, №2, 21 (2018)
**Figure 22 –** The momentum probability density function (pdf)
for Long Range Interactions, α = 0.7, converges to a q-Gaussian with q=1.249

indicating a “weaker” form of chaos as time increases

**Figure 23 –** On the contrary, the pdfs of the momenta for shorter ranges of interactions
with α > 1 (in the above example α = 1.4), are seen

to quickly converge to a pure Gaussian, indicating “strong” chaos

To test these ideas further, we proceeded to carry additional studies to see e.g. how the q index of the pdfs behaves in these very interesting LRI regimes and obtained results of the kind shown in Figure 23.

We also observed, however, that in all cases we studied, if the time interval of our integrations increased, a critical time tc was always reached where the q values started to decrease, exhibiting a tendency to go back to q =1 where strong chaos and real Gaussians prevail.

Another important effect was also observed in the LRI regime, when the number of particles N was increased (see Figure 25(b)): The q index of the pdfs

was also seen to increase as N becomes larger and larger, suggesting that weakly chaotic dynamics may also be found in the thermodynamic limit where N and the total energy E become larger and larger while E/N is kept constant.

Thus, scaling our parameter b by the critical time tc and plotting the pdfs at each tc vs. the value of 1/N we discovered a remarkable “phase diagram”, see Figure 25 below, dividing the parameter plane in two regimes: One characterized by Gaussians, strong chaos and BG statistics and one where pdfs are q – Gaussians and the weakly chaotic motion and Tsallis thermostatistics reign!

International Journal of Mathematics and Physics 9, №2, 21 (2018)
**Figure 24 –** The pdf index q decreases to 1 over the regime of Long Range Interactions,

from α = 0 to just above α = 1, where shorter range interactions take over

(a) (b)

**Figure 25 – **(a): For LRI, α= 0.7, the index q, starts to decrease towards 1 after a time threshold tc ≈ 10^{6},
for these parameter values. (b): Momentum distributions for the system

with b = 10, ε = 9, α = 0.7 for increasing N values.

Note that as N grows the pdfs are described by a q–Gaussian whose index q increases from 1.17 for N = 512 until 1.25 for N = 8192

Finally, we completed our study by plotting the q parameter as a function of 1/log(N) as N grows to higher and higher values. Based in the results plotted in Figure 27 we thus proposed the following formula

���� �� � �_{�}��� �_{����}^{����}. (24)

demonstrating the dependence of q on N in the LRI regime. As a result, we can use this formula to estimate the asymptotic behavior of q → q∞ in the limit N → ∞. Note that as N becomes larger and larger, the above formula allows us to determine the values q∞ (α) that the q-Gaussian pdfs will have in that limit.

International Journal of Mathematics and Physics 9, №2, 21 (2018)
**Figure 26 –** A “phase transition diagram” is obtained, separating BG from
Tsallis thermostatistics, in which the limits *t*→ ∞and *N*→ ∞ do not commute

**Figure 27 –** The index q depends linearly on 1/ \log N for N=4096, 8192, 16384,
as α changes, according the formula (24)

Next, in our more recent work [13], we extended the above results to study of N – dimensional Hamiltonian systems in which the “length” of the linear interactions (represented by quadratic terms in the potential of the Hamiltonian function) are characterized by a different α index than those of the nonlinear interactions, represented by quartic terms in the potential. We also considered the effect od LRI in N – dimensional Hamiltonian lattices, where besides the interparticle interactions each particle possesses an on site potential of its own [14].

It is not the place here, however, to also describe these results, as they are of a more advanced and

specialized character. We thus prefer to encourage the interested reader to consult the corresponding references [13, 14] and limit ourselves to summarizing these findings in the Conclusion section that follows.

**Conclusions **

In this work, we have reviewed a number of results of our research team in the broader field of study that may be called following Complex Dynamics and Statistics of Hamiltonian systems.

Although the work presented here has focused

International Journal of Mathematics and Physics 9, №2, 21 (2018) primarily on 1-D lattices of the so called β – FPU

models involving quartic and quadratic interactions, we believe that they are much more general and can be found also be bound in similar systems with other type of nonlinearities. Indeed, weakly and strongly chaotic motions appear in conservative systems of low dimensionality, even 2 – dimensional area preserving maps like the McMillan model analyzed in section 5.2.

Our main conclusions can, therefore, be summarized as follows:

1. We have demonstrated the importance of Nonlinear Normal Modes in exploring “weak” and

“strong” chaos, depending on the energy density Ec/N > 0 at which they first become unstable. At energies where they have just destabilized it is possible to find in their vicinity weakly chaotic motions that “stick” to the associated saddle points for very long times.

2. We mentioned the significance of Lyapunov spectra in quantifying strong chaos, and introduced the GALIk spectrum of indices k=2,3,4,…, which are best suited for identifying chaos, when they vanish exponentially. We also stressed the fact that in the case of quasiperiodic motion, where the GALIs decay as power laws, they offer the most convenient strategy known to date by which the dimensionality of the torus can be determined.

3. The GALI indices can also be used to study the breakdown of localization in 1-dimensional lattices: (i) In modal space connected to FPU recurrences and (ii) in position space, occurring in the form of discrete breathers, for which we can predict their breakdown long before it can be detected by other methods.

4. When Long Range Interactions are imposed (LRI) on the nonlinear forces (quartic terms in the potential of our FPU models) – for any range of linear interactions – we obtain weakly chaotic motion characterized by q-Gausian pdfs with q>1 (Tsallis thermostatistics).

5. More specifically, in the LRI regime, we
find a new “phase transition diagram”, separating
BG from Tsallis thermostatistics, in which the limits* *
*t *→ ∞and *N *→ ∞ do not commute.

6. When we introduce LRI *only on the linear *
*forces* (quadratic part of the potential) we obtain
strongly chaotic motion demonstrated by pure
Gausian pdfs with q = 1 (Boltzmann Gibbs
thermodynamics).

7. When LRI are imposed on the nonlinear forces, we find for long times limiting values q∞ > 1 as N → ∞ , showing that the system remains weakly

chaotic (with Tsallis and not Boltzmann Gibbs thermostatistics) in the thermodynamic limit.

8. Finally, when we include in our potentials, besides the interparticle interactions, terms associated with local potentials at the site of each separate particle, LRI again yields evidence of highly regular dynamics, as single--site excitations lead to special low--dimensional solutions, that may be well described by a 2 – degree of freedom Duffing oscillator. On the other hand, the behavior of the maximal Lyapunov exponent, suggests in that case an approach to “quasi integrable” behavior in the thermodynamic limit, characterized by non- Gaussian momentum distributions.

**References **

1. M. Lieberman, A. Lichtenberg. “Regular and Chaotic Dynamics.” Springer Verlag, New York (1986).

2. Bountis T. and Skokos H. “Complex Hamiltonian Dynamics.” Springer Synergetics series, Springer, Berlin (2012).

3. S. Flach, M.V. Ivanchenko and O.I.

Kanakov “q-Breathers and the Fermi-Pasta-Ulam Problem.” Phys. Rev. Lett. 95, 064102 (2005).

4. Skokos, Ch., Bountis, T., Antonopoulos, C.

“Geometrical properties of local dynamics in Hamiltonian systems: The generalized alignment index (GALI) method.” Physica D231, (2007): 30- 54.

5. Skokos Ch., Bountis T., Antonopoulos C.

“Detecting chaos, determining the dimensions of tori and predicting slow diffusion in Fermi Pasta- Ulam lattices by the generalized alignment index method.” Eur. Phys. J. Sp. Top. 165, (2008):

5-14.

6. Christodoulidi H., Efthymiopoulos C., Bountis T. “Energy localization on q-tori, long-term stability, and the interpretation of Fermi-Pasta-Ulam recurrences.” Phys. Rev. E 81, 016210 (2010).

7. J. Bergamin, T. Bountis and C. Jung “A Method for Locating Symmetric Homoclinic Orbits Using Symbolic Dynamics.” J. Phys. A: Math. Gen.

33, (2000): 8059-8070.

8. J. Bergamin, T. Bountis and M. Vrahatis,

“Homoclinic Orbits of Invertible Maps.”

Nonlinearity 15, (2002): 1603–1619.

9. C. Tsallis “Introduction to Nonextensive Statistical Mechanics – Approaching a Complex World.” (Springer, New York, 2009).

10. Antonopoulos C., Bountis T., Basios V.

“Quasi-stationary chaotic states of multidimensional

International Journal of Mathematics and Physics 9, №2, 21 (2018) Hamiltonian systems.” Physica A 390, (2011):

3290-3307.

11. Ruiz Lopez G., Bountis T. and Tsallis C. “Time – Evolving Statistics of Chaotic Orbits in Conservative Maps in the Context of the Central Limit Theorem.”

Intern. J. Bifurc. Chaos 22 (9), (2012): 12502.

12. Christodoulidi H., Tsallis C., Bountis T.

“Fermi-Pasta-Ulam model with long-range interactions: Dynamics and thermostatistics.” EPJ Letters 108, 40006 (2014).

13. Christodoulidi H. Bountis T., Tsallis C., Drossos L. “Chaotic Behavior of the Fermi- Pasta-Ulam Model with Different Ranges of Particle interactions.” J. Stat. Mech. 12 (12), 123206 (2016).

14. Christodoulidi H., A. Bountis A., Drossos L. “The Effect of Long-Range Interactions on the Dynamics and Statistics of 1-D Hamiltonian Lattices with On-Site Potential.” EPJST 227 (5, 6), (2018): 563.