horizontal dashed line in Fig. 4.14 (d)]. Then, at later stages of the
evolution ($\log_{10}t>7$), the $\alpha_{R}(t)$ asymptotically converges to
values $\alpha_{R}\approx 0.24$ seen in the weak chaos regime.

This transitional behavior is well in agreement with previous works [see e.g.
Flach, 2010; Laptyeva et al., 2010; Skokos and Flach, 2010; Bodyfelt et al.,
2011a]. Indeed, the fact that the $R$ values reached in the weak chaos regime
($R\approx 10^{2}$) are smaller than the ones observed in the strong chaos
case ($R\approx 10^{3}$) at the final time of our numerical integration is
likely explained by the fact that the faster spreading in the strong chaos
regime implies faster and wider wandering of chaotic seeds as a more rapid
randomization and thermalization are required. In addition, because the strong
chaos is a transient regime, an asymptotic crossover of the wave packet
dynamics to the weak chaos regime should take place. This is what we see in
Fig. 4.14(d) for example since the rate of growth of $R(t)\propto
t^{\alpha_{R}}$ i.e. $\alpha_{R}$ is always larger for the strong chaos regime
and converge to the $\alpha_{R}$ values seen in the weak chaos dynamical
behavior at later stages. This numerically demonstrates that the nature of the
chaotic dynamics cross over well before the wave packet dynamics itself for
which up to now a complete crossover of a strong chaos initial condition has
never been reported in the literature [see e.g. Laptyeva et al., 2010]. It is
worth noticing that the same behavior of the DVD in the strong chaos regime
was retrieved also for the 1D DKG model in [Senyange et al., 2018].

Figure 4.14:  Results of the time dependence of the second moment $m_{2}^{D}$
(3.5) in panel (a) and the participation number $P^{D}$ (3.6) in panel (b) of
the DVD $\xi_{l}^{D}$ (4.8). The region visited by the distribution’s center
$R$ (4.10) is in panel (c) and the power exponent $\alpha_{R}$ in panel (d)
when the values of $R(t)$ are fitted with the relation $t^{\alpha_{R}}$ for
the strong chaos regime. The used initial wave packet parameters are for 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 black lines in the panels (a), (b) guide
the eye for slope $0.0$. The horizontal dashed line in panel (d) indicates the
value $\alpha_{R}=0.24$. The curves presented in panels (a), (b) and (c) show
average results over $100$ disorder realizations.

###  4.5  Discussion and summary

Let us now discuss the main implications of our results in this chapter. The
outputs of the computation of the ftMLE $\Lambda\propto t^{-0.25}$ for the
weak and $\Lambda\propto t^{-0.3}$ for the strong chaos spreading regimes are
quite interesting. In studies of finite size multidimensional dynamical
systems the chaotic dynamics show a plateau of the $\Lambda(t)$ values for
chaotic trajectories (see e.g. [Benettin et al., 1980; Mithun et al., 2018;
Benettin et al., 2018; Hillebrand et al., 2019; Mithun et al., 2018;
Hillebrand et al., 2020]), while quasi-periodic trajectories reveal a
$\Lambda(t)\propto t^{-1}$ tendency [Benettin et al., 1980; Skokos, 2010]. The
results we are obtaining here can be interpreted in the following way. As the
constant total norm $\mathcal{S}_{1D}$ (3.1) (energy $\mathcal{H}_{1K}$ (2.24)
for the 1D DKG model) is shared among more and more thermalized sites as
additional degrees of freedom get activated, the average norm per site $s$
(2.68) of the excited part decreases in time. Therefore both the single
oscillator nonlinear frequency shift $\delta\propto\beta s$ [used for example
in the the Chirikov criterion 2.34] and the effective nonlinear strength
$\propto\beta s^{2}$ [present for example in the equations of motion (2.18)
and (2.105)] decay in time, taking the system closer to the integrable limit.
While one may expect a cross over to quasi-periodic (regular) motion as
speculated in [Johansson et al., 2010; Aubry, 2011], this is never the case,
as the exponents of the power law $\Lambda\propto t^{\alpha_{\Lambda}}$ with
$\alpha_{\Lambda}<0$ are $\alpha_{\Lambda}=-0.25$ (weak chaos) and
$\alpha_{\Lambda}=-0.3$ (strong chaos) which are different from the value
$\alpha_{\Lambda}=-1$ seen in the case of regular motion [see for example Fig.
4.6(b)]. Thus, our results suggest that the system always remains chaotic
although its chaoticity decreases in time.

Further, the fact that the exponent of the strong chaos regime
$\alpha_{\Lambda}=-0.3$ is smaller than the one in the weak chaos case
$\alpha_{\Lambda}=-0.25$ is explained by the fact that the wave packet spreads
faster in the strong chaos regime where $m_{2}(t)\propto t^{1/2}$ is observed
than in the weak chaos region where $m_{2}(t)\propto t^{1/3}$, so that the
average norm per site $s$ within the wave packet excited part and the related
effective nonlinearity $\beta s^{2}$ of the system decrease faster in the
strong chaos case. Thus orbits in the parameter spaces of the 1D DDNLS model
(2.17) in the strong chaos region move faster toward the integrable limit than
in the weak chaos region. This means that the value of the exponent in the
power law $\Lambda\propto t^{\alpha_{\Lambda}}$ in the strong chaos case is
closer to the $-1$ value observed in the integrable limit compared to the weak
chaos case (indeed we see that $-1<-0.3<-0.25$).

Since the wave packet chaoticity is decreasing, an important aspect is to know
whether the wave packet chaotic dynamics will always be able to support the
subdiffusive spreading. In order to analyze that point, we compare the
system’s timescales. This was first done by [Skokos et al., 2013] for the weak
chaos regime of the 1D DKG model (2.24). Here, we extend this analysis to the
weak and strong chaos regimes of the 1D DDNLS system (2.17).

A characteristic timescale for the system’s chaotic dynamics is the Lyapunov
time (3.24)

$T_{\Lambda}=\frac{1}{\Lambda},$ (4.15)

which gives an estimation of the time it takes for a dynamical system to
become chaotic. On the other hand, if we assume that the wave packet dynamics
is mainly characterized by its second moment $m_{2}$ (3.5), then

$T_{M}\sim\frac{1}{D},$ (4.16)

describes a timescale above which the subdiffusive spreading effectively
starts. Note that

$D\sim\frac{dm_{2}}{dt},$ (4.17)

as the behavior $m_{2}\sim Dt$ is assumed by definition of the diffusion
coefficient $D$ (see also Sec. 2.5.1). Alternatively, if we assume that it is
the participation number which mainly measures the wave packet subdiffusive
spreading, we could define in a similar way the spreading timescale

$T_{P}\sim\frac{1}{\dot{P}},$ (4.18)

where $\dot{P}=dP/dt$. In the case of the spreading regime mentioned here, we
have $m_{2}\sim t^{\alpha_{m}}$ and $P\sim t^{\alpha_{m}/2}$. We also saw that
$\Lambda\sim t^{\alpha_{\Lambda}}$. Consequently,

$\frac{T_{M}}{T_{L}}\sim
t^{1+\alpha_{\Lambda}-\alpha_{m}},\qquad\frac{T_{P}}{T_{L}}\sim
t^{1+\alpha_{\Lambda}-\alpha_{m}/2},$ (4.19)

which reduce to

$\frac{T_{M}}{T_{L}}\sim t^{5/12},\qquad\frac{T_{P}}{T_{L}}\sim t^{7/12},$
(4.20)

for the weak chaos spreading regime in which $\alpha_{m}=1/3$,
$\alpha_{\Lambda}=-0.25$ was found, and to

$\frac{T_{M}}{T_{L}}\sim t^{1/5},\qquad\frac{T_{P}}{T_{L}}\sim t^{9/20}$
(4.21)

in the strong chaos case where we obtained $\alpha_{m}=1/2$ and
$\alpha_{\Lambda}=-0.3$. Consequently the Lyapunov time $T_{L}$ (4.15) always
remains smaller than the spreading timescales defined by $T_{M}$ (4.16) and
$T_{P}$ (4.18), so that the wave packet chaotization is always faster than its
spreading to new lattice’s sites.

The fact that the DV needed for the computation of the MLE is aligned to the
most chaotic direction is helpful to visualize the motion of chaotic seeds. We
identified chaotic seeds as regions in the lattice which are highly sensitive
to initial perturbations. Because the region where the NMs overlap and
exchange the energy in a deterministic chaotic manner is finite, we expected
to see chaotic spots (seeds) which are localized in time. In addition, these
chaotic seeds randomly meander within the wave packet interior due to the
presence of a stochastic force $f(t)$ in Eq. (2.108) acting on the modes when
the system is within a thermalized state. This process is often related to the
homogenization of chaos and promote the wave packet thermalization of more
lattice sites.

In order to test this hypothesis, we computed the DVD’s $\xi_{l}^{D}$ (3.25)
second moment $m_{2}^{D}$ (3.5), and participation number $P^{D}$ (3.6) for
both the weak and strong chaos regimes. Our numerical results suggest that the
DVD remains localized within the wave packet with size of about a mode
localization volume. This saturation is reached faster in the strong chaos
case compared to the weak chaos where increase of the DVD extent, following
$m_{2}^{D}\propto t^{0.14}$ and $P^{D}\propto t^{0.04}$ was still seen up to
the largest simulation time. In addition, the center of the DVD exhibits
random fluctuations with increasing amplitudes inside the wave packet width.

We computed the size of the region visited by the DVD center and saw that its
extent monotonically increases with time, reaching values of $R\approx 10^{3}$
for the strong chaos case at the final time of integration ($\log_{10}t\approx
10^{8}$). This is a clear indication that the DVD is visiting modes (and
sites) far outside from the initially excited one. A power law $R\propto
t^{0.24}$, was noticed for the weak chaos, while for the strong chaos a faster
growth $R\propto t^{\alpha_{R}}$ is reported with $\alpha_{R}>0.24$ for
several decades, which asymptotically tends to the same $\alpha_{R}\approx
0.24$ values observed in the weak chaos regime. This is related to the fact
that strong chaos is a transient regime which asymptotically evolves to the
weak chaos one.

All these results suggest that a non-zero chaotic strength coupled to random
fluctuations of the chaotic seeds are the two ingredients for the subdiffusive
spreading of the wave packet in the 1D DDNLS model. Thus, as long as they are
present in the system no crossover of the extended wave packet to regular
dynamics is expected to take place. This fact is general as it has been also
observed for the 1D DKG model Senyange et al. [2018].

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

###  5.1 Introduction

In this chapter, we will investigate the characteristics of wave packet
spreading in the weak and strong chaos regimes conjectured in Secs. 2.3 and
2.5.2 for the 2D DDNLS model (2.20). The main obstacle facing this study is
related to the long-time numerical integration of the 2D DDNLS model (2.20)
which is a computationally hard task. This is seen for instance in the
literature where up to date only a few works on the subject are available.
More specifically, [Garcia-Mata and Shepelyansky, 2009] studied the wave
packet spreading in the 2D DDNLS model (2.20) only for a single weak chaos
case which requires smaller lattice sizes up to a final time of $t_{f}=10^{6}$
time units, averaging over only $10$ disorder realizations. Further [Sales et
al., 2018] studied the diffusion of localized excitations in a similar model
of the 2D DDNLS system, with non-diagonal nonlinear terms, integrating the
phase space trajectories up to $t_{f}=10^{5}$, and averaging over few disorder
realizations for each case. On the other hand, [Laptyeva et al., 2012]
performed an analysis of over $400$ disorder realizations of the spreading in
the weak and strong chaos regimes for the generalized 2D DKG system
$\mathcal{H}_{2K}$ (2.27) whose equations of motion are much easier to
simulate numerically. Thus the numerical investigation of wave packet
spreading in the 2D DDNLS lattice (2.20) we consider in this thesis remains
largely unexplored. The novelty of this work is two fold as we not only
compute the characteristics of the wave packet spreading, but also we
investigate the system’s chaotic dynamics in detail, to the best of our
knowledge for the first time.

The chapter is structured as follows: In Sec. 5.2 we present additional
numerical considerations for the integration of the variational equations of
the 2D lattice. Section 5.3 is devoted to presenting the results of our
computations of the MLE for several cases of localized excitations in the weak
and strong chaos regimes. In Sec. 5.4 we investigate the local nature of
chaos, analyzing the dynamical behavior of the DVD in every spreading regime.
The discussion of our findings and the summary of our work are presented in
Sec. 5.5. The outputs of this chapter are heavily based on the research paper
[Many Manda et al., 2020].

###  5.2 Some aspects related to our numerical simulations

####  5.2.1 Choice of the initial deviation vector

We start our investigation by first analyzing the dynamical behavior of the
tangent dynamics of the 2D DDNLS model (2.20) and more specifically its
dependence on the type of the used initial DV when a central localized
excitation is set as initial state of the lattice. To that end, we consider a
typical initial condition consisting of exciting only the central site of a
$201\times 201$ lattice with a total norm (2.23) $\mathcal{S}_{2D}=0.375$. The
total energy (2.20) is $\mathcal{H}_{2D}\approx 0.7$ for the particular
realization of disorder we choose, with the other parameters set to $W=10$ and
$\beta=10$. We integrate the equations of motion and the variational equations
of the system up to a final time $t_{f}=10^{6}$, using the $s11\mathcal{ABC}6$
SI (3.84). The integration time step $\tau=0.24$ used, leads to relative
energy (3.92) $E_{r}$ and norm (3.93) $S_{r}$ errors always remaining below
$10^{-3}$. For the initial condition of the variational equations, we
implement the following choices of DVs:

  * Case DV1. We give random values to the $\delta q_{l,m}$ and $\delta p_{l,m}$ coordinates of a square with size equal to about the average localization volume of a NM (see Table 5.1 below) i.e. $L\times L=5\times 5$ at the center of the lattice. The related spatial distribution of the DVD $\xi_{l,m}^{D}$ (5.2) is shown in Fig. 5.1(b).

  * Case DV2. Similar to the case of DV1 but the size of the initially excited square is approximately $100$ times larger than the average localization volume of a NM i.e. $L\times L=50\times 50$, leading to the spatial arrangement of the DVD $\xi_{l,m}^{D}$ (5.2) seen in Fig. 5.2(b).

  * Case DV3. Similar to case DV2 but the size of the initially excited square is $4$ times larger, corresponding to a spatial distribution of the DVD $\xi_{l,m}^{D}$ (5.2) in Fig. 5.3(b).

  * Case DV4. We randomly excite the $\delta q_{l,m}$ and $\delta p_{l,m}$ coordinates of only the edge layer of the 2D lattice, leading to the spatial arrangement of the DVD $\xi_{l,m}^{D}$ (5.2) in Fig. 5.4(b).

The results of our computations for each considered DV are shown as follows:
for DV1 in Fig. 5.1, DV2 in Fig. 5.2, DV3 in Fig. 5.3 and DV4 in Fig. 5.4. In
each of these figures, panels (a), (c) and (e) present snapshots of the
temporal evolution of the spatial extent of the norm density distribution
(3.1)

$\xi_{l,m}(t)=\frac{s_{l,m}(t)}{\mathcal{S}_{2D}},\quad
s_{l,m}(t)=\frac{1}{2}\left(q_{l,m}^{2}(t)+p_{l,m}^{2}(t)\right),$ (5.1)

while panels (b), (d) and (f) show the related DVD (3.25)

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

at times $\log_{10}t=0$ [panels (a) and (b)], $\log_{10}t\approx 2.45$ [panels
(c) and (d)] and $\log_{10}t\approx 5.95$ [panels (e) and (f)]. It is worth
noticing that the $s_{l,m}$ in Eq. (5.1) is the local norm at site $(l,m)$ of
the 2D DDNLS lattice (as in Chap. 4, we will often use $s$ instead of
$s_{l,m}$ when referring to its average value over the entire excited part of
the lattice). Evidently, the time evolution of the $\xi_{l,m}$ (5.1) is the
same in all simulations [panels (a), (c) and (e)], indicating a chaotic wave
packet whose extent increases in time, as this state belongs to the thermal
region of the energy-norm parameter space of Fig. 2.6. On the other hand, as
the initial DVs are different in each computation, we expect the tangent
dynamics to depict qualitatively different outputs at least in the initial
parts of the evolution of the variational equations.

Figure 5.1: 3D spatial distributions and the associated 2D color maps (upper
sides) snapshots of the norm density $\xi_{l,m}$ (5.1) [panels (a), (c) and
(e)] and the DVD $\xi_{l,m}^{D}$ (5.2) [panels (b), (d) and (f)] profiles for
a representative case of localized initial excitation in the 2D DDNLS lattice
(2.20). We excite the central site of a $201\times 201$ lattice with a norm
(2.23) $\mathcal{S}_{2D}=0.375$ at $W=10$ and $\beta=10$. The system’s total
energy is $\mathcal{H}_{2D}\approx 0.7$ (2.20). The initial condition for the
variational equations is DV1 (see text for details). Panels (a) and (b)
correspond to the configurations of the $\xi_{l,m}$ and $\xi_{l,m}^{D}$ at
time $\log_{10}t=0$, while panels (c) and (d) at $\log_{10}t\approx 2.45$ and
panels (e) and (f) at $\log_{10}t\approx 5.95$. The red curves on the sides
are the distribution’s projections along the $l$ and $m$ directions. The
coloring of the distributions is done with respect to the $\log_{10}\xi_{l,m}$
and $\log_{10}\xi_{l,m}^{D}$ values according to the horizontal color scales
presented at the top of the figure.  Figure 5.2: Similar to Fig. 5.1, but for
the initial condition DV2 (see text for details) for the variational equations
of the Hamiltonian function of the 2D DDNLS model $\mathcal{H}_{2D}$ (2.20).
Figure 5.3: Similar to Fig. 5.1, but for the initial condition DV3 (see text
for details) for the variational equations of the Hamiltonian function of the
2D DDNLS model $\mathcal{H}_{2D}$ (2.20). Figure 5.4: Similar to Fig. 5.1, but
for the initial condition DV4 (see text for details) for the variational
equations of the Hamiltonian function of the 2D DDNLS model $\mathcal{H}_{2D}$
(2.20). Note that as only the last sites of the lattice’s edges initially
receive the $\xi_{l,m}^{D}$ (5.2) values, it is visible only on the 3D heatmap
of panel (b). The size of the DVD’s excited region is too small compared to
the lattice’s dimension to be visible in the upper color map of panel (b).
Figure 5.5: Time evolution of the ftMLE $\Lambda(t)$ (3.21) for the
representative initial condition of the 2D DDNLS lattice (2.20) used in Figs.
5.1, 5.2, 5.3 and 5.4. The curves correspond to results obtained for the DV1
[(b)lue], DV2 [(g)reen], DV3 [(r)ed] and DV4 [(c)yan] (see text for details).

