PHYSICAL REVIEW E 84, 016205 (2011)

Nonlinear waves in disordered chains: Probing the limits of chaos and spreading

J. D. Bodyfelt,1 T. V. Laptyeva,1 Ch. Skokos,1,2,* D. O. Krimer,1,3 and S. Flach1

1
Max Planck Institute for the Physics of Complex Systems, Nöthnitzer Str. 38, D-01187 Dresden, Germany
2
Center for Research and Applications of Nonlinear Systems, University of Patras, GR-26500 Patras, Greece
3
Theoretische Physik, Universität Tübingen, D-72076 Tübingen, Germany
(Received 19 March 2011; published 11 July 2011)
We probe the limits of nonlinear wave spreading in disordered chains which are known to localize linear waves.
We particularly extend recent studies on the regimes of strong and weak chaos during subdiffusive spreading
of wave packets [Europhys. Lett. 91, 30001 (2010)] and consider strong disorder, which favors Anderson
localization. We probe the limit of infinite disorder strength and study Fröhlich-Spencer-Wayne models. We find
that the assumption of chaotic wave packet dynamics and its impact on spreading is in accord with all studied
cases. Spreading appears to be asymptotic, without any observable slowing down. We also consider chains with
spatially inhomogeneous nonlinearity, which give further support to our findings and conclusions.

DOI: 10.1103/PhysRevE.84.016205 PACS number(s): 05.45.−a, [Link], [Link]

I. INTRODUCTION observed, last forever or will it stop at certain (though probably

very large) time? (iii) Is the shape of the initial wave packet
Anderson localization [1] was discovered 50 years ago in
crucial for the details of spreading? We will mainly address
disordered crystals as an accumulation of single-particle elec-
question (ii) here.
tronic wave functions and can be interpreted as an interference
Johansson et al. [12] conjectured that spreading must
effect between multiple scatterings of the electron by random
eventually stop and dynamics will become close to regular,
defects of the potential. As a consequence eigenstates are
assuming that in these limits the Kolmogorov-Arnold-Moser
no longer spatially extended but are exponentially localized.
(KAM) theorem is applicable, i.e., that for small wave density
Anderson localization is a universal phenomenon of wave
regular nonergodic phase space structures predominate and the
physics, unrestricted to quantum mechanics. Experimental
dynamics develops along KAM tori. Other attempts consist in
observations were made in noninteracting Bose-Einstein con-
a numerical scaling analysis, in order to predict and extend
densates (BEC) expanding in random optical potentials [2,3],
results beyond computational ability [13]. Analytical studies
light propagation in spatially random nonlinear optical media
perform perturbation theory to higher order by treating the
[4,5], and microwave cavities filled with randomly distributed
strength of nonlinearity as a small parameter [14], conflicting
scatterers [6]. Anderson localization is a linear wave effect, i.e.,
with the explosive growth of secular terms in higher orders
it is well established for wave equations which are linear in
of perturbation theory. This theory states that for the disor-
the wave amplitude. However, in many cases one is confronted
dered discrete nonlinear Schrödinger model with nonlinearity
with a nonlinear response of the wave-carrying medium; for
strength exceeding a finite threshold, any initial localized wave
instance, high light intensities induce a nonlinear response
packet cannot fully spread to zero amplitudes at infinite time.
of the optical medium. Electron-electron and electron-phonon
In this case, a part of the excitation is self-trapped as a result of
interactions also result in substantial deviations from Anderson
nonlinearity-induced frequency shifts, which tune a localized
localization in solids. In experiments of Bose-Einstein conden-
excitation out of resonance with its surrounding nonexcited
sates the interatomic interactions are always present, although
linear modes. However, even in the case of strong nonlinearity,
they can be diminished by either decreasing atomic densities
subdiffusion of the non-self-trapped part is observed [9].
or by exploiting magnetically tunable Feshbach resonances.
When strong nonlinearities are avoided numerical studies
From a mathematical perspective, a linear wave equation
showed a rather universal asymptotic subdiffusive spreading of
is integrable, with each normal mode evolving independently
initial single-site excitations [8,9,15], which is characterized
in time. A localized wave packet in the presence of Anderson
by a growth of the second moment of the wave packet as
localization will therefore stay localized as time evolves. Non-
t α , with α < 1 [9,15]. The self-trapping theorem [16] holds
linearity will usually destroy the integrability of a system and
irrespective of the strength of disorder; therefore it is reflecting
induce mode-mode interactions. It was observed numerically
the properties of a strongly nonlinear lattice wave equation,
that wave packets in such nonlinear disordered wave equations
rather than peculiarities of waves propagating in disordered
delocalize in time without respecting Anderson localization
media. Also the self-trapping theorem crucially depends on
limits [7–11]. Thus, there are several intriguing questions
the presence of at least two integrals of motion, and it fails
which have attracted much attention during the last few years:
for most nonlinear wave equations with only one integral of
(i) Will Anderson localization be destroyed by arbitrary small
motion.
strength of nonlinearity or is there a threshold below which
In Ref. [9] the observed wave packet spreading was
the localization is restored? (ii) Will wave packet spreading, if
assumed to be due to an incoherent excitation of the wave
packet exterior, induced by the chaotic dynamics of the
wave packet interior. The number of resonant modes in the
*
hskokos@[Link] packet was estimated by considering quadruplet and triplet

1539-3755/2011/84(1)/016205(10) 016205-1 ©2011 American Physical Society

BODYFELT, LAPTYEVA, SKOKOS, KRIMER, AND FLACH PHYSICAL REVIEW E 84, 016205 (2011)

