approximately how much time the system needs to become chaotic [Skokos, 2010;
Pikovsky and Politi, 2016],

$T_{L}=\frac{1}{\Lambda_{1}},$ (3.24)

refered to as Lyapunov time.

####  3.4.2 Deviation vector distribution

As already mentioned, the ftMLE measures the average exponential rate of
growth (or shrinking) of perturbations of an orbit and it is influenced by the
global dynamics of the system. Therefore is a macroscopic observable,
influenced by the overall dynamics. Trying to define a quantity which measures
the local contribution of perturbations to the chaotic dynamics of the system,
we compute what is called DVD [Skokos et al., 2013; Senyange et al., 2018;
Hillebrand et al., 2019; Miranda Filho et al., 2019; Ngapasare et al., 2019;
Many Manda et al., 2020] defined as the distribution of quantities

$\xi_{l}^{D}=\frac{\delta q_{l}^{2}+\delta p_{l}^{2}}{\sum_{l}\delta
q_{l}^{2}+\delta p_{l}^{2}},$ (3.25)

where $\delta q_{l}$ and $\delta p_{l}$ are the components of the deviation
vector $\bm{w}$ and $l$ is the lattice’s site. This distribution attains its
highest values at sites which are more sensitive to small changes of initial
conditions i.e. more chaotic. Consequently, it is a very useful tool to
visualize the location of chaotic contributions (see the papers cited above).

Similarly to what we did for the normalized norm/energy density $\xi$, during
our simulations we also calculate the statistical properties of the DVD
$\xi^{D}$. We compute the first $\overline{l}^{D}$ and second $m_{2}^{D}$
moments which measure respectively the center’s coordinates and the extent of
the distribution. In addition, we also evaluate the participation number
$P^{D}$ of the DVD which gives the number of highly chaotic sites. It is worth
noting that $P^{D}$ has already been used to study the localization of chaos
in the FPUT lattice [Dauxois et al., 1997] and DNA models [Barré and Dauxois,
2001; Hillebrand et al., 2019].

###  3.5 Numerical determination of power law exponents

As we will see in the two next chapters in several places we will need to
accurately determine the power law exponent of the time evolution of dynamical
quantities. Since the numerical data will not be smoothly located on a curve,
we need a reliable numerical approach to determine such exponents in expected
behavior of the form $Q(t)=Q_{0}t^{\alpha_{Q}}$. The methodology we briefly
describe here has been implemented successfully in many studies [Laptyeva et
al., 2010; Bodyfelt et al., 2011a, b; Laptyeva et al., 2012; Laptyeva, 2013;
Skokos et al., 2013; Laptyeva et al., 2014; Yusipov et al., 2017; Senyange et
al., 2018; Vakulchyk et al., 2019; Many Manda et al., 2020].

Let us present the whole process by considering that our observable $Q$ grows
as $Q(t)=t^{2/5}$. As mentioned, the raw data obtained from the numerical
simulations are not smooth. In order to mimic this behavior in this simple
example, we create a time series of data by introducing a random noise to the
exact values obtained by the law $t^{2/5}$ whose magnitude depends on the
expected $Q$ values.

Figure 3.1: Randomly created values [(b)lue line connected points] around an
artificial function $Q(t)=t^{2/5}$ [(r)ed curves]. The amplitude of the
randomness is $\frac{1}{50}Q(t)$ (a) and $\frac{1}{25}Q(t)$ (b).

Fig. 3.1 presents the evolution of such noisy outcome of $Q(t)$ for two
different amplitudes of the noise. For example, in Fig. 3.1 we see such a data
set created by random additions (or subtractions) to the $Q(t)=t^{2/5}$ values
when the magnitude of this randomness is $\frac{1}{50}Q(t)$ [panel (a)] and
$\frac{1}{25}Q(t)$ [panel (b)]. Now, the question is how can we obtain the
correct exponent ($2/5$ in our simple example) from the randomly scattered
points.

More generally, if we want to estimate the exponent $\alpha_{Q}$ of the
expression $Q(t)=Q_{0}t^{\alpha_{Q}}$, it is better to see this law in the
$\log-\log$ scale as

$\log_{10}Q=\alpha_{Q}\log_{10}t+\log_{10}Q_{0},$ (3.26)

and determine $\alpha_{Q}$ as

$\alpha_{Q}=\frac{d\log_{10}Q}{d\log_{10}t}.$ (3.27)

Since by line connecting the noisy data of Fig. 3.1 we do not get a
differentiable curve, we first need to smooth the raw data. We do so by
applying a locally weighted regression method developed in [Cleveland and
Devlin, 1988], then use finite difference methods to evaluate the slope (3.27)
[Fornberg, 1988, 1998]. More specifically in this work, we adapt the computer
codes developed by Hilboll [Hilboll, 2007] for the regression process and by
Xavier [Xavier, 2019] for the finite difference techniques.

We present in Fig. 3.2 the results of the smoothing procedure for the data of
Fig. 3.1 [(b)lue line connected points] along with original function
$Q=t^{2/5}$ [(r)ed curves]. There is a good agreement between the two, further
supported by the value of the numerically computed exponent $\alpha_{Q}$ of
the smoothed data presented in the insets of Fig. 3.2 whose values are very
close to $2/5$. Obviously, the deviation from the ‘actual’ value
$\alpha_{Q}=2/5=0.4$ is larger in the case of the more noisy points (Fig.
3.2(b)).

Figure 3.2: The smoothed data [(b)lue line connected points] of Fig. 3.1 along
with the function $Q=t^{2/5}$ [(g)reen curves]. The insets show the time
dependence in semi-log scale of the value of the numerically obtained exponent
$\alpha_{Q}$ (3.27) for the smoothed data [(b)lue curves]. The thick dashed
lines denote the value $0.4$.

###  3.6 General purpose integration schemes

Let us now briefly present some possible ways to numerically integrate the
system of ODEs (3.20). We will first discuss some general purpose integrators.

####  3.6.1 The explicit Euler method

Let us first describe the so-called explicit Euler method which is a basic
numerical technique for solving system (3.20). So, starting with the initial
condition, $\bm{Z}(t)$ at $t\geq 0$, we propagate the Hamiltonian flow and the
variational equations to $t+\tau$ using a linear interpolation

$\bm{Z}(t+\tau)=\bm{Z}(t)+\tau\bm{f}(\bm{Z}(t),t),$ (3.28)

where $\tau$ is the integration time increment. In what follows we will
occasionally use $\bm{Z}^{\prime}$ to denote $\bm{Z}(t+\tau)$ and $\bm{Z}$ for
$\bm{Z}(t)$. Following the process of [Yoshida, 1993], in order to illustrate
the applicability and the limitations of this method, let us use as a ‘toy
model’ the Hamiltonian of a single harmonic oscillator

$\mathcal{H}=\frac{\varepsilon}{2}\left(q^{2}+p^{2}\right),$ (3.29)

where $q$ and $p$ are the canonical generalized position and momentum
respectively and $\varepsilon$ is the proper frequency of the oscillator. Its
tangent dynamics is described by the Hamiltonian function of the form [Skokos
and Gerlach, 2010]

$\mathcal{H}_{V}=\varepsilon\left(q\delta q+p\delta p\right),$ (3.30)

with canonical coordinates $\delta p$ and $\delta q$. $\mathcal{H}_{V}$
depends explicitly on time $t$ through the presence of $q(t)$ and $p(t)$. In
this case Eq. (3.20) take the form

$\begin{pmatrix}\dot{q}\\\ \dot{p}\\\ \dot{\delta q}\\\ \dot{\delta
p}\end{pmatrix}=\bm{f}(q,p,\delta q,\delta p)=\begin{pmatrix}\varepsilon p\\\
-\varepsilon q\\\ \varepsilon\delta p\\\ -\varepsilon\delta q\end{pmatrix}.$
(3.31)

Consequently, the explicit Euler scheme (3.28) evolves Hamiltonians (3.29) and
(3.30) through

$\begin{pmatrix}q^{\prime}\\\ p^{\prime}\\\ \delta q^{\prime}\\\ \delta
p^{\prime}\end{pmatrix}=\begin{pmatrix}1&\varepsilon\tau&0&0\\\
-\varepsilon\tau&1&0&0\\\ 0&0&1&\varepsilon\tau\\\ 0&0&-\varepsilon\tau&1\\\
\end{pmatrix}\begin{pmatrix}q\\\ p\\\ \delta q\\\ \delta p\end{pmatrix}.$
(3.32)

Let us now, compute these Hamiltonian functions after one integration step
$\tau$, getting

$\mathcal{H}^{\prime}=\frac{\varepsilon}{2}\left(q^{\prime\leavevmode\nobreak\
2}+p^{\prime\leavevmode\nobreak\
2}\right)=\frac{\varepsilon}{2}\left(1+\varepsilon^{2}\tau^{2}\right)\left(q^{2}+p^{2}\right),$
(3.33)

and

$\mathcal{H}_{V}^{\prime}=\varepsilon\left(q^{\prime}\delta
q\prime+p^{\prime}\delta
p^{\prime}\right)=\varepsilon\left(1+\varepsilon^{2}\tau^{2}\right)\left(q\delta
q+p\delta p\right).$ (3.34)

We see that the new Hamiltonians for both the phase space orbit and the
related deviation vector have some additional terms as they now have the form
$\mathcal{H}^{\prime}=\mathcal{H}+\varepsilon^{2}\tau^{2}\mathcal{H}$ and
$\mathcal{H}_{V}^{\prime}=\mathcal{H}_{V}+\varepsilon^{2}\tau^{2}\mathcal{H}_{V}$.
This means that, an artificial drifting force proportional to $\tau^{2}$ has
been introduced introduced to the system whose origin is purely numerical.
Furthermore at each iteration, the cumulative drift from the initial
Hamiltonian will become larger up to the point where it could not be anymore
neglected in comparison to the original system. This deviation could lead to
the appearance of an artificial chaotic behavior. A common practice to try to
address this problem is the use of small integration time steps $\tau$ along
with the development of time adaptative step size integrators, which is not
very efficient for very long time integrations as the ones needed for our
study.

####  3.6.2 The Taylor/Lie series method

In general, the solution of Eq. (3.20) $\bm{Z}(t+\tau)$ at time $t+\tau$ can
be expressed as a Taylor series of $\bm{Z}(t)$

$\bm{Z}(t+\tau)=\bm{Z}(t)+\tau\frac{d\bm{Z}(t)}{dt}+\frac{\tau^{2}}{2}\frac{d^{2}\bm{Z}(t)}{dt^{2}}+\ldots+\frac{\tau^{n}}{n!}\frac{d^{n}\bm{Z}(t)}{dt^{n}}+\mathcal{O}\left(\frac{\tau^{n+1}}{(n+1)!}\frac{d^{n+1}\bm{Z}(t)}{dt^{n+1}}\right).$
(3.35)

If we consider the first two terms, we obtain the explicit Euler method
presented in Sec. 3.6.1. Thus an artificial drift from the primary system’s
dynamics is also present. Obviously, increasing the number of terms considered
in the expansion (3.35) this artificial drift reduces. The scheme presented in
Eq. (3.35) where the first $n$ terms of the expansion are used is said to be
of order $n$. Further information about this method is presented in [Hairer et
al., 1993, Sec. I.8] and references therein, while a detailed discussion of
the performances of this technique in the present context can be found for
example in [Gerlach et al., 2012; Danieli et al., 2019a].

Following [Gröbner, 1967; Hanslmeier and Dvorak, 1984; Eggl and Dvorak, 2010;
Blanes and Casas, 2016], we will discuss an elegant way to present the Taylor
series using Lie algebra. Indeed, Eq. (3.20) can be written in terms of Lie
operators as

$\frac{d\bm{Z}(t)}{dt}=\bm{L}_{\mathcal{H}_{V}}\bm{Z}(t),$ (3.36)

where the Lie operator $\bm{L}_{\mathcal{H}_{V}}$ is defined as

$\bm{L}_{\mathcal{H}_{V}}=\sum_{i=1}^{2N}f_{i}\frac{\partial}{\partial
z_{i}}+f_{i+2N}\frac{\partial}{\partial\delta z_{i}}.$ (3.37)

Then the solution of Eq. (3.36) can be formally written as

$\bm{Z}(t+\tau)=e^{\tau\bm{L}_{\mathcal{H}_{V}}}\bm{Z}(t),$ (3.38)

where the exponential differential operator $e^{\tau\bm{L}_{\mathcal{H}_{V}}}$
plays the role of the propagator of the solution as its application on the
initial state $\bm{Z}(t)$ advances the solution by time step $\tau$ i.e. to
$\bm{Z}(t+\tau)$. In the case of integrable systems, such operators represent
the exact solution of the equations of motion and the variational equations.
However in the case of non-integrable systems, the Lie series integrator
approximates the solution through a series expansion of the operator
$\bm{L}_{\mathcal{H}_{V}}$

$\bm{Z}(t+\tau)=\bm{L}_{\mathcal{H}_{V}}^{0}\bm{Z}(t)+\tau\bm{L}_{\mathcal{H}_{V}}\bm{Z}(t)+\frac{\tau^{2}}{2}\bm{L}_{\mathcal{H}_{V}}^{2}\bm{Z}(t)+\ldots+\frac{\tau^{n}}{n!}\bm{L}_{\mathcal{H}_{V}}^{n}\bm{Z}(t)+\mathcal{O}\left(\frac{\tau^{n+1}}{(n+1)!}\bm{L}_{\mathcal{H}_{V}}^{n+1}\bm{Z}(t)\right).$
(3.39)

The next step in order to practically implement Eq. (3.39) is to find an
analytical expression for the action of $\bm{L}_{\mathcal{H}_{V}}$ and its
powers $\bm{L}_{\mathcal{H}_{V}}^{n}$ on the system state $\bm{Z}$ at time
$t$. Equation (3.39) corresponds to the Lie series method of order $n$ and is
equivalent to expression (3.35) by considering
$\bm{L}_{\mathcal{H}_{V}}^{0}\bm{Z}=\bm{Z}$,
$\bm{L}_{\mathcal{H}_{V}}\bm{Z}=\frac{d\bm{Z}}{dt}$,
$\bm{L}_{\mathcal{H}_{V}}^{2}\bm{Z}=\frac{d^{2}\bm{Z}}{dt^{2}}$ etc.

Actually, the determination of an analytical closed form of
$\bm{L}_{\mathcal{H}_{V}}^{n}\bm{Z}$ is a very challenging task even for
systems with few degrees of freedom [Eggl and Dvorak, 2010], and becoming
almost impossible for multidimensional systems. Nevertheless the application
of modern computer algebra softwares, can automate successive differentiations
of the force field vector $\bm{f}$ (see e.g. [Barrio, 2005]). This approach
has been followed for instance using the software Mathematica111Wolfram
Mathematica: Modern technical computing. https://www.wolfram.com/mathematica/.
in the freely available package called $TIDES$ [Abad et al., 2014]. The
$TIDES$ package has already been used in several investigations of dynamical
systems (see e.g. [Barrio and Shilnikov, 2011; Gerlach et al., 2012; Barrio et
al., 2015; Wilczak et al., 2016; Danieli et al., 2019a]). It requires as input
the set of differential equations to be integrated along with a set of initial
conditions, as well as the final time of integration. It creates internally a
Fortran77 or C code which is immediately compiled producing an executable
program which generates the numerical solution. The accuracy of the performed
integration does not only depend on the integration time step $\tau$ but also
on a parameter (denoted as $\delta_{a}$) defining the integrator’s single step
precision. As a consequence, the $TIDES$ package creates an adaptive step size
integrator.

####  3.6.3 The General purpose DOP853 Runge-Kutta integrator

The Runge-Kutta (RK) numerical integration methods are some of the most used
numerical schemes. Their simplicity and robustness make them widely used in
several scientific fields such as physics, astrophysics, population dynamics,
etc. Details on this family of integrators can be found in [Hairer et al.,
1993, 2006; Eggl and Dvorak, 2010]. A p-stage RK method for Eq. (3.20), is
described by the scheme

$\displaystyle\bm{Z}(t+\tau)$
$\displaystyle=\bm{Z}(t)+\tau\sum_{i=1}^{p}b_{i}\bm{k}_{i},\quad\text{with}$
(3.40) $\displaystyle\bm{k}_{i}$
$\displaystyle=\bm{f}\left(t+c_{i}\tau,\bm{Z}(t)+\tau\sum_{j=1}^{p-1}a_{ij}\bm{k}_{j}\right),$
(3.41)

where the values of the coefficients $a_{i,j}$, $b_{i}$ and $c_{i}$ depend on
the desired accuracy. In order to obtain numerical values for these
coefficients, one has to perform Taylor series expansions and deduce order
conditions which are algebraic equations on $a_{ij}$, $b_{i}$ and $c_{i}$. As
seen previously, Taylor series expansions introduce an artificial drift which
is also expected to appear for this method.

In the rest of this thesis, unless otherwise state, we implement a $12$-stage
explicit RK algorithm called DOP$853$ which is freely available online
[Hairer, 2009] and is based on the Dormand and Prince method [Hairer et al.,
1993]. The DOP$853$ has already been extensively used in lattice dynamics (see
e.g. [Mulansky et al., 2009; Sales et al., 2018]).

###  3.7 Symplectic integrators

Based on the nature of the physical phenomena they describe, Hamiltonian
models should conserve certain macroscopic observables such as mass, energy
and number of particles. These constants are called integrals of motion and,
in principle, described by analytical functions whose values remain constant
throughout the system’s evolution. Since the basic principle of finding
accurate numerical solutions is to get as close as possible to the properties
of the real analytical solutions, finding numerical integration schemes which
preserve the integrals of motion is highly desirable. SIs are numerical
integration schemes specifically designed for autonomous Hamiltonian systems
as they preserve by default the system’s phase space volume
$d\mathfrak{V}=d\bm{q}\wedge d\bm{p}$. The numerical schemes presented in
Secs. 3.6.2 and 3.6.3 were losing the preservation of symplecticity once we
considered that the action of the operator $e^{L_{\mathcal{H}_{V}}}$ was
approximated by a truncated Taylor series expansion (3.39) (see for example
Sec. 3.6.1). Let us consider a particular class of Hamiltonian functions which
is encountered in many applications, systems whose Hamiltonian can be
expressed as a sum of a kinetic energy $A(\bm{p})$ depending only on the
momentum variables $\bm{p}$ and a potential energy $B(\bm{q})$ term depending
only on the position variables $\bm{q}$

$\mathcal{H}(\bm{q},\bm{p})=A(\bm{p})+B(\bm{q}).$ (3.42)

The formal solution of Eq (3.20) (see (3.38)) becomes

$\bm{Z}(t+\tau)=e^{\tau(\bm{L}_{A_{V}}+\bm{L}_{B_{V}})}\bm{Z}(t).$ (3.43)