The results of the computation of the DVD $\xi_{l,m}^{D}$ (5.2) up to times
around $\log_{10}t\approx 2.45$ in panels (d) of Figs. 5.1, 5.2, 5.3 and 5.4,
respectively for the cases DV1, DV2, DV3 and DV4 confirm our predictions, as
we see a completely different behavior of the DVs in each case. For the DV1, a
clear expansion of the DVD takes place with the distribution remaining well
inside the wave packet extent [see panels (c) and (d) of Fig. 5.1], where the
deterministic chaotic dynamics is concentrated in the 2D lattice. For the DV2,
the DVD $\xi_{l,m}^{D}$ spreads over the whole lattice with the unexcited
background [from the point of view of $\xi_{l,m}$ in Fig. 5.2(c)] having
smaller values compared to the chaotic part of the lattice [panel (d) of Fig.
5.2]. In the case of DV3, a similar behavior to the one of DV2 is observed.
However, a much stronger contribution of the unexcited background is seen, as
the values of the $\xi_{l,m}^{D}$ there are much larger when compared to the
DV2 case. The DV4 choice depicts a completely different behavior to the three
previous cases as only the unexcited part of the lattice possesses
$\xi_{l,m}^{D}$ values during the transient phase of the dynamics. As we did
for the 1D DDNLS model (2.17) in Sec. 4.2.1, we correlate the dynamics of the
DVD to the behavior of the evolution of the ftMLEs (3.21) plotted in Fig. 5.5.
Indeed, at short times, the DV1 leads to the largest ftMLE [see (b)lue curve
in Fig. 5.5], as its evolution is contained within the chaotic domain of the
lattice from the earlier stages of the evolution. On the other hand, the
contribution of the wide unexcited backgrounds whose dynamics are regular in
the cases DV2 and DV3, lowers the value of their ftMLEs compared to the DV1
case, with the DV2 displaying a larger ftMLE than DV3 [respectively the
(g)reen and (r)ed curves in Fig. 5.5]. With the DV4 evolving solely on the
initially unexcited part of the lattice, we expect its ftMLE values to be very
small, at least compared to the other ftMLEs at the first stages of the
dynamics, something which is clearly seen on the evolution of the (c)yan curve
in Fig. 5.5.

During the last two decades of the evolution of the ftMLEs in Fig. 5.5, all
the DV cases present similar values for their ftMLEs, as a result of the fact
that all DVs align to the same most chaotic direction related to the system’s
reference orbit, whose exponential rate of changes measures its MLE.
Consequently, all the spatial distributions of the DVDs tend to be practically
the same as shown in the panels (f) of Figs. 5.1, 5.2, 5.3 and 5.4 whose
snapshots were taken at time $\log_{10}t\approx 5.95$. Note that the most
chaotic direction possesses coordinates $\delta q_{l,m}$ and $\delta p_{l,m}$
only at the lattice’s excited sites where the chaotic dynamics of the system
takes place.

Based on these numerical experiments, we conclude that initial DVs with random
values of the $\delta q_{l,m}$ and $\delta p_{l,m}$ coordinates on a square
with area smaller or of the order of the NM average localization volume in the
region of the lattice’s initially excited sites are good choices. Otherwise
the DV dynamics has a tendency to initially spread over the unexcited part of
the lattice and consequently converging much more slowly to the most chaotic
direction of the system’s computed trajectory.

####  5.2.2 The dynamic lattice expansion

Let us now discuss how we implement the dynamic lattice technique of Sec.
4.2.2 in the 2D system, as a way of saving CPU time. This is in contrast to
using a fixed large $N\times M$ lattice from the beginning of the dynamical
evolution of orbits, so that the lattice edges are practically not reached by
the wave packet at the end of the integration. In the dynamic lattice approach
we start with a much smaller lattice around the initially excited sites and
expand it whenever it is required by $\mathcal{N}$ sites at the edges of all
lattice’s sides. As for the 1D DDNLS model (2.17) this method requires a
criterion which checks if the wave packet’s front has hit the boundaries. For
the 2D DDNLS system (2.20), we use the criterion

$\max_{l\in\mathcal{D}}s_{l,m}>\delta_{\mathcal{D}},$ (5.3)

where $\mathcal{D}$ represents a ribbon at the outermost layer of the lattice
and $\delta_{\mathcal{D}}$ is a threshold number. Notice that this condition
is different to the one used in the 1D DDNLS system (4.1)

$\frac{1}{\mathcal{D}}\sum_{l,m}s_{l,m}>\delta_{\mathcal{D}}.$ (5.4)

Nevertheless, we numerically checked that results obtained from the two
conditions (5.4) and (5.3) are similar. In our study we prefer using the
condition (5.3) as it can be easily applied in the case of 2D lattice.

Figure 5.6: Spatiotemporal dependence of the projection of the norm density
distribution $\xi_{l,m}$ (5.1) and its related DVD $\xi_{l,m}^{D}$ (5.2) for
$7\times 7$ central square excitation in a $401\times 401$ 2D DDNLS lattice
(2.20). Each excited site is given a norm $s=0.72$ such that the total norm
$\mathcal{S}_{2D}=35.28$ (2.23) with $W=10$ and $\beta=1$. For the particular
disorder realization, the total energy (2.20) is
$\mathcal{H}_{2D}\approx-4.714$. The color scales at the top of the figure are
used for coloring the lattice sites according to their $\log_{10}\xi_{l,m}$
values. In both panels, the white curves represent the motion of the
distribution projections’ centers. Panels (c) and (d) contain snapshots of the
spatial dependence of the distribution’s projections respectively at times
$\log_{10}t\approx 1.8$ [(b)lue curves], $\log_{10}t\approx 3.6$ [(g)reen
curves] and $\log_{10}t\approx 4.8$ [(r)ed curves]. These times are also
indicated using the same colors for the dashed horizontal lines in panels (a)
and (b).

The results of a numerical comparison between the fixed and dynamic lattice
simulations are shown in Fig. 5.6. The initial excitation considered there
corresponds to all sites in a square of size $L\times L=7\times 7$ at the
center of a $401\times 401$ lattice (the center of the square coincides with
the lattice center) having the same norm per site $s_{l,m}=0.72$. This means
that the total norm of the system (2.23) is $\mathcal{S}_{2D}=35.28$, while
the value of the Hamiltonian (2.20) is $\mathcal{H}_{2D}\approx-4.174$ for the
disorder realization we choose, corresponding to a thermal state in the
energy-norm parameter space in Fig. 2.6 of the system. In addition, we set
$W=10$ and $\beta=1$. We integrate the system using the $s11\mathcal{ABC}6$ SI
(3.84) of order $6$ with an integration time step $\tau=0.135$ which kept the
relative energy $E_{r}$ (3.92) and norm $S_{r}$ (3.93) errors below $10^{-4}$
up to the final time of integration $t_{f}\approx 10^{5}$.

Figure 5.7: Time evolution of the second moment $m_{2}(t)$ (5.5) [panel (a)],
the ftMLE $\Lambda(t)$ (3.21) [inset of panel (a)] and the CPU time $T_{C}$
[panel (b)] for a particular setup of the initial excitation of the 2D DDNLS
system (2.20) when the dynamic lattice approach is used in comparison to the
fixed lattice computation. The curve’s coloring corresponds to
$\delta_{\mathcal{D}}=10^{-12}$ [(b)lue curves],
$\delta_{\mathcal{D}}=10^{-10}$ [(g)reen curves],
$\delta_{\mathcal{D}}=10^{-8}$ [(r)ed curves] and the fixed lattice [(c)yan
curves] computations.

The initial size of the active dynamic lattice was set to $101\times 101$
which is much larger than the average localization volume of the NMs
($\overline{V}\approx 35$ for $W=10$ as seen in Table 5.1 below). In addition,
we use $\delta_{\mathcal{D}}=10^{-12}$ and $\mathcal{N}=10$.

The spatiotemporal evolution of the projection along the $m$ axis of the
system of the norm density distribution $\xi_{l,m}$ (5.1) in the case of the
fixed lattice computation is presented in Fig. 5.6(a) along with snapshots in
Fig. 5.6(c), at times indicated by the dashed lines in Fig. 5.6(a). A residual
norm $s\lesssim 10^{-12}$ spreads over the entire lattice due to numerical
fluctuations [see lighter areas in Fig. 5.6(a)]. The same results obtained
using the dynamic lattice simulation after $11$ expansions of the initial
lattice is shown in Fig. 5.6(b) with its related snapshots taken at the same
times [in Fig. 5.6(d)] as for the fixed lattice run. Note that the white areas
in Fig. 5.6(b) correspond to unexcited regions i.e. these sites possess zero
norm ($s_{l,m}=0$). Clearly, the two simulations give practically the same
wave packet evolution.

We repeat the same simulations, using different thresholds
$\delta_{\mathcal{D}}$ [see Eq. (5.3)] for the dynamic lattice simulations. In
Fig. 5.7 we show the result of such numerical computations where the values
$\delta_{\mathcal{D}}=10^{-12}$ [(b)lue curves],
$\delta_{\mathcal{D}}=10^{-10}$ [(g)reen curves],
$\delta_{\mathcal{D}}=10^{-8}$ [(r)ed curves] and the fixed lattice [(c)yan
curves] computations are used. We plot the evolution of the second moment of
the norm density distribution $\xi_{l,m}$ (5.1)

$m_{2}=\sum_{l,m}\left[(l-\overline{l})^{2}+(m-\overline{m})^{2}\right]\xi_{l,m},\quad\overline{l}=\sum_{l,m}l\xi_{l,m},\quad\overline{m}=\sum_{l,m}m\xi_{l,m},$
(5.5)

in Fig. 5.7(a), and the evolution of the associated ftMLEs $\Lambda$ (3.21)
for all cases in the inset of Fig. 5.7(a). A good correspondence is visible
between all the cases. Similarly to what we did for the 1D DDNLS model (2.17)
in Sec. 4.2.1, our choice of the $\delta_{\mathcal{D}}$ depends on the final
values of the norm per site attained in our widest wave packet. In the 2D
DDNLS model (2.20), we are reaching wave packet width $m_{2,f}\propto 10^{3}$
sites. Since the constant total norm (2.23) $\mathcal{S}_{2D}\approx
10^{0}-10^{2}$ is shared among all the excited sites, the final average norm
per site

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

such that using

$\delta_{\mathcal{D}}=10^{-8},$ (5.7)

we capture the wave packet spreading with a rather good precision as shown in
Fig. 5.7(a) up to the final time of integration $t_{f}\approx 2\times 10^{6}$.

In Fig. 5.7(b) we plot the time evolution of the cumulative CPU time
$T_{C}(t)$ necessary for the simulations. As in the 1D case [Figs. 4.4(d) and
4.5(b)], the fixed lattice simulation takes longer to complete than the
dynamic lattice ones. The dynamic lattice approach with a threshold (5.3)
$\delta_{\mathcal{D}}=10^{-8}$ gives the smallest $T_{C}$ value as expected.
Overall, a gain (4.2) $g\approx 36\%$, $40\%$ and $43\%$ is seen for the
$\delta_{\mathcal{D}}=10^{-12}$, $\delta_{\mathcal{D}}=10^{-10}$ and
$\delta_{\mathcal{D}}=10^{-8}$ threshold values respectively. In the remaining
part of this chapter, we will use $\delta_{\mathcal{D}}=10^{-8}$ for our
simulations.

###  5.3 The maximum Lyapunov exponent

In order to avoid repetition of similar information in the rest of this
chapter, let us note that we perform simulations of lattice up to sizes
$451\times 451$ in order to avoid finite size effects, integrating the
equations of motion and the variational equations up to final time
$t_{f}\approx 2\times 10^{6}$ by the $s11\mathcal{ABC}6$ SI (3.84). The choice
of the integration time step, $\tau=0.1-0.25$, leads to relative energy
$E_{r}$ (3.92) and norm $S_{r}$ (3.93) errors always remaining below $10^{-3}$
and $10^{-2}$ respectively. We excite all sites in a square of side $L$ at the
center of the lattice with the same norm per site $s$ (or the same rescaled
norm $x=\beta s$), averaging our observables over $50$ disorder realizations.
Notice that similar to what was done in the 1D DDNLS system, when $\beta>0$ in
$\mathcal{H}_{2D}$ (2.20), we set the disorder parameter $\epsilon_{l,m}$ of
lattice’s central site to zero. In that way, single site excitation of all
disorder realizations with the same total norm $\mathcal{S}_{2D}$ (2.23) has
the same total energy value (2.20) $\mathcal{H}_{2D}=\beta\mathcal{S}_{2D}/2$,
which always belongs to the thermal phase of the energy-norm space in Fig.
2.6. On the other hand, for multi-site wave packet excitations, all disorder
realizations have the same value of the total norm $\mathcal{S}_{2D}$ (2.23)
and we ensure that the values of the total energy $\mathcal{H}_{2D}$ (2.20)
remain bounded within the thermal phase of the energy-norm space of the 2D
DDNLS system. As an additional remark, let us note that we are using the
energy-norm space of Fig. 2.6 of the 1D DDNLS model (2.17) for the 2D case
(2.20), as they are very similar to each other. In fact, the only difference
is that the extra dimension of the 2D system renormalizes the ground states of
Fig. 2.6 in the thermodynamic limit, which is now located well below that of
the 1D system. Nevertheless, as we are not interested to the dynamics of the
system close to the ground states, it makes Fig. 2.6 sufficient for exploring
the thermal states of the 2D DDNLS Hamiltonian model (2.20).

####  5.3.1 The weak chaos regime

Let us now investigate the chaotic dynamics of wave packet spreading in the
weak chaos regime of the 2D DDNLS system (2.20). We consider the following
cases:

  1. $\bullet$

Case W$1_{2D}$. $W=10$, $\beta=0.15$, $L=2$, $x=0.15$.

  2. $\bullet$

Case W$2_{2D}$. $W=10$, $\beta=0.92$, $L=1$, $x=0.92$.

  3. $\bullet$

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

  4. $\bullet$

Case W$4_{2D}$. $W=12$, $\beta=1.75$, $L=1$, $x=1.75$.

Let us explain why these cases theoretically belong to the weak chaos regime
of the 2D DDNLS model of Figs. 2.5 and 2.6. In the case of W$1_{2D}$, which is
a multi-site initial excitation, we apply the first condition of Eq. (2.56),

$x\frac{L^{2}}{\overline{V}}<d,$ (5.8)

where $x=\beta s$ is the initial frequency shift, $\overline{V}$ is the
average localization volume, $d$ is the average spacing between NMs in a NM
volume and $L^{2}$ is the area of the initial excitation (note that for the 1D
DDNLS model, we have $L$ instead of $L^{2}$). The estimations of
$\overline{V}$ (2.61), and $d$ (2.16) are listed in Table 5.1 for the values
of the disorder strength $W$ we use. It is worth mentioning that the
$\overline{V}$ is not obtained through the analysis of the NM themselves, but
rather by evolving a single site excitation within the 2D DDLS model (i.e.
$\mathcal{H}_{2D}$ (2.20) with $\beta=0$), up to a time when a halt in the
evolution is observed. Then, the average participation number $P$ (3.6) over
several disorder realizations of the norm density distribution $\xi_{l,m}$
(5.1) is a rough measure of the NM localization volume. Thus, using the values
$x=0.15$, $L=2$, $\overline{V}=35$ and $d=0.51$ we obtain the inequality
$0.017<0.51$ from Eq. (5.8) which is true.

2D NM properties | $W=10$ | $W=10.5$ | $W=11$ | $W=12$ | $W=14$
---|---|---|---|---|---
$\overline{V}$ | $35$ | $28$ | $23$ | $16$ | $10$
$\Delta$ | $18$ | $18.5$ | $19$ | $20$ | $22$
$d=\Delta/\overline{V}$ | $0.51$ | $0.66$ | $0.83$ | $1.25$ | $2.2$
Table 5.1: Average NM localization volume $\overline{V}$, frequency bandwidth
$\Delta=W+8$ and average spacing between NM $d=\Delta/\overline{V}$ (2.16) in
the case of a central multi-site excitation of width
$L^{2}\approx\overline{V}$ sites at the center of the 2D DDNLS lattice system
(2.20). Note that the values of $\overline{V}$ are obtained by averaging the
participation number $P$ (3.6) of the excited part at the end of the
integration of the 2D DDLS model (i.e. $\beta=0$ in $\mathcal{H}_{2D}$ (2.20))
up to final times $t_{f}\approx 10^{5}-10^{6}$ over $10$ disorder
realizations.

The remaining initial conditions of cases W$2_{2D}$, W$3_{2D}$ and W$4_{2D}$
correspond to single site excitations. Consequently, we apply the first
condition of Eq. (2.57)

$x<\Delta,$ (5.9)

such that the $x=0.92$, $x=1$ and $x=1.75$ of cases W$2_{2D}$, W$3_{2D}$ and
W$4_{2D}$ are respectively always smaller than the frequency bandwidths
$\Delta=18$, $18$ and $20$ of $W=10$, $10$ and $12$ seen in Table 5.1.

In Fig. 5.8(a), we plot the time evolution of the second moment $m_{2}(t)$
(5.5) of the norm distribution $\xi_{l,m}$ (5.1). For all cases, a diffusive
behavior of the wave packet is seen, which is well fitted by the monotonic
power law $m_{2}\propto t^{\alpha_{m}}$ with $0<\alpha_{m}<1$ in the last
decades of the evolution. The time dependence of the power law exponents
$\alpha_{m}(t)$ is shown in Fig. 5.8(b) for each studied case. A clear
saturation of the $\alpha_{m}(t)$ values to $\alpha_{m}\approx 0.2$ is
visible, in agreement with the predictions of [Flach, 2010; Laptyeva et al.,
2012] (see also Sec. 2.5). This power exponent value is numerically retrieved
from our data, averaging all the $\alpha_{m}(t)$ values in Fig. 5.8(b) after
the $6^{\text{th}}$ decade where a clear convergence is seen. This power law
exponent was also numerically obtained by [Laptyeva et al., 2012] for the 2D
DKG system (2.27).

Figure 5.8: Results of the time evolution of the second moment $m_{2}(t)$
(5.5) [panel (a)] of the norm density distribution $\xi_{l,m}$ (5.1) and of
the related ftMLEs $\Lambda(t)$ (3.21) [panel (c)] for the weak chaos cases
W$1_{2D}$ [(b)lue curves], W$2_{2D}$ [(g)reen curves], W$3_{2D}$ [(r)ed
curves] and W$4_{2D}$ [(c)yan curves] (see text for details). In panels (b)
and (d) we see the associated power law exponents when the data are fitted
with the relations $m_{2}(t)\propto t^{\alpha_{m}}$ and $\Lambda(t)\sim
t^{\alpha_{\Lambda}}$ respectively. The black dashed lines guide the eye for
slopes $0.2$ [panel (a)] and $-0.37$ [panel (c)]. The horizontal dashed lines
indicate the values $0.2$ in panel (b) and $-0.37$ in panel (d). Each case has
been averaged over $50$ disorder realizations.

For all these weak chaos cases, we evaluate the wave packet’s chaotic
strength, through the computation of their ftMLEs $\Lambda(t)$ (3.21) depicted
in Fig. 5.8 (c). A clear decay of the $\Lambda(t)$ values is visible for all
the cases following a power law $\Lambda\propto t^{\alpha_{\Lambda}}$
($\alpha_{\Lambda}<0$) similarly to what was observed in the 1D case [Sec.
4.3.2]. The values of the time evolution of $\alpha_{\Lambda}(t)$ are
presented in Fig. 5.8(d), showing a tendency to saturate around
$\alpha_{\Lambda}\approx-0.37$. This result of the wave packet chaoticity is
not particular to the 2D DDNLS system (2.20), as it was also found for the 2D
DKG model (2.27) [Many Manda et al., 2020].

####  5.3.2 The strong chaos regime

We also study the strong chaos regime of the 2D DDNLS system (2.20) of Figs.
2.5 and 2.6, by considering the following cases:

  1. $\bullet$

Case S$1_{2D}$. $W=10.5$, $\beta=0.145$, $L=21$, $x=0.145$.

  2. $\bullet$

Case S$2_{2D}$. $W=11$, $\beta=0.68$, $L=10$, $x=0.68$.

  3. $\bullet$

Case S$3_{2D}$. $W=12$, $\beta=0.35$, $L=21$, $x=0.35$.

  4. $\bullet$

Case S$4_{2D}$. $W=14$, $\beta=6$, $L=15$, $x=0.72$.

Figure 5.9: Similar to Fig. 5.8, but for the strong chaos cases S$1_{2D}$
[(b)lue curves], S$2_{2D}$ [(g)reen curves], S$3_{2D}$ [(r)ed curves] and
S$4_{2D}$ [(c)yan curves] (see text for details). The dashed lines in panels
(a) and (c) guide the eye for slopes $0.33$ and $-0.46$ respectively. The
horizontal dashed lines indicate the value $0.33$ in panel (b) and $-0.46$ in
panel (d). The presented results are averaged over $50$ disorder realizations.