mode-mode interactions [17]. A generalization to higher di- For β = 0 and ψl = Al exp(−iλt), Eq. (2) reduces to the
mensions D and different nonlinearity powers was performed. linear eigenvalue problem
This led to a quantitative prediction for the subdiffusive
wave packet spreading characteristic α [9]. Its validity was λAl = l Al − Al−1 − Al+1 . (5)
confirmed numerically in [9,15,18]. Recently, it has been pre-  2
The normalized eigenvectors Aν,l ( l Aν,l = 1) are the
dicted theoretically [19,20] and verified numerically [19] that corresponding normal modes (NMs), and the eigenvalues
a potentially long-lasting strong chaos regime induces faster λν are the frequencies of these NMs. The width of the
(though still subdiffusive) spreading, which is followed by eigenfrequency spectrum λν in Eq. (5) is D = W + 4 with
the asymptotic and slower weak chaos subdiffusive spreading. λν ∈ [−2 − W2 ,2 + W2 ]. The coefficient 1/(2W ) in Eq. (3) was
Notably, published numerical data did not reveal a further chosen so that the linear parts of the Hamiltonians Eqs. (1) and
slowing down of spreading, when starting from the weak chaos (3) would correspond to the same eigenvalue problem. In the
regime. limit H → 0 (in practice by neglecting the nonlinear term
In this paper we present results of extensive numerical u4l /4) the KG model of Eq. (3)—with ul = Al exp(iωt)—is
studies of wave packet spreading in various models of reduced to the same linear eigenvalue problem of Eq. (5), under
disordered nonlinear one-dimensional lattices. In particular, the substitutions λ = W ω2 − W − 2 and l = W (˜l − 1). The
we consider different initial excitations and scan the parameter width of the squared frequency ων2 spectrum is K = 1 + W4
space of disorder strength and nonlinearity over a wide region.
with ων2 ∈ [ 12 , 32 + W4 ]. Note that D = W K . As in the case
The main aim is to test the applicability of previously derived
of DNLS, W determines the disorder strength.
spreading laws and to search for indications of a continuation
The asymptotic spatial decay of an eigenvector is given
of the weak chaos spreading, or for indications of a slowing
by Aν,l ∼ e−l/ξ (λν ) , where ξ (λν ) is the localization length. In
down, as conjectured by others.
the case of weak disorder, W → 0, the localization length
ν )  ξ (0) ≈ 100/W . The NM
2
is approximated [17,21] as ξ (λ
II. MODELS participation number pν = 1/ l Aν,l characterizes the spatial
4

A. Discrete nonlinear Schrödinger and Klein-Gordon chains extend of the NM. An average measure of this extent is the
localization volume V , which is of the order of 3.3ξ (0) ≈
In our study we consider various one-dimensional lattice 330/W 2 for weak disorder and unity in the limit of strong
models. The first one is the disordered discrete nonlinear disorder, W → ∞ [17]. The average spacing of eigenvalues
Schrödinger equation (DNLS) described by the Hamiltonian of NMs within the range of a localization volume is then
function d ≈ /V , with being the spectrum width. The two
 β frequency scales d  determine the packet evolution details
HD = l |ψl |2 + |ψl |4 − (ψl+1 ψl∗ + ψl+1
∗
ψl ), (1) in the presence of nonlinearity.
l
2
In order to write the equations of motion of Hamiltonian
in which ψl are complex variables, l are the lattice site indices,  in the normal mode space2 of the system we insert ψl =
(1)
and β  0 is the nonlinearity strength. The random on-site ν Aν,l φν in (2), with |φν | denoting the time-dependent
energies l are chosen uniformly from the interval [− W2 , W2 ], amplitude of the νth NM. Then, using Eq. (5) and the
with W denoting the disorder strength. The equations of orthogonality of NMs the equations of motion (2) read
motion are generated by ψ̇l = ∂HD /∂(iψl ): 
i φ̇ν = λν φν + β Iν,ν1 ,ν2 ,ν3 φν∗1 φν2 φν3 (6)
i ψ̇l = l ψl + β|ψl |2 ψl − ψl+1 − ψl−1 . (2) ν1 ,ν2 ,ν3

This set of equations conserves both the energy of Eq. (1) and with the overlap integral
 
the norm S = l |ψl |2 . Iν,ν1 ,ν2 ,ν3 = Aν,l Aν1 ,l Aν2 ,l Aν3 ,l . (7)
The second model we consider is the quartic Klein-Gordon
l
(KG) lattice, given as
The frequency shift of a single-site oscillator induced by
 p2 ˜l 1 1 the nonlinearity is δl = β|ψl |2 for the DNLS model. The
HK = l
+ u2l + u4l + (ul+1 − ul )2 , (3)
2 2 4 2W squared frequency shift of a single-site oscillator induced
l
by the nonlinearity for the KG system is δl = (3El )/(2˜l ),
where ul and pl are, respectively, the generalized coordinates with El being the energy of the oscillator. Since all NMs are
and momenta on site l, and ˜l are chosen uniformly from the exponentially localized in space, each effectively couples to a
interval [ 21 , 32 ]. The equations of motion are ül = −∂HK /∂ul finite number of neighbor modes. The nonlinear interactions
and yield are thus of finite range; however, the strength of this coupling
is proportional to the norm (energy) density for the DNLS
1
ül = −˜l ul − u3l + (ul+1 + ul−1 − 2ul ). (4) (KG) model. If the packet spreads far enough, we can
W generally define two norm (energy) densities: one in real
This set of equations only conserves the energy of Eq. (3). space, nl = |ψl |2 (El ), and the other in NM space, nν = |φν |2
The scalar measure of energy resulting from Eq. (3) we shall (Eν ). By averaging over realizations, no strong difference is
henceforth label as H . This scalar value H  0 serves as a seen between the two, and therefore we treat them generally
control parameter of nonlinearity, similar to β for the DNLS as some characteristic norm (n) or energy (E) density. The
case. frequency shift due to nonlinearity is then δD ∼ βn for the

016205-2
NONLINEAR WAVES IN DISORDERED CHAINS: PROBING . . . PHYSICAL REVIEW E 84, 016205 (2011)

TABLE I. Characteristic quantities of the DNLS (1) and the KG (3) models. The
dependence on the strength of disorder W of both the localization volume V and of the
average spacing d of NM eigenvalues within the range of V is given for the limiting cases of
weak W → 0 and strong disorder W → ∞. Note that n (E) represents a general characteristic
norm (energy) of wave packets of the DNLS (KG) model.

DNLS KG
 W W 1 3
On-site energies l ∈ − 2 , 2 ˜l ∈ 2 , 2
   
Spectrum λν ∈ −2 − W2 ,2 + W
2
ων2 ∈ 12 , 23 + 4
W
W +4
Spectrum width D =W +4 K =
W
 330
W →0 V =
Localization volume V W2
W →∞
V ∼1

W →0 dD ∼ W 2 dK ∼ W
Average spacing d
W →∞ dD ∼ W dK ∼ const.
3
Nonlinear energy shift δ δD ∼ βn δK ∼ E
2

DNLS model, while the squared frequency shift is δK ∼ 3E/2 interactions between neighboring NMs in order to rescale
for the KG lattice. The basic characteristics of both models are time [22]:
summarized in Table I.  p2
We order the NMs in space by ν 2 1
 increasing value of the H= + u + (uν+1 − uν )4 ,
ν
(8)
center-of-norm coordinate Xν = l lA2ν,l [9,15,18,19]. For ν
2 2 ν 4
DNLS we  follow normalized norm density distributions
zν ≡ |φν |2 / μ |φμ |2 , while for KG we follow normalized where the NMs are equivalent to the single-site oscillators.
 The NM eigenvalues ν are considered to be uncorrelated,
energy density distributions zν ≡ Eν / μ Eμ with Eν =
also for nearest neighbors. This is different from the DNLS
Ȧ2ν /2 + ων2 A2ν /2, where Aν is the amplitude of the νth NM and KG models. Also the FSW chain has only pair interactions
and ων2 its squared frequency. Wemeasure the second moment between NMs (sites). Note also that the nonlinear part of the
m2 = ν (ν − ν̄)2 zν (with ν̄ = ν νzν ), which quantifies the FSW Hamiltonian is invariant under any shift uν → uν + a,
wave packet’s
 spreading width; the participation number as opposed to the KG model.
P = 1/ ν zν2 , i.e., the number of the strongest excited modes
in zν ; and the compactness index ζ = P 2 /m2 , which quantifies
the inhomogeneity of a wave packet. Thermalized distributions C. Models with spatially inhomogeneous nonlinearity
have ζ ≈ 3, while ζ 3 indicates very inhomogeneous pack- We also consider two variants of DNLS and KG models
ets, e.g., sparse (with many holes) or partially self-trapped ones with spatially inhomogeneous nonlinearity. The first type of
(see [15] for more details). In addition, following Anderson’s lattice is composed of linear coupled oscillators except for a
definition of localization [1], we measure the fraction SV central region of length L where nonlinearities are present.
(HV ) of the wave packet norm (energy) in a localization We refer to this type as the LNL (linear-nonlinear-linear)
volume V around the initially excited state in real space. model. The second type is called NLN (nonlinear-linear-
For a localized state this fraction asymptotically tends to a nonlinear) and is the exact counterpart of the previous one,
constant nonzero value, while it goes to zero in the case of since the linear part of the lattice is located at the central
delocalization. L sites.

III. WAVE PACKET EVOLUTION

B. Fröhlich-Spencer-Wayne chain A. Theoretical predictions

In the limit of strong disorder (W → ∞) the DNLS 1. DNLS and KG

and KG models suffer from increasing computational times The evolution of wave packets in nonlinear disordered
needed to observe any nontrivial dynamics. This is because chains can be expected to be self-trapped for strong nonlin-
the eigenvectors tend to single-site profiles; i.e., the overlap earities, or show no self-trapping for weaker nonlinearities.
integrals become very small. Fröhlich, Spencer, and Wayne The existence of the self-trapping regime was theoretically
(FSW) suggested considering a modified Hamiltonian, which predicted for the DNLS model in [16] (see also [15] for
operates directly in normal mode space for the strong more details). According to the theorem stated in [16], for
disorder limit, but considers artificial rescaled anharmonic large enough nonlinearities (δD > D ) single-site excitations

016205-3
BODYFELT, LAPTYEVA, SKOKOS, KRIMER, AND FLACH PHYSICAL REVIEW E 84, 016205 (2011)

cannot uniformly spread over the entire lattice. Consequently, emulating the dynamics of the latter in NM space. However,
a part of the wave packet will remain localized, although the in the KG and DNLS case, triplet interactions between NMs
theorem does not prove that the location of this inhomogeneity are present and necessary in order to allow for finite (though
is constant in time. small) resonance probabilities for small (but finite) energy and
If the nonlinear shift δ moves the frequencies of some of the norm densities [17,20]. Pair interactions will cease to produce
initially excited oscillators out of the linear spectrum, it tunes NM resonances for sufficiently small densities due to level
them out of resonance and part of the wave packet will be self- repulsion within one localization volume [15,17,20]. The FSW
trapped. In our study we consider initial “block” wave packets, chain keeps only pair interactions. At the same time, level
where L central oscillators of the lattice are excited having repulsion between neighboring (interacting) NMs is absent in
the same norm (energy). Since we consider many random the FSW case. Therefore, a theoretical analysis similar to the
disorder realizations (of the order of a few hundred) we expect DNLS and KG case [20] appears to be possible. Its details
that, on average, the linear frequencies of the initially excited will be considered in a future publication. The width of the
lattice sites l (˜l ) cover the whole range of permitted values linear spectrum is FSW = 1. The average spacing of nearest
[− W2 , W2 ] ([ 21 , 32 ]). Thus, some of these frequencies are tuned neighbor eigenvalues is dFSW ∼ FSW . Therefore, we expect
out of resonance if δD  2 (δK  1/W ). These conditions for that only the asymptotic regime of weak chaos and the regime
the possible appearance of self-trapping are less strict than of self-trapping can be expected.
the theoretically defined ones [15,16] and are, in general, in
good agreement with numerical simulations. In particular, the 3. LNL and NLN chains
self-trapping regime was numerically observed for single-site
excitations [9,15,18] and for extended excitations [19], both The LNL chain can be expected to start in a chaotic wave
for the DNLS and the KG models, despite the fact that the KG packet spreading regime as long as the wave packet is confined
system conserves only the total energy H , and the self-trapping mainly to the finite-size N (nonlinear) part. However, the more
theorem cannot be applied there. the wave packet spreads, the more it extends into the infinitely
When self-trapping is avoided for δD < 2 (δK  1/W ), two extended L (linear) parts. Resonances and chaos are therefore
different spreading regimes were predicted, having different confined to the finite N part. Since distant NMs in the L part
dynamical characteristics: an asymptotic weak chaos regime are exponentially weakly interacting with the chaotic NMs in
and a potential intermediate strong chaos one [20]. Numerical the N part, their excitation—if present—will take times which
verifications of the existence of these two regimes were increase exponentially with growing distance. Therefore, the
presented in [18,19]. In the weak chaos regime, for L  V wave packet will spread (if at all) slower than any power law.
and δ < d most of the NMs are weakly interacting with each Thus, the LNL model is the only model we consider here,
other. Then the subdiffusive spreading of the wave packet where almost trivial slowing down of spreading is expected.
is characterized by m2 ∼ t 1/3 . If the nonlinearity is weak Recently, the dynamics of a similar system was theoret-
enough to avoid self-trapping, yet strong enough to ensure ically investigated in [24], where a disordered subsystem
δ > d, the strong chaos regime is realized. Wave packets in of coupled anharmonic oscillators, linearly coupled to an
this regime initially spread faster than in the case of weak infinite harmonic system, was considered. In this work the
chaos, with m2 ∼ t 1/2 . Since the norm density drops with conditions which permit the persistence of the discrete breather
further spreading, δ is dropping in time as well, and eventually of the isolated anharmonic system for small but nonvanishing
the wave packet enters the weak chaos regime, where its couplings to the harmonic lattice were derived, and cases
evolution is characterized by slower spreading with m2 ∼ t 1/3 . characterized by energy transfer to the harmonic system were
The wave packet evolution in both the weak and the strong also discussed.
chaos spreading regimes is also expected to be characterized The NLN chain is expected to behave differently. As long
by an increase of the participation number as P ∼ t α/2 when as the wave packet is confined mainly to the L region, the
m2 ∼ t α . dynamics is regular, and no spreading should occur. For large
Let as now discuss the spreading of wave packets when enough time, some part of the packet will leak out into the N
L < V . The packet will initially spread over the localization regions. Therefore, spreading of the wave packet should finally
volume V during a time interval τin ∼ 2π/d, even in the occur.
absence of nonlinearities [19,20]. The initial average norm
(energy) density nin (Ein ) of the wave packet is then lowered B. Numerical results
to n(τin ) ≈ nin L/V [E(τin ) ≈ Ein L/V ]. The further spreading
of the wave packet in the presence of nonlinearities is then 1. Methods
determined by these reduced densities. Note that for single- We consider compact DNLS wave packets at t = 0 span-
site excitations (L = 1), the strong chaos regime therefore ning a width L centered in the lattice, such that within L there
completely disappears and the wave packet evolves either in is a constant initial norm of nin = 1 and a random phase at
the weak chaos or in the self-trapping regime [8,9,15,23]. each site, while outside the width L the norm density is zero.
In the KG case, we excite each site in the width L with √ the
same energy, E = H /L, i.e., initial momenta of pl = ± 2E
2. FSW chain with randomly assigned signs.
There are no existing theoretical predictions for wave We use symplectic integration schemes of the SABA family
packet dynamics in the FSW chain. The FSW case can of integrators [15,25,26] for the integration of Eqs. (2) and
be considered as a strong disorder limit of the KG model, (4). The particular symplectic scheme used for the DNLS

016205-4
NONLINEAR WAVES IN DISORDERED CHAINS: PROBING . . . PHYSICAL REVIEW E 84, 016205 (2011)

5 b 4
(a) br (b) b bl (c)
g r
〈 log10 m2 〉

4 g 3

〈 log10 P 〉
b 2
m

〈 〉
r br 2
3 r
m 1.5 g
m 1
2 bl bl br
0.8
(d) 0.25 (e) bl (f)
b 0.8
0.6 g br 0.2 r
0.6 m

〈 SV 〉
b m
αm

0.15
αP
0.4 r
0.1 g 0.4 br
0.2 m br r
bl 0.05 0.2 g
bl b
0 0 0
0 2 4 6 0 2 4 6 0 2 4 6
log10 t log10 t log10 t
FIG. 1. (Color online) DNLS, W = 4: Evolution of (a) log10 m2 (t) , (b) log10 P (t) , (c) ζ (t) , (d) finite-difference derivative αm (t) for
the smoothed m2 data of panel (a), (e) finite-difference derivative αP (t) for the smoothed P data of panel (b), and (f) SV (t) for the spreading
of wave packets with initial width L = 21 and β = 0.012, 0.04, 0.18, 0.72, 3.6, 8.4 [(bl) black; (m) magenta; (r) red; (b) blue; (g) green; (br)
brown]. In panels (a), (b), (d), and (e) straight lines correspond to theoretically predicted power laws m2 ∼ t α , P ∼ t α/2 with α = 1/3 (dashed
lines) and α = 1/2 (dash-dotted lines). Error bars in panels (a) and (b) denote representative one-standard-deviation errors.

model is described in the Appendix. The number of lattice The weak chaos dynamics is observed in Fig. 1 for β =
sites, N, and the integration time step τ varied between 0.012 (black curves) and β = 0.04 (magenta curves). Initially,
N = 1000 to N = 2000 and τ = 0.01 to τ = 0.1, in order the wave packets remain localized and all quantities of Fig. 1
to exclude finite-size effects in the wave packet evolution, are constant with αm , αP being practically zero. After some
and in order to reach long integration times up to 107 –109 detrapping time td the wave packets start to subdiffusively
time units with feasible CPU times. In all our simulations, the spread with m2 ∼ t 1/3 and P ∼ t 1/6 [Figs. 1(a), 1(b), 1(d), and
relative energy and norm errors are kept smaller than 10−3 . For 1(e)]. In addition, the compactness index ζ ≈ 3 [Fig. 1(c)],
each parameter set we averaged our data over 1000 different indicating that wave packets are well thermalized inside. The
disorder realizations, unless otherwise stated, and denote this tendency toward complete delocalization of wave packets is
by . . . . In particular, we compute m2 and P , and we smooth clearly depicted in the evolution of the averaged norm fraction
log10 m2 and log10 P with a locally weighted regression SV , which remains at the L = 21 initially excited sites
algorithm [27], and then apply a central finite difference to [Fig. 1(f)]. After the detrapping time td SV decreases
calculate the local derivatives continuously up to the final integration time t = 107 .
By increasing the value to β = 0.18, the initial spreading
d log10 m2 d log10 P
αm = , αP = . (9) dynamics enters the crossover between weak and strong chaos,
d log10 t d log10 t and a faster spreading is observed. Spreading sets in earlier, the
compactness index again indicates thermalized wave packets,
and the local derivatives αm and αP increase up to 0.4 and 0.2,
2. Weak disorder respectively, with a possibly very slow decrease at even larger
In this section we considerably extend the reports on times. SV again continuously decreases to zero, indicating
the observation of weak chaos, strong chaos, the crossover complete delocalization.
between both and the self-trapping regime in Ref. [19]. In For β = 0.72, we fully enter the strong chaos regime. Most
Fig. 1 we show results for the DNLS model with W = 4 importantly we observe a saturation of the local exponent αm
and L = V = 21, for six different values of the nonlinearity around the theoretical value 1/2, with a subsequent decay,
strength β. The time evolution of log10 m2 (t) and log10 P (t) again as predicted by theory [20] and first observed in [19].
is plotted in Figs. 1(a) and 1(b), respectively. The evolution of Finally, for β = 3.6 and β = 8.4 (green and brown curves,
the compactness index ζ (t) is shown in Fig. 1(c). In Figs. respectively, in Fig. 1) the dynamics enters the self-trapping
1(d) and 1(e) we plot the time dependence of the numerically regime. We observe that a part of a wave packet remains
computed derivatives αm (t) and αP (t) (9) obtained from the localized, while the remainder spreads. The spreading portion
smoothed curves of Figs. 1(a) and 1(b), respectively. Finally, results in a continuous increase of m2 [Fig. 1(a)], which ini-
in Fig. 1(f) the values of SV (t) are plotted. tially is characterized by large values of αm > 1/2 [Fig. 1(d)].

016205-5
BODYFELT, LAPTYEVA, SKOKOS, KRIMER, AND FLACH PHYSICAL REVIEW E 84, 016205 (2011)

FIG. 2. (Color online) DNLS, W = 4: Time evolution of average norm density distributions zl in real space for (a) β = 0.04, (b) β = 0.72,
and (c) β = 3.6. The color scales shown on top of panels (a)–(c) are used for coloring each lattice site according to its log10 zl value.

For larger time αm decreases below 1/2. The evolution appears wave packets remain compact as they spread since ζ ≈ 3
to be rather complex and is not captured by the theoretical [Fig. 3(b)], and the fraction HV of the energy of the initially
considerations in [20]. The large values of αm may be due to L = 21 excited sites decreases [Fig. 3(d)]. For E = 0.04
temporal trapping and detrapping processes in this strongly we enter the crossover region between the weak and the
nonequilibrium dynamics of the wave packet, with a final strong chaos regimes, with characteristics similar to the DNLS
trapped packet fraction remaining. Therefore, we observe that case. For E = 0.2 (red curves in Fig. 3) we observe the
log10 P starts to level off [Fig. 1(b)], αP tends to very small typical behavior of the strong chaos scenario: spreading is
values [Fig. 1(e)], and ζ tends to zero [Fig. 1(c)]. In Fig. 1(f), characterized by a saturated αm ≈ 1/2 [Fig. 3(c)] for about
we observe that the values of SV saturate to higher values two decades (log10 t ≈ 3.5–5.5), followed by a crossover to
for β = 8.4. the weak chaos dynamics with αm decreasing. Getting closer
For these typical cases of weak chaos (β = 0.04), strong to the self-trapping regime for E = 0.75 (blue curves in Fig. 3),
chaos (β = 0.72), and self-trapping (β = 3.6) we present in or being deep inside it for E = 3 (green curves in Fig. 3) the
Figs. 2(a)–2(c) the time evolution of the averaged norm density
distributions zl in real space. All simulations presented in
Fig. 2 started from the same initial profile with size L = V ; 5 (a) g (b)
3 r m
therefore the width of the localization volume set by the linear b
〈 log10 m2〉

4 r
case is the width of the distributions at the shortest times in the bl
〈 〉

2
plots. In the weak chaos regime [Fig. 2(a)], the wave packets m
3
remain close to their initial configuration for some times, as the bl b
1
high zl values in the region of the initial excitation indicate, 2 g
followed by delocalization at larger times. Thus, at t = 107
the averaged wave packet has spread to about 100 sites with 0.8 (c) g 1 g (d)
zl > 10−5 . Therefore, the wave packet spreads continuously b 0.8
0.6 bl
over distances which are an order of magnitude larger than r
〈 HV〉

0.6
αm

the limits set by the linear theory and destroy Anderson 0.4 m
0.4
localization. In the case of strong chaos [Fig. 2(b)], spreading bl r
0.2
is even faster, leading to more extended profiles at t = 107 : 0.2 b m
there are about 400 sites with zl > 10−5 . In the self-trapping 0 0
2 4 6 8 2 4 6 8
regime [Fig. 2(c)], the spreading part of the wave packet covers log10 t log10 t
700 sites, with another clearly visible part staying self-trapped
at the initial excitation region. The curved fronts in the density FIG. 3. (Color online) KG, W = 4: Evolution of (a) log10 m2 (t) ,
plots in Fig. 2 follow from the theoretical
√ prediction m2 ∼ t α , (b) ζ (t) , (c) αm (t), and (d) HV (t) vs log10 t for the spreading
which leads to a packet width N ∼ m2 and N ∼ e(α log10 t)/2 . of initially compact wave packets of width L = 21 with E =
For the KG model (3) with W = 4 we present in Fig. 3 0.01, 0.04, 0.2, 0.75, 3 [(bl) black; (m) magenta; (r) red; (b) blue;
similar results to the ones for the DNLS model. For small (g) green]. In panels (a) and (c) straight lines correspond to the
values of the initial energy density E = 0.01, the charac- theoretically predicted power laws m2 ∼ t α with α = 1/3 (dashed
teristics of the weak chaos regime are observed: m2 ∼ t 1/3 lines) and α = 1/2 (dotted lines). Error bars in panel (a) denote
[black curve in Fig. 3(a)] after a detrapping time td ≈ 105 , representative one-standard-deviation errors.

016205-6
NONLINEAR WAVES IN DISORDERED CHAINS: PROBING . . . PHYSICAL REVIEW E 84, 016205 (2011)

characteristics of the self-trapping behavior appear, since ζ spreading after long time intervals. For W = 15 and β = 0.5
decreases [Fig. 3(b)], and HV tends to stabilize to nonzero we find td ≈ 103 [Fig. 4(a)]. The local derivative αm increases
small values [Fig. 3(d)]. Similar to the β = 3.6 and β = 8.4 from zero, showing a tendency to approach the theoretically
cases of the DNLS model, the evolution of m2 [Fig. 3(a)] predicted value α = 1/3 [Fig. 4(c)], and both ζ [Fig. 4(b)]
shows an initial phase of fast growth with αm > 1/2 [Fig. 3(c)] and SL [Fig. 4(d)] start to decrease. We note that since V ∼ 1
followed by a lowering in the values of αm . for large W , we measure the time evolution of the fraction SL (t)
Our numerical results are in accord with the predictions of of the norm density of the L = 10 initially excited sites.
weak and strong chaos regimes, as well as of the crossover The detrapping time td increases as W increases. This is
from strong to weak chaos. Self-trapping is observed as well, seen from the results for W = 40, β = 1 (magenta curves in
with less understood strongly nonequilibrium dynamics of the Fig. 5). In this case, we have to wait at least up to td = 105
trapped and spreading packet parts. We vehemently stress in order to get some evidence that spreading starts, since after
that in all our simulations we never observed any evidence that time log10 m2 starts to slightly grow [Fig. 5(a)], while ζ
of a wave packet transition from the weak chaos regime, [Fig. 5(b)] and SL [Fig. 5(d)] start to decrease. This increase
characterized by αm = 1/3, to a subsequent slowing down of td happens despite the fact that the nonlinearity strength β
of spreading, which would lead to αm < 1/3. also increased by a factor of 2 as compared to the W = 15
case. Nevertheless, even in this case of large W = 40 we are
3. Strong disorder able to numerically observe the onset of spreading in the weak
chaos regime. Increasing W to even higher values pushes td
To search for potential deviations from the predicted to values larger than the final integration time t = 107 used in
spreading laws, we turn to large values of W for the DNLS our simulations.
system. In all our simulations we had L  V . In particular, With increasing β we observe self-trapping, but a part of
we set L = 10 and considered the cases with W = 15 and the wave packet spreads and m2 increases [green and brown
W = 40. For large values of W we expect only the weak curves in Figs. 4(a) and 5(a)]. The average compactness index
chaos and the self-trapping dynamical regimes to be observed ζ decreases, which is a clear indication that a part of the wave
[19]. For each value of W three different values of β were packet remains localized, and reaches smaller final values for
considered, one being in the weak chaos regime and the other larger β [Figs. 4(b) and 5(b)]. The self-trapping of the wave
two in the self-trapping regime. In particular, we considered packets is also clearly seen from the evolution of SL . In
β = 0.5, β = 9, and β = 30 for W = 15 and β = 1, β = 25, Fig. 4(d) we see that for W = 15, β = 9 (green curve) and
and β = 100 for W = 40. The obtained results are shown in β = 30 (brown curve) SL decreases due to the spreading
Fig. 4 for W = 15 and in Fig. 5 for W = 40. of a part of wave packets, while, finally, it shows a tendency
For large W , the localization volume V decreases drasti- to level off to a positive value, indicating that part of the
cally, so that eventually the overlap integrals become small as wave packets remains localized. Similar behaviors of SL are
well. Therefore, the values of β at which spreading becomes observed in Fig. 5(d) for the W = 40 case with β = 25 (green
visible will increase. In the weak chaos spreading regime curve) and β = 100 (brown curve), although the plateauing of
the detrapping time td becomes large and we start observing

8 1.8 (a) 10
(a) br (b)
2.5 m
g 6 m 1.6 8
〈 log10 m2 〉

〈 log10 m2〉

2
〈 〉
〈 〉

4 1.4 6 br
g
1.5 1.2 g
2 4
m g br
br 1 m
1 (b)
2
(c) (d) (c) m
0.5 br
0.8 m 0.4 0.9
br
0.4 0.8
0.6 0.3 br
〈 SL〉

g
〈 SL〉

0.3
αm

αm

br 0.7
0.4 0.2 g
0.2 0.6
m g
0.2 g 0.1
0.1 0.5
m (d)
0 0 0 0.4
0 2 4 6 0 2 4 6 0 2 4 6 0 2 4 6
log10 t log10 t log10 t log10 t

FIG. 4. (Color online) DNLS, W = 15: Evolution of (a) FIG. 5. (Color online) DNLS, W = 40: Evolution of (a)
log10 m2 (t) , (b) ζ (t) , (c) αm (t), and (d) SL (t) vs log10 t for the log10 m2 (t) , (b) ζ (t) , (c) αm (t), and (d) SL (t) vs log10 t for the
spreading of initially compact wave packets of width L = 10 with spreading of initially compact wave packets of width L = 10 with
β = 0.5, 9, 30 [(m) magenta; (g) green; (br) brown]. Mean values are β = 1, 25, 100 [(m) magenta; (g) green; (br) brown]. Mean values are
averaged quantities over 100 disorder realizations. In panels (a) and averaged quantities over 100 disorder realizations. In panels (a) and
(c) straight lines correspond to theoretically predicted weak chaos (c) straight lines correspond to theoretically predicted weak chaos
behavior m2 ∼ t 1/3 . behavior m2 ∼ t 1/3 .

016205-7
BODYFELT, LAPTYEVA, SKOKOS, KRIMER, AND FLACH PHYSICAL REVIEW E 84, 016205 (2011)

(a)
6
(b)
Therefore, even for strong disorder, the dynamics of wave
3 packets evolves according to the theoretical predictions. Most
〈 log10 P, m2 〉

5 importantly, we do not observe a slowing down of the wave

〈 〉
4
packet below the limits set by the weak chaos regime.
2
r
3
b 4. FSW model
1
0.3 (c) 1 (d) To further probe a possible slowing down of wave packet
0.8
spreading beyond the limits of weak chaos, we turn to the
0.2 FSW model (8). In Fig. 6 we present results with initial
αp , αm

〈 HL 〉
0.6 compact wave packets of width L = 21, with energy density
0.1
0.4 E = 0.05, similarly to the KG model (3). We observe that
r
0 0.2 also for this model subdiffusive spreading occurs, because the
b second moment and the participation number [red and blue
0
2 4 6 8 2 4 6 8 curves, respectively, in Fig. 6(a)] increase continuously, and
log10 t log10 t the corresponding exponents αm and αP [Fig. 6(c)] tend to
eventual constant nonzero values. The fraction HL of energy
FIG. 6. (Color online) FSW: Evolution of (a) log10 m2 (t) [(r) remaining in the L = 21 initially excited sites [Fig. 6(d)]
red curve] and log10 P (t) [(b) blue curve], (b) ζ (t) , (c) αm (t) [(r) decreases as time increases, indicating the delocalization of
red curve] and αP (t) [(b) blue curve], and (d) HL (t) vs log10 t for
the wave packets. The compactness index [Fig. 6(b)] has a
initially compact wave packets of width L = 21 with E = 0.05.
different behavior with respect to what we have seen in the
rest of our simulations, as it decreases slowly up to t ≈ 107 ,
with a subsequent increase. Therefore, the wave packet is
highly inhomogeneous for almost all the integration time,
SL is not as clear as in Fig. 4(d). The numerically computed violating the assumptions which are used in the theoretical
exponents αm exhibit the typical behavior of self-trapping seen considerations of weak chaos [20]. The observed subdiffusive
in Fig. 1(d). For W = 15 they increase, reaching values larger spreading may still not be in its final asymptotic range. Still,
than 1/3, and afterward decrease toward αm = 1/3 [Fig. 4(c)]. we again do not see any signature of a slowing down of this
For W = 40 a similar behavior is observed for β = 100, while subdiffusive process, as was reported in [28] where a similar
for β = 25 αm seems to approach the theoretically predicted model was considered. Clearly the FSW model calls for a
value 1/3 from below. thorough and independent study.

3.5 (a) 2 (b) (c)

5
3
〈 log10 m2〉

〈 log10 P 〉

1.75 m
m 4
〈 〉

m
2.5 r r
1.5 g
g g 3
2
1.25 r
1.5 2
0.5 (d) g 0.25 (e) g 1 (f)
0.4 m 0.2 m 0.8
〈 HV〉

0.3 0.15 0.6 g

αm

αP

0.2 0.1 0.4 r

r
0.1 0.05 r 0.2 m
0 0 0
2 4 6 8 2 4 6 8 2 4 6 8
log10 t log10 t log10 t

FIG. 7. (Color online) Evolution of (a) log10 m2 (t) , (b) log10 P (t) , (c) ζ (t) , (d) αm (t), (e) αP (t), and (f) HV (t) vs log10 t for the
spreading of initially compact wave packets of width L = 21 with W = 4 and E = 0.02 in the KG [(m) magenta curves], the NLN KG [(g)
green curves] and the LNL KG [(r) red curves] models. In panels (a), (b), (d), and (e) straight lines correspond to theoretically predicted power
laws m2 ∼ t 1/3 , P ∼ t 1/6 of the weak chaos regime. Error bars in panels (a) and (b) denote representative one-standard-deviation errors.

016205-8
NONLINEAR WAVES IN DISORDERED CHAINS: PROBING . . . PHYSICAL REVIEW E 84, 016205 (2011)

5. LNL and NLN models take the disorder strength to its infinite limits, are also
In Fig. 7 we present results for the KG, LNL, and NLN showing subdiffusive growth. Most remarkably, in none of
models for W = 4, L = 21, and E = 0.02. For comparison we our studies (except the artificial LNL case) did we encounter a
also include the results for the KG model (3) (magenta curves) slowing down of spreading beyond the limits set by the weak
with E = 0.02, for which subdiffusive spreading in the weak chaos predictions. Therefore, our numerical data support the
chaos regime is observed. The NLN KG model (green curves conjecture that the wave packets, once they spread, will do so
in Fig. 7) exhibits a similar behavior, since both the second up to infinite times in a subdiffusive way, bypassing Anderson
moment [Fig. 7(a)] and the participation number [Fig. 7(b)] localization of the linear wave equations.
start to grow after some detrapping time td ≈ 105 . This time The only cases where spreading shows a tendency to stop
is larger than the detrapping time of the KG model (td ≈ 104 ), are the LNL models, for which nonlinearities are absent
because a wave packet in the NLN KG model initially evolves everywhere except inside a finite-size central region where the
in an almost linear system and only after some large time, initial wave packet is launched. In these models, when wave
when it has spread significantly to the nonlinear part of the packets have spread substantially, their chaotic component
lattice, does spreading takes on characteristics of the purely in the central region of the lattice becomes weak, and
nonlinear model. distant normal modes in the linear parts of the system are
On the other hand the evolution of all quantities of Fig. 7 exponentially weakly coupled to the central nonlinear region.
for the LNL KG system (red curves) follows the KG model When the nonlinearity strength tends to smaller values,
until t ≈ 104 , because initially the wave packets evolve in the waiting (detrapping) times for wave packet spreading of
same nonlinear system. Later on the wave packet enters the L compact initial excitations increase beyond the detection
(linear) parts of the system. Thus, spreading starts to retard, and capabilities of our computational tools. The corresponding
both log10 m2 (t) [Fig. 7(a)] and log10 P (t) [Fig. 7(b)] show question of whether a KAM regime can be entered at finite
a characteristic slowing down in the exponents αm [Fig. 7(d)] nonlinearity strength was addressed in [12] and is analyzed in
and αP [Fig. 7(e)]. In addition, SV [Fig. 7(f)] saturates at detail in a forthcoming work [29].
finite nonzero values, indicating that wave packets tend to
localize again. For all three KG models, the values of ζ ACKNOWLEDGMENTS
[Fig. 7(c)] show that wave packets do not become sparse and
We thank S. Aubry, S. Fishman, M. Mulansky, A. Pikovsky,
inhomogeneous in the course of time. We obtained similar
R. Schilling, and D. Shepelyansky for useful discussions. Ch.S.
results for the LNL and NLN DNLS models for W = 4, L =
was partly supported by the European research project “Com-
21, and β = 0.04.
plex Matter,” funded by the GSRT of the Ministry Education
Thus, spreading is also observed for the LNL and NLN
of Greece under the ERA-Network Complexity Program.
models, and only in the case of the LNL models have we
observed a slowing down of the spreading, as expected.
APPENDIX: SYMPLECTIC INTEGRATION
OF THE DNLS EQUATIONS
IV. SUMMARY AND CONCLUSIONS We discuss a novel method which we designed to integrate
We considered several models of disordered nonlinear the DNLS equations locally, which we shall call the ‘PQ’
one-dimensional lattices and performed extensive numerical method. Previously used methods employ a transformation of
simulations of norm (energy) propagations. Since we focused the wave function from real into Fourier space and back, at
on the dynamical spreading of fronts, we prepared initial block each integration step. These transformations induce small but
wave packet profiles, having widths equal to or larger than the observable corrections in the tails of the wave packet, which
average localization volume defined by the linear problem. We slowly but steadily grow in time. In such a case we will have
would expect similar behaviors for initial Gaussian profiles to stop the integration once this noisy background reaches
(although calculations were not performed), again where the a substantial level. The PQ method avoids the generation
width (for Gaussians, say the standard deviation) is on par with of this background by simply not performing the Fourier
the average localization volume. transformation. Instead the PQ method integrates the DNLS
We carefully studied statistical properties of the dynamics, equations in real space.
by varying the values of disorder and nonlinearity strengths The canonical transformation
over a wide interval, and by averaging results over many 1
ψl = √ (ql + ipl ) (A1)
disorder realizations. Our results agree quite well with our 2
theoretical expectations for the existence of the weak and
strong chaos regimes. of the complex variable ψl in Eq. (1) transforms (1) into
The main outcome of our study is that in the presence  l   β 2
of nonlinearities we always observe subdiffusive spreading, so HD = ql2 +pl2 + ql2 +pl2 −(ql+1 ql +pl+1 pl ), (A2)
2 8
that the second moment grows initially as m2 ∼ t α with α < 1, l

showing signs of a crossover to the asymptotic m2 ∼ t 1/3 law where ql and pl are generalized coordinates and momenta,
at larger times. Remarkably, subdiffusive spreading is also respectively.
observed for large disorder strengths, when the localization If a Hamiltonian function can be split into two integrable
volume (which defines the number of interacting partner parts, then a symplectic integration scheme can be used for the
modes) tends to one. Fröhlich-Spencer-Wayne models, which integration of its equations of motion. One possible splitting

016205-9
BODYFELT, LAPTYEVA, SKOKOS, KRIMER, AND FLACH PHYSICAL REVIEW E 84, 016205 (2011)

of the DNLS Hamiltonian (A2) into two separate Hamiltonian and

functions A and B is 
 l   β 2 ql = ql ,
A= ql2 + pl2 + ql2 + pl2 , e τ LQ : (A7)
2 8 pl = pl + (ql−1 + ql+1 )τ.
l
 (A3)
B=− (ql+1 ql + pl+1 pl ). This technique of splitting the Hamiltonian into multiple
l parts has been used in different applications of symplectic
Hamiltonian A is integrable and the operator eτ LA which integrators (see, for example [30]).
propagates the set of initial conditions (ql ,pl ) at time t to In our simulations we successively apply the SBAB2
their final values (ql ,pl ) at time t + τ is symplectic integrator [15,25,26] twice: first for the split
 (HD = A + B) of the DNLS Hamiltonian Eq. (A2) and
ql = ql cos(αl τ ) + pl sin(αl τ ), second for the split B = P + Q in Eq. (A3). The so-
e τ LA : (A4)
pl = pl cos(αl τ ) − ql sin(αl τ ), lution for the equations of motion from the Hamilto-
nian Eq. (A2) is thus approximated by the application
with αl = l + β(ql2 + pl2 )/2. Hamiltonian B of Eq. (A3) is of 13 simple operators on an initial condition (ql ,pl ),
not integrable, and thus the operator eτ LB cannot be written since
explicitly. If we consider B as a separate Hamiltonian function
and again split it as B = P + Q, the component parts eτ HD = eτ (A+B) ≈ ed1 τ LA ec2 τ LB ed2 τ LA ec2 τ LB ed1 τ LA
  2 2
P =− pl+1 pl , Q = − ql+1 ql (A5) ≈ ed1 τ LA ed1 c2 τ LP ec2 τ LQ ec2 d2 τ LP ec2 τ LQ ed1 c2 τ LP ed2 τ LA
2 2
l l ×ed1 c2 τ LP ec2 τ LQ ec2 d2 τ LP ec2 τ LQ ed1 c2 τ LP ed1 τ LA ,
are integrable, under the corresponding operators (A8)

pl = p l ,
τ LP
e : (A6) with the SBAB2 coefficients [25] of d1 = 1/6, d2 = 2/3, and
ql = ql − (pl−1 + pl+1 )τ c2 = 1/2.

[1] P. W. Anderson, Phys. Rev. 109, 1492 (1958). [14] S. Fishman, Y. Krivolapov, and A. Soffer, Nonlinearity 22, 2861
[2] J. Billy, V. Josse, Z. Zuo, A. Bernard, B. Hambrecht, P. Lugan, (2009).
D. Clément, L. Sanchez-Palencia, P. Bouyer, and A. Aspect, [15] Ch. Skokos, D. O. Krimer, S. Komineas, and S. Flach, Phys.
Nature (London) 453, 891 (2008). Rev. E 79, 056211 (2009).
[3] G. Roati, C. D’Errico, L. Fallani, M. Fattori, C. Fort, M. Zaccanti, [16] G. Kopidakis, S. Komineas, S. Flach, and S. Aubry, Phys. Rev.
G. Modugno, M. Modugno, and M. Inguscio, Nature (London) Lett. 100, 084103 (2008).
453, 895 (2008). [17] D. O. Krimer and S. Flach, Phys. Rev. E 82, 046221
[4] T. Schwartz, G. Bartal, S. Fishman, and M. Segev, Nature (2010).
(London) 446, 52 (2007). [18] Ch. Skokos and S. Flach, Phys. Rev. E 82, 016208 (2010).
[5] Y. Lahini, A. Avidan, F. Pozzi, M. Sorel, R. Morandotti, D. N. [19] T. V. Laptyeva, J. D. Bodyfelt, D. O. Krimer, Ch. Skokos, and
Christodoulides, and Y. Silberberg, Phys. Rev. Lett. 100, 013906 S. Flach, Europhys. Lett. 91, 30001 (2010).
(2008). [20] S. Flach, Chem. Phys. 375, 548 (2010).
[6] R. Dalichaouch, J. P. Armstrong, S. Schultz, P. M. Platzman, and [21] B. Kramer and A. MacKinnon, Rep. Prog. Phys. 56, 1469
S. L. Mccall, Nature (London) 354, 53 (1991); C. Dembowski, (1993).
H.-D. Gräf, R. Hofferbert, H. Rehfeld, A. Richter, and T. [22] J. Fröhlich, T. Spencer, and C. E. Wayne, J. Stat. Phys. 42, 247
Weiland, Phys. Rev. E 60, 3942 (1999); J. D. Bodyfelt, M. C. (1986).
Zheng, T. Kottos, U. Kuhl, and H.-J. Stöckmann, Phys. Rev. [23] H. Veksler, Y. Krivolapov, and S. Fishman, Phys. Rev. E 80,
Lett. 102, 253901 (2009). 037201 (2009).
[7] M. I. Molina, Phys. Rev. B 58, 12547 (1998). [24] S. Aubry and R. Schilling, Physica D 283, 2045 (2009).
[8] A. S. Pikovsky and D. L. Shepelyansky, Phys. Rev. Lett. 100, [25] J. Laskar and P. Robutel, Celest. Mech. Dyn. Astron. 80, 39
094101 (2008). (2001).
[9] S. Flach, D. O. Krimer, and Ch. Skokos, Phys. Rev. Lett. 102, [26] Ch. Skokos and E. Gerlach, Phys. Rev. E 82, 036704 (2010).
024101 (2009). [27] W. S. Cleveland and S. J. Devlin, J. Am. Stat. Assoc. 83, 596
[10] D. M. Basko, Ann. Phys. 36, 1577 (2011). (1988).
[11] I. Březinová, L. A. Collins, K. Ludwig, B. I. Schneider, and J. [28] M. Mulansky, K. Ahnert, and A. Pikovsky, Phys. Rev. E 83,
Burgdörfer, Phys. Rev. A 83, 043611 (2011). 026205 (2011).
[12] M. Johansson, G. Kopidakis, and S. Aubry, Europhys. Lett. 91, [29] M. V. Ivanchenko, T. V. Laptyeva, and S. Flach (in preparation).
50001 (2010). [30] K. Goździewski, S. Breiter, and W. Borczyk, Mon. Not. R.
[13] A. Pikovsky and S. Fishman, Phys. Rev. E 83, 025201(R) (2011). Astron. Soc. 383, 989 (2008).

016205-10