Since each part $A(\bm{p})$ and $B(\bm{q})$ possess $N$ cyclic coordinates
which act as integrals of motion, the analytical forms of the operators
$e^{\tau\bm{L}_{A_{V}}}$ and $e^{\tau\bm{L}_{B_{V}}}$ are explicitly known.
These operators represent the exact solutions of the equations of motion and
the related variational equations of the respectively simple Hamiltonian
functions $A(\bm{p})$ and $B(\bm{q})$

Figure 3.3: Geometric representations of the symplectic explicit Euler (3.28)
(a) and leapfrog (LF) (3.55) (b) methods. The bold arrows represent the paths
followed by the successive steps within the numerical scheme and the dashed
arrows joins the initial and final integration points. In (a)
$\bm{Z}_{1}=e^{\tau\bm{L}_{A_{V}}}\bm{Z}$ and
$\bm{Z}^{\prime}=e^{\tau\bm{L}_{B_{V}}}\bm{Z}_{1}$ respectively, giving a $2$
step SI. In (b) $\bm{Z}_{1}=e^{\tau\bm{L}_{A_{V}}/2}\bm{Z}$,
$\bm{Z}_{2}=e^{\tau\bm{L}_{B_{V}}}\bm{Z}_{1}$ and
$\bm{Z}^{\prime}=e^{\tau\bm{L}_{A_{V}}/2}\bm{Z}_{2}$ such that the LF
corresponds to a $3$ step SI.

$e^{\tau\bm{L}_{A_{V}}}\colon\begin{dcases}\bm{q}^{\prime}&=\bm{q}+\tau\left.\frac{\partial
A}{\partial\bm{p}}\right\rvert_{\begin{subarray}{c}\bm{p}=\bm{p}\\\
\delta\bm{p}=\delta\bm{p}\end{subarray}}\\\
\delta\bm{q}^{\prime}&=\delta\bm{q}+\tau\left[\delta\bm{p}\cdot\bm{D}^{2}A\right]_{\begin{subarray}{c}\bm{p}=\bm{p}\\\
\delta\bm{p}=\delta\bm{p}\end{subarray}}\end{dcases},\quad
e^{\tau\bm{L}_{B_{V}}}\colon\begin{dcases}\bm{p}^{\prime}&=\bm{p}-\tau\left.\frac{\partial
B}{\partial\bm{q}}\right\rvert_{\begin{subarray}{c}\bm{q}=\bm{q}\\\
\delta\bm{q}=\delta\bm{q}\end{subarray}}\\\
\delta\bm{p}^{\prime}&=\delta\bm{p}-\tau\left[\delta
q\cdot\bm{D}^{2}B\right]_{\begin{subarray}{c}\bm{q}=\bm{q}\\\
\delta\bm{q}=\delta\bm{q}\end{subarray}}\end{dcases}.$ (3.44)

Here $\bm{D}^{2}$ represents a multi-variable derivative operator of order two
with respect to the momenta ($\bm{D}^{2}A$ in $e^{\tau\bm{L}_{A_{V}}}$) or the
positions ($\bm{D}^{2}B$ in $e^{\tau\bm{L}_{B_{V}}}$). Further assuming that
$\bm{L}_{A_{V}}$ and $\bm{L}_{B_{V}}$ commute (something which is generally
not true), we can write

$e^{\tau\bm{L}_{\mathcal{H}_{V}}}=e^{\tau\bm{L}_{B_{V}}}e^{\tau\bm{L}_{A_{V}}},$
(3.45)

such that

$\bm{Z}(t+\tau)=e^{\tau\bm{L}_{B_{V}}}e^{\tau\bm{L}_{A_{V}}}\bm{Z}(t).$ (3.46)

The geometrical representation of Eq. (3.46) is presented in Fig. 3.3(a). We
see that $\bm{Z}$ and $\bm{Z}^{\prime}$ are related through the application of
a two step process involving two perpendicular mappings. Thus using Eq.
(3.45), we approximate the exact solution described by
$e^{\tau\bm{L}_{\mathcal{H}_{V}}}$, by the composition of two operators
$e^{\tau\bm{L}_{A_{V}}}$ and $e^{\tau\bm{L}_{B_{V}}}$. Since each operator is
symplectic, so is the composition. Thus we have constructed a very simple SI

$S_{1}(\tau)=e^{\tau\bm{L}_{B_{V}}}e^{\tau\bm{L}_{A_{V}}},$ (3.47)

which is called the symplectic explicit Euler method.

For the sake of completeness, let us apply the SI (3.47) to the single
harmonic oscillator (3.29) as done in [Yoshida, 1993]. Indeed, Hamiltonain
(3.29) can be splitted in two parts $A(p)=\varepsilon p^{2}/2$ and
$B(q)=\varepsilon q^{2}/2$. Thus Eq. (3.44), gives