Note that here, we use extended localized initial excitations as the value of
$L^{2}$ is much greater than our estimation of the NM average localization
volume shown in Table 5.1. In Sec. 2.6, we discussed a particular condition
(2.119) for setting up the initial wave packet within the strong chaos region
of Figs. 2.5 and 2.6, which for the 2D DDNLS system becomes

$L^{2}\gtrsim\frac{W}{x}.$ (5.10)

The advantage of this approach is that it does not depend on the
characteristics of the NMs of the 2D DDLS system (i.e. $\beta=0$ in
$\mathcal{H}_{2D}$ (2.20)), which have not yet been extensively studied.
Consequently, using the parameter values of the cases S$1_{2D}$, S$2_{2D}$,
S$3_{2D}$ and S$4_{2D}$ we compute $L^{2}$ to respectively be
$441,\leavevmode\nobreak\ 100,441$ and $225$. These values are always much
greater than $W/x\approx 72$, $16$, $34$ and $19$ respectively so that
condition (5.10) is satisfied.

Let us note that, both the conditions of strong chaos in Eq. (2.56)

$x\frac{L^{2}}{\overline{V}}>d,$ (5.11)

which depends on the NM characteristics ($\overline{V}$ (2.61) and $d$
(2.16)), and the one of Eq. (5.10) are equivalent. Indeed, using the following
simple algebraic manipulations

$x\frac{L^{2}}{\overline{V}}>d\Rightarrow
x\frac{L^{2}}{\overline{V}}>\frac{\Delta}{\overline{V}}\Rightarrow
xL^{2}>\Delta\Rightarrow L^{2}\gtrsim\frac{W}{x},$ (5.12)

[where we assumed $\Delta\sim W$ for the 2D system (2.20)] we can relate the
two inequalities. In order to be completely certain about our argumentation,
we have also checked using Eq. (5.11) that our cases qualitatively belong to
the strong chaos region of the 2D DDNLS model (2.20). Indeed, using the
parameter values for the cases S$1_{2D}$, S$2_{2D}$, S$3_{2D}$ and S$4_{2D}$
and the NM average localization volume $\overline{V}$ in Table 5.1, we compute
$xL^{2}/\overline{V}$ to respectively get $0.145\cdot 21^{2}/28\approx 2.28$,
$0.68\cdot 10^{2}/23\approx 2.96$, $0.35\cdot 21^{2}/16\approx 9.65$ and
$0.72\cdot 15^{2}/10\approx 16.2$. These numbers are always much higher than
the $d$ values $0.66$, $0.83$, $1.25$ and $2.2$ for the disorder strengths
$W=10.5$, $11$, $12$ and $14$.

The result of the computation of the second moment $m_{2}(t)$ (5.5) of the
norm distribution $\xi_{l,m}$ (5.1) is shown in Fig. 5.9 (a) for all the
strong chaos cases. Again, a spreading of wave packets appears for all the
cases as the $m_{2}(t)$ values monotonically increase for all times. These
increases are well fitted by a power law $m_{2}(t)\propto t^{\alpha_{m}}$,
with $0<\alpha_{m}<1$. The time evolution of the power law exponents
$\alpha_{m}(t)$ are shown in Fig. 5.9(d) for all these cases. Clearly, the
values of $\alpha_{m}(t)$ reach $\alpha_{m}\approx 0.33$ at some point of
their evolution, in agreement with the theory [Flach, 2010; Laptyeva et al.,
2012] (see also Sec. 2.5). Results about of the time dependence of the ftMLEs
$\Lambda(t)$ (3.21) for all the strong chaos cases are shown in Fig. 5.9(c).
Again, a monotonic decay is evident for the ftMLEs, well fitted by
$\Lambda(t)\propto t^{\alpha_{\Lambda}}$ with $\alpha_{\Lambda}<0$. The values
of the power law exponents converge toward $\alpha_{\Lambda}\approx-0.46$
[Fig. 5.9(d)]. Both estimations $\alpha_{m}=0.33$ and
$\alpha_{\Lambda}=-0.46$, are obtained from averaging results from all the
cases studied in the interval of time where the saturation is seen. These
behaviors of the wave packet’s chaoticity in the strong chaos regime were also
found for the 2D DKG model (2.27) in [Many Manda et al., 2020].

###  5.4 The deviation vector distribution

The computation of the spatiotemporal evolution of the DVD obtained during the
computation of the ftMLEs was used in the case of 1D DDNLS model (2.17) to
visualize the motion of chaotic seeds in the lattice’s excited part. We
perform a similar dynamical investigation for the DVD (5.2) in the 2D DDNLS
model (2.20).

####  5.4.1 The weak chaos regime

Let us first study the behavior of the DVD $\xi_{l,m}^{D}$ (5.2) in the weak
chaos regime. In Fig. 5.10 we show the time evolution of the norm density
distribution $\xi_{l,m}$ (5.1) [panels (a), (c) and (e)] and the DVD
$\xi_{l,m}^{D}$ (5.2) [panels (b), (d) and (f)] for a representative disorder
realization of the $W1_{2D}$ case. The snapshots of the distributions are
taken at times $\log_{10}t\approx 3.15$ [panels (a) and (b)],
$\log_{10}t\approx 5.04$ [panels (c) and (d)] and $\log_{10}t\approx 6.09$
[panels (e) and (f)].

Figure 5.10: 3D spatial distribution and the associated 2D color map (upper
sides) snapshots of the norm density $\xi_{l,m}$ (5.1) [panels (a), (c) and
(e)] and the DVD $\xi_{l,m}^{D}$ (5.2) for a representative realization of the
case W$1_{2D}$ (see text for details). Panels (a) and (b) correspond to the
configurations of the $\xi_{l,m}$ and $\xi_{l,m}^{D}$ at time
$\log_{10}t\approx 3.15$, panels (c) and (d) at $\log_{10}t\approx 5.04$ and
panels (e) and (f) at $\log_{10}t\approx 6.09$. The red curves on the sides
are the distributions’ projections along the $l$ and $m$ directions. The white
areas in the upper 2D color map represent the motion of the distribution’s
center. The coloring of the distributions is done with respect to the
$\log_{10}\xi_{l,m}$ and $\log_{10}\xi_{l,m}^{D}$ values according to the
horizontal color scales presented at the top of the figure.

The results of Figs. 5.10(a), (c) and (e) clearly show that the norm
distribution is spreading as its extent (represented for instance by the red
curves of the side projections of $\xi_{l,m}$) grows rather smoothly and
symmetrically about the wave packet center which remains practically fixed
[see white region at the center of the 2D color map on the upper sides of
Figs. 5.10(a), (c) and (e)]. On the other hand, the DVD $\xi_{l,m}^{D}$ (5.2)
keeps a rather localized pointy shape (although its extent slightly increases
in time), which remains well inside the wave packet width [see for example the
red curves of the side projections in Figs. 5.10(b), (d) and (f)]. Its center
is meandering inside the lattice’s excited part with increasing amplitudes
[see white areas at the upper sides of the 2D color maps in Figs. 4.11(b), (d)
and (f)].

Figure 5.11: Time evolution of observables of the tangent dynamics in the weak
chaos regime. (a) The second moment $m_{2}^{D}$ (5.5) and (b) the
participation number $P^{D}$ (3.6) of the DVD $\xi_{l,m}^{D}$ (5.2). The
dashed lines in panels (a) and (b) guide the eye for slopes $0.12$ and $0.0$
respectively. The curve’s colors correspond to the ones used in Fig. 5.8.
Figure 5.12: Time dependence of the area $A$ (5.13) covered by the DVD’s
center $\bm{r}^{D}=\left(\overline{l}^{D},\overline{m}^{D}\right)$ in cases
belonging to the weak chaos regime of the 2D DDNLS system (2.20). The dashed
line indicates slope $0.5$. The curve’s colors correspond to the ones used in
Fig. 5.8.

In order to quantify the observations made about the DVD’s behavior in Fig.
5.10, we compute the second moment $m_{2}^{D}(t)$ (5.5) and participation
number $P^{D}(t)$ (3.6) of $\xi_{l,m}^{D}$ (5.2). The obtained results are
shown in Figs. 5.11(a) and (b) respectively. A steady increase of the
$m_{2}^{D}(t)$ values [Fig. 5.11(a)] is observed for all the weak chaos cases,
which is well fitted by $m_{2}^{D}(t)\propto t^{0.12}$, exhibiting a much
slower growth than the $m_{2}(t)\propto t^{0.2}$ behavior of the wave packet
expansion observed in Figs. 5.8(a) and (b). In addition, a clear saturation of
the $P^{D}(t)$ values is visible at the later stage of the evolution, which
tend to fluctuate around numbers of the order of the NM average localization
volume. For instance, for cases W$1_{2D}$ [(b)lue curve], W$2_{2D}$ [(g)reen
curve] and W$3_{2D}$ [(r)ed curve] in Fig. 5.8(b), the DVD’s participation
number attains values $P^{D}\approx 10^{1.5}\approx 31$ similar to the
$\overline{V}\approx 35$ observed for $W=10$ in Table 5.1. This behavior is
due to the localized pointy shapes of the $\xi_{l,m}^{D}$ (5.2) seen in Figs.
5.10(b), (d) and (f) whose extents approximately span the average volume of a
NM. Note that larger fluctuations of the $P^{D}(t)$ values in Figs. 5.8(b)
compared to the 1D case in Fig. 4.12(b) can be attributed to the smaller
number of disorder realizations used in the 2D case and to the fact that the
geometrical characteristics of the NMs in 2D lattices is more complex compared
to the 1D case, leading for instance to larger discrepancies in the
measurement of the NM localization length.

We also investigate the motion of the DVD’s center
$\overline{\bm{r}}^{D}=\left(\overline{l}^{D},\overline{m}^{D}\right)$, by
computing the quantity

$A(t)=R_{x}(t)\cdot R_{y}(t),$ (5.13)

at each time $t>0$ of the evolution. Here,

$R_{x}(t)=\max_{[0,t]}\left\\{\overline{l}^{D}(t)\right\\}-\min_{[0,t]}\left\\{\overline{l}^{D}(t)\right\\},\qquad
R_{y}(t)=\max_{[0,t]}\left\\{\overline{m}^{D}(t)\right\\}-\min_{[0,t]}\left\\{\overline{m}^{D}(t)\right\\}.$
(5.14)

similar to the quantity $R$ (4.10) used in the 1D case. The $A$ quantity
(5.13) tries to quantify the area visited by the DVD (whose extent possesses
fuzzy borders), by measuring the surface of the smallest rectangle containing
the white area in panels (b), (d) and (f) of Figs. 5.10 and 5.13.

The fate of the wave packet can be conjectured based on the behavior of the
area $A$ (5.13) in a similar way to the one followed in the 1D case for the
length $R$ (4.10). Indeed, if

$A\lesssim\overline{V},\quad A(t)\sim t^{0},$ (5.15)

the DVD does not visit sites outside the initially excited central ones. In
this case, no spreading of the wave packet will be observed. On the other
hand, when

$A(t)>\overline{V},$ (5.16)

after a certain amount of time, two scenarios are possible. If the values of
$A$ saturate i.e $A(t)\sim t^{0}$ after some time, then the DVD’s center has
only reached sites in close proximity of the central ones. Thus, a
subdiffusion of the wave packet will be visible, followed by a slowing down of
the wave packet spreading which will eventually halt forever. This corresponds
to the situation discussed in [Johansson et al., 2010; Aubry, 2011; Pikovsky
and Fishman, 2011]. On the other hand, if the $A(t)$ values show a monotonic
increase such that at large times $A(t)\gg V$, then the DVD center is visiting
sites far away from the initially excited ones, leading to the case of wave
packet which is spreading indefinitely as described in [Pikovsky and
Shepelyansky, 2008; Garcia-Mata and Shepelyansky, 2009; Flach, 2010; Laptyeva
et al., 2012].

The time evolution of the $A(t)$ (5.14) values is depicted in Fig. 5.12 for
all the weak chaos cases of Fig. 5.8. A clear monotonic increase of $A(t)$ is
visible. The values of $A(t)$ are well fitted by a power law $A(t)\propto
t^{\alpha_{A}}$ with $\alpha_{A}\approx 0.5$ [black dashed line in Fig. 5.12].
The slight slowing down observed for cases W$3_{2D}$ and W$4_{2D}$ [(r)ed and
(c)yan curves in Fig. 5.12] at the end of the integration is likely due to the
fact that the integration time is not enough to display a firm asymptotic
behavior of the chaotic seeds. This is supported by the fact that for all
cases, the second moment $m_{2}(t)$ (5.5) in Figs. 5.8(a) and (b) showed that
the system has truly reached the weak chaos regime around the last decade of
the evolution. Nevertheless, the fact that the $A(t)$ values monotonically
increase up to $A(t)\propto 10^{2.5}\approx 300$ sites at the final stages of
the simulations, which correspond to areas much larger than the estimates of
the NM localization volume is a strong evidence that the chaotic seeds are
visiting sites far away from the lattice’s center.

####  5.4.2 The strong chaos regime

We also investigate the behavior of the DVD (5.2) in the strong chaos regime.

Figure 5.13: Similar to Fig. 5.10 but for a representative realization of the
S$1_{2D}$ strong chaos case. Panels (a) and (b) correspond to the
configurations of the $\xi_{l,m}$ (5.1) and $\xi_{l,m}^{D}$ (5.2) at
$\log_{10}t\approx 3$, panels (c) and (d) at $\log_{10}t\approx 4.8$ and
panels (e) and (f) at $\log_{10}t\approx 5.90$.

In Fig. 5.13, we show the evolution of the norm density $\xi_{l,m}$ (5.1) and
its related DVD $\xi_{l,m}^{D}$ (5.2) at three different times
$\log_{10}t\approx 3.0$ [panels (a) and (b)], $\log_{10}t\approx 4.8$ [panels
(c) and (d)] and $\log_{10}t\approx 5.9$ [panels (e) and (f)] for a
representative disorder realization of the $S1_{2D}$ case. Similarly to the
weak chaos case [Figs. 5.10(a), (c) and (e)], in Figs. 5.13(a), (c) and (e) we
see that the norm density $\xi_{l,m}$ (5.1) spreads more or less evenly around
its center [white region in the color map of the upper sides of Figs. 5.13(a),
(c) and (e)] which remains practically fixed. On the other hand the DVD
$\xi_{l,m}^{D}$ (5.2) maintains a rather localized pointy shape. In addition,
its center meanders inside the wave packet interior with increasing amplitude,
as indicated by the always growing white regions at the center of the upper
sides of Fig. 5.13 (b), (d) and (f). These oscillations are much larger
compared to what was observed in the weak chaos case [see Figs. 5.10(b), (d)
and (f)].

In order to substantiate these observations, in Figs. 4.14(a) and (b) we
respectively show the time evolution of the DVD’s second moment $m_{2}^{D}(t)$
(5.5), and participation number $P^{D}(t)$ (3.6). A slow increase of the
values of the $m_{2}^{D}(t)$, which are well fitted by the power law
$m_{2}^{D}(t)\propto t^{0.17}$ is observed, much smaller than the wave
packet’s second moment evolution $m_{2}(t)\propto t^{0.33}$ in Figs. 5.9(a)
and (b). The saturation of the $P^{D}$ values is clearly seen in Fig. 5.14(b),
where fluctuations in the interval $P^{D}\approx 10$ for S$4_{2D}$ [(c)yan
curve in Fig. 5.14(b)] up to $P^{D}\approx 10^{1.4}\approx 25$ in the case of
S$1_{2D}$ [(b)lue curve in Fig. 5.14(b)] are observed. In both cases the
$P^{D}$ values are of the order of the NM localization volume
$\overline{V}\approx 10$ and $28$, respectively at $W=14$ and $10.5$ [Table
5.1]. Consequently, the DVD $\xi_{l,m}^{D}$ (5.2) remains localized during the
system’s evolution, with its extent reaching approximately the NM’s average
localization volume. Similarly to the weak chaos case, we understand the
larger fluctuations in the calculations of the $m_{2}^{D}(t)$ (5.5) and
$P^{D}(t)$ (3.6) based on the rather small number of considered disorder
realizations and the complex geometry of 2D NMs.

Figure 5.14: Similar to Fig. 5.11, but for the strong chaos regime of the 2D
DDNLS system (2.20). The dashed lines in panels (a) and (b) guide the eye for
slopes $0.17$ and $0.0$ respectively. The curve’s colors correspond to the
ones used in Fig. 5.9.  Figure 5.15: Similar to Fig. 5.12, but for the strong
chaos regime of the 2D DDNLS system (2.20). The dashed line indicates slope
$0.55$. The curve’s colors correspond to the ones used in Fig. 5.9.

The time dependence of the area $A(t)$ (5.14) visited by the DVD center is
shown in Fig. 5.15 for all the strong chaos cases of Fig. 5.14. As in the case
of weak chaos, a monotonic increase of the $A$ values is observed for all
cases. At the final integration time ($\log_{10}t=10^{6}$), this quantity has
reached values $A\approx 10^{3.5}\approx 3162$ for all cases, which are much
larger than any NM average localization volume [see Table 5.1]. In addition,
the values of $A(t)$ are well fitted by a power law $A\propto t^{\alpha_{A}}$
with $\alpha_{A}\approx 0.55$. Consequently, this provides a strong indication
that the chaotic seeds are consistently visiting sites far from the central
part of the lattice. It is worth mentioning that similar results was also
obtained for the 2D DKG system (2.27) in [Many Manda et al., 2020].

###  5.5 Discussion and summary

Let us now summarize the main features of the chaotic dynamics of the wave
packet spreading in the 2D DDNLS lattice (2.20). The chaos strength of the
wave packet spreading was performed through the computation of the system’s
ftMLE $\Lambda(t)$ (3.21). More specifically, we found that the
$\Lambda(t)\propto t^{-0.37}$ in the weak and $\Lambda(t)\propto t^{-0.46}$
for the strong chaos regimes. These behavior are clearly different from what
is observed in the case of regular motion ($\Lambda\propto t^{-1}$),
confirming that the wave packet remains chaotic although its chaos strength
decreases in time. As to the question of whether wave packet chaotization is
fast enough to support the wave packet subdiffusion, we can compare the wave
packet’s spreading timescale $T_{M}=1/D$, with $D=dm_{2}/dt$ to its Lyapunov
time $T_{L}=1/\Lambda$ (4.15). By doing that, we find that since $m_{2}\sim
t^{\alpha_{m}}$ and $\Lambda\sim t^{\alpha_{\Lambda}}$

$\frac{T_{M}}{T_{L}}\sim t^{1-\alpha_{\Lambda}-\alpha_{m}},$ (5.17)

which leads to

$\frac{T_{M}}{T_{L}}\sim t^{0.43},$ (5.18)

for the weak chaos regime where $\alpha_{m}=1/5$ and $\alpha_{\Lambda}=-0.37$,
and to

$\frac{T_{M}}{T_{L}}\sim t^{0.21},$ (5.19)

for the strong chaos regime where $\alpha_{m}=1/3$ and
$\alpha_{\Lambda}=-0.46$. Thus, $T_{L}$ is always smaller than $T_{M}$, such
that the wave packet chaotization happens faster than its spreading to more
sites.

By examining the ratios $T_{M}/T_{L}$ obtained in all spreading dynamical
regimes for both the 1D and 2D DDNLS models i.e.

$\frac{T_{M}}{T_{L}}\sim t^{0.43}\quad\mbox{(weak chaos),
}\quad\frac{T_{M}}{T_{L}}\sim t^{0.21}\quad\mbox{(strong chaos)}$ (5.20)

obtained in Eq. (5.18) and (5.19) for the 2D DDNLS model (2.20) and

$\frac{T_{M}}{T_{L}}\sim t^{0.42}\quad\mbox{(weak chaos),
}\quad\frac{T_{M}}{T_{L}}\sim t^{0.20}\quad\mbox{(strong chaos)},$ (5.21)

observed for the 1D DDNLS system (2.17) in Sec. 4.5, we see a remarkable
closeness of exponents in Eqs. (5.20) and (5.21). This is a strong indication
that nonlinear interactions of the same nature are responsible for the wave
packet spreading in the cases of both spatial dimensions within the thermal
region of the energy-norm parameter space. This is in agreement with the
theoretical calculations performed in Sec. 2.5. Indeed, for both the 1D and 2D
DDNLS systems, we ended up with the same reduced equation of the asymptotic
dynamics

$i\frac{\partial\phi_{\nu}}{\partial t}\sim
E_{\nu}\phi_{\nu}+\beta\mathcal{P}\left(\beta\phi_{\nu_{1}}^{2}\right)\phi_{\nu_{1}}^{3}f(t).$
(5.22)

Here, the nonlinear term
$\beta\mathcal{P}\left(\beta\phi_{\nu_{1}}^{2}\right)\phi_{\nu_{1}}^{3}f(t)$
is responsible of the thermalization of more sites in both lattices. Note that
in Eq. (5.22) $\nu$ and $\nu_{1}$ represent respectively the non-thermalized
and thermalized modes. Thus, only the geometry of the problem affects the
spreading and the chaotic dynamics of the wave packet.

Notice that equating the $T_{M}/T_{L}$ ratios of Eqs. (5.20) and (5.21) for
each dynamical regime, we can write

$\frac{\Lambda^{1}}{m_{2}^{1}}=\frac{\Lambda}{m_{2}},$ (5.23)

where the subscript ‘1’ refers to the 1D lattice. Equation (5.23) allows the
prediction of the characteristics of chaos of the wave packet dynamics. For
example, using the relations describing the time evolution of the wave packet
second moment $m_{2}$ in Chap. 2 and the results obtained for the computation
of the ftMLEs of the 1D DDNLS system in Chap. 4, we get from Eq. (5.23)

$\Lambda\sim t^{\alpha_{m}-\alpha_{m}^{1}+\alpha_{\Lambda}^{1}}.$ (5.24)

Then, we obtain for the behavior of the ftMLE in the case of the 2D lattice

$\Lambda\sim t^{-0.38},$ (5.25)

for the weak chaos regime where $\alpha_{m}^{1}=1/3$, $\alpha_{m}=1/5$ and
$\alpha_{\Lambda}=-0.25$, and

$\Lambda\sim t^{-0.47},$ (5.26)

in the strong chaos case where $\alpha_{m}^{1}=1/2$, $\alpha_{m}=1/3$ and
$\alpha_{\Lambda}=-0.3$. Both laws (5.25) and (5.26) are extremely similar to
the $\Lambda\sim t^{-0.37}$ and $\Lambda\sim t^{-0.46}$ found numerically in
Sec. 5.3.

Similarly to what was seen in the 1D case, the evolution of the DVD revealed
the mechanism of wave packet chaotization: Localized chaotic seeds randomly
oscillate in the wave packet’s interior, homogenizing chaos for the lattice’s
excited region and thermalizing neighboring initially unexcited sites. Indeed,
our numerical computations showed a DVD’s second moment $m_{2}^{D}(t)$ (5.5)
which slowly increases in time following a $m_{2}(t)\propto t^{0.12}-t^{0.17}$
law. In addition, the DVD’s participation number $P^{D}$ saturated at values
of the order of the NM average localization volume. Thus, the DVD always
retains a rather localized and pointy shape, which randomly meanders inside
the wave packet’s interior in both the weak and strong chaos regimes.

The numerical computations of the area $A$ (5.13) visited by the the DVD’s
center showed a monotonically increasing relation with time, reaching values
$A\approx 300$ sites in the weak and $A\approx 3000$ sites in the strong chaos
regimes for final integration time $t_{f}\approx 2\times 10^{6}$. Fitting the
variations of $A$ with the power law $t^{\alpha_{R}}$, we found
$\alpha_{A}\approx 0.5$ (weak chaos) and $\alpha_{A}\approx 0.55$ (strong
chaos). The fact that we observed larger values of $A$ and larger increase
rates in the strong chaos case compared to the weak chaos one is in agreement
with the results of Chap. 2 [see e.g. Flach, 2010]. Indeed, as in the 1D DDNLS
lattice, the faster spreading in the strong chaos case implies that chaotic
seeds visit larger regions of the lattice as faster homogenization of chaos
and thermalization of more lattice sites is required.

The fact that all these behaviors have been found in both the 1D and 2D
lattices of the DDNLS and DKG models is a clear indication of the generality
of our findings.

##  Chapter 6 Summary and outlook

This work presented results regarding the contribution of chaotic dynamics to
wave packet spreading in disordered multidimensional Hamiltonian models. In
particular, we focussed on the well known one (1D) and two dimensional (2D)
disordered discrete nonlinear Schrödinger (DDNLS) models, which are generic in
modelling physical systems such as the conductivity of metals, the propagation
of light in optical media or the expansion of trapped Bose-Einstein
condensates. The aim of this work was to give more insights into the long
standing question of whether a wave packet spreading in a disordered medium
will expand forever or halt after a certain amount of time, based on the
computation of observables of the system’s tangent dynamics such as the
maximum Lyapunov exponent (MLE) and its related deviation vector distribution
(DVD). This question is of utmost importance as it links up to several
physical aspects like the ergodic, chaotic and statistical properties of
multidimensional nonlinear systems which are still widely debated. In the next
paragraphs, we recapitulate our work and give some of the possible open
avenues along which this study could be taken forward.

As a preliminary step to reaching our goal, in Chap. 2 we presented the theory
behind the wave packet spreading in the 1D and 2D DDNLS systems. Firstly, we
outlined the so called Anderson localization (AL), whereby in a disordered
linear model, localized excitations remain in that state forever, not
diffusing through the lattice. Then, we discussed what happens to the AL when
nonlinearity is introduced into the system. Upon examining the dynamical
behavior of normal modes (NMs), we found that when AL prevails, the NMs
possess a discrete frequency spectrum coupled to a spatially localized nature
which gets more pronounced when the disorder becomes stronger in the system.
On the other hand, when nonlinearity is switched on, the NMs overlap and
interact, allowing the energy to move from mode to mode in a deterministic
chaotic manner. Consequently, the AL does not survive, leading to an avenue
for a subdiffusive spreading of the initially localized excitation. Depending
on the strength of the overlapping part, which is expressed as a function of
the norm $\mathcal{S}$, the total energy $\mathcal{H}$, the magnitude of the
nonlinear coefficient $\beta$ and the disorder strength $W$, we discussed how
two different subdiffusive behaviors of the localized excitation can take
place. In particular, the so-called weak chaos regime where the wave packet
extent, measured through the second moment $m_{2}$ (3.5) [participation number
$P$ (3.6)], increases according to $m_{2}\sim t^{\alpha_{m}}$ [$P\sim
t^{\alpha_{P}}$] with $\alpha_{m}=1/3$ [$\alpha_{P}=1/6$] for the 1D DDNLS
model and $\alpha_{m}=1/5$ [$\alpha_{P}=1/5$] in the case of the 2D system. In
addition, the so-called strong chaos regime, with a faster wave packet
expansion where the power law exponents for the $m_{2}$ [$P$] are
$\alpha_{m}=1/2$ [$\alpha_{P}=1/4$] in the case of the 1D and $\alpha_{m}=1/3$
[$\alpha_{P}=1/3$] for the 2D DDNLS system can also appear.

In order to investigate the chaotic dynamics of the spreading regimes
mentioned above, we mainly relied on numerical experiments. This is motivated
by the difficulty involved with an analytical treatment of the stochasticity
appearing both in the system’s equations of motion and tangent dynamics. Due
to the complexity of integrating the equations of motion of the DDNLS models
compared to other multidimensional lattice models like the disordered Klein-
Gordon (DKG) and Fermi-Pasta-Ulam-Tsingou (FPUT) chains of nonlinear
oscillators, no particular effort in designing efficient numerical methods for
integrating the variational equations of the former were undertaken to date.
Thus, a large part of Chap. 3 is dedicated to designing efficient, accurate
and fast numerical integration techniques for the long time integration of the
system’s equations of motion and the variational equations required in this
study. We exploited symplectic integrators (SIs), which are numerical schemes
specifically built for integrating the equations of motion of Hamiltonian
models, along with the tangent map method [see e. g. Skokos and Gerlach, 2010]
in order to evolve the variational equations. Note that SIs keep the computed
energy of the system constant up to a certain level of accuracy which can be
associated with the precision of the integrator. We presented methods based on
the two part (only for the 1D DDNLS model) and three part split approaches of
the Hamiltonian functions of the 1D and 2D DDNLS systems and compared their
performances based on standard schemes like the $DOP853$ [Hairer et al.,
1993], a Runge-Kutta (RK) integrator of order $8$ based on the Dormand and
Prince method and the TIDES algorithm based on Taylor/Lie series expansions of
the system’s field equations [Abad et al., 2012, 2014].

Then, we studied the performance of the designed SIs in Sec. 3.8 of Chap. 3
based on their ability to accurately reproduce the system’s dynamical
behavior. For instance, this was quantified by checking the constancy of the
total energy $\mathcal{H}$ and norm $\mathcal{S}$ through the computation of
the relative energy and norm errors, respectively denoted as $E_{r}$ and
$S_{r}$. We also calculated the characteristics of the wave packet spreading,
i.e. the time evolution of the second moment $m_{2}$ and the participation
number $P$, and evaluated the precision of the evolution of the variational
equations based on the computed finite time maximum Lyapunov exponent (ftMLE)
$\Lambda$ for several initial conditions. We focussed on two energy accuracy
levels: moderate accuracy with $E_{r}\lesssim 10^{-5}$, which is often
considered to be a good precision in nonlinear lattice studies [see e.g.
Laptyeva et al., 2010; Skokos et al., 2013], as well as high accuracy where
the threshold values $E_{r}\lesssim 10^{-9}$ were used [see e.g. Sales et al.,
2018]. All the tested numerical integrators produced similar results for both
the evolution of system’s phase space and tangent dynamics. Although the
$DOP853$ and the TIDES are the most accurate of all the integrators we tested,
the fact that they do not preserve the values of the computed total energy and
norm ($E_{r}$ and $S_{r}$ are not bounded from above) makes them a not ideal
choice for the long time evolution of the system’s dynamical equations. On the
other hand, we found that SIs in general take less CPU time than the $DOP853$
and TIDES in order to reproduce the same system dynamics with the three part
split SIs performing better than the two part split SIs for large lattice
sizes, i.e. when the number of degrees of freedom $N\gtrsim 70$. In particular
the $s11\mathcal{ABC}6$ and $s9\mathcal{ABC}6$ SIs of order $6$ along with the
$s17\mathcal{ABC}8$ and $s19\mathcal{ABC}8$ SIs of order $8$ proved to be the
best integrators respectively for moderate ($E_{r}\lesssim 10^{-5}$) and high
($E_{r}\lesssim 10^{-9}$) accuracies. The results of these investigations were
presented in a recent research paper [Danieli et al., 2019a].

In Chap. 4, we computed the time evolution of the finite time MLE (ftMLE) and
the spatiotemporal evolution of its related DVD for the 1D DDNLS model in the
weak and strong chaos spreading regimes. We considered a total of $8$ cases,
$4$ cases in the weak chaos and $4$ in the strong chaos regimes. For each
case, we averaged the observables of the system’s phase space and tangent
dynamics over $100$ different realizations i.e. different values of the
lattice’s disorder parameters $\epsilon_{l}$, for norm and energy values
located in the same thermal regime. In all cases, a decay of the ftMLE
$\Lambda$ following the power law $\Lambda\sim t^{\alpha\Lambda}$ was seen.
For the weak chaos case the power law exponent was fitted to
$\alpha_{\Lambda}\approx-0.25$ while in the strong chaos regime, a smaller
exponent value $\alpha_{\Lambda}\approx-0.3$ was found. As both exponents are
different from $\alpha_{\Lambda}=-1$ observed for regular motion, we concluded
that these wave packet dynamics remain always chaotic up to the largest times
of our simulations ($t_{f}\approx 10^{8}$). We correlated the different values
of $\alpha_{\Lambda}$ of the weak and strong chaos regimes to the wave packet
spreading dynamics. Indeed, as the wave packet spreading is faster in the
strong chaos, the constant total norm is shared to a larger number of
thermalized lattice sites resulting in a faster decay of the norm per site $s$
and consequently of the effective nonlinearity, which is approximately $\beta
s^{2}$. As a result, the system is driven to the integrable limit quicker in
the strong chaos regime, and the $\alpha_{\Lambda}=-0.3$ value is closer to
the regular case ($\alpha_{\Lambda}=-1$) compared to $\alpha_{\Lambda}=-0.25$
seen in the weak chaos case.

The motion of the DVD related to the computation of the MLE showed random
fluctuations of chaotic spots (called chaotic seeds) which are known to
homogenize chaos inside the excited part of the lattice and contribute to the
wave packet thermalization. The computation of the second moment $m_{2}^{D}$
and the participation number $P^{D}$ of the DVD in the weak chaos regime
showed that the DVD extent increases slowly in time as $m_{2}^{D}\sim
t^{0.14}$ and $P^{D}\sim t^{0.04}$, with $P^{D}$ showing a tendency to
saturate to values of the order of the NM’s localization volume. This process
is faster in the strong chaos case where a clear constancy of the $m_{2}^{D}$
and $P^{D}$ was visible from time $t\gtrsim 10^{5}$. In addition, we computed
the length $R$ of the region visited by the DVD center and found that this
quantity monotonically increases in time for all cases, to values $R\approx
10^{2}$ and $R\approx 10^{3}$ at the largest times of our simulations in the
weak and strong chaos regimes respectively. These values of $R$ are much
larger than the NM localization volume and therefore constitute clear
numerical evidence that modes other than the initially excited ones also
contribute to the wave packet thermalization. The obtainment of larger $R$
values in the strong chaos case is likely related to the fact that a faster
spreading of the wave packet is seen in that regime, so that chaotic seeds
have to travel to larger excited regions of the lattice. In addition, we found
that the time dependence of $R$ is well fitted by $R\sim t^{\alpha_{R}}$,
where $\alpha\approx 0.24$ for the weak chaos case, while the strong chaos
case, exhibits higher values of $\alpha_{R}$ which decrease in time toward
$\alpha_{R}\approx 0.24$. This result can be explained by taking into account
the fact that the strong chaos regime is a transient one, meaning that an
asymptotic crossover of the chaotic dynamics of the system to the weak chaos
behavior has to be seen after a finite amount of time. A large part of the
results of Chap. 4 have been reported in [Senyange et al., 2018].

Chapter 5, focussed to the same questions as in Chap. 4 but this time for the
2D DDNLS model. Notice that, to the best of our knowledge, this is the first
time that such a detailed investigation has been performed for this model
mainly due to the fact that the integration of the equations of motion and
variational equations of this system is a computationally very difficult task.
In our study, we considered $4$ initial conditions in the weak chaos and $4$
in the strong chaos regions. For each case, we averaged the observables of the
system’s phase space and tangent dynamics over $50$ disorder realizations. Our
numerical computations confirmed the theoretical expected values of the wave
packet spreading characteristics i.e. $m_{2}\sim t^{\alpha_{m}}$ with
$\alpha_{m}=1/5$ in the weak and $\alpha_{m}=1/3$ for the strong chaos regimes
respectively. Furthermore, the computation of $\Lambda$ showed a decaying
chaos strength, whose time dependence was well fitted by $\Lambda\sim
t^{\alpha_{\Lambda}}$ with $\alpha_{\Lambda}<0$ being similar to the 1D
system. Actually, the values $\alpha_{\Lambda}\approx-0.37$ and
$\alpha_{\Lambda}\approx-0.46$ were found respectively for the weak and strong
chaos cases. It is worth noticing again that these values are different from
$\alpha_{\Lambda}=-1$ which is seen in the case of regular motion. Thus the
wave packet always remains chaotic in the 2D DDNLS model as well. Furthermore
we found a dimension independent scaling between the wave packet dynamics
(measured through the $m_{2}$) and its degree of chaoticity (measured through
the $\Lambda$) between the 1D and 2D models. It will be interesting to check
if this relation also holds for three dimensional models as well.

Chaos is not the only ingredient needed for wave packet spreading, as the
other necessary component is the random oscillations of chaotic seeds within
the lattice’s excited part. The evolution of the DVD obtained during the
computation of the ftMLEs in the 2D DDNLS system also revealed the existence
of localized chaotic spots which randomly wander inside the excited part of
the lattice. Indeed, in the weak chaos case a timid spatial expansion of the
DVD was observed over time as its second moment grew as $m_{2}^{D}\sim
t^{0.12}$ and the participation number $P^{D}$ showed a tendency to saturate
to values similar to the 2D NM’s localization volume at later stages of the
evolution. The saturation of these observables was again much faster for the
strong chaos cases as a relation of the form $P^{D}\sim t^{0}$ could be fitted
to the data over a larger period of time compared to the weak chaos case. In
addition, we computed the area $A$ ($A\sim R^{2}$) visited by the DVD center
and found a monotonic increase, leading to values much larger than the NM’s
localization volume. Specifically, the values $A\approx 10^{2}$ and $A\approx
10^{4}$ are found for the weak and strong chaos regimes at the final time of
integration, suggesting that the DVD visited sites far from the lattice’s
center. Further, the temporal dependence of $A$ was well fitted by $A\sim
t^{\alpha_{A}}$, where the values $\alpha_{A}\approx 0.5$ and
$\alpha_{A}\approx 0.55$ were obtained for the weak and strong chaos regimes
respectively. Similarly to the 1D DDNLS model, the higher values of $A$ and
stronger variations of $\alpha_{A}$ observed in the strong chaos regime can be
related to the faster spreading of the wave packet in that regime, meaning
that the DVD has to cover larger areas of the excited part of the lattice than
in the weak chaos case. Consequently, the mechanisms of chaos responsible of
the wave packet spreading in the 1D and 2D DDNLS systems are similar. Most of
findings of this chapter have been reported in a recent research paper [Many
Manda et al., 2020].

It is worth noting that similar results to the ones presented in Chaps. 4 and
5 of this thesis have been obtained for the 1D and 2D DKG models in [Skokos et
al., 2013; Senyange et al., 2018; Many Manda et al., 2020], emphasizing the
generality of our findings.

Although our study contributed in comprehending several aspects of the
dynamical behavior of localized initial excitations in disordered systems,
several issues still need to be investigated. For instance the energy-norm
density spaces of the 1D and 2D DDNLS systems discussed in Sec. 2.4 revealed a
more complex partition between the wave packet spreading regions and the
thermal and non-thermal phases, also underlined in a recent work [Kati et al.,
2020]. Thus, the techniques used in this work to study the weak and strong
chaos spreading regimes could be used in further investigating the connection
between the wave packet spreading and the ergodic properties of these two
models. Another important aspect is the theoretical determination of the
observed power law exponents of the ftMLE seen in this work for both the 1D
and 2D DDNLS systems. Achieving this goal could help, for instance in
predicting the system’s dynamical characteristics for nonlinear strengths and
model’s dimensionalities far from the reach of current computer simulations.
In addition, we pinned the motion of the DVD as random and therefore the
evolution of the DVD center could be approximated by the spatial evolution of
a random walker evolving in a field of one or two spatial dimensions. This
means that, the total length/area covered represents only a single possible
measurement among many others of the characteristics of the DVD as a wide
theory exists around random walks. Furthermore, as the DVD generated during
the computation of the MLE shows the location of only the chaotic spot with
the largest instantaneous chaoticity, another important challenge is related
to the construction of a method to investigate most, if not all, of the
chaotic spots, as they are central in understanding the NM interactions.
Pursuing methods to extend all these ideas could potentially reveal more
subtle changes in the wave packet spreading dynamical behavior, allowing a
better in depth understanding of chaos properties in multidimensional systems.