$e^{\bm{L}_{A_{V}}}\bm{Z}(t)\colon\left\\{\begin{array}[]{ll}q^{\prime}&=q+\tau\varepsilon
p\\\ p^{\prime}&=p\\\ \delta q^{\prime}&=\delta q+\tau\varepsilon\delta p\\\
\delta p^{\prime}&=\delta p\\\ \end{array}\right.,\quad\mbox{and}\quad
e^{\bm{L}_{B_{V}}}\bm{Z}(t)\colon\left\\{\begin{array}[]{ll}q^{\prime}&=q\\\
p^{\prime}&=p-\tau\varepsilon p\\\ \delta q^{\prime}&=\delta q\\\ \delta
p^{\prime}&=\delta p-\tau\varepsilon\delta p\\\ \end{array}\right.,$ (3.48)

such that the numerical scheme $S_{1}(\tau)$ (3.47) gives

$\begin{pmatrix}q^{\prime}\\\ p^{\prime}\\\ \delta q^{\prime}\\\ \delta
p^{\prime}\end{pmatrix}=\begin{pmatrix}1&\varepsilon\tau&0&0\\\
-\varepsilon\tau&1-\varepsilon^{2}\tau^{2}&0&0\\\ 0&0&1&\varepsilon\tau\\\
0&0&-\varepsilon\tau&1-\varepsilon^{2}\tau^{2}\\\
\end{pmatrix}\begin{pmatrix}q\\\ p\\\ \delta q\\\ \delta p\end{pmatrix}.$
(3.49)

The computation of the energy (3.29) after one time step $\tau$

$\mathcal{H}^{\prime}=\mathcal{H}+\tau^{2}\frac{\varepsilon^{3}}{2}\left(q^{2}-p^{2}\right)+\tau^{3}\frac{\varepsilon^{4}}{2}qp+\tau^{4}\frac{\varepsilon^{5}}{2}p^{2},$
(3.50)

depicts a non-conserved total energy of the system. Nevertheless, the presence
of a negative term in $\tau^{2}$, ensures that there is no monotonic increase
of the energy error which is bounded by $\mathcal{O}\left(\tau\right)$. This
observation is further explained by the fact the mapping (3.49) has an
integral of motion [Yoshida, 1993]

$\frac{\varepsilon}{2}\left(q^{2}+p^{2}\right)+\frac{\tau}{2}\varepsilon
qp=const.$ (3.51)

Indeed, if we start with a natural initial condition $(0,1)$, all the points
generated iterating the mapping (3.49) lie on the conic section
$q^{2}+p^{2}+\tau qp=1$, which deviates from the exact orbit $q^{2}+p^{2}=1$
by a term of the order of $\tau$. Consequently the error in energy in (3.50)
cannot grow in time. Further, for small enough values of $\tau$, the orbit of
the system’s perturbed trajectory also takes the form (3.51), meaning that no
exponenital growth of the DV is possible. Thus, the system’s MLE will always
exhibit a regular behavior $\Lambda\propto t^{-1}$ as expected from the
harmonic oscillator (3.29). Obviously, the computer precision also has to be
taken into account on the above argumentation.

In general, any composition of the operators $e^{\tau\bm{L}_{A_{V}}}$ and
$e^{\tau\bm{L}_{B_{V}}}$ leads to a symplectic numerical scheme. A SI of order
$n$ is

$S_{n}(\tau)=\prod_{i=1}^{k}e^{\tau c_{i}\bm{L}_{A_{V}}}e^{\tau
d_{i}\bm{L}_{B_{V}}}+\mathcal{O}(\tau^{n+1}).$ (3.52)

We note that the SI (3.47) is of order $1$ as it is indicated by the subscript
‘$1$’ in $S_{1}$. The problem of constructing the SIs in (3.52) is reduced in
finding the coefficients $c_{i}$ and $d_{i}$ with $i=1,2,\ldots,k\leq n$ such
that Eq. (3.52) approximates $e^{\tau(A_{V}+\bm{L}_{B_{V}})}$ up to a certain
order of precision. This problem, was initially set and tackled in several
pioneering papers [Neri, preprint, 1988; Suzuki, 1990; Yoshida, 1990; Suzuki,
1991].

From Eqs (3.47) and (3.52) we see that the symplectic explicit Euler method
corresponds to $n=1$ and $k=1$ which leads to the trivial values of the
coefficients $c_{1}=1$ and $d_{1}=1$. For $n=2$ and $k=2$ we get the second
order SI [Hairer et al., 1993, 2006; Blanes and Casas, 2016]

$S_{2}(\tau)=e^{\tau c_{1}\bm{L}_{A_{V}}}e^{\tau d_{1}\bm{L}_{B_{V}}}e^{\tau
c_{2}\bm{L}_{A_{V}}}e^{\tau d_{2}\bm{L}_{B_{V}}}.$ (3.53)

For this scheme the coefficients $c_{1}$, $c_{2}$, $d_{1}$ and $d_{2}$ should
satisfy the order conditions

$\left\\{\begin{matrix}c_{1}+c_{2}=1\\\ d_{1}+d_{2}=1\\\ \end{matrix}\right.,$
(3.54)

whose solution is not uniquely determined. For example the choice

$c_{1}=c_{2}=\frac{1}{2},\quad d_{1}=1\quad\mbox{and}\quad d_{2}=0,$

satisfies the order conditions (3.54), leading to the SI

$S_{2}^{[1]}(\tau)=e^{\tau\frac{1}{2}\bm{L}_{A_{V}}}e^{\tau\bm{L}_{B_{V}}}e^{\tau\frac{1}{2}\bm{L}_{A_{V}}}+\mathcal{O}(\tau^{3}).$
(3.55)

This integrator is the symplectic leapfrog or Störmer-Verlet scheme [Hairer et
al., 1993]. The geometrical interpretation of the leapfrog scheme (3.55) is
depicted in Fig. 3.3 (b). We see that $\bm{Z}^{\prime}$ and $\bm{Z}$ are
related through $3$ different steps each conserving the local constants of
motion. It is worth mentioning that not all choices of the parameters $c_{1}$,
$c_{2}$, $d_{1}$ and $d_{2}$ satisfying the order conditions (3.54) lead to
SIs. For example, the choice

$c_{1}=c_{2}=d_{1}=d_{2}=\frac{1}{2},$

leads to the numerical scheme

$S_{2}^{[2]}(\tau)=e^{\tau\frac{1}{2}\bm{L}_{A_{V}}}e^{\tau\frac{1}{2}\bm{L}_{B_{V}}}e^{\tau\frac{1}{2}\bm{L}_{A_{V}}}e^{\tau\frac{1}{2}\bm{L}_{B_{V}}}+\mathcal{O}(\tau^{3}).$
(3.56)

Although this integrator is also symplectic, it is not time reversible as is
the leapfrog method $S_{2}^{[1]}$ (3.55). This means that making a time step
$\tau$ with this integrator followed by another step $-\tau$ i.e. going back
and forth (and vise versa) in time for one time step will not lead us to our
initial state:

$S_{2}^{[2]}(\tau)S_{2}^{[2]}(-\tau)\bm{Z}(t)\neq
S_{2}^{[2]}(-\tau)S_{2}^{[2]}(\tau)\bm{Z}(t)\neq\bm{Z}(t).$ (3.57)

Since the property of time reversibility is quite natural, we will restrict
our study only to SIs which have this property.

####  3.7.1 The Baker-Campbell-Hausdorff formula

Let us now see how we can use the Baker-Campbell-Hausdorff (BCH) formula to
obtain SIs. The BCH formula gives the solution for operator $W_{1}$ in the
equation

$e^{X}e^{Y}=e^{W_{1}},$ (3.58)

when $X$ and $Y$ are non-commutative elements. In particular, we get

$\displaystyle W_{1}$ $\displaystyle=$ $\displaystyle
X+Y+\frac{1}{2}[X,Y]+\frac{1}{12}([X,X,Y]+[Y,Y,X])+\frac{1}{24}[X,Y,Y,X]+\ldots$
(3.59)

where $[\cdot,\cdot]$ and $[\cdot,\cdot,\cdot]$ are commutators such that
$[X,Y]=XY-YX$, $[X,X,Y]=[X,[X,Y]]$ etc.

Setting $X=\tau\bm{L}_{A_{V}}$ and $Y=\tau\bm{L}_{B_{V}}$ in (3.58) and
(3.59), we obtain

$\displaystyle W_{1}$
$\displaystyle=\tau\bm{L}_{A_{V}}+\tau\bm{L}_{B_{V}}+\frac{1}{2}\left[\tau\bm{L}_{A_{V}},\tau\bm{L}_{B_{V}}\right]+\frac{1}{12}\left(\left[\tau\bm{L}_{A_{V}},\tau\bm{L}_{A_{V}},\tau\bm{L}_{B_{V}}\right]+\left[\tau\bm{L}_{B_{V}},\tau\bm{L}_{B_{V}},\tau\bm{L}_{A_{V}}\right]\right)+\ldots,$
$\displaystyle=\tau\bm{L}_{\mathcal{H}_{V}}+\frac{\tau^{2}}{2}\left[\bm{L}_{A_{V}},\bm{L}_{B_{V}}\right]+\frac{\tau^{3}}{12}\left(\left[\bm{L}_{A_{V}},\bm{L}_{A_{V}},\bm{L}_{B_{V}}\right]+\left[\bm{L}_{B_{V}},\bm{L}_{B_{V}},\bm{L}_{A_{V}}\right]\right)+\ldots.$
(3.60)

So neglecting the terms of order $2$ and above in $\tau$, we obtain the
symplectic explicit Euler integrator of Eq. (3.47)

$S_{1}\left(\tau\right)=e^{\tau\bm{L}_{A_{V}}}e^{\tau\bm{L}_{B_{V}}}\approx
e^{\tau\bm{L}_{\mathcal{H}_{V}}}.$ (3.61)

Let us now implement the BCH formula in order to obtain the leapfrog
integrator of Eq. (3.55), setting

$e^{X}e^{Y}e^{X}=e^{W_{1}}e^{X}=e^{W_{2}},$ (3.62)

and applying again (3.59), we get

$\displaystyle W_{2}$
$\displaystyle=2\frac{\tau}{2}\bm{L}_{A_{V}}+\tau\bm{L}_{B_{V}}+\frac{1}{6}\left(\left[\tau\bm{L}_{B_{V}},\tau\bm{L}_{B_{V}},\frac{\tau}{2}\bm{L}_{A_{V}}\right]-\left[\frac{\tau}{2}\bm{L}_{A_{V}},\frac{\tau}{2}\bm{L}_{A_{V}},\tau\bm{L}_{B_{V}}\right]\right)+\ldots,$
$\displaystyle=\tau\bm{L}_{\mathcal{H}_{V}}+\frac{\tau^{3}}{12}\left(\left[\bm{L}_{B_{V}},\bm{L}_{B_{V}},\bm{L}_{A_{V}}\right]-\frac{1}{2}\left[\bm{L}_{A_{V}},\bm{L}_{A_{V}},\bm{L}_{B_{V}}\right]\right)+\ldots.$
(3.63)

So neglecting terms of order $3$ and above in $\tau$ we obtain the leapfrog
integrator of Eq. (3.55)

$S_{2}(\tau)=e^{\frac{1}{2}\tau\bm{L}_{A_{V}}}e^{\tau\bm{L}_{B_{V}}}e^{\frac{1}{2}\tau\bm{L}_{A_{V}}}\approx
e^{\tau\bm{L}_{\mathcal{H}_{V}}}.$ (3.64)

This process can be followed in order to build higher order SIs. The error of
such integrator is defined by the power of the integration time step $\tau$
encountered in the first nonlinear term. For Eq. (3.60), we see that this term
is $\tau^{2}$ for the symplectic explicit Euler method, which according to
(3.52) means that this integrator is of order $1$, while Eq. (3.63) denotes
that the leapfrog integrator is of order $2$. The approach describes here is
based on Lie algebra concepts and is more complex as we go to higher order
integrators. The interested reader can get a very concise and complete
overview on these method in [Hairer et al., 2006; Blanes and Casas, 2016].

####  3.7.2 The Yoshida composition technique

Let us now see how we can use these ideas in order to build higher order
methods in a repeatable way. Let us rewrite Eq. (3.63) as

$W_{2}=\tau\alpha_{1}+\tau^{3}\alpha_{3}+\tau^{5}\alpha_{5}+\ldots,$ (3.65)

with $\alpha_{1}=\bm{L}_{A_{V}}+\bm{L}_{B_{V}}$,
$\alpha_{3}=\frac{1}{12}[\bm{L}_{B_{V}},\bm{L}_{B_{V}},\bm{L}_{A_{V}}]-\frac{1}{24}[\bm{L}_{A_{V}},\bm{L}_{A_{V}},\bm{L}_{B_{V}}]$,
etc. Expression (3.65) stands at the basis of what is called Yoshida
composition technique [Yoshida, 1990]. According to that technique, we can
build higher order symmetric SIs using palindromic composition of lower order
SIs. For instance, let us try to construct a SI of order $4$ of the form

$S_{4}(\tau)=S_{2}\left(x_{1}\tau\right)S_{2}\left(x_{0}\tau\right)S_{2}\left(x_{1}\tau\right),$
(3.66)

where $x_{0}$, $x_{1}$ are real coefficients we have to define. We get from
(3.65)

$W_{2}=\tau\left(x_{0}+2x_{1}\right)\alpha_{1}+\tau^{3}\left(x_{0}^{3}+2x_{1}^{3}\right)\alpha_{3}+\tau^{5}\left(x_{0}^{5}+2x_{1}^{5}\right)\alpha_{5}+\ldots.$
(3.67)

In order to achieve the desired order $4$, we want the first nonlinear term in
$\tau$ we encounter in (3.67) to be $\tau^{5}$. This means that we want to
kill the order $3$ term by setting $x_{0}^{3}+2x_{1}^{3}=0$ and at the same
time keeping the $\alpha_{1}=\bm{L}_{\mathcal{H}_{V}}$ term by setting
$x_{0}+2x_{1}=1$. Thus we obtain the conditions

$\left\\{\begin{matrix}x_{0}+2x_{1}&=1\\\
x_{0}^{3}+2x_{1}^{3}&=0\end{matrix}\right.,$ (3.68)

which can be solved analytically to give [Yoshida, 1990]

$x_{0}=-\frac{2^{1/3}}{2-2^{1/3}},\quad x_{1}=\frac{1}{2-2^{1/3}}.$ (3.69)

In fact this process can be repeated over and over again in order to create
higher order SIs. Its generalization states that from any symmetric $n$th
order SI we build a $(n+2)$th scheme using the composition (see also [Iserles
and Quispel, 2018] and references therein)

$S_{n+2}(\tau)=S_{n}(x_{1}\tau)S_{n}(x_{0}\tau)S_{n}(x_{1}\tau),$ (3.70)

with coefficients

$x_{0}=-\frac{2^{1/(2n+1)}}{2-2^{1/(2n+1)}},\quad
x_{1}=\frac{1}{2-2^{1/(2n+1)}}.$ (3.71)

Equation (3.70) is usually referred to as the Yoshida composition. Let us make
an additional remark here. The palindromic composition as in (3.70) ensures
that the SIs built using the Yoshida composition technique are always time
reversible. It is worth noting that there also exist other composition schemes
which follow the same basic principle of Yoshida’s method having some
differences on the way the Lie algebra is applied (see e.g. [Suzuki, 1991]),
based sometimes on the particular type of physical problem they were initially
created to solve (see e.g. [Ruth, 1983; Laskar and Robutel, 2001; Laskar and
Vaillant, 2019]).

####  3.7.3 Two-part split symplectic integrators

In this section we present several two-part split SIs i.e. schemes which can
be used in case the Hamiltonian function can be split into two integrable
parts like in Eq. (3.42), which have been recently used in studies of
multidimensional Hamiltonian models. These studies considered for example
models such as the DKG chain of anharmonic oscillators [Skokos et al., 2013;
Danieli et al., 2017; Senyange et al., 2018; Many Manda et al., 2020], the
DDNLS system in one [Gerlach et al., 2016; Mithun et al., 2018; Mallick et
al., 2020] and two [Many Manda et al., 2020] spatial dimensions, the Josephson
junction chain [Mithun et al., 2019], the FPUT system [Danieli et al., 2017,
2019a], the Peyrard-Bishop-Dauxois model of the DNA [Hillebrand et al., 2019],
as well as graphene models [Bountis et al., 2020; Hillebrand et al., 2020] to
name a few.

##### Symplectic integrators of order two.

  1. The leapfrog (LF) method. We explicitly presented this scheme in Sec. 3.7.1. More details about the LF scheme can be found in [Ruth, 1983] and [Hairer et al., 2006, Sect. VI.3]. This integrator has $3$ individual steps i.e. the number of applications of simple operators $e^{\tau\bm{L}_{A_{V}}}$, $e^{\tau\bm{L}_{B_{V}}}$ (see also Fig. 3.3(b)) and is given by the relation

$LF(\tau)=e^{a_{1}\tau\bm{L}_{A_{V}}}e^{b_{1}\tau\bm{L}_{B_{V}}}e^{a_{1}\tau\bm{L}_{A_{V}}},$
(3.72)

with $a_{1}=\frac{1}{2}$ and $b_{1}=1$.

  2. The $SABA_{2}/SBAB_{2}$ integrators. These SIs possess $5$ individual steps and have the following expressions

$SABA_{2}(\tau)=e^{a_{1}\tau\bm{L}_{A_{V}}}e^{b_{1}\tau\bm{L}_{B_{V}}}e^{a_{2}\tau\bm{L}_{A_{V}}}e^{b_{1}\tau\bm{L}_{B_{V}}}e^{a_{1}\tau\bm{L}_{A_{V}}},$
(3.73)

with $a_{1}=\frac{1}{2}-\frac{1}{2\sqrt{3}}$, $a_{2}=\frac{1}{\sqrt{3}}$ and
$b_{1}=\frac{1}{2}$, and

$SBAB_{2}(\tau)=e^{b_{1}\tau\bm{L}_{B_{V}}}e^{a_{1}\tau\bm{L}_{A_{V}}}e^{b_{2}\tau\bm{L}_{B_{V}}}e^{a_{1}\tau\bm{L}_{A_{V}}}e^{b_{1}\tau\bm{L}_{B_{V}}},$
(3.74)

with $a_{1}=\frac{1}{2}$, $b_{1}=\frac{1}{6}$ and $b_{2}=\frac{2}{3}$. This
class of schemes was studied in [McLachlan, 1995], where they were referred to
as the $(4,2)$ methods. In addition, they have been employed in [Laskar and
Robutel, 2001] to study Hamiltonian functions of the form $A(\bm{p})+\varrho
B(\bm{q})$ or $A(\bm{q},\bm{p})+\varrho B(\bm{q},\bm{p})$, where the
contribution $\varrho B$ was considered small since $\varrho\ll 1$.
Nevertheless these integrators have shown good performance also for cases
$\varrho=1$ (see e.g. [Skokos et al., 2009; Veksler et al., 2009; Danieli et
al., 2017; Mithun et al., 2018; Kati et al., 2020; Mallick et al., 2020]).

  3. The $ABA82$ scheme. This integrator is given by [McLachlan, 1995; Farrés et al., 2013]

$ABA82(\tau)=e^{a_{1}\tau\bm{L}_{A_{V}}}e^{b_{1}\tau\bm{L}_{B_{V}}}e^{a_{2}\tau\bm{L}_{A_{V}}}e^{b_{2}\tau\bm{L}_{B_{V}}}e^{a_{3}\tau\bm{L}_{A_{V}}}e^{b_{2}\tau\bm{L}_{B_{V}}}e^{a_{2}\tau\bm{L}_{A_{V}}}e^{b_{1}\tau\bm{L}_{B_{V}}}e^{a_{1}\tau\bm{L}_{A_{V}}},$
(3.75)

and has $9$ individual steps. The values of the coefficients $a_{i},b_{i}$,
$i=1,2,3$ in Eq. (3.75) are given in Table 2 of [Farrés et al., 2013]. It is
worth mentioning that the $ABA82$ SI was also used in [Laskar and Robutel,
2001] where it was named $SABA_{4}$.

##### Symplectic integrators of order four.

  1. The $ABA864/ABAH864$ integrators. These SIs were introduced in [Blanes et al., 2013; Farrés et al., 2013]. The first scheme has the form

$\displaystyle ABA864(\tau)=$ $\displaystyle
e^{a_{1}\tau\bm{L}_{A_{V}}}e^{b_{1}\tau\bm{L}_{B_{V}}}e^{a_{2}\tau\bm{L}_{A_{V}}}e^{b_{2}\tau\bm{L}_{B_{V}}}e^{a_{3}\tau\bm{L}_{A_{V}}}e^{b_{3}\tau\bm{L}_{B_{V}}}e^{a_{4}\tau\bm{L}_{A_{V}}}$
$\displaystyle\times
e^{b_{4}\tau\bm{L}_{B_{V}}}e^{a_{4}\tau\bm{L}_{A_{V}}}e^{b_{3}\tau\bm{L}_{B_{V}}}e^{a_{3}\tau\bm{L}_{A_{V}}}e^{b_{2}\tau\bm{L}_{B_{V}}}e^{a_{2}\tau\bm{L}_{A_{V}}}e^{b_{1}\tau\bm{L}_{B_{V}}}e^{a_{1}\tau\bm{L}_{A_{V}}},$
(3.76)

and includes $15$ individual steps. The coefficients
$a_{i},b_{i}\leavevmode\nobreak\ i=1,2,3,4$ of Eq. (3.76) are listed in Table
3 of [Blanes et al., 2013]. The form of the second SI is

$\displaystyle ABAH864(\tau)=$ $\displaystyle
e^{a_{1}\tau\bm{L}_{A_{V}}}e^{b_{1}\tau\bm{L}_{B_{V}}}e^{a_{2}\tau\bm{L}_{A_{V}}}e^{b_{2}\tau\bm{L}_{B_{V}}}e^{a_{3}\tau\bm{L}_{A_{V}}}e^{b_{3}\tau\bm{L}_{B_{V}}}e^{a_{4}\tau\bm{L}_{A_{V}}}e^{b_{4}\tau\bm{L}_{B_{V}}}e^{a_{5}\tau\bm{L}_{A_{V}}}$
$\displaystyle\times
e^{b_{4}\tau\bm{L}_{B_{V}}}e^{a_{4}\tau\bm{L}_{A_{V}}}e^{b_{3}\tau\bm{L}_{B_{V}}}e^{a_{3}\tau\bm{L}_{A_{V}}}e^{b_{2}\tau\bm{L}_{B_{V}}}e^{a_{2}\tau\bm{L}_{A_{V}}}e^{b_{1}\tau\bm{L}_{B_{V}}}e^{a_{1}\tau\bm{L}_{A_{V}}}.$
(3.77)

This integrator has $17$ steps and its constants
$a_{i},b_{i},\leavevmode\nobreak\ i=1,2,3,4,5$ are presented in Table 4 of
[Blanes et al., 2013]. The main difference between the two schemes (3.76) and
(3.77) is the fact that the $ABAH864$ integrator was constructed for
Hamiltonians for which one of their parts $A$ or $B$ is not directly
integrable and its action can be approximated by an order two SI. Thus it can
be considered as a two, or a three-part split SI.

##### Symplectic integrators of order six.

  1. The $SABA_{2}Y6$ integrator. As mentioned earlier in Sec. 3.7.2, apart from the Yoshida composition technique (3.70), there exist other composition techniques leading to higher order SIs. A composition method for creating a sixth order SI, starting from a second order scheme $S_{2}$ was presented in [Yoshida, 1990] having the form

$S_{6}(\tau)=S_{2}(w_{3}\tau)S_{2}(w_{2}\tau)S_{2}(w_{1}\tau)S_{2}(w_{0}\tau)S_{2}(w_{1}\tau)S_{2}(w_{2}\tau)S_{2}(w_{3}\tau).$
(3.78)

It is worth mentioning that the numerical scheme (3.78) was also presented in
[Kahan and Li, 1997] where it was called the $s6odr6$ method. The constants
$w_{i},\leavevmode\nobreak\ i=0,1,2,3$ of (3.78) are listed in Table 1 of
[Yoshida, 1990] corresponding to the so-called ‘solution A’ of that table.
Using Eq. (3.78) with the $SABA_{2}$ (3.73) SI in the place of $S_{2}(\tau)$,
we construct the order 6 SI $SABA_{2}Y6$ having $29$ steps.

  2. The $ABA864Y6$ scheme. Implementing the composition method (3.70) for $n=4$, with the $ABA864$ (3.77) in the role of $S_{4}(\tau)$, we construct the order $6$ integrator $ABA864Y6$, possessing $43$ steps. This integrator was used in [Senyange and Skokos, 2018].

  3. The $s9SABA_{2}6$ integrator. We also consider the composition technique [Kahan and Li, 1997]

$s9odr6b(\tau)=S_{2}(\delta_{1}\tau)S_{2}(\delta_{2}\tau)S_{2}(\delta_{3}\tau)S_{2}(\delta_{4}\tau)S_{2}(\delta_{5}\tau)S_{2}(\delta_{4}\tau)S_{2}(\delta_{3}\tau)S_{2}(\delta_{2}\tau)S_{2}(\delta_{1}\tau),$
(3.79)

with the values of the coefficients $\delta_{i}$, $i=1,2,3,4,5$ given in the
Appendix A of [Kahan and Li, 1997]. Implementing the composition method (3.79)
using the SI $SABA_{2}$ (3.73) in place of $S_{2}(\tau)$, we build the
$s9SABA_{2}6$ order 6 integrator, having $37$ steps.

####  3.7.4 Three-part split symplectic integrators

For some models such as the 1D and 2D DDNLS systems [Senyange et al., 2018;
Many Manda et al., 2020; Sun and Yu, 2020] studied here, it is preferable (or
even necessary) to split the Hamiltonian in more than two integrable parts
[Skokos et al., 2014; Gerlach et al., 2016; Laskar and Vaillant, 2019; Casas
and Escorihuela-Tomàs, 2020]. Here, we construct some basic SIs for such
systems, i.e. for Hamiltonian functions which split in three integrable parts
as

$\mathcal{H}=\mathcal{A}+\mathcal{B}+\mathcal{C}.$ (3.80)

In this case, we are looking for appropriate coefficients
$c_{i},\leavevmode\nobreak\ d_{i}$ and $e_{i}$ in order to obtain the
approximation

$e^{\tau(L_{\mathcal{A}_{V}}+L_{\mathcal{B}_{V}}+L_{\mathcal{C}_{V}})}=\prod_{i=1}^{k}e^{\tau
a_{i}L_{\mathcal{A}_{V}}}e^{\tau b_{i}L_{\mathcal{B}_{V}}}e^{\tau
c_{i}L_{\mathcal{C}_{V}}}+\mathcal{O}(\tau^{n+1}).$ (3.81)

Before proceeding, it is worth mentioning that some additional well-performing
three-part split schemes were proposed recently in [Casas and Escorihuela-
Tomàs, 2020] which were also implemented for the 1D DDNLS model.

##### Symplectic integrators of order two.

  1. The $\mathcal{A}\mathcal{B}\mathcal{C}2$ integrator. The simplest possible scheme of this type is

$\mathcal{ABC}2(\tau)=e^{a_{1}\tau\bm{L}_{\mathcal{A}_{V}}}e^{b_{1}\tau\bm{L}_{\mathcal{B}_{V}}}e^{c_{1}\tau\bm{L}_{\mathcal{C}_{V}}}e^{b_{1}\tau\bm{L}_{\mathcal{B}_{V}}}e^{a_{1}\tau\bm{L}_{\mathcal{A}_{V}}},$
(3.82)

with $a_{1}=b_{1}=1/2$ and $c_{1}=1$, which is in some sense the equivalent to
the LF (3.72) we saw in the case of two-part split. This integrator has $5$
steps and is of order $2$. Equation (3.82) can be considered as the outcome of
two successive applications of the LF (3.72), first to the Hamiltonian
$\mathcal{H}=\mathcal{A}+B$ and then to the Hamiltonian
$B=\mathcal{B}+\mathcal{C}$, considering it to be split into two parts. This
scheme has already been used in several studies (see for example [Koseleff,
1996; Chambers, 1999; Goździewski et al., 2008; Chambers, 2010; Quinn et al.,
2010; Sun and Yu, 2020]). As in the case of two-part split SIs we can
construct higher order schemes based on composition methods. Some integrators
constructed this way are reported in the rest of this section.

##### Symplectic integrators of order four.

  1. The $\mathcal{ABC}Y4$ scheme. Implementing the Yoshida composition technique (3.70) for $n=2$ with $\mathcal{ABC}2$ (3.82) in the place of the second order SI $S_{2}(\tau)$, we obtain the SI $\mathcal{ABC}Y4\left(\tau\right)$ [Koseleff, 1996] having $13$ steps. This integrator has already been used by [Gerlach et al., 2016] for studying multidimensional Hamiltonian systems. In that work, the integrator was named $\mathcal{ABC}^{4}_{[Y]}$.

  2. The $\mathcal{ABC}S4$ integrator. Implementing the composition technique proposed in [Suzuki, 1990], we build the fourth order SI

$\mathcal{ABC}S4\left(\tau\right)=\mathcal{ABC}2\left(p_{2}\tau\right)\mathcal{ABC}2\left(p_{2}\tau\right))\mathcal{ABC}2\left((1-4p_{2})\tau\right)\mathcal{ABC}2\left(p_{2}\tau\right)\mathcal{ABC}2\left(p_{2}\tau\right),$
(3.83)

with $p_{2}=\frac{1}{4-4^{1/3}}$ and $1-4p_{2}=\frac{-4^{1/3}}{4-4^{1/3}}$.
This integrator has $21$ individual steps. The $\mathcal{ABC}S4$ SI was also
considered in [Kahan and Li, 1997] where it was denoted as the $s5odr4$
method. In addition, it was also used in [Gerlach et al., 2016] where it was
called $\mathcal{ABC}^{4}_{[S]}$.

  3. The $SS864S$ scheme. As we mentioned in Sec. 3.7.3, the $ABAH864$ (3.77) SI was initially designed for Hamiltonians which can be split in three integrable parts and for which the numerical evolution of two of these parts (say $\mathcal{B}+\mathcal{C}$ in (3.80)) is approximated by the application of a second order two-part split SI. Using the $SABA_{2}$ (3.73) SI for the integration of the $\mathcal{B}+\mathcal{C}$ (3.80) part, the application of the $ABAH864$ (3.77) SI results in a method having $49$ steps which we name $SS864S$. This integrator was called $SS_{864}^{4}$ in [Skokos et al., 2014].

##### Symplectic integrators of order six.

  1. The $\mathcal{ABC}Y4Y6/\mathcal{ABC}S4Y6$ integrators. Implementing the Yoshida composition technique (3.70) for $n=4$, with the $\mathcal{ABC}Y4$ and $\mathcal{ABC}S4$ (3.83) methods in the place of $S_{4}(\tau)$, we respectively construct the $\mathcal{ABC}Y4Y6$ and $\mathcal{ABC}S4Y6$ schemes having $37$ and $49$ steps.

  2. The $\mathcal{ABC}Y6\\_A$ scheme. Applying the composition technique (3.78) with $\mathcal{A}\mathcal{B}\mathcal{C}2$ (3.82) in the role of $S_{2}(\tau)$, we create the $\mathcal{ABC}Y6\\_A(\tau)$ integrator with $29$ steps. This SI has already been implemented in [Skokos et al., 2014; Gerlach et al., 2016], where it was refereed to as $\mathcal{ABC}^{6}_{[Y]}$.

  3. The $s9\mathcal{ABC}6$ integrator. Implementing the composition technique (3.79), with $\mathcal{ABC}2$ (3.82) in the place of $S_{2}(\tau)$, we build the SI $s9\mathcal{ABC}6(\tau)$ which has $37$ steps. This integrator was also studied in [Skokos et al., 2014; Gerlach et al., 2016] where it was named $\mathcal{ABC}^{6}_{[KL]}$.

  4. The $s11\mathcal{ABC}6$ scheme. The application of the composition technique [Sofroniou and Spaletta, 2005]

$\begin{split}s11odr6(\tau)=S_{2}(\gamma_{1}\tau)S_{2}(\gamma_{2}\tau)&S_{2}(\gamma_{3}\tau)S_{2}(\gamma_{4}\tau)S_{2}(\gamma_{5}\tau)S_{2}(\gamma_{6}\tau)S_{2}(\gamma_{5}\tau)\\\
&\times
S_{2}(\gamma_{4}\tau)S_{2}(\gamma_{3}\tau)S_{2}(\gamma_{2}\tau)S_{2}(\gamma_{1}\tau),\end{split}$
(3.84)

with $\mathcal{ABC}2$ (3.82) in the role of $S_{2}(\tau)$, leads to the scheme
$s11\mathcal{ABC}6(\tau)$ which has $45$ steps. The values of the coefficients
$\gamma_{i}$, $i=1,2,\ldots,6$ of (3.84) are given in section 4.2 of
[Sofroniou and Spaletta, 2005]. The $s11\mathcal{ABC}6$ SI was called
$\mathcal{ABC}^{6}_{[SS]}$ in [Skokos et al., 2014; Gerlach et al., 2016].

##### Symplectic integrators of order eight.

  1. The $\mathcal{ABC}Y8\\_A/\mathcal{ABC}Y8\\_D$ schemes. Implementing the composition technique [Yoshida, 1990]

$\displaystyle S_{8}(\tau)=S_{2}(w_{7}\tau)S_{2}(w_{6}\tau)S_{2}(w_{5}\tau)$
$\displaystyle
S_{2}(w_{4}\tau)S_{2}(w_{3}\tau)S_{2}(w_{2}\tau)S_{2}(w_{1}\tau)S_{2}(w_{0}\tau)S_{2}(w_{1}\tau)S_{2}(w_{2}\tau)$
$\displaystyle\times
S_{2}(w_{3}\tau)S_{4}(w_{4}\tau)S_{2}(w_{5}\tau)S_{2}(w_{6}\tau)S_{2}(w_{7}\tau),$
(3.85)

with the $\mathcal{ABC}2$ (3.82) in the role of $S_{2}(\tau)$, we build an
order $8$ SI $\mathcal{ABC}Y8$ which has $61$ steps. For this scheme, we use
two sets of coefficient values $w_{i}$, $i=0,\ldots,7$ corresponding to the
so-called ‘solution A’ and ‘solution D’ in [Yoshida, 1990]. The particular
sets of these values can be found in Table 2 of [Yoshida, 1990]. Here, we
generate the SI $\mathcal{ABC}Y8\\_A(\tau)$ [$\mathcal{ABC}Y8\\_D(\tau)$]
using the coefficient values of ‘solution A’ [‘solution D’].

  2. The $s17\mathcal{ABC}8$ integrator. Implementing the composition technique [Kahan and Li, 1997]

$\displaystyle
s17odr8b(\tau)=S_{2}(\delta_{1}\tau)S_{2}(\delta_{2}\tau)\times\dots\times
S_{2}(\delta_{8}\tau)S_{2}(\delta_{9}\tau)S_{2}(\delta_{8}\tau)\dots\times
S_{2}(\delta_{2}\tau)S_{2}(\delta_{1}\tau),$ (3.86)

with the $\mathcal{ABC}2$ (3.82) SI playing the role of the second order SI
$S_{2}(\tau)$, we create the scheme $s17\mathcal{ABC}8$ having $69$ steps. The
coefficient values $\delta_{i}$, $i=1,2,\ldots,9$ are given in Appendix A of
[Kahan and Li, 1997]. The $s17\mathcal{ABC}8$ integrator was named
$\mathcal{ABC}^{8}_{[KL]}$ in [Skokos et al., 2014; Gerlach et al., 2016].

  3. The $s19\mathcal{ABC}8$ scheme. Implementing the composition technique [Sofroniou and Spaletta, 2005]

$\displaystyle
s19odr8b(\tau)=S_{2}(\gamma_{1}\tau)\mathcal{ABC}2(\gamma_{2}\tau)\times\dots\times
S_{2}(\gamma_{8}\tau)S_{2}(\gamma_{9}\tau)S_{2}(\gamma_{8}\tau)\dots\times
S_{2}(\gamma_{2}\tau)S_{2}(\gamma_{1}\tau),$ (3.87)

with the $\mathcal{ABC}2$ (3.82) integrator being the required second order SI
$S_{2}(\tau)$, we obtain the scheme $s19\mathcal{ABC}8$ having $72$ steps. In
[Sofroniou and Spaletta, 2005], the values of the coefficients $\gamma_{i}$,
$i=1,2,\ldots,9$ are listed in equation ($13$). The $s19\mathcal{ABC}8$ SI was
called $\mathcal{ABC}^{8}_{[SS]}$ in [Skokos et al., 2014; Gerlach et al.,
2016].

###  3.8 Efficient integration of the 1D and 2D DDNLS systems

We will compare the numerical efficiency of several of the integration schemes
presented in Secs. 3.7.3 and 3.7.4 for studying the DDNLS model in one (2.17)
and two (2.20) spatial dimensions. In particular, we will compute the time
evolution of initial excitations by solving the system’s equations of motion
and variational equations through numerical integration. The efficiency of the
integration methods will be judged by their ability to correctly reproduce the
system’s dynamical behaviors in the least possible needed central processing
unit (CPU) time. Since our focus is on the implementation of SIs for these
tasks, our results will be presented mainly in relation to such integration
techniques.

####  3.8.1 Two-part split symplectic integrators

Two-part split SIs can be used for the Hamiltonian functions
$\mathcal{H}_{1D}$ (2.17) and $\mathcal{H}_{2D}$ (2.20) of the 1D and 2D DDNLS
models, although as we will disucss below this is not the best approach in
integrating such systems if the number of sites is large. Let us discuss in
more detail the case of the 1D DDNLS model (2.17). Its Hamiltonian function
splits into two integrable parts $\mathcal{H}_{1D}=A+B$ as

$A=\sum_{l}\epsilon_{l}\lvert\psi_{l}\rvert^{2}+\frac{\beta}{2}\lvert\psi_{l}\rvert^{4},$
(3.88)

possessing $N$ integrals of motion, namely the quantities
$s_{l}=\lvert\psi_{l}\rvert^{2}$, and

$B=-\sum_{l}\psi_{l+1}^{\star}\psi_{l}+\psi_{l+1}\psi_{l}^{\star},$ (3.89)

where $\psi_{l}$ is a complex valued quantity representing the wave packet’s
state at site $l$. Following [Skokos et al., 2014; Gerlach et al., 2016], we
introduce real phase space variables

$\psi_{l}=\frac{1}{\sqrt{2}}\left(q_{l}+ip_{l}\right),\quad\psi_{l}^{\star}=\frac{1}{\sqrt{2}}\left(q_{l}-ip_{l}\right),$
(3.90)

with $q_{l}$ and $p_{l}$ respectively being the generalized position and
momentum of the system at site $l$. This transformation brings the A (3.88)
and B (3.89) parts of the $\mathcal{H}_{1D}$ (2.17) Hamiltonian to

$A=\sum_{l}\frac{\epsilon_{l}}{2}\left(q_{l}^{2}+p_{l}^{2}\right)+\frac{\beta}{8}\left(q_{l}^{2}+p_{l}^{2}\right)^{4},\quad
B=\sum_{l}-p_{l+1}p_{l}-q_{l+1}q_{l}.$ (3.91)

In [Gerlach et al., 2016], the efficiency of several two-part split SIs (in
particular the $LF$ (3.72), $SABA_{2}$ (3.73), $ABA864$ (3.76) and
$SABA_{2}Y6$ (3.78)) were tested mainly in comparison with three-part split
schemes. In that work, the ability of the integrator to accurately reproduce
the correct system’s dynamics in the least CPU time was investigated primarily
through the computation of the wave packet’s second moment (3.5) and
participation number (3.6). As the 1D DDNLS possesses two integrals of motion,
namely the total norm $\mathcal{S}_{1D}$ (2.19) and energy $\mathcal{H}_{1D}$
(2.17), the accuracy of the integration of the equations of motion was checked
via the computation of the relative energy error

$E_{r}(t)=\left\lvert\frac{\mathcal{H}_{1D}(t)-\mathcal{H}_{1D}(0)}{\mathcal{H}_{1D}(0)}\right\rvert,$
(3.92)

and relative norm error

$S_{r}(t)=\left\lvert\frac{\mathcal{S}_{1D}(t)-\mathcal{S}_{1D}(0)}{\mathcal{S}_{1D}(0)}\right\rvert.$
(3.93)

The main outcome of that study was that for large 1D DDNLS systems having
$N\gtrsim 70$ lattice sites, three-part split SIs (whose efficiency will be
discussed below) have to be preferred over two-part split schemes. Although
the investigation in [Gerlach et al., 2016] was limited to the numerical
integration of the system’s equations of motion it is reasonable to assume
that the simultaneous integration of the variational equations (needed for the
computation of chaos indicators like the MLE) will not change the superiority
of three-part split methods for large systems. Since in our work, we always
consider systems with $N>70$ sites, typically $N$ is of the order of several
hundreds, we will not use two-part split SIs. Nevertheless, for the sake of
completeness, we present in Sec. D.1 of Appendix D the explicit expression of
the operators appearing in SI implementations for solving both the equations
of motion and the variational equations of the 1D DDNLS system (2.17). Due to
the inefficiency of two-part split SIs for the 1D DDNLS system, such SI were
not also considered for the 2D DDNLS model.

As a final remark, let us note that the two-part split SIs are adequate for
computing the dynamics (integration of the equations of motion) and tangent
dynamics (integration of the variational equations) for several many body
systems like the 1D (2.24) and 2D (2.27) DKG models (see [Skokos et al., 2013;
Danieli et al., 2017]), the FPUT model and the XY model of Josephson junction
array [Danieli et al., 2017; Mithun et al., 2019]. The performance of several
such schemes were studied in detail in [Senyange and Skokos, 2018; Danieli et
al., 2019a]. There, it was found that the integrators $ABA864$ (3.76),
$ABAH864$ (3.77) of order $4$ and the $SABA_{2}Y6$ (3.78), $s9SABA_{2}6$
(3.79) and $ABA864Y6$ (3.70) SIs of order $6$ exhibited the best performance
when a moderate energy conservation ($E_{r}\approx 10^{-5})$ was required.

####  3.8.2 Three-part split symplectic integrators

Let us now consider the case of SIs for Hamiltonian systems which can be split
in three integrable parts. Example cases of such Hamiltonian functions include
the 1D (2.17) and 2D (2.20) DDNLS models [Skokos et al., 2014; Gerlach et al.,
2016; Danieli et al., 2019a], system’s describing the motion of a rotating
body [Laskar and Vaillant, 2019] and of charged particles under Lorentz force
fields [He et al., 2015, 2016; Casas and Escorihuela-Tomàs, 2020].

The 1D DDNLS Hamiltonian (2.17) can be split in three integrable parts

$\mathcal{A}=\sum_{l}\frac{\epsilon_{l}}{2}\left(q_{l}^{2}+p_{l}^{2}\right)+\frac{\beta}{8}\left(q_{l}^{2}+p_{l}^{2}\right)^{2},\quad\mathcal{B}=-\sum_{l}p_{l+1}p_{l},\quad\mathcal{C}=-\sum_{l}q_{l+1}q_{l},$
(3.94)

where obviously $\mathcal{A}=A$ (3.88). Similarly, the Hamiltonian function of
the 2D DDNLS system (2.20) can be split as

$\begin{split}\mathsf{A}=\sum_{l,m}\frac{\epsilon_{l,m}}{2}\left(q_{l,m}^{2}+p_{l,m}^{2}\right)+\frac{\beta}{8}&\left(q_{l,m}^{2}+p_{l,m}^{2}\right)^{2},\quad\mathsf{B}=-\sum_{l,m}p_{l+1,m}p_{l,m}+p_{l,m+1}p_{l,m},\\\
&\mathsf{C}=-\sum_{l,m}q_{l+1,m}q_{l,m}+q_{l,m+1}q_{l,m}.\end{split}$ (3.95)

The explicit expressions of the operators needed for the application of three-
part split SIs for these two systems are respectively given in sections D.2.1
and D.2.2 of Appendix D.

Let us first focus on the 1D DDNLS system (3.94) and check the ability of the
various integration schemes to accurately reproduce the system’s evolution in
its phase space, as well as its tangent dynamics. To that end, we monitor the
time evolution of the energy $E_{r}(t)$ (3.92) and norm $S_{r}(t)$ (3.93)
errors along with a basic observable like the wave packet’s second moment
$m_{2}(t)$ (3.5). For monitoring the accurate description of the system’s the
tangent dynamics we compute the motion’s ftMLE (3.21).

Figure 3.4: Results for the numerical integration of case $\text{A}1$ (see
text for details) of the 1D DDNLS system (2.17) by the second order SI
$\mathcal{ABC}2$ (3.82) for $\tau=0.0002$ [(r)ed curves], the fourth order SI
$\mathcal{ABC}Y4$ (3.7.4) for $\tau=0.0125$ [(b)lue curves] and the non-
symplectic schemes $DOP853$ (Sec. 3.6.3) [(g)reen curves] and $TIDES$ (Sec.
3.6.2) [(br)own curves]. The time evolution of (a) the relative energy error
$E_{r}(t)$ (3.92), (b) the relative norm error $S_{r}(t)$ (3.93) and (c) the
second moment $m_{2}(t)$ (3.5) of the wave packet. In (d) the norm density
distribution $\xi_{l}$ (3.1) at time $t_{f}=10^{6}$ is shown. The dashed line
in (c) guides the eye for slope $1/2$ (from [Danieli et al., 2019a]).

We consider two cases belonging respectively to the strong and weak chaos
spreading regimes, which were also investigated in [Senyange et al., 2018]:

  1. $\bullet$

Case $\text{A}1$. We excite $L=21$ central sites, with a constant norm per
site $s=1$. In addition, the parameters $W=3.5$ and $\beta=0.62$ were chosen
for the $\mathcal{H}_{1D}$ (2.17) such that its value is approximately
$0.0212$ for the disorder realization we are using.

  2. $\bullet$

Case $\text{A}2$. Here, we consider a similar setup as in case $\text{A}1$,
but the parameters of the $\mathcal{H}_{1D}$ (2.17) are $W=3$ and
$\beta=0.03$, giving for the disorder realization we selected
$\mathcal{H}_{1D}\approx 3.444$.

In both cases, the total size of the lattice is $N=1024$ sites and the initial
DV needed for the computation of the MLE is chosen such that only one site at
the center of the lattice has non zero $\delta q_{l}$ and $\delta p_{l}$
values, while its other components are set to zero. As usual, the initial DV
was normalized so that $\lVert\bm{w}(0)\rVert=1$.

Some of the result of our analysis for the case $\text{A}1$ are presented in
Fig. 3.4. In particular, we depict the evolution in time of the energy error
$E_{r}(t)$ (3.92) in Fig. 3.4(a), the norm error $S_{r}(t)$ (3.93) in Fig.
3.4(b), and the $m_{2}(t)$ (3.5) in Fig. 3.4(c), while in Fig. 3.4(d) we plot
the norm profile at the end of our integration. We see that the used SIs:
namely the $\mathcal{ABC}2$ [(r)ed curves] and $\mathcal{ABC}Y4$ [(b)lue
curves] schemes conserve quite well the value of the computed energy as the
relevant error is bounded from above by $E_{r}\approx 10^{-5}$ [Fig. 3.4(a)].
On the other hand, the general purpose integrators (non-symplectic schemes):
$DOP853$ [(g)reen curves] and $TIDES$ [(br)own curves] integrators although
having $E_{r}(t)$ values much smaller than that of the SIs, their errors are
increasing in time. From the results of Fig. 3.4 we see that the value of
$S_{r}(t)$ grows in time for all the used integrators, with $TIDES$ showing
the best behavior as its $S_{r}$ values are the smallest and show a mild
increase up to the maximum considered time $t_{f}=10^{6}$. We see that the SIs
do not conserve the second integral of motion i.e. the system’s norm, although
the increase of the related error $S_{r}(t)$ is in general slower than the
ones observed for the non-symplectic methods. We emphasize that all the used
schemes accurately reproduce the system’s dynamics as the evolution of the
wave packet’s $m_{2}$ [Fig. 3.4(c)] and its profile [Fig. 3.4(d)] are
correctly captured by all methods. We also note that similar behaviors are
seen for case A$2$, although we do not report these results here in order to
avoid repetition of similar outcomes.

The performance of seven three-part split SIs along with non-symplectic
schemes for the integration of the equations of motion of the 1D DDNLS system
was also investigated in [Skokos et al., 2014; Gerlach et al., 2016]. Here, we
want to go one step further and investigate the ability of three-part split
SIs to efficiently integrate the system’s variational equations through the
application of the ‘tangent map’ method. In Figs. 3.5(a) and (c), we present
the evolution of the ftMLEs, computed through the numerical solution of the
variational equations and see that for both sets of conditions considered
here, all schemes practically reproduce similar results. In Fig. 3.5(b) and
(d) we show the required CPU time $T_{C}$ for these computations and see that
the best performing integrator among the ones presented in Fig. 3.4 is the
$ABCY4$ scheme.

Figure 3.5: Results obtained by the numerical integration of the variational
equations of the 1D DDNLS Hamiltonian (2.17) for the initial conditions
described in cases $\text{A}1$ [panels (a) and (b)] and $\text{A}2$ [panels
(c) and (d)] (see text for details). The time evolution of, (a) and (c) the
ftMLE $\Lambda(t)$ (3.21), and (b) and (d) the required CPU time $T_{C}$ in
seconds. The integrators and the curves’ colors are the ones used in Fig. 3.4.

In Figs. 3.4 and 3.5, we showed some representative results of only four
integration schemes, although we performed a more detailed study by repeating
the same simulations with several integrators. The obtained results for cases
A$1$ and A$2$ are respectively presented in Tables 3.1 and 3.2. In these
Tables accurate simulations with $E_{r}\lesssim 10^{-9}$ are also considered.
We see that overall, both tables give qualitatively the same information. The
best performing integrators for moderate accuracy ($E_{r}\lesssim 10^{-5}$)
are the $s11\mathcal{ABC}6$ and $s9\mathcal{ABC}6$ SIs, while for higher
accuracies ($E_{r}\lesssim 10^{-9}$), the $s19\mathcal{ABC}8$ and
$s17\mathcal{ABC}8$ schemes are behaving better. Consequently, we conclude
that the $s11\mathcal{ABC}6$ and $s9\mathcal{ABC}6$ SIs are particularly well
suited for the integration of the 1D DDNLS for moderate accuracies that are
typically required in dynamical studies of this model [Skokos et al., 2009;
Bodyfelt et al., 2011a; Danieli et al., 2017, 2019b; Senyange et al., 2018;
Many Manda et al., 2020].

Table 3.1: Information on the performance of the numerical shcmes used for the integration of the equations of motion and the variational equations of the 1D DDNLS system (2.17) up to final time of $t_{f}=10^{6}$ time units for the case A$1$ (see text for details). The order $n$, the number of ‘Steps’ of each SI along with the required CPU time $T_{C}$ in seconds are reported. The integration time step $\tau$ used for each integrator fixes the relative energy error $E_{r}(t)\lesssim 10^{-5}$ and $E_{r}(t)\lesssim 10^{-9}$. All simulations were performed on a workstation using $3.00$ GHz Intel Xeon E$5-2623$ processors. All codes were written in Fortran90 language and were compiled by using the gfortran compiler (https://gcc.gnu.org/) with $-$O$3$ optimization flag. No advanced vectorization mode has been implemented.  $E_{r}\approx 10^{-5}$ | $E_{r}\approx 10^{-9}$
---|---
Integrator | $n$ | Steps | $\tau$ | $T_{\text{C}}$ | Integrator | $n$ | Steps | $\tau$ | $T_{\text{C}}$
$s11\mathcal{ABC}6$ | $6$ | $45$ | $0.115$ | $3395$ | $s19\mathcal{ABC}8$ | $8$ | $77$ | $0.09$ | $7242$
$s9\mathcal{ABC}6$ | $6$ | $37$ | $0.095$ | $3425$ | $s17\mathcal{ABC}8$ | $8$ | $69$ | $0.08$ | $7301$
$\mathcal{ABC}Y6\\_A$ | $6$ | $29$ | $0.07$ | $3720$ | $s11\mathcal{ABC}6$ | $6$ | $45$ | $0.025$ | $15692$
$SS864S$ | $4$ | $17$ | $0.05$ | $6432$ | $s9\mathcal{ABC}6$ | $6$ | $37$ | $0.02$ | $16098$
$\mathcal{ABC}Y4$ | $4$ | $13$ | $0.0125$ | $10317$ | $DOP853$ | $8$ | $\delta_{a}=10^{-16}$ | $0.05$ | $18408$
$\mathcal{ABC}S4Y6$ | $6$ | $49$ | $0.015$ | $35417$ | $\mathcal{ABC}Y8\\_D$ | $8$ | $61$ | $0.002$ | $258891$
$\mathcal{ABC}Y4Y6$ | $6$ | $37$ | $0.008$ | $40109$ | $TIDES$ | $-$ | $\delta_{a}=10^{-16}$ | $0.05$ | $419958$
$\mathcal{ABC}S4$ | $4$ | $21$ | $0.00085$ | $267911$ |  |  |  |  |
$\mathcal{ABC}2$ | $2$ | $5$ | $0.0002$ | $320581$ |  |  |  |  |
Table 3.2: Similar to Table 3.1 but for case $\text{A}2$ (see text for details) of the 1D DDNLS system (2.17) $E_{r}\approx 10^{-5}$ | $E_{r}\approx 10^{-9}$
---|---
Integrator | $n$ | Steps | $\tau$ | $T_{\text{C}}$ | Integrator | $n$ | Steps | $\tau$ | $T_{\text{C}}$
$s11\mathcal{ABC}6$ | $6$ | $45$ | $0.4$ | $1132$ | $s19\mathcal{ABC}8$ | $8$ | $77$ | $0.3$ | $2184$
$s9\mathcal{ABC}6$ | $6$ | $37$ | $0.285$ | $1147$ | $s17\mathcal{ABC}8$ | $8$ | $69$ | $0.225$ | $2632$
$\mathcal{ABC}Y6\\_A$ | $6$ | $29$ | $0.2$ | $1308$ | $s11\mathcal{ABC}6$ | $6$ | $45$ | $0.1$ | $4137$
$SS864S$ | $4$ | $17$ | $0.265$ | $1365$ | $s9\mathcal{ABC}6$ | $6$ | $37$ | $0.075$ | $4462$
$\mathcal{ABC}Y4$ | $4$ | $13$ | $0.055$ | $2354$ | $\mathcal{ABC}Y8\\_D$ | $8$ | $61$ | $0.065$ | $8528$
$\mathcal{ABC}S4Y6$ | $6$ | $49$ | $0.105$ | $4965$ | $DOP853$ | $8$ | $\delta_{a}=10^{-16}$ | $0.05$ | $14998$
$\mathcal{ABC}Y4Y6$ | $6$ | $37$ | $0.04$ | $8091$ | $TIDES$ | $-$ | $\delta_{a}=10^{-16}$ | $0.05$ | $420050$
$\mathcal{ABC}S4$ | $4$ | $21$ | $0.02$ | $9774$ |  |  |  |  |
$\mathcal{ABC}2$ | $2$ | $5$ | $0.0055$ | $11700$ |  |  |  |  |

Let us now perform similar computations for the 2D DDNLS system (2.20) and in
particular for two typical particular setups.

  1. $\bullet$

Case $\text{B}1$. The initial wave packet is located at $L\times L=7\times 7$
central sites having a uniformly distributed norm per site $s_{l,m}=1/6$. The
disorder strength is set at $W=15$, the nonlinear coefficient $\beta=6$ and
the specific disorder realization gives $\mathcal{H}_{2D}\approx 1.96$ (2.20).

  2. $\bullet$

Case $\text{B}2$. An initial wave packet of size $L\times L=1$ at the central
site with a unit norm i.e. $s_{l,m}=1$ per site for $W=16$ and $\beta=1.25$.
The specific realization of disorder parameters sets up
$\mathcal{H}_{2D}\approx 0.625$.

In both cases, a lattice of $200\times 200$ sites was considered and the
initial DV needed for the variational equations has non zero components chosen
from a random uniform distribution for only the initially excited sites and
its norm was set to unity i.e. $\lVert\bm{w}(0)\rVert=1$. The obtained
classifications for the performance of several integration schemes are shown
in Tables 3.3 and 3.4 respectively for cases B$1$ and B$2$. From these tables,
we see again that for moderate energy accuracy ($E_{r}\lesssim 10^{-5}$) the
$s9\mathcal{ABC}6$ and $s11\mathcal{ABC}6$ SIs have the best performance,
while at higher energy accuracy ($E_{r}\lesssim 10^{-9}$) both the
$s19\mathcal{ABC}8$ and $s17\mathcal{ABC}8$ SIs exhibit the best performance.

Table 3.3: Similar to Table 3.1 but for the case B$1$ of the 2D DDNLS model (2.20). $E_{r}\approx 10^{-5}$ | $E_{r}\approx 10^{-9}$
---|---
Integrator | $n$ | Steps | $\tau$ | $T_{\text{C}}$ | Integrator | $n$ | Steps | $\tau$ | $T_{\text{C}}$
$s9\mathcal{ABC}6$ | $6$ | $45$ | $0.105$ | $13914$ | $s17\mathcal{ABC}8$ | $8$ | $77$ | $0.075$ | $36528$
$s11\mathcal{ABC}6$ | $6$ | $37$ | $0.125$ | $14000$ | $s19\mathcal{ABC}8$ | $8$ | $69$ | $0.08$ | $38270$
$\mathcal{ABC}Y6\\_A$ | $6$ | $29$ | $0.08$ | $15344$ | $s9\mathcal{ABC}6$ | $6$ | $45$ | $0.0235$ | $65287$
$\mathcal{ABC}Y4$ | $4$ | $13$ | $0.025$ | $23030$ | $s11\mathcal{ABC}6$ | $6$ | $37$ | $0.0275$ | $67314$
$SS864S$ | $4$ | $17$ | $0.085$ | $23887$ | $\mathcal{ABC}Y8\\_D$ | $8$ | $61$ | $0.008$ | $140506$
$\mathcal{ABC}S4Y6$ | $6$ | $49$ | $0.03$ | $77424$ | $DOP853$ | $8$ | $\delta_{a}=10^{-16}$ | $0.05$ | $218704$
$\mathcal{ABC}Y4Y6$ | $6$ | $37$ | $0.0165$ | $87902$ |  |  |  |  |
$\mathcal{ABC}S4$ | $4$ | $21$ | $0.0065$ | $132713$ |  |  |  |  |
$\mathcal{ABC}2$ | $2$ | $5$ | $0.005$ | $157694$ |  |  |  |  |
Table 3.4: Similar to Table 3.1 but for the case B$2$ of the 2D DDNLS model (2.20) $E_{r}\approx 10^{-5}$ | $E_{r}\approx 10^{-9}$
---|---
Integrator | $n$ | Steps | $\tau$ | $T_{\text{C}}$ | Integrator | $n$ | Steps | $\tau$ | $T_{\text{C}}$
$s11\mathcal{ABC}6$ | $6$ | $45$ | $0.1515$ | $11443$ | $s19\mathcal{ABC}8$ | $8$ | $77$ | $0.135$ | $20952$
$s9\mathcal{ABC}6$ | $6$ | $37$ | $0.11$ | $13408$ | $s17\mathcal{ABC}8$ | $8$ | $69$ | $0.0875$ | $28966$
$\mathcal{ABC}Y6\\_A$ | $6$ | $29$ | $0.0775$ | $14607$ | $s11\mathcal{ABC}6$ | $6$ | $45$ | $0.0335$ | $50301$
$SS864S$ | $4$ | $17$ | $0.0915$ | $15564$ | $s9\mathcal{ABC}6$ | $6$ | $37$ | $0.024$ | $58187$
$\mathcal{ABC}Y4$ | $4$ | $13$ | $0.0215$ | $25898$ | $\mathcal{ABC}Y8\\_D$ | $8$ | $61$ | $0.009$ | $150045$
$\mathcal{ABC}S4Y6$ | $6$ | $49$ | $0.035$ | $64423$ | $DOP853$ | $8$ | $\delta_{a}=10^{-16}$ | $0.05$ | $166145$
$\mathcal{ABC}Y4Y6$ | $6$ | $37$ | $0.01375$ | $102580$ |  |  |  |  |
$\mathcal{ABC}S4$ | $4$ | $21$ | $0.005$ | $185615$ |  |  |  |  |
$\mathcal{ABC}2$ | $2$ | $5$ | $0.00155$ | $198534$ |  |  |  |  |

Let us also stress a point which we consider important. Although the results
presented in this section suggest that the $DOP853$ and $TIDES$ schemes are
not well suited for very long time integrations as they require a lot of CPU
time, they are nevertheless extremely versatile and perform quite well for
several tasks we expect an integrator to do. This fact is emphasized in Fig.
4.2 of [Hairer et al., 1993] and in the related discussion where the $DOP853$
scheme proved to be the best performing integrator among the ones considered
there. In the same spirit, the advantages of the $TIDES$ package were
discussed in [Barrio, 2005; Abad et al., 2012].

###  3.9 Summary

In this chapter, we presented the numerical methods we use in our study.
Particular emphasis was given to a plethora of SIs developed for the
integration of the equations of motion and the variational equations of the 1D
and 2D DDNLS systems. For the 1D DDNLS model, we considered symplectic schemes
based on splitting the system’s Hamiltonian into two and three integrable
parts, while for the 2D DDNLS model, only the three-part split of the
Hamiltonian function was explored. We studied the performance of the obtained
numerical SIs by comparing their efficiency to some commonly used integrators:
The RK scheme of order $8$ based on the Dormand-Prince method $DOP853$ and the
$TIDES$ package which is extracted from Taylors/Lie series expansions of the
system’s vector field. Our results suggest that SIs should be preferred for
very long time integrations as they need less CPU time to achieve a certain
level of accuracy. In the case of three-part split methods, the
$s11\mathcal{ABC}6$ and $s9\mathcal{ABC}6$ as well as the $s19\mathcal{ABC}8$
and $s17\mathcal{ABC}8$ SIs proved to be the most efficient ones for moderate
and high accuracy respectively. In addition, the three-part split approach
leads to SIs which are faster than those produced by the two-part split method
for large system sizes. This result justifies our choice of the former
approach as our used integration techniques.

The proper selection of initial excitations of the DDNLS models has been
discussed in [Kati et al., 2020] as one has to take into account the values of
both integrals of motion i.e. the total energy and total norm. Therefore,
numerical approaches for designing initial conditions for these systems in
order to obtain specific energy and norm values were also presented in detail.
We further outlined other numerical techniques we use in our study, such as
the determination of the power law exponents from non smooth numerical data
sets, the numerical characterization of the wave packet as well as the
numerical calculations of its extent through definitions of the first and
second moments, and the participation number. Furthermore, we introduced the
notion of tangent space of a trajectory in the phase space and presented
numerical techniques for the computation of observables of the orbit’s
dynamics like the MLE and the DVD.

So far we have discussed and presented the analytical approaches, as well as
numerical techniques we are using in our work. The two next chapters will be
devoted on the application of these techniques in order to investigate in
detail the characteristics of chaos of the 1D and 2D DDNLS systems.

##  Chapter 4 Chaotic dynamics of wave packet spreading in the one-dimensional
disordered nonlinear Schrödinger equation

###  4.1 Introduction

In this chapter we turn the focus to one of the main objectives of this
thesis, namely the investigation of the chaotic dynamics of wave packet
spreading in the 1D DDNLS model (2.17). Although this subject has been
considered in the past [see e.g. Mulansky et al., 2009; Skokos et al., 2013;
Antonopoulos et al., 2014, 2016, 2017], the detailed investigation of the
asymptotic behavior of chaos for the different dynamical regimes in the 1D
DDNLS system is still incomplete. Due to both the presence of disorder in the
equations of motion (2.18) and the stochastic nature of the tangent dynamics,
visible for example if we consider a small perturbation in the reduced
nonlinear stochastic ODE (2.105) approximating the system’s dynamics in the
limit of small amplitude wave packet spreading, the analytical study of the
system chaotic dynamics is very difficult. Thus, we rely on numerical
simulations for our investigations.

The chapter is organized as follows: In Sec. 4.2, we discuss some of the
computational issues mainly related to the integration of the model’s
variational equations. Section 4.3 is devoted to presenting our computations
of the MLE for the AL, weak chaos, and strong chaos dynamical regimes. In Sec.
4.4 we analyse the local characteristics of chaos. We investigate how the DVD
behaves in the linear regime and discuss the idea of the presence of chaotic
hotspots within the wave packet interior as a mechanism for promoting
spreading in the weak and strong chaos regimes. In Sec. 4.5 we discuss our
findings. More specifically, we discuss the implication of the wave packet
chaoticity for its spreading dynamics in the weak and strong chaos cases. The
results of this chapter are heavily based on findings in [Senyange et al.,
2018].

###  4.2 Some aspects related to our numerical simulations

####  4.2.1 Choice of the initial deviation vector

Let us now look more closely at the numerical aspect of the DV in our
simulations. In general the choice of the DVs is of practical importance for
the computation of the Lyapunov exponents [Skokos, 2010]. The most commonly
suggested option being initial DVs having randomly selected elements in order
to avoid any linear dependence between the coordinates of the DVs. As the
computation of the MLE only requires a single DV the practice above is not
absolutely necessary. However, most papers suggest it as good choice of DV
[see e.g. Skokos, 2010]. For low degrees of freedom Hamiltonian models such as
the ones seen in pioneering fields of chaos theory like celestial mechanics
[Hénon and Heiles, 1964; Froeschlé, 1972; Skokos, 2001] or plasma physics
[Chirikov, 1971; Wingen et al., 2005; Kominis et al., 2008; White, 2013], the
dynamics of the system takes place in a rather ‘small’ phase space such that
if the system’s orbit is chaotic, all choices of DVs for the computation of
the ftMLE converge to the MLE practically at the same rate. In systems with
many more degrees of freedom, especially the ones we are working in this
thesis which have in theory an infinite phase space dimension, the system’s
dynamics in the phase space take place in ‘tiny’ regions compared to the
overall size of the system’s phase space. Thus, one has to carefully select
the type of DVs in order to correctly capture the dynamics.

To better understand this let us consider the example in Figs. 4.1 and 4.2(a),
where we initially excite in a lattice of $N=1021$ sites the $q_{l}$ and
$p_{l}$ coordinates of $L=21$ central sites with a uniformly distributed norm
per site $s=0.3$, such that the total norm (2.19) is $\mathcal{S}_{1D}=6.3$.
Moreover, the nonlinear coefficient $\beta=1$, the disordered strength $W=4$
and the random values $\epsilon_{l}$ are chosen such that the total energy
(2.17) $\mathcal{H}_{1D}\approx 1.05$ is uniformly distributed among all the
excited sites with an energy per site $h\approx 0.05$. This initial wave
packet belongs to the weak chaos regime of the 1D DDNLS model (2.17) according
to the parameter spaces of Figs. 2.4 and 2.6. We integrate the system’s
equations of motion and variational equations using the $s11\mathcal{ABC}6$
scheme (3.84) up to $t_{f}=10^{6}$, fixing the integration time step
$\tau=0.15$ which leads to a maximal relative energy error (3.92)
$E_{r}(t)\lesssim 10^{-6}$ and relative norm error (3.92) $S_{r}(t)\lesssim
10^{-4}$. For the initial DV, we use four cases:

  1. Case DV1. We give random values to $\delta q_{l}$ and $\delta p_{l}$ of only the central site of the lattice. This leads to the spatial distribution of the DVD $\xi_{l}^{D}$ (3.25) as shown in upper left panel of Fig. 4.2(a).

  2. Case DV2. We randomly excite the $\delta q_{l}$ and $\delta p_{l}$ coordinates of $L=21$ central sites i.e. the DVD $\xi_{l}^{D}$ (3.25) has non-zero elements only in the lattice’s initially excited region as shown in the upper right panel of Fig. 4.2(a).

  3. Case DV3. We give to the $\delta q_{l}$ and $\delta p_{l}$ coordinates of the $L=21$ sites at either end of the lattice, randomly chosen values, which corresponds to the DVD $\xi_{l}^{D}$ (3.25) depicted in the lower left panel of Fig. 4.2(a)

  4. Case DV4. We excite all the elements of the DV with values drawn from a random uniform distribution. This results to a DVD $\xi_{l}^{D}$ (3.25) having non-zero elements at all lattice sites as shown in the lower right panel of Fig. 4.2(a).

In Fig. 4.1 we present the spatiotemporal evolution of the DVD, $\xi_{l}^{D}$
(3.25), for the cases DV1 in Fig. 4.1(a), DV2 in Fig.4.1(b), DV3 in Fig.
4.1(c) and DV4 in Fig. 4.1(d). A transient time where the dynamical behaviors
of each DVD differs is clearly visible. For DV1 and DV2, a timid expansion of
the DVD is observed which do not spread over the whole lattice, but always
remain within the area of the wave packet extent where the dynamics of the
system take place. This is further supported by the observations of the
spatiotemporal evolution of DV3 and DV4, where the tangent dynamics has a
tendency to enhance the magnitude of the coordinates of the DV inside the
region where the motion of the wave packet occurs and the values of the
coordinates of the DV outside the wave packet die out. In the case of DV3,
this transition happens around $\log_{10}t\approx 3.5$ [see Fig. 4.1(c)],
while for the case DV4 it happens at $\log_{10}t\approx 2$ [see Fig. 4.1(d)].

During these transient times, the ftMLEs in each DV case show different
evolutions as seen in Fig. 4.2(b). Indeed, the ftMLEs obtained with DV1 and
DV2 display larger values compared to the other DVs (DV3 and DV4). This is due
to the fact that these DVs evolve within the active domain, whose evolution is
controlled by nonlinear deterministic equations (2.18). Consequently, a
deterministic chaotic dynamics is visible from the beginning of their
evolutions which explain the large values of the ftMLEs. Next is the ftMLE of
DV4 which is at its peak two order of magnitudes smaller than the values
observed for the ftMLEs of DV1 and DV2. The reason for such a measurement is
the following: For time $t$ such that $\log_{10}t\lessapprox 2$ in the early
evolution of DV4, a large part of the DV coordinates evolve in the regular
background outside the system’s excited region. Thus the contribution of the
regular part is non-negligible compared to the contribution of the excited
region of the lattice which is chaotic. This leads to the smaller value of the
ftMLE observed in Fig. 4.2(b) compared to the ones obtained from DV1 and
DV$2$. In the case of DV3, when $\log_{10}t\lessapprox 3.5$, all the non-zero
coordinates of the DV evolve only into the regular part of the lattice which
therefore solely contribute to the value of the ftMLE which indicates
$\Lambda(t)\approx 10^{-10}$ at its lowest pick i.e. an almost zero ftMLE.

For longer times, well above the largest transient time of DV3, all the DVs
converge to a consistent behavior as all the asymptotic spatiotemporal
evolution of the DVDs $\xi_{l}^{D}$ (3.25) are practically the same in Fig.
4.2. According to the theory of chaos [see e.g. Skokos, 2010; Skokos and
Manos, 2016], this asymptotic behavior of the DVs in Fig. 4.1 corresponds to
the spatial alignment of the tangent dynamics along the most chaotic direction
of the system’s reference orbit. Thus, the average exponential rate of growth
(or shrinking) of the magnitude of this DV measures the system’s MLE. This is
what is observed in Fig. 4.2(b), where for time $t$ above the transient time,
all the ftMLEs practically indicate the same value.

Figure 4.1: Spatiotemporal evolution of the DVD $\xi_{l}^{D}$ (3.25) profile
for different initial DVs. The panel (a) corresponds to DV1, the panel (b) to
DV2, the panel (c) to DV3 and the panel (d) to DV4 (see text for details). The
color scales at the top of the figure are used for coloring the lattice sites
according to their $\log_{10}\xi_{l}^{D}$ values.

Figure 4.2: (a) Initial profiles $\xi_{l}^{D}$ (3.1) of the DVDs used for
obtaining the results of Fig. 4.1. The positions of each panel corresponds to
the arrangement seen in Fig. 4.1. (b) Time evolution of the ftMLE,
$\Lambda(t)$ (3.21) for the DVs of panel (a). The colors of the curves
correspond to the ones in panel (a).

Based on this experiment, we can decide which type of DVs is the best from a
practical point of view. Indeed, our objective is to be able to retrieve the
asymptotic behavior of the tangent dynamics as fast as possible since we are
constrained by computational time limitations. Thus, DV3 and DV4, which
possess the largest transient times are bad choices for our purpose. This is
worsened by the fact that for lattice size $N>1021$ (with the size of the
excited region $L$ almost constant), the regular background outside the
excited central part of the lattice becomes larger, therefore longer transient
times are expected for the DVs similar to DV3 and DV4 as the elements of the
tangent dynamics outside the wave packet evolution take larger times to vanish
(only the elements of the DV inside the lattice’s excited part contribute to
the MLE). On the other hand the cases DV1 and DV2 are good choices since they
are located inside the central excited part of the lattice where the dynamics
take place. Consequently, they take less time to converge to the most chaotic
direction whose non zero elements of the DV are always located inside the wave
packet as indicated by the evolution of the DVDs in Fig. 4.1. For the rest of
this chapter, unless otherwise stated, our default initial DV is the DV1.

####  4.2.2 The dynamic lattice expansion

Because the initially excited sites are located at the central part of a
lattice of large size, in our simulations we can consider a smaller subset of
the lattice at the beginning of our simulations and expand it when ever it is
necessary. We refer to this approach as the dynamic lattice simulation which
has been successfully implemented for both the 1D (2.17) and 2D (2.20) DDNLS
models and DKG systems (Eqs. (2.24) and (2.27)) in [Senyange et al., 2018;
Many Manda et al., 2020] in contrast to the fixed lattice method used in
several studies on the subject [Flach et al., 2009; Skokos and Flach, 2010;
Laptyeva et al., 2010, 2012; Kati et al., 2020]. Let us note that this is not
particularly unusual as even in the theory (see Chap. 2 and references
therein), the infinite lattice size i.e. $N\rightarrow\infty$ in
$\mathcal{H}_{1D}$ (2.17) is often substituted by a large but finite value of
$N$. Then, we have to define a condition which controls our lattice expansion.
For the 1D DDNLS model (2.17), we monitor the value of the average norm per
site $\sum_{l\in\mathcal{D}}s_{l}/\mathcal{D}$ of the last $\mathcal{D}$ sites
at both edges of the ‘active’ (dynamic) lattice with size
$N_{\mathcal{D},\delta_{\mathcal{D}}}$ sites. Then, whenever

$\frac{1}{\mathcal{D}_{L}}\sum_{l\in\mathcal{D}_{L}}s_{l}>\delta_{\mathcal{D}_{L}},\mbox{
or
}\frac{1}{\mathcal{D}_{R}}\sum_{l\in\mathcal{D}_{R}}s_{l}>\delta_{\mathcal{D}_{R}},$
(4.1)

the lattice is expanded at both ends (in order to maintain the symmetry with
respect to the initial central excitation) by $\mathcal{N}$ sites. Note that
in Eq. (4.1), $\mathcal{D}_{L}$ and $\mathcal{D}_{R}$ describe the number of
sites at the left and right edges respectively (we often assume
$\mathcal{D}_{L}=\mathcal{D}_{R}=\mathcal{D}$). In addition,
$\delta_{\mathcal{D}_{L}}=\delta_{\mathcal{D}_{R}}=\delta_{\mathcal{D}}$ is
the threshold above which we expand our current dynamic lattice at the edges.
It is worth mentioning that $\mathcal{N}>\mathcal{D}$.

Figure 4.3: Heatmaps of the spatiotemporal evolutions of the normalized norm
distribution $\xi_{l}$ (3.25) for an initial condition belonging to the strong
chaos spreading regime (see text for details). We excited $L=21$ central sites
on a lattice of $N=1021$ sites with the same norm per site $s=1$ so that the
total norm (2.19) $\mathcal{S}_{1D}=21$. The random values $\epsilon_{l}$,
$1\leq l\leq N$, give a total energy (2.17) $\mathcal{H}_{1D}\approx 3.21$,
using the parameters $W=4$ and $\beta=0.72$. The panel (a) corresponds to the
fixed lattice size as we see the residual norm i.e. $s\sim 10^{-12}$ (light
yellow in the heatmap) spreads over the entire lattice. The panel (b)
corresponds to the lattice expansion approach as the residual norm is
concentrated around the wave packet (darker region in the heatmap). The white
area, corresponds to (empty) lattice sites with no norm $s$. The color scales
at the top of the figure are used for coloring the lattice sites according to
the magnitude of their $\log_{10}\xi_{l}$ values.

To be more specific, in Fig. 4.3, we plot the evolution of the normalized norm
distribution $\xi_{l}$ (3.1) on a lattice of $N=1021$ sites using at time
$t=0$, $L=21$ centrally excited sites such that each site receives the same
constant norm per site $s=1$, leading to the total norm (2.19)
$\mathcal{S}_{1D}=21$. Further, the values of $\epsilon_{l}$ are chosen in a
way that the Hamiltonian (2.17) $\mathcal{H}_{1D}\approx 3.21$, with
parameters $W=4$ and $\beta=0.72$. This initial condition belongs to the
strong chaos spreading regime according to Figs. 2.4 and 2.6 of Chap. 2. We
integrate the system’s phase space and tangent dynamics using the
$s11\mathcal{ABC}6$ SI (3.84), with time step $\tau=0.2$ up to the final
integration time $t_{f}=10^{7}$. In doing so, we obtain a relative energy
error (3.92) $E_{r}(t)\lesssim 10^{-5}$ and a relative norm error (3.93)
$S_{r}(t)\lesssim 10^{-3}$. In Fig. 4.3(a), we plot the spatiotemporal wave
packet evolution $\xi_{l}(t)$ for the fixed lattice size system, along with
the one of the dynamic lattice implementation [see Fig. 4.3(b)], starting from
$N_{\mathcal{D},\delta_{\mathcal{D}}}=321$ at time $t=0$. In addition, we used
as threshold $\delta_{\mathcal{D}}=10^{-12}$. We monitor $\mathcal{D}=10$
points at the edges and increase each edges by $\mathcal{N}=50$ sites each
time that the wave packet has reached the borders of the lattice i.e. whether
the condition (4.1) is verified.

In Fig. 4.3(b), the dynamic lattice size
$N_{\mathcal{D},\delta_{\mathcal{D}}}=321$ sites at time $t=0$, was increased
by $2\mathcal{N}=100$ at time $\log_{10}t\approx 3.64$ to
$N_{\mathcal{D},\delta_{\mathcal{D}}}=421$ sites, at $\log_{10}t\approx 4.72$
to $N_{\mathcal{D},\delta_{\mathcal{D}}}=521$ sites, at $\log_{10}t\approx
5.57$ to $N_{\mathcal{D},\delta_{\mathcal{D}}}=621$ sites, at
$\log_{10}t\approx 5.75$ to $N_{\mathcal{D},\delta_{\mathcal{D}}}=721$ sites,
at $\log_{10}t\approx 6.34$ to $N_{\mathcal{D},\delta_{\mathcal{D}}}=821$
sites, and at time $\log_{10}t\approx 6.86$ reaching a total lattice size at
the final time $t_{f}=10^{7}$ of the integration
$N_{\mathcal{D},\delta_{\mathcal{D}}}=921$ sites. Consequently, after
increasing the dynamic lattice $6$ times, we obtain the results shown in Fig.
4.3(b) where the white areas at the edges indicate empty space i.e. lattice
sites with norm per site $s=0$. Further, it clearly appears from Figs. 4.3(a)
and (b) that the spatiotemporal evolutions of the normalized norm distribution
$\xi_{l}(t)$ (3.1) are practically the same. It is worth mentioning that the
initial lattice size $N_{\mathcal{D},\delta_{\mathcal{D}}}=321$ is much larger
than the NM average localization volume (2.61) $\overline{V}\approx 21$ sites
in this context.

In order to quantify the qualitative observations of the evolution of the norm
distribution $\xi_{l}$ (3.1) in Fig. 4.3, we plot the time evolution of the
$m_{2}(t)$ (3.5) [panel (a) of Fig. 4.4] and $P(t)$ (3.6) [inset of panel (a)
in Fig. 4.4] of $\xi_{l}$ (3.1) for the fixed [(g)reen curve] and dynamic
[(b)lue curves] lattice simulations. Both the $m_{2}$ and $P$ of $\xi_{l}$ are
practically overlapping up to the largest simulation time $t_{f}=10^{7}$. As
in this work we are also interested in the evolution of the variational
equations, we depict the time evolution of the second moment $m_{2}^{D}(t)$
(3.5) and participation number $P^{D}(t)$ (3.6) of the DVD $\xi_{l}^{D}$
(3.25) as well as the values of the ftMLE $\Lambda(t)$ (3.21) in Figs.
4.4(b-c) for both the fixed [(g)reen curves] and dynamic [(b)lue curves]
lattices computations of Fig. 4.3. We notice a good correspondence in
magnitude of the values of the $m_{2}^{D}$ and $P^{D}$ between the two lattice
types as they are randomly oscillating in the same range of $m_{2}^{D}$ and
$P^{D}$ values. Further, the chaoticity of the reference orbit depicted in the
initial condition used in Fig. 4.3 is well reported by the two type of
simulations as the curves of the ftMLEs practically overlap. Finally, we
analyze the CPU time $T_{C}$ necessary to reach the final time of integration
$t_{f}=10^{7}$ [Fig. 4.4(d)]. Clearly, the dynamic lattice appears
advantageous as the final CPU time is noticeably smaller on the logarithmic
scale. More specifically, we calculate the gain factor

$g=\frac{T_{C}^{(F)}-T_{C}^{(D)}}{T_{C}^{(F)}},$ (4.2)

with $T_{C}^{(F)}$ and $T_{C}^{(D)}$ being respectively the CPU time needed
for the fixed and dynamics lattice simulations. In the example above, we
obtain an improvement of approximately $20\%$, using $T_{C}^{(F)}\approx
5.5\leavevmode\nobreak\ \text{hours}$ and $T_{C}^{(D)}\approx
4.43\leavevmode\nobreak\ \text{hours}$.

Figure 4.4: Results of the time evolution of the observables of the phase
space and tangent dynamics of the case used in Fig. 4.3. The panel (a)
presents the second moment $m_{2}(t)$ (3.5) and the participation number
$P(t)$ (3.6) [inset of panel (a)] of the norm distribution $\xi_{l}$ (3.1).
The panel (b) shows the second moment $m_{2}^{D}(t)$ and participation number
$P^{D}(t)$ [inset of panel (b)] of the DVD $\xi_{l}^{D}$ (3.25). In panel (c)
we see the ftMLEs $\Lambda(t)$ (3.21) and in panel (d) the required CPU time.
The (g)reen curves are for the fixed and the (b)lue curves the dynamic lattice
simulations.

It is worth commenting a bit more on the bent nature of the curve of the CPU
time for the dynamic lattice simulation. For the fixed lattice simulation, a
fairly straight line is seen on the dependence of $\log_{10}T_{C}$ on the
logarithmic time i.e. $\log_{10}t$. This is due to the fact that the lattice
size remains fixed along the integration of the system’s orbit. Thus the
number of basic computer operations (addition, multiplication, and so on)
necessary to evolve the equations of motion and the variational equations, as
well as to compute all the observables of the phase space and tangent
dynamics, remains constant for each integration time step $\tau$. So the
constant CPU time necessary for a single integration time step $\tau$
evolution adds up to the cumulative CPU time $\log_{10}T_{C}$ in the constant
rate (or slope) seen on the (g)reen curve of Fig. 4.4(d). This is not the case
for the dynamic lattice simulation. Indeed, the number of lattice sites is
constantly increased, leading to an augmentation of the number of basic
computer operations for the time evolution of the equations of motion and
variational equations along with the calculation of observables of the phase
space and tangent dynamics. Consequently, the CPU time needed for a single
integration time step $\tau$ increases in time leading to a constant increase
of the slope of the (b)lue curve in Fig. 4.4(d) which create the bending
effects of the $\log_{10}T_{C}$ curve. Ultimately, there exists a time $t$
where the CPU time of both lattice simulations overlap and no gain from using
the dynamic lattice is visible.

Figure 4.5: Similar to Fig. 4.4, but only the $m_{2}(t)$ (3.5) [panel (a)],
the ftMLE $\Lambda(t)$ (3.21) [inset of panel (a)] and the CPU time $T_{C}(t)$
[panel (b)] are reported. The curve coloring corresponds to
$\delta_{\mathcal{D}}=10^{-8}$ [(b)lue], $\delta_{\mathcal{D}}=10^{-10}$
[(g)reen], $\delta_{\mathcal{D}}=10^{-12}$ [(r)ed] and the fixed lattice
simulation [(c)yan]. See text for details.

Let us now discuss what happens when we change the value of the threshold
$\delta_{\mathcal{D}}$ in Eq. (4.1). In Fig. 4.5, we perform simulations with
the initial condition used in the example of Fig. 4.3. We use different
parameter values of the norm per site threshold at the edge points
$\delta_{\mathcal{D}}$ ranging from $10^{-8}$ to $10^{-12}$. The result is
depicted in Fig. 4.5. There, we plot the time evolution of $m_{2}(t)$ (3.5) of
$\xi_{l}$ (3.1) as observable of the evolution of the equations of motion [see
Fig. 4.5(a)] and the ftMLE $\Lambda(t)$ (3.21) [see inset of Fig. 4.5(a)] for
of the tangent dynamics with $\delta_{\mathcal{D}}=10^{-8}$ [(b)lue curves],
$\delta_{\mathcal{D}}=10^{-10}$ [(g)reen curves],
$\delta_{\mathcal{D}}=10^{-12}$ [(r)ed curves] and the fixed lattice
computation [(c)yan curves]. The evolutions of the $m_{2}(t)$ (3.5) and the
$\Lambda(t)$ (3.21) practically overlap in all used cases, indicating a
tendency to follow the same power law behaviors. In Fig. 4.5 (b), we depict
the time evolution of the CPU time $T_{C}$ needed for these simulations. Once
again the fixed lattice takes a longer time to complete the integration
compared to the dynamic lattice. In addition, dynamic lattice simulations with
larger norm per site threshold $\delta_{\mathcal{D}}$ for the control of the
endpoints take less time. Overall, a gain of $g\approx 27\%$, $24\%$, $15\%$
respectively when $\delta_{\mathcal{D}}=10^{-8}$, $10^{-10}$ and $10^{-12}$
are observed. The closeness of the gains in the case
$\delta_{\mathcal{D}}=10^{-8}$ and $\delta_{\mathcal{D}}=10^{-10}$ is
explained by the fact that the norm flows rather quickly at the boundaries, as
the case used here belongs to the strong chaos spreading regime. A larger
difference appears between the different $\delta_{\mathcal{D}}$ values when
initial conditions are taken from the weak chaos spreading regime.

From what has been observed above, the overall idea of the dynamic lattices is
to find the optimal lattice size for the numerical evolution of the
Hamiltonian models (in the present context $\mathcal{H}_{1D}$ (2.17) and
$\mathcal{H}_{2D}$ (2.20)) for a specific set of parameter values $W$,
$\beta$, $\mathcal{H}_{1D}$ (2.17), and $\mathcal{S}_{1D}$ (3.1) realizing a
specific level of accuracy $\delta_{\mathcal{D}}$ in capturing the evolution
of a central initial wave packet excitation. Of course, the smaller
$\delta_{\mathcal{D}}$ is the larger lattice size we are going to need. In
practice, our choice of $\delta_{\mathcal{D}}$ depends on the final values of
the average norm per site attained in our longest simulations. In the case of
the 1D DDNLS model, we are reaching a wave packet width $P_{w}\propto 10^{3}$
sites. Since the constant total norm $\mathcal{S}_{1D}\propto 10$ is shared
among all the excited sites, the final avarage norm per site

$s_{f}\propto\frac{\mathcal{S}_{1D}}{P_{w}}\sim 10^{-2}-10^{-3},$ (4.3)

such that using

$\delta_{\mathcal{D}}\sim 10^{-8},$ (4.4)

we capture the wave packet evolution with a rather good accuracy as shown in
Fig. (4.4) up to our largest simulation time $t_{f}=10^{8}$.

In order to avoid repetition, in all our numerical simulations, we consider
for the 1D DDNLS system (2.17) initially localized wave packet excitations of
size $L$ at the center of a lattice whose maximal size reach in some cases
$N\approx 5000$ sites in order to exclude effects from the lattice edges (i.e.
in order to secure that the evolved wave packet does not reach the fixed
edges). We give to each excited site a single oscillator frequency shift (or
rescaled norm per site) $x=\beta s$ (2.68), where $s$ is the constant norm at
each site so that the total norm (2.19) is $\mathcal{S}_{1D}=\sum_{l\in
L}s=Ls$ (or $\mathcal{S}_{1D}=Lx/\beta$, $\beta$ being the nonlinear
coefficient of the Hamiltonian $\mathcal{H}_{1D}$ (2.17)). As mentioned in
Sec. (4.2.1), for the initial DV $\bm{w}(0)$, we choose to excite only the
$\delta q_{l}$ and $\delta p_{l}$ values of the site at the centre the lattice
(it corresponds to DV1 in Sec. 4.2.1) and ensure that the norm
$\lVert\bm{w}(0)\rVert=1$. Further we integrate the equations of motion and
the variational equations using the SI $s11\mathcal{ABC}6$ (3.84) of order $6$
with an integration time step $\tau\approx 0.15-0.25$, always ensuring that
the maximal relative energy (3.92) $E_{r}$ and norm (3.1) $S_{r}$ errors
remain below $10^{-3}$. Our typical final integration time is $t_{f}\approx
10^{7}-10^{8}$ time units. For all the considered realizations with $\beta>0$
in the $\mathcal{H}_{1D}$ (2.17), we set the same random parameter
$\epsilon_{l}=0$ at the lattice’s centre. This ensures that all simulations of
single (the central) site excitation ($L=1$) results in the same energy value
for the same initial norm per site $s_{l}$. On the other hand, in case of
multi-site ($L>1$) initial excitation, all the active sites receive the same
value of the norm per site $s_{l}$ (2.84) and we take particular care so the
value of the corresponding $\mathcal{H}_{1D}$ (2.17) is in a region belonging
to the system’s thermal phase [Flach, 2016; Kati et al., 2020] (see also Sec.
3.3).

###  4.3  The maximum Lyapunov exponent

####  4.3.1 The linear regime

Here, we look at the linear regime in the 1D DDNLS model (2.17). Its
Hamiltonian corresponds to $\mathcal{H}_{L}$ (2.1), or equivalently
$\mathcal{H}_{1D}$ (2.17) with the nonlinear coefficient $\beta=0$. The theory
in Sec. 2.1 tells us that no diffusion of initially localized excitations is
present as AL prevails. In this situation, a saturation of the $m_{2}$ and $P$
is observed once the wave packet extent reaches the border of the localization
volume $\overline{V}$ (2.14). This corresponds to equilibrium non-ergodic
(regular) dynamics in the system’s phase space where the trajectories belong
to a torus associated with energy level $E$ of the periodic motion
$\psi_{l}=A_{l}e^{-iEt}$ at each site $l$ (2.4) of the chain of oscillators.
Thus no chaos is present and one expects the ftMLE $\Lambda(t)\propto t^{-1}$.

Figure 4.6: The time evolution of the second moment $m_{2}(t)$ (3.5) [(b)lue
curve in panel (a)], the participation number $P(t)$ (3.6) [(g)reen curve in
panel (a)] of the norm density distribution $\xi_{l}$ (3.1). The related ftMLE
$\Lambda(t)$ is shown in panel (b) for a representative of case A$1_{1D}$ (see
text for details). The black lines in panel (a) guide the eye for zero slopes
along the observable’s mean values $\approx 32$ for the $m_{2}$ [upper black
line] and $\approx 8$ [lower black line] for $P$ respectively. In addition,
the one in panel (b) guides the eye for slope $-1$.

In Fig. 4.6, we plot the results of the time evolution of the $m_{2}(t)$ (3.5)
and $P(t)$ (3.6) of the norm distribution $\xi_{l}$ (3.1) as observables of
the evolution of the equations of motion, and the ftMLE $\Lambda(t)$ (3.21)
for the variational equations for a single site $L=1$ excitation (we label
this as case A$1_{1D}$ for future reference) with a unit norm per site i.e.
$s=1$ at the center of a lattice of size $N=1021$ sites. The energy
$\mathcal{H}_{1D}=1$ is obtained from the disorder realization we used, with
parameters $W=4$ and $\beta=0$. In panel (a) of Fig. 4.6, we depict the
obtained evolution of $m_{2}(t)$ (3.5) [(blue) curve], and $P(t)$ (3.6)
[(g)reen curve]. We see that the $m_{2}(t)$ and $P(t)$ saturate as their
values randomly oscillate around their the means $\overline{m}_{2}\approx 32$
[upper black line in Fig. 4.6(a)] and $\overline{P}\approx 8$ [lower black
line in Fig. 4.6(a)] rather quickly ($\log_{10}t\approx 1$). In panel (b) of
Fig. 4.6, we show the evolution of the ftMLE $\Lambda(t)$ (3.21). A clear
behavior $\Lambda(t)\propto t^{-1}$ is visible as the (b)lue curve in Fig.
4.6(b) has the tendency to lay along the black straight line whose slope is
$-1$. This was also reported for the DKG of harmonic oscillators [see e. g.
Skokos et al., 2013]. Next, we inch up the nonlinearity $\beta$ in
$\mathcal{H}_{1D}$ (2.17) and numerically investigate the fate of the wave
packet’s localization.

####  4.3.2 The weak chaos regime

So far, the results predicted in Chap. 2 regarding the spreading of a
nonlinear localized excitation in the 1D DDNLS model (2.17) appear to be
accurately reproduced by numerical simulations. Indeed, the first numerical
experiments performed in [Flach et al., 2009] were a complete success in
confirming these predictions. Furthermore, the observation of two distinct
power laws of wave packet expansion during the spreading was extensively
studied in [Skokos et al., 2009; Flach, 2010; Laptyeva et al., 2010; Bodyfelt
et al., 2011a, b; Laptyeva, 2013] with the crossover from the strong to weak
chaos spreading characteristics also reported in [Laptyeva et al., 2010]. Here
we present the results of the numerical investigation of the chaotic dynamics
of the wave packet spreading within the weak chaos regime. In particular we
focus on the following cases:

  1. $\bullet$

Case $\text{W}1_{\text{1D}}$. $W=3$, $\beta=0.03$, $L=21$ and $x=0.03$.

  2. $\bullet$

Case $\text{W}2_{\text{1D}}$. $W=3$, $\beta=0.6$, $L=1$ and $x=0.6$.

  3. $\bullet$

Case $\text{W}3_{\text{1D}}$. $W=4$, $\beta=1$, $L=1$ and $x=1$.

  4. $\bullet$

Case $\text{W}4_{\text{1D}}$. $W=4$, $\beta=0.04$, $L=21$ and $x=0.04$.

We can justify the choice of these parameters based on the first condition of
Eq. (2.56). Indeed, a multisite excitation is expected to evolve in the the
weak chaos regime whenever

$x\frac{L}{\overline{V}}<d.$ (4.5)

where $x=\beta s$ (2.68) is the rescaled average norm per site, $\overline{V}$
(2.14) is the average localization volume and $d$ (2.16) is the average
spacing between frequencies in a NM volume. In Table 4.1, we give the various
values of $d$ (2.16), $\overline{V}$ (2.14), and $\Delta$ (2.15) the frequency
bandwidth of NMs for the different disorder strength $W$ values of the 1D
DDNLS Hamiltonian model $\mathcal{H}_{1D}$ (2.17) we are using in this
chapter. Applying the condition (4.5) in case W$1_{1D}$, we see that
$xL/\overline{V}=0.03\cdot 21/37\approx 0.017$ which is one order of magnitude
smaller than the value of $d=0.2$ indicated in Table 4.1. For the case
W$4_{1D}$, the same calculation leads to the number $0.04\cdot 21/21$ which is
again much smaller than $d=0.39$ of $W=4$ in Table 4.1. In the case of single
excitation wave packet, it is the first condition of inequality (2.57)

$x<\Delta,$ (4.6)

which governs the weak chaos spreading dynamics appearance. For the case
W$2_{1D}$ the inequality (4.6) gives $0.6<7$, while in the case W$3_{1D}$ we
obtain $1<8$, both of which are true. It is worth mentioning that the cases
$\text{W}1_{\text{1D}}$, $\text{W}2_{\text{1D}}$, $\text{W}3_{\text{1D}}$ and
$\text{W}4_{\text{1D}}$ have already been used in several works [Laptyeva et
al., 2010; Bodyfelt et al., 2011a; Senyange and Skokos, 2018; Senyange et al.,
2018]. In all the cases, we are averaging our observables of the phase space
and tangent dynamics over $100$ disorder realizations in order to obtain a
robust representation of the system’s dynamical behavior which takes into
account several trajectories of the same total norm $\mathcal{S}_{1D}$ (3.1).

NM properties | $W=3$ | $W=3.5$ | $W=4$
---|---|---|---
$\overline{V}$ | $37$ | $27$ | $21$
$\Delta$ | $7$ | $7.5$ | $8$
$d=\Delta/\overline{V}$ | $0.2$ | $0.28$ | $0.39$
Table 4.1: Average localization volume $\overline{V}$ (2.14), frequency
bandwidth $\Delta$ (2.15) and average spacing between NM $d$ (2.16) in case of
a multisite excitation of width $L=\overline{V}$ sites at the center of the 1D
DDNLS lattice system (2.17).

Figure 4.7: Time dependence of the second moment $m_{2}(t)$ (3.5) [panel (a)],
the participation number $P(t)$ (3.6) [panel (b)] of $\xi_{l}$ (3.1), the
ftMLE $\Lambda(t)$ (3.21) (c) and $\alpha_{\Lambda}(t)$ (3.27) [panel (d)] for
the weak chaos cases W$1_{1D}$ [(b)lue], W$2_{1D}$ [(g)reen], W$3_{1D}$
[(r)ed] and W$4_{1D}$ [(c)yan] (see text for details). The black dashed lines
guide the eye for slopes $1/3$ [panel (a)], $1/6$ [panel (b)] and $-0.25$
[panel (c)]. The horizontal dashed black line in panel (d) indicates the value
$\alpha_{\Lambda}=-0.25$.

The fact that all these cases belong to weak chaos becomes evident in Fig.
4.7(a) and (b), where we plot the time evolution of the second moment
$m_{2}(t)$ (3.5) and the participation number $P(t)$ (3.6) of the norm density
distribution $\xi_{l}$ (3.1) respectively. For both quantities, a subdiffusive
spreading is noticeably well fitted by the monotonic asymptotic power laws
$m_{2}(t)\propto t^{\alpha_{m}}$ [see black line in Fig. 4.7(a)] with
$\alpha_{m}=\frac{1}{3}$ and $P(t)\propto t^{\alpha_{P}}$ [see black line in
Fig. 4.7(b)] where $\alpha_{P}=\frac{1}{6}$ in agreement with previous works
[Flach et al., 2009; Skokos et al., 2009; Laptyeva et al., 2010; Flach, 2010;
Bodyfelt et al., 2011a]. Further, in Fig. 4.7(c) we present the time evolution
of the ftMLE, $\Lambda(t)$ (3.21) which is the primary quantity of our
interest. An asymptotic decay following the power law $\Lambda\propto
t^{\alpha_{\Lambda}}$ ($\alpha_{\Lambda}<0$) is observed for all cases. In
Fig. 4.7(d), we present the time evolution of the related power law exponent
$\alpha_{\Lambda}(t)$ computed following the techniques presented in Sec. 3.5.
The dependence on time of the power law exponent values $\alpha_{\Lambda}(t)$
shows a clear tendency to saturate at $\alpha_{\Lambda}\approx-0.25$ in
agreement with [Skokos et al., 2013; Senyange et al., 2018]. This number is
retrieved in this work taking the average value of all the cases in the last
decade of the evolution of $\alpha_{\Lambda}(t)$. These results are not only
specific to the 1D DDNLS model, as they have also been found for the 1D DKG
model of anharmonic oscillators (2.24) in [Skokos et al., 2013; Senyange et
al., 2018].

####  4.3.3 The strong chaos regime

We also examine the strong chaos regime of the 1D DDNLS system (2.17). We use
the following initial localized wave packet excitations:

  * $\bullet$

Case $\text{S}1_{\text{1D}}$. $W=3$, $\beta=1$, $L=37$ and $x=0.3$.

  * $\bullet$

Case $\text{S}2_{\text{1D}}$. $W=3.5$, $\beta=0.62$, $L=21$ and $x=0.62$.

  * $\bullet$

Case $\text{S}3_{\text{1D}}$. $W=3.5$, $\beta=0.72$, $L=21$ and $x=0.72$.

  * $\bullet$

Case $\text{S}4_{\text{1D}}$. $W=3.5$, $\beta=1$, $L=27$ and $x=0.525$.

Again, we justify the clustering of these initial parameter in the strong
chaos spreading regime applying the second condition in Eq. (2.56). It states
that if

$x\frac{L}{\overline{V}}>d,$ (4.7)

the initial wave packet excitation will evolve for a finite amont of time
within the strong chaos region of Fig. 2.6. As we did in the case of weak
chaos, we numerically check the condition (4.7) on the parameter values above.
For the case S$1_{1D}$ the left hand side of inequality (4.7) results in
$0.3\cdot 37/37=0.3$, which is larger than $d=0.2$ of $W=3$. For the cases
S$2_{1D}$, S$3_{1D}$ and S$4_{1D}$ which all use the value $W=3.5$, we
calculate the left hand side of Eq. (4.7) as $0.62\cdot 21/27\approx 0.482$,
$0.72\cdot 21/27=0.56$ and $0.525\cdot 27/27=0.525$ respectively which are all
greater than $d=0.28$ in Table 4.1. Thus all the cases above belong to the
strong chaos regimes according to the theory of Chap. 2. Note that the cases
S$2_{1D}$ and S$3_{1D}$ were used in [Senyange et al., 2018], while we
introduce the cases S$1_{1D}$ and S$4_{1D}$ in order to investigate other
initial conditions than the ones already presented in the literature [Laptyeva
et al., 2010; Bodyfelt et al., 2011a; Senyange et al., 2018]. Like for the
weak chaos regime, we average our observables of the phases space and tangent
dynamics over $100$ disorder realizations in order to ensure the mixing of
several different trajectories with constant $\mathcal{S}_{1D}$ values in the
system’s phase space in order to ensure the robustness of our computations.

Figure 4.8: Similar to Fig. 4.7 but for strong chaos regime. The coloring
corresponds to the cases S$1_{1D}$ [(b)lue], S$2_{1D}$ [(g)reen], S$3_{1D}$
[(r)ed] and S$4_{1D}$ [(c)yan] (see text for details). The straight lines
guide the eye for slopes $\frac{1}{2}$ [panel (a)], $\frac{1}{4}$ [panel (b)]
and $-0.3$ [panel (c)]. The horizontal dashed line in panel (d) indicates the
value $\alpha_{\Lambda}=-0.3$.

The panels of Figs. 4.8(a) and (b) depict the evolution of the wave packet
second moment $m_{2}(t)$ (3.5) and participation number $P(t)$ (3.6)
respectively. Again, a subdiffusive spreading of the wave packet initial
excitation following the power laws $m_{2}(t)\propto t^{\alpha_{m}}$ and
$P(t)\propto t^{\alpha_{P}}$ with $\alpha_{m},\alpha_{P}<1$ is clearly
visible. This time the growth rates are well fitted by
$\alpha_{m}=\frac{1}{2}$ [see black curve in Fig. 4.8(a)] and
$\alpha_{P}=\frac{1}{4}$ [see black curve in Fig. 4.8(b)] in a region of
around the early stages of the time evolution i.e. $\log_{10}t\approx
10^{3}-10^{6}$. This numerically confirms our choice of these cases as strong
chaos ones. We can therefore compute the ftMLEs of all these strong chaos
cases in Fig. 4.8(c). There, a similar tendency $\Lambda(t)\propto
t^{\alpha_{\Lambda}}$ ($\alpha_{\Lambda}<0$) to the one noticed for the weak
chaos cases is seen for all the used cases. Nevertheless, the
$\alpha_{\Lambda}(t)$ values tend to saturate at $\alpha_{\Lambda}\approx-0.3$
as shown in Fig. 4.8(d). Again, this value is obtained averaging the data of
the time evolution of $\alpha_{\Lambda}(t)$ in the last decade of the
evolution where the saturation is observed. It is worth mentioning that once
again the same results have also been found for the 1D DKG model (2.24)
[Senyange et al., 2018] underlining the generality of this power exponent.

###  4.4 The deviation vector distribution

There is a need of a second ingredient, a non-correlated spatiotemporal force
$f(t)$ acting on the chaotic fraction of the system’s NMs evolution as seen in
Eq. (2.105) in order to obtain the chaotic spreading of the wave packet. These
random fluctuations of chaotic seeds are known to support the subdiffusive
spreading of the lattice’s excited part [Flach et al., 2009; Skokos et al.,
2013]. We have been stressing in Sec. 3.4.2 and in the first part of this
chapter that the DV used for the computation of the MLE aligns with the most
chaotic direction intrinsic to a system’s orbit in its phase space. Therefore,
the DVD (3.25)

$\xi_{l}^{D}(t)=\frac{\delta q_{l}^{2}(t)+\delta p_{l}^{2}(t)}{\sum_{l}\delta
q_{l}^{2}(t)+\delta p_{l}^{2}(t)},$ (4.8)

can be used to represent sites which are more sensitive to perturbation at
time $t>0$. Consequently, it can be used to visualize the motion of chaotic
seeds [Oganesyan et al., 2009; Basko, 2011, 2012; Mulansky, 2014; Skokos et
al., 2013; Senyange et al., 2018; Many Manda et al., 2020]. In this section,
we analyse in depth the dynamical behavior of the DVD in the AL, weak chaos
and strong chaos regimes of the initially localized wave packet excitation in
the 1D DDNLS lattice system (2.17).

####  4.4.1 The linear regime

Let us first look more closely to the behavior of the norm density
distribution $\xi_{l}(t)=s_{l}(t)/\sum_{l}s_{l}(t)$ (3.1) with
$s_{l}=(q_{l}^{2}+p_{l}^{2})/2$ and the DVD $\xi_{l}^{D}$ (3.25) generated
during the computation of the ftMLE $\Lambda(t)$ (3.22) when nonlinearity is
absent from the system’s dynamical equations i.e. $\beta=0$ in
$\mathcal{H}_{1D}$ (2.17). We use the case A$1_{1D}$ of Fig. 4.6. For this
case, we compute the spatiotemporal evolution of the norm distribution
$\xi_{l}(t)$ and the DVD $\xi_{l}^{D}(t)$ in Figs. 4.9(a) and (b)
respectively. For both distributions, snapshots at specific times
$\log_{10}t\approx 5.6$ [(b)lue dashed line in Figs. 4.9(a) and (b)],
$\log_{10}t\approx 6.65$ [(g)reen dashed line in Figs. 4.9(a) and (b)] and
$\log_{10}t\approx 7$ [(r)ed dashed line in Figs. 4.9(a) and (b)] are
displayed in Figs. 4.9(c) and (d) respectively for the $\xi_{l}$ and
$\xi_{l}^{D}$.

With $\xi_{l}(t)$ (3.1), a clear spatial saturation of the wave packet extents
appears in Figs. 4.9(a) and (c), in agreement with the saturation of the wave
packet $m_{2}(t)$ and $P(t)$ noticed in Fig. 4.6(a). For $\xi_{l}^{D}$ (3.25)
in Figs. 4.9(b) and (d) a saturation of its spatial evolution is also present.
Further, a striking resemblance between the distributions $\xi_{l}$ and
$\xi_{l}^{D}$ in Fig. 4.9 is visible. This is explained by the fact that if
one computes the variational equations of the 1D DDLS model (2.1), we obtain

$i\dot{\delta\psi}_{l}=\epsilon_{l}\delta\psi_{l}-\left(\delta\psi_{l-1}+\delta\psi_{l+1}\right),$
(4.9)

whose solutions are oscillatory motion of the form $\delta\psi_{l}=\delta
A_{l}\exp{iEt}$ which oscillates in the same manner as the system’s reference
orbit $\psi_{l}=A_{l}\exp iEt$ which resolves the equations of motion of
$\mathcal{H}_{L}$ (2.1) having the same form to (4.9).

Figure 4.9: Spatiotemporal dependence of the norm density distribution
$\xi_{l}$ (3.1) and DVD $\xi_{l}^{D}$ (4.8) for a representative case
A$1_{1D}$ (see text for details) of the linear regime. The color scales at the
top of the figure are used for coloring the lattice sites according to their
$\log_{10}\xi_{l}$ [panel (a)] and $\log_{10}\xi_{l}^{D}$ [panel (b)] values.
In both panels, the white curves represent the motion of the distribution’s
centers. The panels (c) and (d) show snapshots of the spatial dependence of
the distributions $\xi_{l}$ and $\xi_{l}^{D}$ respectively at times
$\log_{10}t\approx 5.6$ [(b)lue curves], $\log_{10}t\approx 6.65$ [(g)reen
curves] and $\log_{10}t\approx 7$ [(r)ed curves]. These times are also
indicated using the same coloring for the dashed horizontal lines in panels
(a) and (b).

Consequently, both the $\xi_{l}(t)$ and $\xi_{l}^{D}(t)$, are expected to
depict qualitatively the same spatiotemporal dynamics. These properties are
expressions of the regularity (non-ergodicity) of the system’s phase space in
the linear regime and tell us that an initial condition taken a bit further
away from the system’s reference orbit will also experience AL as seen in Sec.
2.1.

Figure 4.10: Time dependence of the second moment $m_{2}^{D}(t)$ (3.5) in
panel (a), participation number $P^{D}(t)$ (3.6) in panel (b) of $\xi_{l}^{D}$
(3.25). In panel (c) we plot the length covered by the distribution’s center
$R(t)$ (4.10), while in panel (d) we show the related power law exponent
$\alpha_{R}$ when $R$ is fitted by $R\propto t^{\alpha_{R}}$ at every point of
its evolution. The (b)lue curves correspond to case A$1_{1D}$ and the (g)reen
ones to case A$2_{1D}$ of the linear regime (see text for details). The
horizontal lines in panels (a), (b) and (c) guide the eye for slope $0$.
Similarly, the horizontal dashed line in panel (d) indicates the value
$\alpha_{R}=0$. The results of panels (a), (b) and (c) are averaged over $48$
disorder realizations.

Our observations, on the spatial evolution of the DVD (3.25) are confirmed in
Fig. 4.10(a) and (b) where we plot the time evolution of the DVD’s second
moment $m_{2}^{D}(t)$ [see Fig. 4.10(a)] and participation number $P^{D}(t)$
[see Fig. 4.10(b)] of $\xi_{l}^{D}$ (3.25). Indeed, a clear saturation of the
$m_{2}^{D}$ [panel (a) of Fig. 4.10] and the $P^{D}$ [panel (b) of Fig. 4.10]
is seen, indicative of a halt of the DVD $\xi_{l}^{D}(t)$ (4.8). Note that in
Fig. 4.10 the (b)lue curves corresponds to the case A$1_{1D}$ of Figs. 4.6 and
4.9, while the (g)reen one (Case A$2_{1D}$) is the same condition but this
time obtained with $W=3$ in order to stress the general nature of these
findings. The curves represent averages over $48$ disorder realizations in
order to obtain a robust representation of the behavior of the tangent
dynamics in case of regular motion (we do this because of the stochasticity of
the random disorder parameters).

In [Skokos et al., 2013], it was shown that the DVD oscillates inside the wave
packet extent with increasing amplitudes. In order to capture this dynamical
behavior, in this work we also investigate the spatiotemporal evolution of the
center $\overline{l}^{D}(t)=\sum_{l}l\xi_{l}^{D}(t)$ of the DVD (3.25) by
computing

$R(t)=\max_{[0,t]}\left\\{\overline{l}^{D}(t)\right\\}-\min_{[0,t]}\left\\{\overline{l}^{D}(t)\right\\},$
(4.10)

which measures the distance between the two extreme positions reached by
$\overline{l}^{D}(t)$ at time $t>0$. Thus, $R(t)$ (4.10) measures the total
size of the region visited by the DVD’s center along its evolution. If the
$\overline{l}^{D}$ is fixed in time then
$\max_{\left[0,t\right]}\left\\{\overline{l}^{D}(t)\right\\}=\min_{\left[0,t\right]}\left\\{\overline{l}^{D}(t)\right\\}$,
so that

$R(t)=0.$ (4.11)

This situation is highly improbable as small oscillations of the DVD center
can appear due to numerical fluctuations. In the present case, we say that
whether the $R(t)$ has a tendency to saturate to values smaller or comparable
to the localization volume $\overline{V}$ (2.14) i.e.

$R\lesssim\overline{V},\quad R(t)\sim t^{0},$ (4.12)

we conclude that the DVD center $\overline{l}^{D}(t)$ is practically fixed in
time. Notice that Eq. (4.12) implies that no motion of the DVD center outside
the initially excited region is observed, leading to no thermalization of more
lattice’s sites (or modes). On the other hand, if

$R(t)>\overline{V},$ (4.13)

after a certain time $t$, then we can distinguish between two cases. (i) A
saturation of the values of $R$ take place i.e. $R(t)\sim t^{0}$, or (ii) we
observe a consistent monotonic increasing tendency of $R(t)$ such that after a
certain time $R\gg\overline{V}$. In the case (i) the DVD visited the
neighboring modes of the initially excited one, therefore a wave packet
expansion takes place in one of the regimes mentionned above and then slow
down up to a complete halt after a certain amount of time, restoring the AL
into the system. This situation, corresponds to the predictions of [Johansson
et al., 2010; Aubry, 2011; Mulansky et al., 2011; Pikovsky and Fishman, 2011].
On another hand, the case (ii) means that the DVD’s center is visiting sites
far from the initially excited ones. This is the indication of an infinite
spreading of the wave packet into the system and corresponds to situations
metioned in [Shepelyansky, 1993; Pikovsky and Shepelyansky, 2008; Flach, 2010;
Laptyeva et al., 2010, 2012].

Figure 4.10(c) shows the evolution of the $R(t)$ (4.10) of the center of the
DVD $\overline{l}^{D}$ for the initial conditions A$1_{1D}$ [(b)lue curve] and
A$2_{1D}$ [(g)reen curve] used above. A clear tendency to saturate a value
$R(t)\approx 3$ for the case A$1_{1D}$ and $R(t)\approx 5$ in case of
A$2_{1D}$ is seen. These values are much smaller than the localization volume
$\overline{V}=21$ and $37$ respectively for the values $W=4$ and $W=3$.
Further, we assumed that the behavior of $R(t)$ can be fitted with a
functional form

$R(t)\propto t^{\alpha_{R}},$ (4.14)

and plotted the time evolution of the power law exponent $\alpha_{R}(t)$ in
Fig. 4.10(d), where a clear tendency to converge to values $\alpha_{R}\approx
0$ is seen, supporting the fact that the center of the DVD [see white curve in
Fig. 4.9(b)] remains practically fixed in Figs. 4.9(b) and (d). Therefore, no
modes’ chaotic heating take place as expected in the linear regime of
disordered systems. Note that the value $\alpha_{R}\approx 0$ was calculated
averaging the slopes $\alpha_{R}(t)$ of the cases A$1_{1D}$ and A$2_{1D}$ over
the last decade of evolution.

####  4.4.2 The weak chaos regime

Let us now look at the dynamical behavior of the DVD (3.25) in the case
$\beta>0$ for the 1D DDNLS model with Hamiltonian $\mathcal{H}_{1D}$ (2.17).
In Fig. 4.11(a), we sketch the time evolution of the norm density distribution
$\xi_{l}$ (3.1) along with its related DVD $\xi_{l}^{D}$ (4.8) in Fig. 4.11(b)
for an individual realization of the weak chaos case $\text{W}4_{\text{1D}}$.
Snapshots of these distributions are presented at times $\log_{10}t\approx
4.8$ [(b)lue dashed horizontal line] $\log_{10}t\approx 6.82$ [(g)reen
horizontal dashed line] and $\log_{10}t\approx 7.94$ [(r)ed dashed horizontal
line] in Figs. 4.11(c) and (d).

From Figs. 4.11(a) and (c), we see that the norm density $\xi_{l}$ is
expanding smoothly and in a more or less symmetric way about the distribution
center which remains (almost) fix at the position of the initial excitation
[see white curve in Fig. 2.2(a)]. On the other hand, the DVD $\xi_{l}^{D}$
(4.8) in Figs. 4.11(b) and (d) retains a rather localized character,
especially compared to the norm extent. In addition, its centre
$\overline{l}^{D}$ [see the white curve in Figs. 4.11(b) and (d)] is randomly
wandering inside the lattice excited part as suggested by the fuzziness of the
white curve in Fig. 4.11(b). These random fluctuations increase in amplitude
with time which becomes evident for $\log_{10}t\gtrsim 6$ where large
excursions of the $\xi_{l}^{D}$ centre are clearly visible.

Figure 4.11: Spatiotemporal dependence of the norm density distribution
$\xi_{l}$ (3.1) and DVD $\xi_{l}^{D}$ (4.8) for the representative case
W$4_{1D}$ (see text for details) of the weak chaos regime. The color scales at
the top of the figure are used for coloring the lattice sites according to
their $\log_{10}\xi_{l}$ [panel (a)] and $\log_{10}\xi_{l}^{D}$ [panel (b)]
values. In both panels, the white curves represent the motion of the
distribution’s center. The panels (c) and (d) show snapshots of the spatial
dependence of the distributions $\xi_{l}$ and $\xi_{l}^{D}$ respectively at
times $\log_{10}t\approx 4.8$ [(b)lue curves], $\log_{10}t\approx 6.82$
[(g)reen curves] and $\log_{10}t\approx 7.94$ [(r)ed curves]. These times are
also indicated using similar colors for the dashed horizontal lines in panels
(a) and (b).

In order to validate our observations, we compute the evolution of the DVD’s
(3.25) second moment $m_{2}^{D}(t)$ (3.5) (see Fig. 4.12(a)) and participation
number $P^{D}(t)$ (3.6) (see Fig. 4.12(b)) for all the weak chaos cases of
Sec. 4.3.2. A growth of the $m_{2}^{D}(t)$ well fitted by the power law
$m_{2}^{D}\propto t^{0.14}$ (see black line in Fig. 4.12(a)) is observed
smaller than the $m_{2}(t)\propto t^{1/3}$ in the case of the $\xi_{l}$ (3.1).
In addition, the $P^{D}(t)$ slowly increases in all the studied cases
following a timid power law $P^{D}\propto t^{0.04}$ in Fig. 4.12(b) and
reaching values similar to the mode’s localization volume $\overline{V}$
(2.14) (in the cases W$1_{1D}$ [(b)lue] and W$2_{1D}$ [(g)reen] $P^{D}\approx
10^{1.3}\approx 20$) before showing a tendency to saturate at the last stages
of the evolution in accordance to the localized character of the DVD observed
in Figs. 4.11(b) and (d).

Figure 4.12: Results of the temporal dependence of the DVD’s $\xi_{l}^{D}$
(4.8) second moment $m_{2}^{D}$ (3.5) in panel (a) and participation number
$P^{D}$ (3.6) in panel (b). Panel (c) shows the region $R$ (4.10) visited by
the distribution’s center and in panel (d) we depict the power exponent
$\alpha_{R}$ when the values of $R(t)$ are fitted with $t^{\alpha_{R}}$. All
results corresponds to the weak chaos regime of the 1D DDNLS system (2.17).
The used initial wave packet parameters are for cases W$1_{1D}$ [(b)lue],
W$2_{1D}$ [(g)reen], W$3_{1D}$ [(r)ed] and W$4_{1D}$ [(c)yan] (see text for
details). The black lines in the panels (a), (b) and (c) guide the eye along
the slopes $0.14$, $0.04$ and $0.24$. The horizontal line in panel (d)
indicates the value $\alpha_{R}=0.24$ obtained averaging the two last decades
of the evolution of $\alpha_{R}(t)$ for all the cases. The curves presented in
panels (a), (b) and (c) show average results over $100$ disorder realizations.

We also compute the size $R(t)$ (4.10) of the region visited by the centre of
the DVD $\overline{l}^{D}(t)=\sum_{l}l\xi_{l}^{D}(t)$. The result is depicted
in Fig. 4.12(c), where a monotonic increase of $R$ with time is seen which
attains values $R(t)\approx 10^{2}$, much larger than the NM localization
volume. This increase is well fitted by an asymptotic power law $R\propto
t^{\alpha_{R}}$ ($\alpha_{R}>0$) where the evolution of $\alpha_{R}(t)$ is
shown in Fig. 4.12(d). There, all the curves of the power law exponent
$\alpha_{R}(t)$ exhibit a tendency to saturate to an estimated exponent
$\alpha_{R}\approx 0.24$ (we averaged the $\alpha_{R}(t)$ of all the studied
cases in the two last decades of their evolutions to obtain this value). This
asymptotic value of the variations is reached at an early stage
($\log_{10}t\approx 2$) of the evolution and persists for the remaining $6$
decades of our numerical integration of the phase space and tangent space
trajectories. Consequently, the fluctuations of the chaotic seeds exist from
near the beginning of the evolution and sustain their growth along with the
wave packet expansion. It is worth mentioning that the same results were found
also for the 1D DKG model of anharmonic oscillators (2.24) in [Senyange et
al., 2018].

####  4.4.3 The strong chaos regime

We also investigate the dynamical behavior of the DVD $\xi_{l}^{D}$ (4.8) in
the strong chaos region of Fig. 2.4 and 2.6. Here, we use a realization of the
$\text{S}3_{\text{1D}}$ case as a representative of the strong chaos region
dynamical behavior in Fig. 4.13. We plot the evolution of the norm density
$\xi_{l}$ (3.1) (see Fig. 4.13(a)), and its related DVD $\xi_{l}^{D}$ (see
Fig. 4.13(b)) along with snapshots of these distributions at specific times
(see Figs. 4.13(c) and (d)) namely $\log_{10}t\approx 6$ [(b)lue horizontal
lines], $\log_{10}t\approx 6.54$ [(g)reen horizontal lines] and
$\log_{10}t\approx 7.24$ [(r)ed horizontal lines] in Figs. 4.13(a) and (b).
Similar to what was observed for the weak chaos case in Fig. 4.11, the norm
distribution expands rather smoothly and symmetrically around its center [the
white curve in Fig. 4.13(a)]. The $\xi_{l}$ spreads more widely compared to
the weak chaos case in Fig. 4.11(a) as in this regime the wave packet spreads
faster. On the other hand, the $\xi_{l}^{D}$ exhibits a localized character
[see Fig. 4.13(b)], with the distribution profile retaining a confined
‘pointy’ shape [see Fig. 4.13(d)]. In addition, we see that the centre of
$\xi_{l}^{D}$ randomly oscillates inside the wave packet width [see the
fuzziness of the white curve of Fig. 4.13(b) and the profile positions in (d)]
with increasing amplitudes, moving faster and wider than in the weak chaos
regime [see Figs. 4.11(b) and (d)].

Figure 4.13: Spatiotemporal dependence of the norm density distribution
$\xi_{l}$ (3.1) and DVD $\xi_{l}^{D}$ (4.8) for the representative case
S$3_{1D}$ (see text for details) of the strong chaos regime. The color scales
at the top of the figure are used for coloring the lattice sites according to
their $\log_{10}\xi_{l}$ [panel (a)] and $\log_{10}\xi_{l}^{D}$ [panel (b)]
values. In both panels, the white curves represent the motion of the center’s
distributions. The panels (c) and (d) show snapshots of the spatial dependence
of the distributions $\xi_{l}$ and $\xi_{l}^{D}$ respectively at times
$\log_{10}t\approx 6.0$ [(b)lue curves], $\log_{10}t\approx 6.53$ [(g)reen
curves] and $\log_{10}t\approx 7.23$ [(r)ed curves]. These times are also
indicated using similar colors for the dashed horizontal lines in panels (a)
and (b).

We quantify these observations by computing the DVD’s $\xi_{l}^{D}$ (4.8)
second moment $m_{2}^{D}(t)$ (3.5) [see panel (a) of Fig. 4.14] and
participation number $P^{D}(t)$ [see panel (b) of Fig. 4.14] for all the
strong chaos cases of Sec. 4.3.3. After a fast growth of the $m_{2}^{D}$ and
$P^{D}$ values at the early stage of the evolution ($\log_{10}t\lesssim 5$), a
tendency to saturate to a constant value is clearly visible for both
quantities. For instance, we have fitted the last decades of evolution of the
$P^{D}(t)$ with a power law $P^{D}(t)\propto t^{a}$ and obtained $a\approx 0$.
It is worth mentioning that the number of highly excited sites reached by the
$\xi_{l}^{D}$ is $P^{D}\approx\overline{V}$ (for example W$1_{1D}$ gives
$P^{D}\approx 10^{1.4}\approx 25$) which suggests that the DVD always remains
confined in a single mode’s localization volume at a time while moving around
from site to site. These results support the pointy and concentrated profile
seen in Figs.4.13(b) and (d).

The evolution of the size $R(t)$ (4.10) of the region visited by the DVD
center $\overline{l}^{D}$ is depicted in Fig. 4.14(c). A clear monotonic
increase is visible from the early stage of the evolution with $R(t)$ values
reaching up to $R\approx 10^{2.8}\approx 700$ sites at the final time of the
integration, several times larger than the average localization volume
$\overline{V}$ of a NM. The dependence of $R(t)$ on time can be fitted with a
power law $R\propto t^{\alpha_{R}}$ with $\alpha_{R}>0$. In Fig. 4.14(d), the
evolution of the power law exponent $\alpha_{R}(t)$ is shown. We see that, in
a large part of the dynamical evolution, the value of $\alpha_{R}(t)$ is
always well above the ones observed for the weak chaos regime [black