## List of abbreviations

The following abbreviations are used throughout the content of this thesis.

  * •

1D – One-Dimensional/One-Dimension.

  * •

2D – Two-Dimensional/Two-Dimension.

  * •

AL – Anderson Localization.

  * •

BCH – Baker-Campbell-Hausdorff.

  * •

CPU – Central Processing Unit.

  * •

DDLS – Disordered Discrete Linear Schrödinger.

  * •

DDLSE – Disordered Discrete Linear Schrödinger Equation.

  * •

DDNLS – Disordered Discrete NonLinear Schrödinger.

  * •

DKG – Disordered Klein-Gordon.

  * •

DOP – Dormand and Prince.

  * •

DV – Deviation Vector.

  * •

DVD – Deviation Vector Distribution.

  * •

FPUT – Fermi-Pasta-Ulam-Tsingou.

  * •

ftMLE – finite-time Maximum Lyapunov Exponent.

  * •

KAM – Kolmogorov-Arnol’d-Moser.

  * •

MLE – Maximum Lyapunov Exponent.

  * •

NM – Normal Mode.

  * •

ODE – Ordinary Differential Equation.

  * •

RK – Runge-Kutta.

  * •

SI – Symplectic Integration/Symplectic Integrator.

##  Appendix A Frequency spectrum of the 1D DDLS model and the 1D DKG system
of harmonic oscillators

###  A.1 The 1D DDLS model

Here we analytically determine the characteristics of the bandwidth of the
frequency spectrum of the DDLS model. Its Hamiltonian takes the form [see Eq.
(2.1)]

$\mathcal{H}_{L}=\sum_{l}\epsilon_{l}\lvert\psi_{l}\rvert^{2}-\left(\psi_{l}^{\star}\psi_{l+1}+\psi_{l}\psi_{l+1}^{\star}\right),$
(A.1)

where $\psi_{l}\in\mathbb{C}$ is the wave function and $\epsilon_{l}$ the
disorder parameter at site $l$. The values of ${\epsilon_{l}}$ in (A.1) are
chosen from a uniform probability distribution $\mathcal{P}_{W}=1/W$ on the
interval $[-W/2,W/2]$, with $W$ denoting the system’s disorder strength. As
has already been mentioned many times in this thesis, plane waves describe the
solution of Hamiltonian (A.1). Consequently we write the wave function at site
$l$ with norm $s_{l}$ and phase $\theta_{l}$ (which modulates the value of the
energy $h_{l}$), as

$\psi_{l}=\sqrt{s_{l}}\exp\left(i\theta_{l}\right).$ (A.2)

In this thesis we have been working (without loss of generality) with
homogeneous initial wave packets so that $s=s_{l}$, for all the initially
excited sites. In this situation, the Hamiltonian (A.1) becomes

$\mathcal{H}_{L}=\sum_{l}\epsilon_{l}s-2s\cos\left(\theta_{l+1}-\theta_{l}\right).$
(A.3)

Since the cosine function is bounded, the value of $\mathcal{H}_{L}$ can also
be bounded as

$\sum_{l}\left(\epsilon_{l}-2\right)s\leq\sum_{l}\epsilon_{l}s-2s\cos\left(\theta_{l+1}-\theta_{l}\right)\leq\sum_{l}\left(\epsilon_{l}+2\right)s.$
(A.4)

In particular, the extreme case for the lower bound of (A.4) leads to a
minimum value when all sites have $\epsilon_{l}^{-}=-W/2$. In the same way,
the scenario corresponding to the maximum possible value of the upper bound in
(A.4) is all sites having $\epsilon_{l}^{+}=W/2$. Thus,

$\sum_{l}\left(-\frac{W}{2}-2\right)s\leq\sum_{l}\epsilon_{l}s-2s\cos\left(\theta_{l+1}-\theta_{l}\right)\leq\sum_{l}\left(\frac{W}{2}+2\right)s.$
(A.5)

Summing up each term of (A.5), we get

$\mathcal{H}_{L}^{-}=E^{-}\mathcal{S}_{1D}\leq\mathcal{H}_{L}\leq
E^{+}\mathcal{S}_{1D}=\mathcal{H}_{L}^{+},$ (A.6)

with

$E^{-}=-2-\frac{W}{2},\quad E^{+}=2+\frac{W}{2},$ (A.7)

and $\mathcal{S}_{1D}=\sum_{l}s$ being the total norm of the system which is a
conserved quantity. It is worth noticing that the Hamiltonian functions
$\mathcal{H}^{\pm}_{L}=E^{\pm}S_{1D}$ model independent harmonic oscillators,
each one of which has the same proper frequency $E^{\pm}$. The inequality in
(A.6) involves solely macroscopic quantities and therefore is valid for all
kinds of initial states of the system.

Applying the linear differentiation operator
$\frac{\partial}{\partial(i\psi_{l}^{\star})}$ in (A.6), we obtain

$\begin{split}\begin{pmatrix}E^{-}&0&\ldots&0&0\\\ 0&E^{-}&\ldots&0&0\\\
\vdots&\vdots&\ddots&\vdots&\vdots\\\ 0&0&\ldots&E^{-}&0\\\
0&0&\ldots&0&E^{-}\end{pmatrix}\begin{pmatrix}\psi_{1}\\\ \psi_{2}\\\
\vdots\\\ \psi_{N-1}\\\
\psi_{N}\end{pmatrix}&\leq\begin{pmatrix}\epsilon_{1}&-1&0&\ldots&0&0\\\
-1&\epsilon_{2}&-1&\ldots&0&0\\\ \vdots&\vdots&\vdots&\ddots&\vdots&\vdots\\\
0&0&0&\ldots&\epsilon_{N-1}&0\\\ 0&0&0&\ldots&-1&\epsilon_{N}\\\
\end{pmatrix}\begin{pmatrix}\psi_{1}\\\ \psi_{2}\\\ \vdots\\\ \psi_{N-1}\\\
\psi_{N}\end{pmatrix}\\\ &\leq\begin{pmatrix}E^{+}&0&\ldots&0&0\\\
0&E^{+}&\ldots&0&0\\\ \vdots&\vdots&\ddots&\vdots&\vdots\\\
0&0&\ldots&E^{+}&0\\\
0&0&\ldots&0&E^{+}\end{pmatrix}\begin{pmatrix}\psi_{1}\\\ \psi_{2}\\\
\vdots\\\ \psi_{N-1}\\\ \psi_{N}\end{pmatrix},\end{split}$ (A.8)

where the ‘$\leq$’ compares the system’s field vectors element by element.
Equation (A.8), can be reduced to a matrix form

$\bm{D}^{-}\bm{\Psi}\leq\tilde{\bm{A}}\bm{\Psi}\leq\bm{D}^{+}\bm{\Psi},$ (A.9)

with

$\bm{D}^{\pm}=\begin{pmatrix}E^{\pm}&0&\ldots&0&0\\\ 0&E^{\pm}&\ldots&0&0\\\
\vdots&\vdots&\ddots&\vdots&\vdots\\\ 0&0&\ldots&E^{\pm}&0\\\
0&0&\ldots&0&E^{\pm}\end{pmatrix},\quad\text{and}\quad\bm{\Psi}=\begin{pmatrix}\psi_{1}\\\
\psi_{2}\\\ \vdots\\\ \psi_{N-1}\\\ \psi_{N}\end{pmatrix}.$ (A.10)

In addition in (A.9) $\bm{\tilde{A}}$ is a tridiagonal matrix with entries

$\tilde{\bm{A}}=\begin{pmatrix}\epsilon_{1}&-1&0&\ldots&0&0\\\
-1&\epsilon_{2}&-1&\ldots&0&0\\\ \vdots&\vdots&\vdots&\ddots&\vdots&\vdots\\\
0&0&0&\ldots&\epsilon_{N-1}&-1\\\ 0&0&0&\ldots&-1&\epsilon_{N}\\\
\end{pmatrix}.$ (A.11)

The matrix $\tilde{\bm{A}}$ is diagonalizable [Gershgorin, 1931; Wolkowicz and
Styan, 1980; Smith, 1985; Losonczi, 1992; Yueh, 2005; Da Fonseca and
Kowalenko, 2019]. Its corresponding diagonal matrix $\bm{D}$ has as diagonal
elements the eigenvalues $E_{\nu}$, $\nu=1,2,\ldots,N$ of $\tilde{\bm{A}}$.
Defining $\bm{P}$ a matrix having as columns the normalized eigenvectors
related to the eigenvalues of $\bm{D}$, we get

$\bm{D}^{-}\bm{\Psi}\leq\bm{P}\bm{D}\bm{P}^{-1}\bm{\Psi}\leq\bm{D}^{+}\bm{\Psi}.$
(A.12)

Multiplying from the left side (A.12) by $\bm{P}^{-1}$ leads to

$\bm{D}^{-}\bm{P}^{-1}\bm{\Psi}\leq\bm{D}\bm{P}^{-1}\bm{\Psi}\leq\bm{D}^{+}\bm{P}{-1}\bm{\Psi},$
(A.13)

where we used $\bm{P}^{-1}\bm{P}=\bm{P}\bm{P}^{-1}=\bm{I}_{N}$ with
$\bm{I}_{N}$ being the $N\times N$ identity matrix. Setting the state vector
$\bm{U}=\bm{P}^{-1}\bm{\Psi}$, we therefore end up with a set of uncoupled
oscillators

$\bm{D}^{-}\bm{U}\leq\bm{D}\bm{U}\leq\bm{D}^{+}\bm{U},$ (A.14)

whose frequency spectra can be compared wave number by wave number in the
matrices of (A.14). What comes out from this comparison is the inequality

$E^{-}\leq E_{\nu}\leq E^{+},$ (A.15)

with $E_{\nu}$, $\nu=1,2,\ldots,N$ representing the eigenvalues of matrix
$\tilde{\bm{A}}$. Therefore, the width of the frequency spectrum is

$\Delta=E^{+}-E^{-}=W+4,$ (A.16)

since

$E_{\nu}\in\left[-2-\frac{W}{2},2+\frac{W}{2}\right].$ (A.17)

###  A.2 The 1D DKG system of harmonic oscillators

We can express the linear characteristics of the DKG model in terms of the
DDLS model. The Hamiltonian function of the linear DKG model of harmonic
oscillators reads [see e.g. Skokos et al., 2009]

$\mathcal{H}_{LK}=\sum_{l}\frac{p_{l}^{2}}{2}+\frac{\tilde{\epsilon}_{l}}{2}q_{l}^{2}+\frac{1}{2W}\left(q_{l+1}-q_{l}\right)^{2},$
(A.18)

and its related equations of motion can be written as

$\ddot{q}_{l}=\tilde{\epsilon}_{l}q_{l}+\frac{1}{W}\left(2q_{l}-q_{l-1}-q_{l+1}\right),$
(A.19)

where $q_{l}$ and $p_{l}$ are respectively the generalized position and
conjugate momentum and $\tilde{\epsilon}_{l}$ denotes disorder parameter at
site $l$. The values of $\epsilon_{l}$ are chosen on the interval $[1/2,3/2]$
following a random uniform distribution. Here $W$ denotes the disorder
strength because it controls the size of the normal mode’s localization volume
$\overline{V}$ [Skokos et al., 2009; Senyange et al., 2020].

The solutions of (A.19) are of the form

$q_{l}=A_{l}\exp\left(i\omega t\right),$ (A.20)

with $A_{l}$ being the amplitude at site $l$ and $\omega$ the frequency of
oscillation. By substituting (A.20) to (A.19) we obtain the following
eigenvalue problem

$\left(\omega^{2}W-W-2\right)A_{l}=W\left(\tilde{\epsilon}_{l}-1\right)A_{l}-A_{l-1}-A_{l+1}.$
(A.21)

The normal modes of (A.21)
$\bm{\Psi}_{\nu}=(A_{\nu,1},A_{\nu,2},\ldots,A_{\nu,N})$, match the ones of
the DDLS system (2.17) given in (2.5) by setting

$\displaystyle E$ $\displaystyle=\omega^{2}W-W-2,$ (A.22)
$\displaystyle\epsilon_{l}$
$\displaystyle=W\left(\tilde{\epsilon}_{l}-1\right).$ (A.23)

Consequently Eq. (A.22) allows us to determine the limits of the frequency
spectrum for the DKG system of linear oscillators (A.18) based on the results
obtained for the DDLS model in Sec. A.1. In particular, using (A.7) the lower
squared frequency $\omega_{-}^{2}$ of the bandwidth is

$E^{-}=\omega_{-}^{2}W-W-2\quad\Rightarrow\quad-2-\frac{W}{2}=\omega_{-}^{2}W-W-2\quad\Rightarrow\quad\omega_{-}^{2}=\frac{1}{2},$
(A.24)

and the upper bound $\omega_{+}^{2}$ is given by

$E^{+}=\omega_{+}^{2}W-W-2\quad\Rightarrow\quad
2+\frac{W}{2}=\omega_{+}^{2}W-W-2\quad\Rightarrow\quad\omega_{+}^{2}=\frac{3}{2}+\frac{4}{W}.$
(A.25)

Thus the size of the squared frequency spectrum of the DKG system of linear
oscillators is

$\Delta_{K}=\frac{3}{2}+\frac{4}{W}-\frac{1}{2}\quad\Rightarrow\quad\Delta_{K}=1+\frac{4}{W}.$
(A.26)

A more general derivation of this relation can be found in [Senyange et al.,
2020].

##  Appendix B Equations of motion of the 1D DDNLS system in the NM space

We recall the equations of motion of the 1D DDNLS system (2.18) are

$i\dot{\psi}_{l}=\epsilon_{l}\psi_{l}+\beta\lvert\psi_{l}\rvert^{2}\psi_{l}-\left(\psi_{l-1}+\psi_{l+1}\right).$
(B.1)

Writing the solution of Eq. (B.1) in the NM space of the system as a weighted
sum over its linear NMs
$\bm{\Psi}_{\nu}\propto\overline{\zeta}^{-1/2}\exp\left[-(l-\overline{l}_{\nu})/\overline{\zeta}\right]$

$\psi_{l}=\sum_{\nu}\phi_{\nu}\bm{\Psi}_{\nu}=\sum_{\nu}\phi_{\nu}A_{\nu,l},$
(B.2)

where $\bm{\Psi}_{\nu}=\left(A_{\nu,1},A_{\nu,2},\ldots,A_{\nu,N}\right)$ is
the NM and $\phi_{\nu}$ the time dependent amplitude at mode $\nu$ with
$\nu=1,2,\ldots,N$ ($N$ being the lattice size), and substituting this
expression in Eq. (B.1), we obtain

$i\sum_{\nu}A_{\nu,l}\dot{\phi}_{\nu}=\sum_{\nu}\epsilon_{l}A_{\nu,l}\phi_{\nu}-\sum_{\nu}\phi_{\nu}A_{\nu,l-1}-\sum_{\nu}\phi_{\nu}A_{\nu,l+1}+\beta\sum_{\nu_{1}\nu_{2}\nu_{3}}A_{\nu_{1},l}^{\star}A_{\nu_{2},l}A_{\nu_{3},l}\phi_{\nu_{1}}^{\star}\phi_{\nu_{2}}\phi_{\nu_{3}}.$
(B.3)

Taking into account Eq. (2.5) i.e. $EA_{l}=\epsilon_{l}A_{l}-A_{l-1}-A_{l+1}$,
we rewrite (B.3) as

$i\sum_{\nu}A_{\nu,l}\dot{\phi}_{\nu}=\sum_{\nu}E_{\nu}A_{\nu,l}\phi_{\nu}+\beta\sum_{\nu_{1}\nu_{2}\nu_{3}}A_{\nu_{1},l}^{\star}A_{\nu_{2},l}A_{\nu_{3},l}\phi_{\nu_{1}}^{\star}\phi_{\nu_{2}}\phi_{\nu_{3}},$
(B.4)

where $E_{\nu}$ is the frequency of the normal mode $\nu$. We then multiply
Eq. (B.4) by $\bm{\Psi}_{\nu^{\prime}}^{\star}$, and get

$i\sum_{\nu}\sum_{l}A_{\nu^{\prime},l}^{\star}A_{\nu,l}\dot{\phi}_{\nu}=\sum_{\nu}\sum_{l}A_{\nu^{\prime},l}^{\star}E_{\nu}A_{\nu,l}\phi_{\nu}+\beta\sum_{\nu_{1}\nu_{2}\nu_{3}}\sum_{l}A_{\nu^{\prime},l}^{\star}A_{\nu_{1},l}^{\star}A_{\nu_{2},l}A_{\nu_{3},l}\phi_{\nu_{1}}^{\star}\phi_{\nu_{2}}\phi_{\nu_{3}}.$
(B.5)

At this point it is useful to emphasize few things:

  1. (a) The orthogonality of the NMs leads to

$\bm{\Psi}_{\nu}^{\star}\bm{\Psi}_{\nu}=\sum_{l}A_{\nu^{\prime},l}^{\star}A_{\nu,l}=\delta_{\nu^{\prime},\nu}.$
(B.6)

  2. (b) In the considered set up, initial conditions are chosen so that the total norm of any initial wave packet is set to unity. This leads to the condition

$\sum_{l}\lvert\psi_{l}\rvert^{2}=\sum_{l}\sum_{\nu^{\prime}}\sum_{\nu}\phi_{\nu^{\prime}}^{\star}\phi_{\nu}\bm{\Psi}_{\nu^{\prime}}^{\star}\bm{\Psi}_{\nu}=\sum_{\nu}\sum_{l}A_{\nu,l}^{\star}A_{\nu,l}\phi_{\nu}^{\star}\phi_{\nu}=\sum_{\nu}\lvert\phi_{\nu}\rvert^{2}=1.$
(B.7)

Substituting Eq. (B.6) into Eq. (B.5) leads to

$i\dot{\phi}_{\nu^{\prime}}=E_{\nu^{\prime}}\phi_{\nu^{\prime}}+\beta\sum_{\nu_{1}\nu_{2}\nu_{3}}I_{\nu^{\prime},\nu_{1},\nu_{2},\nu_{3}}\phi_{\nu_{1}}^{\star}\phi_{\nu_{2}}\phi_{\nu_{3}},$
(B.8)

where

$I_{\nu^{\prime},\nu_{1},\nu_{2},\nu_{3}}=\sum_{l}A_{\nu^{\prime},l}^{\star}A_{\nu_{1},l}^{\star}A_{\nu_{2},l}A_{\nu_{3},l},$
(B.9)

is the so-called overlap-integral. Note that similar calculation can be
performed for the DKG system [Laptyeva, 2013].

##  Appendix C Statistical physics analytical computations

###  C.1 Derivation of Eq. (2.81)

We give here the detailed calculation of the integral

$\tilde{y}_{\beta}\left(\tilde{\nu},\mu_{l}\right)=\int_{0}^{\infty}dxe^{-\tilde{\nu}\mu_{l}x}e^{-\tilde{\nu}\frac{\beta}{2}x^{2}},$
(C.1)

appearing in Eq. (2.80). In (C.1), $\tilde{\nu}=1/T$ with $T$ being the
temperature of the system, $\beta$ is the nonlinear coefficient of the 1D
DDNLS model (2.17) and $\mu_{l}=\mu+\epsilon_{l}$ represents a kind of local
chemical potential at site $l$ of the lattice with $\epsilon_{l}$ being the
disorder parameter of $\mathcal{H}_{1D}$ (2.17). From the Taylor expansion of
the exponential function $e^{x}=\sum_{n=0}^{\infty}x^{n}/n!$ we get

$e^{-\tilde{\beta}\frac{x^{2}}{2}}=1-\frac{\tilde{\beta}x^{2}}{2}+\frac{\tilde{\beta}^{2}x^{4}}{8}+\ldots,$
(C.2)

with $\tilde{\beta}=\beta\tilde{\nu}$. Substituting (C.2) into (C.1) leads to

$\tilde{y}_{\beta}\left(\tilde{\nu},\mu_{l}\right)=\int_{0}^{\infty}dxe^{-\tilde{\nu}\mu_{l}x}\left[1-\frac{\tilde{\beta}x^{2}}{2}+\frac{\tilde{\beta}^{2}x^{4}}{8}+\ldots\right].$
(C.3)

Consequently, we can write

$\tilde{y}_{\beta}\left(\tilde{\nu},\mu_{l}\right)\approx\int_{0}^{\infty}dxe^{-\tilde{\nu}\mu_{l}x}-\int_{0}^{\infty}dx\frac{\tilde{\beta}x^{2}}{2}e^{-\tilde{\nu}\mu_{l}x}+\ldots.$
(C.4)

The first integral in (C.4) can be computed straightforwardly. For the
following terms, we use the fact that [Abramowitz and Stegun, 1964; Johansson
and Rasmussen, 2004]

$\frac{\Gamma(n+1)}{a^{n+1}}=\int_{0}^{\infty}dxx^{n}e^{-ax},$ (C.5)

where $a\in\mathbb{R}$ is a constant and $\Gamma$ is the usual $\Gamma$
function defined such that $\Gamma(n+1)=n\Gamma(n)$ if $n\in\mathbb{N}$.
Therefore

$\tilde{y}_{\beta}\left(\tilde{\nu},\mu_{l}\right)=\frac{1}{\tilde{\nu}\mu_{l}}-\frac{\tilde{\beta}}{2}\frac{\Gamma(3)}{(\tilde{\nu}\mu)^{3}}+\frac{\tilde{\beta}^{2}}{8}\frac{\Gamma(5)}{(\tilde{\nu}\mu_{l})^{5}}+\ldots.$
(C.6)

In the limit of $\tilde{\nu}\rightarrow 0$, because $\beta$ is finite
$\tilde{\beta}\nu\rightarrow 0$ we can neglect higher order terms in
$\tilde{\beta}$ and obtain

$\tilde{y}_{\beta}\left(\tilde{\nu},\mu_{l}\right)\approx\frac{1}{\tilde{\nu}\mu_{l}}\left(1-\frac{\beta\tilde{\nu}}{(\tilde{\nu}\mu_{l})^{2}}\right).$
(C.7)

###  C.2 Derivation of Eq. (2.84)

The total energy $\mathcal{H}_{1D}$ and the total norm $\mathcal{S}_{1D}$ of
the 1D DDNLS system (2.17) can be calculated from the grand partition function
$\mathcal{Z}$ (2.82) [Landau and Lifshitz, 1980; Tong, 2012]. Here we evaluate
both the average norm $s=\mathcal{S}_{1D}/N$ and energy $h=\mathcal{H}_{1D}/N$
per site in terms of the system’s parameters at equilibrium. We recall the
partition function $\ln\mathcal{Z}$ (2.82) containing all the information on
the macroscopic states of the system was presented to be

$\ln\mathcal{Z}\approx N\ln
2\pi-\sum_{l=1}^{N}\left[\ln\left(\tilde{\nu}\mu_{l}\right)+\frac{\tilde{\nu}\beta}{\left(\tilde{\nu}\mu_{l}\right)^{2}}\right],$
(C.8)

where $\tilde{\nu}=1/T$, with $T$ being the system temperature, $\beta$ the
nonlinear coefficient in the expression of $\mathcal{H}_{1D}$ (2.17) and
$\mu_{l}=\mu+\epsilon_{l}$, where $\mu$ is a chemical potential and
$\epsilon_{l}$ the disorder parameter in (2.17). The total norm
$\mathcal{S}_{1D}$ is given by

$\begin{split}\mathcal{S}_{1D}&=-\frac{1}{\tilde{\nu}}\frac{\partial}{\partial\mu_{l}}\ln\mathcal{Z},\\\
&=-\frac{1}{\tilde{\nu}}\left\\{-\sum_{l}\left[\frac{\tilde{\nu}}{\tilde{\nu}\mu_{l}}-\frac{2\tilde{\nu}^{2}\beta}{(\tilde{\nu}\mu_{l})^{3}}\right]\right\\},\\\
&=\sum_{l}\frac{1}{\tilde{\nu}\mu_{l}}-\sum_{l}\frac{2\beta}{(\tilde{\nu}\mu_{l})^{2}\mu_{l}}.\end{split}$
(C.9)

We note that in the infinite temperature limit $\tilde{\nu}\rightarrow 0$,
while the product $\tilde{\nu}\mu$ remains constant by definition [Landau and
Lifshitz, 1980; Rasmussen et al., 2000; Johansson and Rasmussen, 2004; Tong,
2012; Flach, 2016]. Therefore the chemical potential $\mu\sim 1/\tilde{\nu}$.
In addition since $\epsilon_{l}$ is finite, and $\tilde{\nu}$ is large, we
also assume that $\mu_{l}\sim 1/\tilde{\nu}$. Consequently the second term in
(C.9) becomes

$\sum_{l}\frac{2\beta}{(\tilde{\nu}\mu_{l})^{2}\mu_{l}}\sim\tilde{\nu}\sum_{l}\frac{2\beta}{(\tilde{\nu}\mu_{l})^{2}}\rightarrow
0,$ (C.10)

for $\tilde{\nu}\rightarrow 0$. Thus, (C.10) can be used to provide the
approximation

$\mathcal{S}_{1D}\approx\sum_{l}\frac{1}{\tilde{\nu}\mu_{l}}=\sum_{l}\frac{1}{\tilde{\nu}(\mu+\epsilon_{l})}=\sum_{l}\frac{1}{\tilde{\nu}\mu+\tilde{\nu}\epsilon_{l}}\sim\frac{N}{\tilde{\nu}\mu},$
(C.11)

as $\tilde{\nu}\epsilon_{l}\rightarrow 0$ when $\tilde{\nu}\rightarrow 0$.
From (C.11), we conclude that the expression of the average norm per site is

$s=\frac{\mathcal{S}_{1D}}{N}\approx\frac{1}{\tilde{\nu}\mu}.$ (C.12)

Let us now estimate the average energy per site. The total energy
$\mathcal{H}_{1D}$ is given by

$\mathcal{H}_{1D}=\left(\frac{\mu}{\tilde{\nu}}\frac{\partial}{\partial\mu_{l}}-\frac{\partial}{\partial\tilde{\nu}}\right)\ln\mathcal{Z}.$
(C.13)

Then (C.13), using the expression of $\mathcal{Z}$ from (C.8) gives

$\begin{split}\mathcal{H}_{1D}&=\frac{\mu}{\tilde{\nu}}\left(-\sum_{l}\left[\frac{\tilde{\nu}}{\tilde{\nu}\mu_{l}}-\frac{2\tilde{\nu}^{2}\beta}{(\tilde{\nu}\mu_{l})^{3}}\right]\right)+\sum_{l}\left[\frac{\mu_{l}}{\tilde{\nu}\mu_{l}}-\frac{\beta}{(\tilde{\nu}\mu_{l})^{2}}\right],\\\
&=-\sum_{l}\frac{\mu}{\tilde{\nu}\mu_{l}}+\sum_{l}\frac{2\beta\mu}{(\tilde{\nu}\mu_{l})^{2}\mu_{l}}+\sum_{l}\frac{\mu_{l}}{\tilde{\nu}\mu_{l}}-\sum_{l}\frac{\beta}{(\tilde{\nu}\mu_{l})^{2}},\\\
&=\sum_{l}\frac{\epsilon_{l}}{\tilde{\nu}\mu_{l}}-\sum_{l}\frac{\beta\mu_{l}}{(\tilde{\nu}\mu_{l})^{2}\mu_{l}}+\sum_{l}\frac{2\beta\mu}{(\tilde{\nu}\mu_{l})^{2}\mu_{l}},\\\
&=\sum_{l}\frac{\epsilon_{l}}{\tilde{\nu}\mu_{l}}+\sum_{l}\frac{\beta}{(\tilde{\nu}\mu_{l})^{2}}-\sum_{l}\frac{\beta\epsilon_{l}}{(\tilde{\nu}\mu_{l})^{2}\mu_{l}}.\end{split}$
(C.14)

The last term in the final expression of $\mathcal{H}_{1D}$ in (C.14) is
similar to the one in the final expression of $\mathcal{S}_{1D}$ in (C.9),
which converges to zero in the limit of $\tilde{\nu}\rightarrow 0$. Thus

$\mathcal{H}_{1D}\approx\sum_{l}\frac{\epsilon_{l}}{\tilde{\nu}\mu_{l}}+\sum_{l}\frac{\beta}{(\tilde{\nu}\mu)^{2}}.$
(C.15)

Let us look a bit closer each term of $\mathcal{H}_{1D}$ (C.15). The first one
can be written as

$\begin{split}\sum_{l}\frac{\epsilon_{l}}{\tilde{\nu}\mu_{l}}&=\sum_{l}\frac{\epsilon_{l}}{\tilde{\nu}(\mu+\epsilon_{l})}=\frac{1}{\tilde{\nu}\mu}\sum_{l}\left(\frac{\epsilon_{l}}{1+\frac{\epsilon_{l}}{\mu}}\right).\end{split}$
(C.16)

Then, applying the Taylor expansion

$\frac{\epsilon_{l}}{1+\frac{\epsilon_{l}}{\mu}}=\epsilon_{l}\left(1-\frac{\epsilon_{l}}{\mu}+\frac{\epsilon_{l}^{2}}{\mu^{2}}+\ldots\right),$
(C.17)

and assuming that $\epsilon_{l}/\mu\ll 1$ as $\mu\sim 1/\tilde{\nu}\gg 1$, for
$\tilde{\nu}\rightarrow 0$ in (C.16), we obtain

$\sum_{l}\frac{\epsilon_{l}}{\tilde{\nu}(\mu+\epsilon_{l})}\approx\frac{1}{\tilde{\nu}\mu}\left(\sum_{l}\epsilon_{l}-\frac{1}{\mu}\sum_{l}\epsilon_{l}^{2}\right),$
(C.18)

where we considered only the first two leading terms in (C.17) [see also
[Flach, 2016]]. In addition, in order for the value of $\sigma_{\epsilon}$
(see (2.63)) to be finite i.e. $\sup\sigma_{\epsilon}<\infty$, we must have

$\sum_{l}\epsilon_{l}^{2}\sim
N\quad\Rightarrow\quad\sum_{l}\epsilon_{l}\sim\sqrt{N}.$ (C.19)

With these two approximations of (C.19), (C.18) gives

$\sum_{l}\frac{\epsilon_{l}}{\tilde{\nu}(\mu+\epsilon_{l})}\sim\frac{\sqrt{N}}{\tilde{\nu}\mu}-\frac{N\tilde{\nu}}{\left(\tilde{\nu}\mu\right)^{2}}.$
(C.20)

The second term in the right hand side of (C.20) can be neglected when
$\tilde{\nu}\rightarrow 0$. In addition the variation of
$\sqrt{N}/(\tilde{\nu}\mu)$ is too slow to contribute to the final energy-norm
relation at infinite temperature [Flach, 2016]. Finally, only the second term
in the right hand side of (C.15) remains, leading to

$h=\frac{\mathcal{H}_{1D}}{N}\approx\frac{\beta}{(\tilde{\nu}\mu)^{2}}\quad\Rightarrow\quad
h\approx\beta\left(\frac{1}{\tilde{\nu}\mu}\right)^{2}\quad\Rightarrow\quad
h\approx\beta s^{2}.$ (C.21)

##  Appendix D Symplectic integration of the equations of motion and
variational equations of the 1D and 2D DDNLS models

We present here the expression of the various operators needed for the
implementation of SIs for the integration of the equations of motion and the
variational equations for the 1D and 2D DDNLS systems. We note that the
Hamiltonian of the 1D DDNLS model (2.17)

$\mathcal{H}_{1D}=\sum_{l=1}^{N}\epsilon_{l}\lvert\psi_{l}\rvert^{2}+\frac{\beta}{2}\lvert\psi_{l}\rvert^{4}-\left(\psi_{l+1}\psi_{l}^{\star}+\psi_{l+1}^{\star}\psi_{l}\right),$
(D.1)

becomes

$\mathcal{H}_{1D}=\sum_{l}\left\\{\frac{\epsilon_{l}}{2}\left(q_{l}^{2}+p_{l}^{2}\right)+\frac{\beta}{8}\left(q_{l}^{2}+p_{l}^{2}\right)^{4}-p_{l+1}p_{l}-q_{l+1}q_{l}\right\\},$
(D.2)

with $q_{l}$ and $p_{l}$ being real variables, through the implementation of
the transformation (3.90)

$\psi_{l}=\frac{1}{\sqrt{2}}\left(q_{l}+ip_{l}\right).$ (D.3)

In a similar manner, the 2D DDNLS Hamiltonian function (2.20)

$\mathcal{H}_{2D}=\sum_{l=1,m=1}^{N,M}\epsilon_{l,m}\lvert\psi_{l,m}\rvert^{2}+\frac{\beta}{2}\lvert\psi_{l,m}\rvert^{4}-\left(\psi_{l+1,m}\psi_{l,m}^{\star}+\psi_{l,m+1}\psi_{l,m}^{\star}+\psi_{l+1,m}^{\star}\psi_{l,m}+\psi_{l,m+1}^{\star}\psi_{l,m}\right),$
(D.4)

becomes

$\mathcal{H}_{2D}=\sum_{l,m}\left\\{\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}-q_{l,m+1}q_{l,m}-q_{l+1,m}q_{l,m}-p_{l,m+1}p_{l,m}-p_{l+1,m}p_{l,m}\right\\},$
(D.5)

after the implementation of the transformation

$\psi_{l,m}=\frac{1}{\sqrt{2}}\left(q_{l,m}+ip_{l,m}\right).$ (D.6)

We note that for both Hamiltonians (D.2) and (D.5) we use fixed boundary
conditions.

###  D.1 Two-part split symplectic integration of the 1D DDNLS model

As was mentioned in Sec. 3.8.1, the Hamiltonian function of the 1D DDNLS
system can be split into two integrable parts as

$\mathcal{H}_{1D}=A+B,$ (D.7)

with

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

whose solutions can be obtained in closed forms.

####  D.1.1 Analytic expression for the operator $e^{\tau{\bf L}_{A_{V}}}$

Let us first consider the $A$ part of the decomposition (D.8)

$A=\sum_{l}\frac{\epsilon_{l}}{2}(q_{l}^{2}+p_{l}^{2})+\frac{\beta}{8}(q_{l}^{2}+p_{l}^{2})^{2}.$
(D.9)

Its equations of motion are

$\displaystyle\frac{d\bm{z}}{dt}=\bm{L}_{A}\bm{z}\colon\begin{dcases}\dot{q_{l}}&=p_{l}\left(\epsilon_{l}+\beta\frac{q_{l}^{2}+p_{l}^{2}}{2}\right)\\\
\dot{p_{l}}&=-q_{l}\left(\epsilon_{l}+\beta\frac{q_{l}^{2}+p_{l}^{2}}{2}\right)\end{dcases},$
(D.10)

which conserve the value of $A$ (D.9).

Let us express (D.9) and (D.10) in action-angles variables $J$, $\Theta$
through the canonical transformation [Herbert, 1980; Goldstein et al., 2002;
Tong, 2004]

$\Omega_{A}:\left(q_{l},p_{l}\right)\mapsto\left(J_{l},\Theta_{l}\right)\colon\left\\{\begin{array}[]{ll}J_{l}=\frac{1}{2}(q_{l}^{2}+p_{l}^{2})\\\
\Theta_{l}=\arctan\frac{q_{l}}{p_{l}}\\\ \end{array}\right..$ (D.11)

It is easily seen that the inverse of transformation $\Omega_{A}$ (D.11) is
given by

$\Omega_{A}^{-1}:\left(J_{l},\Theta_{l}\right)\mapsto\left(q_{l},p_{l}\right)\colon\left\\{\begin{array}[]{ll}q_{l}=\sqrt{2J_{l}}\sin\Theta_{l}\\\
p_{l}=\sqrt{2J_{l}}\cos\Theta_{l}\end{array}\right..$ (D.12)

In (D.11) and (D.12), the variable $\Theta_{l}$ plays the role of a
generalized position and $J_{l}$ of the generalized momentum. The application
of (D.12) on Hamiltonian $A$ (D.9) brings it to the form

$A=\sum_{l}\epsilon_{l}J_{l}+\frac{\beta}{2}J_{l}^{2}.$ (D.13)

The corresponding equations of motion generated through the Hamilton equations
of motion

$\dot{\Theta}_{l}=\frac{\partial A}{\partial
J_{l}},\quad\dot{J}_{l}=-\frac{\partial A}{\partial\Theta_{l}},$ (D.14)

are

$\displaystyle\begin{dcases}\dot{\Theta}_{l}&=\epsilon_{l}+\beta J_{l}\\\
\dot{J}_{l}&=0\end{dcases}.$ (D.15)

From the second equation of (D.15), we see that the actions $J_{l}$,
$l=1,2,\ldots,N$ are integrals of motion and their values are conserved along
the system’s evolution. Given an initial condition $(\Theta_{l}(0),J_{l}(0))$,
the solution of (D.15) is given by

$\displaystyle\begin{dcases}\Theta_{l}(t)&=(\epsilon_{l}+\beta
J_{l}(0))t+\Theta_{l}(0)\\\ J_{l}(t)&=J_{l}(0)\end{dcases}.$ (D.16)

Using the transformation $\Omega_{A}^{-1}$ (D.12), we rewrite the solution
(D.16) in terms of the $q_{l}$ and $p_{l}$ variables as

$\displaystyle\begin{dcases}q_{l}(t)&=\sqrt{2J_{l}(0)}\sin\left[\left(\epsilon_{l}+\beta
J_{l}(0)\right)t+\Theta_{l}(0)\right]\\\
p_{l}(t)&=\sqrt{2J_{l}(0)}\cos\left[\left(\epsilon_{l}+\beta
J_{l}(0)\right)t+\Theta_{l}(0)\right]\end{dcases},$ (D.17)

with

$J_{l}(0)=\frac{1}{2}\left(q_{l}(0)^{2}+p_{l}(0)^{2}\right),\quad\Theta_{l}(0)=\arctan\frac{q_{l}(0)}{p_{l}(0)}.$
(D.18)

Let us now discuss some particular cases related to (D.17).

  1. Case A: $p_{l}(0)\neq 0$. No problems arise in this case and the substitution of (D.18) onto Eq. (D.17) leads to

$\displaystyle\begin{dcases}q_{l}(t)&=\sqrt{q_{l}(0)^{2}+p_{l}(0)^{2}}\left[\sin\left[(\epsilon_{l}+\beta
J_{l}(0))t\right]\cos\Theta_{l}(0)+\cos\left[(\epsilon_{l}+\beta
J_{l}(0))t\right]\sin\Theta_{l}(0)\right]\\\
p_{l}(t)&=\sqrt{q_{l}(0)^{2}+p_{l}(0)^{2}}\left[\cos\left[(\epsilon_{l}+\beta
J_{l}(0))t\right]\cos\Theta_{l}(0)-\sin\left[(\epsilon_{l}+\beta
J_{l}(0))t\right]\sin\Theta_{l}(0)\right]\end{dcases},$ (D.19)

which can further be rewritten as

$\displaystyle\begin{dcases}q_{l}(t)&=\sqrt{q_{l}(0)^{2}+p_{l}(0)^{2}}\left[\sin\left[(\epsilon_{l}+\beta
J_{l}(0))t\right]\frac{1}{\sqrt{1+\frac{q_{l}(0)^{2}}{p_{l}(0)^{2}}}}+\cos\left[(\epsilon_{l}+\beta
J_{l}(0))t\right]\frac{\frac{q_{l}(0)}{p_{l}(0)}}{\sqrt{1+\frac{q_{l}(0)^{2}}{p_{l}(0)^{2}}}}\right]\\\
p_{l}(t)&=\sqrt{q_{l}(0)^{2}+p_{l}(0)^{2}}\left[\cos\left[(\epsilon_{l}+\beta
J_{l}(0))t\right]\frac{1}{\sqrt{1+\frac{q_{l}(0)^{2}}{p_{l}(0)^{2}}}}-\sin\left[(\epsilon_{l}+\beta
J_{l}(0))t\right]\frac{\frac{q_{l}(0)}{p_{l}(0)}}{\sqrt{1+\frac{q_{l}(0)^{2}}{p_{l}(0)^{2}}}}\right]\end{dcases}.$
(D.20)

We note that in deriving Eq. (D.20), we used the following trigonometric
identities $\sin(a+b)=\sin a\cos b+\cos a\sin b$, $\cos(a+b)=\cos a\cos b-\sin
a\sin b$, $\sin a=\tan a/\sqrt{1+\tan^{2}a}$ and $\cos
a=1/\sqrt{1+\tan^{2}a}$. In this case the solution of the Hamilton equations
of motion of $A$ (D.9) is given by [Skokos et al., 2014; Gerlach et al., 2016;
Danieli et al., 2019a]

$\displaystyle\begin{dcases}q_{l}(t)&=q_{l}(0)\cos\alpha_{l}t+p_{l}(0)\sin\alpha_{l}t\\\
p_{l}(t)&=p_{l}(0)\cos\alpha_{l}t-q_{l}(0)\sin\alpha_{l}t\end{dcases},$ (D.21)

with

$\alpha_{l}=\epsilon_{l}+\beta
J_{l}(0)=\epsilon_{l}+\beta\frac{q_{l}(0)^{2}+p_{l}(0)^{2}}{2},$

being constant coefficients.

  2. Case B: $p_{l}(0)=0$. In this case, we have to consider two different options: $q_{l}(0)\neq 0$ and $q_{l}(0)=0$.

    1. Case B$1$: $q_{l}(0)=0$. The equations of motion (D.10) become

$\displaystyle\begin{dcases}\dot{q}_{l}=0\\\ \dot{p}_{l}=0\end{dcases},$
(D.22)

and can be trivially solved to get $q_{l}(t)=0$ and $p_{l}(t)=0$.

    2. Case B$2$: $q_{l}(0)\neq 0$. In this case, we can have two possible values for $\Theta_{l}(0)$, namely

$\left\\{\begin{array}[]{lcr}\Theta_{l}(0)=&\frac{\pi}{2}&\mbox{ if
$q_{l}(0)>0$}\\\ \Theta_{l}(0)=&-\frac{\pi}{2}&\mbox{ if
$q_{l}(0)<0$}\end{array}\right..$ (D.23)

For the first case in (D.23), the solution (D.17) is written as

$\displaystyle\begin{dcases}q_{l}(t)&=\sqrt{2J_{l}(0)}\sin\left[(\epsilon_{l}+\beta
J_{l}(0))t+\frac{\pi}{2}\right]\\\
p_{l}(t)&=\sqrt{2J_{l}(0)}\cos\left[(\epsilon_{l}+\beta
J_{l}(0))t+\frac{\pi}{2}\right]\end{dcases},$ (D.24)

or simply as

$\displaystyle\begin{dcases}q_{l}(t)&=q_{l}(0)\cos\left[(\epsilon_{l}+\beta
J_{l}(0))t\right]\\\ p_{l}(t)&=-q_{l}(0)\sin\left[(\epsilon_{l}+\beta
J_{l}(0))t\right]\end{dcases}.$ (D.25)

We note that going from (D.24) to (D.25), we used the trigonometric identities
$\sin(a+\pi/2)=\cos a$ and $\cos(a+\pi/2)=-\sin a$, as well as the facts that
$J_{l}(0)=\frac{1}{2}q_{l}(0)^{2}$ and $\lvert
q_{l}(0)\rvert=q_{l}(0)\leavevmode\nobreak\ (q_{l}(0)>0)$. Following a similar
process as in (D.24) and (D.25) in the second case of (D.23), and using the
identities $\sin(a-\pi/2)=-\cos a$ and $\cos(a-\pi/2)=\sin a$, along with
$\lvert q_{l}(0)\rvert=-q_{l}(0)$ on (D.24), we finally obtain the solution

$\displaystyle\begin{dcases}q_{l}(t)&=q_{l}(0)\cos\left[(\epsilon_{l}+\beta
J_{l}(0))t\right]\\\ p_{l}(t)&=-q_{l}(0)\sin\left[(\epsilon_{l}+\beta
J_{l}(0))t\right]\end{dcases}.$ (D.26)

The solutions (D.25) and (D.26) are identical. Consequently, no distinction
between the cases B$1$ and B$2$ needs to be implemented in our computer
program.

Let us now discuss the solution of the variational equations of the
Hamiltonian $A$ (D.9) which have the form

$\frac{d\delta\bm{z}}{dt}=\bm{L}_{A_{V}}\delta\bm{z}\colon\begin{dcases}\dot{\delta
q_{l}}&=\beta q_{l}p_{l}\delta
q_{l}+\left[\epsilon_{l}+\beta\frac{q_{l}^{2}+p_{l}^{2}}{2}+\beta
p_{l}^{2}\right]\delta p_{l}\\\ \dot{\delta p_{l}}&=-\beta q_{l}p_{l}\delta
p_{l}-\left[\epsilon_{l}+\beta\frac{q_{l}^{2}+p_{l}^{2}}{2}+\beta
q_{l}^{2}\right]\delta q_{l}\end{dcases}.$ (D.27)

In the action-angle variables, (D.27) obtains a much simpler form

$\begin{dcases}\dot{\delta\Theta_{l}}&=\beta\delta J_{l}\\\ \dot{\delta
J_{l}}&=0\end{dcases},$ (D.28)

whose analytic solution can be easily found. In particular, the solution of
(D.28) for initial conditions $(\delta\Theta_{l}(0),\delta J_{l}(0))$, is

$\begin{dcases}\delta\Theta_{l}(t)&=\beta\delta
J_{l}(0)t+\delta\Theta_{l}(0)\\\ \delta J_{l}(t)&=\delta
J_{l}(0)\end{dcases}.$ (D.29)

The transformation relating (D.27) to (D.28), can be obtained via total
differentiation of (D.11) and (D.12) and has the form

$\delta\Omega_{A}:\left(\delta q_{l},\delta
p_{l}\right)\mapsto\left(\delta\Theta_{l},\delta
J_{l}\right)\colon\begin{dcases}\delta\Theta_{l}(t)&=\frac{p_{l}(t)}{2J_{l}(t)}\delta
q_{l}(t)-\frac{q_{l}(t)}{2J_{l}(t)}\delta p_{l}(t)\\\ \delta
J_{l}(t)&=q_{l}(t)\delta q_{l}(t)+p_{l}(t)\delta p_{l}(t)\end{dcases},$ (D.30)

while its inverse transformation is given by

$\delta\Omega^{-1}_{A}:\left(\delta\Theta_{l},\delta
J_{l}\right)\mapsto\left(\delta q_{l},\delta
p_{l}\right)\colon\begin{dcases}\delta
q_{l}(t)&=\left[\frac{1}{\sqrt{2J_{l}(t)}}\sin\Theta_{l}(t)\right]\delta
J_{l}(t)+\left[\sqrt{2J_{l}(t)}\cos\Theta_{l}(t)\right]\delta\Theta_{l}(t)\\\
\delta
p_{l}(t)&=\left[\frac{1}{\sqrt{2J_{l}(t)}}\cos\Theta_{l}(t)\right]\delta
J_{l}(t)-\left[\sqrt{2J_{l}(t)}\sin\Theta_{l}(t)\right]\delta\Theta_{l}(t)\end{dcases}.$
(D.31)

The equivalence of solution (D.29) (in action-angle variables) and (D.30) (in
position and momentum variables) can be easily seen by differentiating (D.30).
Let us see that in more detail for the $\delta\Theta_{l}$ variable. Applying
the time derivative to the first equation of (D.30)

$\begin{split}\frac{d\delta\Theta_{l}(t)}{dt}&=\frac{1}{2J_{l}(t)}\left[\dot{p_{l}}(t)\delta
q_{l}(t)+p_{l}(t)\dot{\delta q_{l}}(t)-\dot{q_{l}}(t)\delta
p_{l}(t)-q\dot{\delta p_{l}}(t)\right],\\\
&=\frac{1}{2J_{l}(t)}\left\\{\left[-q_{l}\left(\epsilon_{l}+\beta\frac{q_{l}^{2}+p_{l}^{2}}{2}\right)\right]\delta
q_{l}(t)+p_{l}(t)\left[\beta q_{l}p_{l}\delta
q_{l}+\left(\epsilon_{l}+\beta\frac{q_{l}^{2}+p_{l}^{2}}{2}+\beta
p_{l}^{2}\right)\delta p_{l}\right]\right.\\\
&-\left.\left[p_{l}\left(\epsilon_{l}+\beta\frac{q_{l}^{2}+p_{l}^{2}}{2}\right)\right]\delta
p_{l}(t)-q_{l}(t)\left[-\left(\epsilon_{l}+\beta\frac{q_{l}^{2}+p_{l}^{2}}{2}+\beta
q_{l}^{2}\right)\delta q_{l}-\beta q_{l}p_{l}\delta p_{l}\right]\right\\},\\\
&=\frac{1}{2J_{l}(t)}\left[2J_{l}\beta\left(q_{l}\delta q_{l}+p_{l}\delta
p_{l}\right)\right],\\\ &=\beta\delta J_{l}(t),\end{split}$ (D.32)

which is just the first equation of (D.28). A similar result for $\delta
J_{l}(t)$ (D.30) can also be obtained.

In practice in order to obtain the solution of the system’s variational
equations we do the following: given an initial condition
$\left(q_{l}(0),p_{l}(0),\delta q_{l}(0),\delta p_{l}(0)\right)$, we compute
the quantities

$J_{l}(0)=\frac{1}{2}(q_{l}(0)^{2}+p_{l}(0)^{2}),\quad\alpha_{l}=\epsilon_{l}+\beta
J_{l}(0),\quad l=1,2,\ldots,N,$

and

$\delta J_{l}(0)=q_{l}(0)\delta q_{l}(0)+p_{l}(0)\delta
p_{l}(0),\quad\delta\Theta_{l}(0)=\frac{p_{l}(0)}{2J_{l}(0)}\delta
q_{l}(0)-\frac{q_{l}(0)}{2J_{l}(0)}\delta p_{l}(0),\quad l=1,2,\ldots,N.$

Thus the solution (D.31) of the variational equations takes the form

$\begin{dcases}\delta
q_{l}(t)&=\frac{q_{l}(0)\cos\alpha_{l}t+p_{l}(0)\sin\alpha_{l}t}{2J_{l}(0)}\delta
J_{l}(0)+\left(p_{l}(0)\cos\alpha_{l}t-q_{l}(0)\sin\alpha_{l}t\right)\left(\beta\delta
J_{l}(0)t+\delta\Theta_{l}(0)\right)\\\ \delta
p_{l}(t)&=\frac{p_{l}(0)\cos\alpha_{l}t-q_{l}(0)\sin\alpha_{l}t}{2J_{l}(0)}\delta
J_{l}(0)-\left(q_{l}(0)\cos\alpha_{l}t+p_{l}(0)\sin\alpha_{l}t\right)\left(\beta\delta
J_{l}(0)t+\delta\Theta_{l}(0)\right)\end{dcases},$ (D.33)

using the equalities

$\sin\Theta_{l}(t)=\frac{q_{l}(t)}{\sqrt{2J_{l}(t)}}=\frac{q_{l}(0)\cos\alpha_{l}t+p_{l}(0)\sin\alpha_{l}t}{\sqrt{2J_{l}(0)}},$

and

$\sqrt{2J_{l}(t)}\cos\Theta(t)=p_{l}(t)=p_{l}(0)\cos\alpha_{l}t-q_{l}(0)\sin\alpha_{l}t.$

along with (D.29).

In summary, the expression of the operator $e^{\tau\bm{L}_{A_{V}}}$
propagating $(q_{l},p_{l},\delta q_{l},\delta p_{l})$ at time $t$ to
$(q_{l}^{\prime},p_{l}^{\prime},\delta q_{l}^{\prime},\delta p_{l}^{\prime})$
at time $t+\tau$ is given by

$e^{\tau\bm{L}_{A_{V}}}\colon\begin{dcases}q_{l}^{\prime}&=q_{l}\cos(\tau\alpha_{l})+p_{l}\sin(\tau\alpha_{l})\\\
p_{l}^{\prime}&=p_{l}\cos(\tau\alpha_{l})-q_{l}\sin(\tau\alpha_{l})\\\ \delta
q_{l}^{\prime}&=\frac{q_{l}\cos(\tau\alpha_{l})+p_{l}\sin(\tau\alpha_{l})}{2J_{l}}\delta
J_{l}+\left(p_{l}\cos(\tau\alpha_{l})-q_{l}\sin(\tau\alpha_{l})\right)\left(\beta\delta
J_{l}\tau+\delta\Theta_{l}\right)\\\ \delta
p_{l}^{\prime}&=\frac{p_{l}\cos(\tau\alpha_{l})-q_{l}\sin(\tau\alpha_{l})}{2J_{l}}\delta
J_{l}-\left(q_{l}\cos(\tau\alpha_{l})+p_{l}\sin(\tau\alpha_{l})\right)\left(\beta\delta
J_{l}\tau+\delta\Theta_{l}\right)\end{dcases},$ (D.34)

with

$J_{l}=\frac{1}{2}(q_{l}^{2}+p_{l}^{2}),\quad\alpha_{l}=\epsilon_{l}+\beta
J_{l},$ $\delta J_{l}=q_{l}\delta q_{l}+p_{l}\delta
p_{l},\quad\delta\Theta_{l}=\frac{p_{l}}{2J_{l}}\delta
q_{l}-\frac{q_{l}}{2J_{l}}\delta p_{l}.$

Let us now discuss in detail a particular problematic case appearing in (D.33)
and (D.34), namely the situation where $J_{l}=(q_{l}^{2}+p_{l}^{2})/2=0$. In
that case, $q_{l}=p_{l}=0$. Then, from (D.27), we see that the variational
equations take the form

$\begin{dcases}\dot{\delta q_{l}}=\epsilon_{l}\delta p_{l}&\\\ \dot{\delta
p_{l}}=-\epsilon_{l}\delta q_{l}&\\\ \end{dcases}.$ (D.35)

The solution of this set of ODEs is given by

$\begin{dcases}\delta q_{l}(t)=\delta
q_{l}(0)\cos\left(\epsilon_{l}t\right)+\delta
p_{l}(0)\sin\left(\epsilon_{l}t\right)&\\\ \delta p_{l}(t)=\delta
p_{l}(0)\cos\left(\epsilon_{l}t\right)-\delta
q_{l}(0)\sin\left(\epsilon_{l}t\right)\end{dcases},$ (D.36)

assuming we started from $(\delta q_{l}(0),\delta p_{l}(0))$ at $t=0$.
Consequently, the general form of the operator $e^{\tau\bm{L}_{A_{V}}}$ is
given by (D.34) except for the particular case $J_{l}=0$ where we have

$e^{L_{A_{V}}\tau}:=\begin{dcases}q_{l}^{\prime}=q_{l}=0\\\
p_{l}^{\prime}=p_{l}=0\\\ \delta q_{l}^{\prime}=\delta
q_{l}\cos\left(\epsilon_{l}\tau\right)+\delta
p_{l}\sin\left(\epsilon_{l}\tau\right)\\\ \delta p_{l}^{\prime}=\delta
p_{l}\cos\left(\epsilon_{l}\tau\right)-\delta
q_{l}\sin\left(\epsilon_{l}\tau\right)\end{dcases}.$ (D.37)

####  D.1.2 Analytic expression for the operator $e^{\tau{\bf L}_{B_{V}}}$

Based on the analysis presented in the appendix of [Gerlach et al., 2016], let
us now determine the solution of both the equations of motion and the
variational equations of the $B$ Hamiltonian part (D.8)

$B=-\sum_{l}p_{l}p_{l+1}+q_{l}q_{l+1},$ (D.38)

of the 1D DDNLS system (D.1). More specifically, we will determine the
operator $e^{\tau\bm{L}_{B_{V}}}$ propagating an initial condition $\bm{Z}$ at
time $t$ to $\bm{Z}^{\prime}$ at time $t+\tau$ according to

$e^{\tau L_{B_{V}}}\colon\bm{Z}^{\prime}=\bm{\Gamma}(\tau)\cdot\bm{Z},$ (D.39)

with $\bm{Z}=\left(\bm{z}\leavevmode\nobreak\ \delta\bm{z}\right)^{T}$ such
that $\bm{z}=\left(\bm{q}\leavevmode\nobreak\
\bm{p}\right)^{T}=\left(q_{1}\leavevmode\nobreak\ q_{2}\leavevmode\nobreak\
\ldots\leavevmode\nobreak\ q_{N}\leavevmode\nobreak\ p_{1}\leavevmode\nobreak\
p_{2}\leavevmode\nobreak\ \ldots\leavevmode\nobreak\ p_{N}\right)^{T}$,
$\delta\bm{z}=\bm{w}=\left(\delta\bm{q}\leavevmode\nobreak\
\delta\bm{p}\right)^{T}=\left(\delta q_{1}\leavevmode\nobreak\ \delta
q_{2}\leavevmode\nobreak\ \ldots\leavevmode\nobreak\ \delta
q_{N}\leavevmode\nobreak\ \delta p_{1}\leavevmode\nobreak\ \delta
p_{2}\leavevmode\nobreak\ \ldots\leavevmode\nobreak\ \delta p_{N}\right)^{T}$
and $\bm{\Gamma}$ being a $4N\times 4N$ matrix. Since SIs are used for fixed
integration time steps, matrix $\bm{\Gamma}(\tau)$ in (D.39) is constant (for
fixed $\tau$ value). In addition, Hamiltonian $B$ (D.38) is a quadratic
function in its variables $q_{l}$ and $p_{l}$ which in turn means that the
equations of motion, and the variational equations are linear in $\bm{z}$, and
$\bm{\delta z}$ and independent of each other (see (D.42) and (D.59) below).
Consequently, matrix $\bm{\Gamma}(\tau)$ takes the form

$\bm{\Gamma}\left(\tau\right)=\begin{bmatrix}\bm{C}\left(\tau\right)&\bm{0}_{2N}\\\
\bm{0}_{2N}&\bm{C}^{\prime}\left(\tau\right)\end{bmatrix},$ (D.40)

with $\bm{C}(\tau)$ and $\bm{C}^{\prime}\left(\tau\right)$ being appropriate
$2N\times 2N$ matrices respectively evolving the equations of motion and the
variational equations according to

$\begin{bmatrix}\bm{q}^{\prime}\\\
\bm{p}^{\prime}\end{bmatrix}=\bm{C}\left(\tau\right)\cdot\begin{bmatrix}\bm{q}\\\
\bm{p}\end{bmatrix},\quad\begin{bmatrix}\delta\bm{q}^{\prime}\\\
\delta\bm{p}^{\prime}\end{bmatrix}=\bm{C}^{\prime}\left(\tau\right)\cdot\begin{bmatrix}\delta\bm{q}\\\
\delta\bm{p}\end{bmatrix}.$ (D.41)

More specifically, the equations of motion of Hamiltonian $B$ (D.38) take the
form

$\frac{d\bm{z}}{dt}=\bm{L}_{B}\bm{z}\colon\begin{dcases}\dot{q}_{1}=-p_{2}&\\\
\dot{q}_{l}=-p_{l-1}-p_{l+1}&\quad 2\leq l\leq N-1\\\ \dot{q}_{N}=-p_{N-1}&\\\
\dot{p}_{1}=q_{2}&\\\ \dot{p}_{l}=q_{l-1}+q_{l+1}&\quad 2\leq l\leq N-1\\\
\dot{p}_{N}=q_{N-1}\end{dcases}.$ (D.42)

Equation (D.42) can be written as

$\frac{d\bm{z}}{dt}=\bm{L}_{B}\bm{z}=\begin{pmatrix}\bm{0}&\bm{A}\\\
-\bm{A}&\bm{0}\end{pmatrix}\cdot\bm{z},$ (D.43)

with $\bm{A}$ being a tridiagonal $N\times N$ matrix whose main diagonal
elements are equal to zero ($\bm{A}_{l,l}=0,l=1,2,\ldots,N$), while all the
elements of the first upper and lower diagonals are equal to $-1$
($\bm{A}_{l,l+1}=\bm{A}_{l-1,l}=-1$) i.e.

$\bm{A}=\begin{pmatrix}0&-1&0&\dots&0&0\\\ -1&0&-1&\dots&0&0\\\
\vdots&\vdots&\ddots&\ddots&\vdots&\vdots\\\ 0&0&0&\dots&0&-1\\\
0&0&0&\dots&-1&0\end{pmatrix}.$ (D.44)

The solution of (D.43) for an integration time step $\tau$ can be expressed in
terms of a Lie series expansion (3.39) as

$\bm{z}(t+\tau)=e^{\tau\bm{L}_{B}}\bm{z}(t)=\sum_{k=0}^{\infty}\frac{\tau^{k}}{k!}\bm{L}_{B}^{k}\bm{z}(t)=\bm{C}(\tau)\bm{z}(t),$
(D.45)

with

$\bm{C}\left(\tau\right)=\sum_{k=0}^{\infty}\frac{\tau^{k}}{k!}\bm{L}_{B}^{k}.$
(D.46)

As we mentioned in Sec. 3.6.2, the next step to determine the expression of
$\bm{C}(\tau)$ is to find the form of the operator $\bm{L}_{B}^{k}$ in (D.46).
Luckily, this task can be performed due to the rather simple form of operator
$\bm{L}_{B}$ (D.43). Let us initially note that

$\bm{L}_{B}^{2k}=(-1)^{k}\begin{pmatrix}\bm{A}^{2k}&\bm{0}\\\
\bm{0}&\bm{A}^{2k}\end{pmatrix},\quad\bm{L}_{B}^{2k+1}=(-1)^{k}\begin{pmatrix}\bm{0}&\bm{A}^{2k+1}\\\
-\bm{A}^{2k+1}&\bm{0}\end{pmatrix}.$ (D.47)

Then, using (D.46) and (D.47), we get

$\displaystyle\bm{C}(\tau)$
$\displaystyle=\begin{pmatrix}\sum_{k=0}^{\infty}\frac{(-1)^{k}}{(2k)!}\bm{A}^{2k}\tau^{2k}&\bm{0}\\\
\bm{0}&\sum_{k=0}^{\infty}\frac{(-1)^{k}}{(2k)!}\bm{A}^{2k}\tau^{2k}\end{pmatrix}$
(D.48)
$\displaystyle+\begin{pmatrix}\bm{0}&\sum_{k=0}^{\infty}\frac{(-1)^{k}}{(2k+1)!}\bm{A}^{2k+1}\tau^{2k+1}\\\
-\sum_{k=0}^{\infty}\frac{(-1)^{k}}{(2k+1)!}\bm{A}^{2k+1}\tau^{2k+1}&\bm{0}\end{pmatrix},$
(D.49)

which is equivalent to

$\bm{C}(\tau)=\begin{pmatrix}\cos\left(\bm{A}\tau\right)&\sin\left(\bm{A}\tau\right)\\\
-\sin\left(\bm{A}\tau\right)&\cos\left(\bm{A}\tau\right)\end{pmatrix}.$ (D.50)

The evaluation of the elements of matrices $\cos\left(\bm{A}\tau\right)$ and
$\sin\left(\bm{A}\tau\right)$ can be obtained through the determination of the
eigenvalues and eigenvectors of matrix $\bm{A}$ (D.44) itself [Smith, 1985].
In particular, the eigenvalues $\lambda_{k}$ of $\bm{A}$ are symmetric with
respect to zero and are given by

$\lambda_{k}=-2\cos\left(\frac{k\pi}{N+1}\right),\quad k=1,2,\dots,N,$ (D.51)

while the normalized eigenvectors are orthogonal (since $\bm{A}$ is a
symmetric matrix) and have the form

$\bm{v}_{k}=\frac{1}{V_{k}}\left(\sin\left(\frac{k\pi}{N+1}\right),\sin\left(\frac{2k\pi}{N+1}\right),\dots,\sin\left(\frac{jk\pi}{N+1}\right),\dots,\sin\left(\frac{Nk\pi}{N+1}\right)\right).$
(D.52)

The Euclidean norms $V_{k}$ of these eigenvectors are

$\displaystyle
V_{k}^{2}=\lVert\bm{v}_{k}\rVert^{2}=\sum_{j=1}^{N}\sin^{2}\left(\frac{jk\pi}{N+1}\right)$
$\displaystyle=\sum_{j=1}^{N}\left[\frac{1}{2}-\frac{1}{2}\sum_{j=1}^{N}\cos\left(\frac{2jk\pi}{N+1}\right)\right],$
(D.53)
$\displaystyle=\frac{2N+1}{4}-\frac{\sin\left(\frac{2N+1}{N+1}k\pi\right)}{4\sin\left(\frac{k\pi}{N+1}\right)},$
(D.54) $\displaystyle=\frac{N+1}{2}.$ (D.55)

We note that going from (D.53) to (D.55), we appropriately used the
identities,

$\sum_{j=0}^{N}\cos\left(jx\right)=\frac{1}{2}\left[1+\frac{\sin\left[\left(N+\frac{1}{2}\right)x\right]}{\sin\left(\frac{x}{2}\right)}\right],$
(D.56)

which can be found for example at [Jeffrey and Zwillinger, 2007] and

$\sin\left(\frac{2N+1}{N+1}k\pi\right)=-\sin\left(\frac{k\pi}{N+1}\right).$
(D.57)

The matrix $\bm{P}$ having as columns the eigenvectors (D.52) i.e.
$\bm{P}=(\bm{v}_{1}^{T}\leavevmode\nobreak\ \bm{v}_{2}^{T}\leavevmode\nobreak\
\dots\leavevmode\nobreak\ \bm{v}_{N}^{T})$ can be used to diagonalize $\bm{A}$
(D.44) so that

$\bm{A}=\bm{P}\cdot\bm{D}\cdot\bm{P}^{-1}=\bm{P}\cdot\bm{D}\cdot\bm{P}.$
(D.58)

We note that in (D.58) $\bm{P}$ is a unitary $N\times N$ matrix, i.e.
$\bm{P}=\bm{P}^{-1}=\bm{P}^{T}$, and $\bm{D}$ a diagonal $N\times N$ matrix
having as non-zero elements the eigenvalues (D.51) i.e.
$\bm{D}=diag\left(\lambda_{1},\lambda_{2},\dots,\lambda_{N}\right)$. Thus we
can write the matrices $\cos\left(\bm{A}\tau\right)$ and
$\sin\left(\bm{A}\tau\right)$ defining $\bm{C}(\tau)$ in (D.50) as

$\cos\left(\bm{A}\tau\right)=\bm{P}\cdot\bm{D}_{c}\cdot\bm{P},\quad\sin\left(\bm{A}\tau\right)=\bm{P}\cdot\bm{D}_{s}\cdot\bm{P},$

with
$\bm{D}_{c}=diag\left(\cos\left(\lambda_{1}\tau\right),\cos\left(\lambda_{2}\tau\right),\dots,\cos\left(\lambda_{N}\tau\right)\right)$,
and
$\bm{D}_{s}=diag\left(\sin\left(\lambda_{1}\tau\right),\sin\left(\lambda_{2}\tau\right),\dots,\sin\left(\lambda_{N}\tau\right)\right)$.

A similar treatment can be performed for the variational equations of the
Hamiltonian $B$ (D.38)

$\frac{d\bm{w}}{dt}=\bm{L}_{B_{V}}\cdot\bm{w}\colon\begin{dcases}\dot{\delta
q}_{1}=-\delta p_{2}&\\\ \dot{\delta q}_{l}=-\delta p_{l-1}-\delta
p_{l+1}&\quad 2\leq l\leq N-1\\\ \dot{\delta q}_{N}=-\delta p_{N-1}&\\\
\dot{\delta p}_{1}=\delta q_{2}&\\\ \dot{p}_{l}=\delta q_{l-1}+\delta
q_{l+1}&\quad 2\leq l\leq N-1\\\ \dot{p}_{N}=\delta q_{N-1}\end{dcases},$
(D.59)

which can be rewritten as

$\frac{d\bm{w}}{dt}=\bm{L}_{B_{V}}\bm{w}=\begin{pmatrix}\bm{0}&\bm{A}\\\
\bm{-A}&\bm{0}\end{pmatrix}\cdot\bm{w},$ (D.60)

is the system of ODEs (D.43). Thus,

$\bm{C}^{\prime}\left(\tau\right)=\bm{C}\left(\tau\right),$ (D.61)

in (D.40) and (D.41).

###  D.2 Three-part split symplectic integration of the 1D and 2D DDNLS models

Following [Skokos et al., 2014; Gerlach et al., 2016] and mainly [Danieli et
al., 2019a], we present here the analytic expressions of the operators
involved in the symplectic integration of the 1D and 2D DDNLS systems when
they are split in three integrable parts.

####  D.2.1 The 1D DDNLS system

The Hamiltonian function $\mathcal{H}_{1D}$ (D.2) of the 1D DDNLS system can
be split into three integrable parts as

$\mathcal{H}_{1D}=\mathcal{A}+\mathcal{B}+\mathcal{C},$ (D.62)

with

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

Obviously, the $\mathcal{A}$ part in (D.63) coincides with the $A$ part in
(D.8) so that the corresponding operator $e^{\tau\bm{L}_{\mathcal{A}_{V}}}$ is
given by (D.34) and (D.37). The Hamiltonian functions $\mathcal{B}$ and
$\mathcal{C}$ on the other hand possess $N$ cyclic coordinates as
$\mathcal{B}$ depends only on momenta $p_{l}$ and $\mathcal{C}$ on the
positions $q_{l}$. Therefore it is easy to find the analytic expression of the
operators $e^{\tau\bm{L}_{\mathcal{B}_{V}}}$ and
$e^{\tau\bm{L}_{\mathcal{C}_{V}}}$ for the integration of their equations of
motion and the variational equations.

More specifically, for the $\mathcal{B}$, we have the following system of
equations of motion and variational equations

$\frac{d\bm{Z}}{dt}=\bm{L}_{\mathcal{B}_{V}}\bm{Z}\colon\begin{dcases}\dot{q}_{1}=-p_{2}&\\\
\dot{q}_{l}=-p_{l-1}-p_{l+1}&l=2,3,\dots,N-1\\\ \dot{q}_{N}=-p_{N-1}&\\\
\dot{p}_{l}=0&l=1,2,\dots,N\\\ \dot{\delta q}_{1}=-\delta p_{2}&\\\
\dot{\delta q}_{l}=-\delta p_{l-1}-\delta p_{l+1}&l=2,3\dots,N-1\\\
\dot{\delta q}_{N}=-\delta p_{N-1}&\\\ \dot{\delta p}_{l}=0&l=1,2,\dots,N\\\
\end{dcases},$ (D.64)

which can be easily solved, leading to

$e^{\tau\bm{L}_{\mathcal{B}_{V}}}\colon\begin{dcases}q_{1}^{\prime}=q_{1}-\tau
p_{2}&\\\ q_{l}^{\prime}=q_{l}-\tau(p_{l-1}+p_{l+1})&l=2,3,\dots,N-1\\\
q_{N}^{\prime}=q_{N}-\tau p_{N-1}&\\\ p_{l}^{\prime}=p_{l}&l=1,2,\dots,N\\\
\delta q_{1}^{\prime}=\delta q_{1}-\tau\delta p_{2}&\\\ \delta
q_{l}^{\prime}=\delta q_{l}-\tau(\delta p_{l-1}+\delta
p_{l+1})&l=2,3,\dots,N-1\\\ \delta q_{N}^{\prime}=\delta q_{N}-\tau\delta
p_{N-1}&\\\ \delta p_{l}^{\prime}=\delta p_{l}&l=1,2,\dots,N\\\ \end{dcases}.$
(D.65)

Similarly, the system of equations of motion and variational equations for the
$\mathcal{C}$ Hamiltonian (D.63) is

$\frac{d\bm{Z}}{dt}=\bm{L}_{\mathcal{C}_{V}}\bm{Z}\colon\begin{dcases}\dot{p}_{1}=\delta
q_{2}&\\\ \dot{p}_{l}=\delta q_{l-1}+\delta q_{l+1}&l=2,3,\dots,N-1\\\
\dot{p}_{N}=\delta q_{N-1}&\\\ \dot{q}_{l}=0&l=1,2,\dots,N\\\ \dot{\delta
p}_{1}=\delta q_{2}&\\\ \dot{\delta p}_{l}=\delta q_{l-1}+\delta
q_{l+1}&l=2,3,\dots,N-1\\\ \dot{\delta p}_{N}=\delta q_{N-1}&\\\ \dot{\delta
q}_{l}=0&l=1,2,\dots,N\\\ \end{dcases},$ (D.66)

and its solution is given by

$e^{\tau\bm{L}_{\mathcal{C}_{V}}}\colon\begin{dcases}q_{l}^{\prime}=q_{l}&l=1,2,\dots,N\\\
p_{1}^{\prime}=p_{1}+\tau q_{2}&\\\
p_{l}^{\prime}=p_{l}+\tau(q_{l-1}+q_{l+1})&l=2,3,\dots,N-1\\\
p_{N}^{\prime}=p_{N}+\tau q_{N-1}&\\\ \delta q_{l}^{\prime}=\delta
q_{l}&l=1,2,\dots,N\\\ \delta p_{1}^{\prime}=\delta p_{1}+\tau\delta q_{2}&\\\
\delta p_{l}^{\prime}=\delta p_{l}+\tau(\delta q_{l-1}+\delta
q_{l+1})&l=2,3,\dots,N-1\\\ \delta p_{N}^{\prime}=\delta p_{N}+\tau\delta
q_{N-1}&\\\ \end{dcases}.$ (D.67)

####  D.2.2 The 2D DDNLS system

We can treat the 2D DDNLS Hamiltonian $\mathcal{H}_{2D}$ (D.5) in a similar
manner to the one implemented for the 1D DDNLS system in Sec. D.2.1. In
particular, the Hamiltonian can be split

$\mathcal{H}_{2D}=\mathsf{A}+\mathsf{B}+\mathsf{C},$ (D.68)

with

$\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},$
(D.69)
$\mathsf{B}=\sum_{l,m}-p_{l,m+1}p_{l,m}-p_{l+1,m}p_{l,m},\qquad\mathsf{C}=\sum_{l,m}-q_{l,m+1}q_{l,m}-q_{l+1,m}q_{l,m}.$
(D.70)

For the Hamiltonian $\mathsf{A}$ (D.69), following similar steps to the ones
performed in Sec. D.1.1, we obtain the following set of equations of motion
and variational equations

$\begin{split}\frac{d\bm{Z}}{dt}=\bm{L}_{\mathsf{A}_{V}}\bm{Z}&:\left\\{\begin{array}[]{rl}\dot{q}_{l,m}&=p_{l,m}\dot{\Theta}_{l,m}\\\
\dot{p}_{l,m}&=-q_{l,m}\dot{\Theta}_{l,m},\\\ \dot{\delta
q}_{l,m}&=\left[\dot{\Theta}_{l,m}+\beta p_{l,m}^{2}\right]\delta
p_{l,m}+\beta q_{l,m}p_{l,m}\delta q_{l,m}\\\ \dot{\delta
p}_{l,m}&=-\left[\dot{\Theta}_{l,m}+\beta q_{l,m}^{2}\right]\delta
q_{l,m}-\beta q_{l,m}p_{l,m}\delta p_{l,m}\end{array}\right.,\\\ \end{split}$
(D.71)

where $\dot{\Theta}_{l,m}=\epsilon_{l,m}+\beta(q_{l,m}^{2}+p_{l,m}^{2})/2$ are
constants of motion. Setting

$J_{l,m}=\frac{1}{2}(q_{l,m}^{2}+p_{l,m}^{2}),\quad\alpha_{l,m}=\epsilon_{l,m}+\beta
J_{l,m},$

and

$\delta J_{l,m}=q_{l,m}\delta q_{l,m}+p_{l,m}\delta
p_{l,m},\quad\delta\Theta_{l,m}=\frac{p_{l,m}}{2J_{l,m}}\delta
q_{l,m}-\frac{q_{l,m}}{2J_{l,m}}\delta p_{l,m},$

the solution of system (D.71) for $J_{l,m}\neq 0$ is obtained through the
application of the operator

$\begin{split}e^{\tau\bm{L}_{\mathsf{A}_{V}}}\colon\left\\{\begin{array}[]{rl}q_{l,m}^{\prime}&=q_{l,m}\cos(\tau\alpha_{l,m})+p_{l,m}\sin(\tau\alpha_{l,m})\\\
p_{l,m}^{\prime}&=p_{l,m}\cos(\tau\alpha_{l,m})-q_{l,m}\sin(\tau\alpha_{l,m})\\\
\delta
q_{l,m}^{\prime}&=\frac{q_{l,m}\cos(\tau\alpha_{l,m})+p_{l,m}\sin(\tau\alpha_{l,m})}{2J_{l,m}}\delta
J_{l,m}+\left(p_{l,m}\cos(\tau\alpha_{l,m})-q_{l,m}\sin(\tau\alpha_{l,m})\right)\left(\beta\delta
J_{l,m}\tau+\delta\Theta_{l,m}\right)\\\ \delta
p_{l,m}^{\prime}&=\frac{p_{l,m}\cos(\tau\alpha_{l,m})-q_{l,m}\sin(\tau\alpha_{l,m})}{2J_{l,m}}\delta
J_{l,m}-\left(q_{l,m}\cos(\tau\alpha_{l,m})+p_{l,m}\sin(\tau\alpha_{l,m})\right)\left(\beta\delta
J_{l,m}\tau+\delta\Theta_{l,m}\right)\end{array}\right.,\end{split}$ (D.72)

or when $J_{l,m}=0$ through the operator

$e^{\bm{L}_{\mathsf{A}_{V}}\tau}\colon\begin{dcases}q_{l,m}^{\prime}=q_{l,m}=0\\\
p_{l,m}^{\prime}=q_{l,m}=0\\\ \delta q_{l,m}^{\prime}=\delta
q_{l,m}\cos\left(\epsilon_{l,m}\tau\right)+\delta
p_{l,m}\sin\left(\epsilon_{l,m}\tau\right)\\\ \delta p_{l,m}^{\prime}=\delta
p_{l,m}\cos\left(\epsilon_{l,m}\tau\right)-\delta
q_{l,m}\sin\left(\epsilon_{l,m}\tau\right)\end{dcases}.$ (D.73)

For the Hamiltonian functions $\mathsf{B}$ and $\mathsf{C}$ (D.70) the set of
equations of motion and variational equations are respectively