Computational efficiency of numerical integration methods for the tangent dynamics of many-body Hamiltonian systems in one and two spatial dimensions

1.   Introduction

A huge number of important problems in physics, astronomy, chemistry, etc. are modeled via sets of ordinary differential equations (ODEs) governed by the Hamiltonian formalism. Due to non-integrability, the investigation of the time evolution of these problems, and generally their properties, often rely solely on numerical techniques. As modern research requires numerical simulations to be pushed to their very limit (e.g., large integration times, macroscopic limits), a methodical assessment of the properties of different numerical methods becomes a fundamental issue. Such studies allow to highlight the most suitable scheme for both general purposes and specific tasks, according to criteria of stability, accuracy, simplicity and efficiency.

The beginning of the era of computational physics is considered to be the computer experiment performed in the 1950s by Fermi, Pasta, Ulam and Tsingou (FPUT) ^[1,2,3] to observe energy equipartition due to weak-anharmonic coupling in a set of harmonic oscillators. Indeed, the breaking of integrability in Hamiltonian systems is often performed with the introduction of nonlinear terms in the equations of motion. These additional nonlinear terms describe new physical processes, and led to important questions and significant advancements in condensed matter physics. For instance, the FPUT problem has been used to answer questions related to ergodicity, dynamical thermalization and chaos occurrence (see e.g., ^[4,5,6] and references therein) and led to the observation of solitons ^[7,8], and progress in Hamiltonian chaos ^[9]. The interested reader can find a concise review of numerical results concerning the FPUT system in ^[10].

Another important example concerns disordered media. In 1958, Anderson ^[11] theoretically showed the appearance of energy localization in one-dimensional lattices with uncorrelated disorder (a phenomenon which is now referred to as Anderson localization). This phenomenon was later investigated also for two-dimensional lattices ^[12]. An important question which has attracted wide attention in theory, numerical simulations and experiments is what happens when nonlinearity (which appears naturally in real world problems) is introduced in a disordered system ^{[13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30]}.

Two basic Hamiltonian models are at the center of these investigations: the disordered Klein-Gordon (DKG) chain and the disordered, discrete nonlinear Schrödinger equation (DDNLS). In Refs. ^{[22,23,24,25,26,28,29,30,31,32,33,34]} it was shown that Anderson localization is destroyed by the presence of nonlinearity, resulting to the subdiffusive spreading of energy due to deterministic chaos. Such results rely on the accurate, long time integration of the systems' equations of motion and variational equations. We note that the variational equations are a set of ODEs describing the time evolution of small deviations from a given trajectory, something which is needed for the computation of chaos indicators like the maximum Lyapunov exponent (mLE) ^[35,36,37]. The numerical integration of these sets of ODEs can be performed by any general purpose integrator. For example Runge-Kutta family schemes are still used ^[38,39]. Another category of integrators is the so-called symplectic integrators (SIs), which were particularly designed for Hamiltonian systems (see e.g., ^{[40,41,42,43,44]} and references therein). The main characteristic of SIs is that they conserve the symplectic structure of Hamiltonian systems, so that the numerical approximation they provide corresponds to the exact solution of a system which can be considered as a perturbation of the original one. Consequently the error in the computed value of the Hamiltonian function (usually refereed to as the system's energy) is kept almost constant over long integration times. This almost constant error can be used as an indicator of the accuracy of the used integration scheme.

SIs have been successfully implemented for the long time integration of multidimensional Hamiltonian systems like the $\alpha$- and $\beta$-FPUT systems and FPUT-like models, the DKG, the DDNLS and systems of coupled rotors (see e.g., ^{[22,23,24,25,26,28,29,30,32,33,45,46,47,48,49,50,51,52,53,54]}). In these studies SIs belonging to the so-called SABA family ^[55] were mostly used, since the FPUT, the DKG and systems of coupled rotors can be split into two integrable parts (the system's kinetic and the potential energy), while in some more recent studies ^[33,56] the $ABA864$ ^[57] SI was implemented as it displayed an even better performance. As the DDNLS model is not explicitly decomposed into two integrable parts, the implementation of SABA schemes requires a special treatment involving the application of fast Fourier transforms which are computationally time consuming ^[29]. In ^[49,50] it was shown that for the DDNLS model, the split of the Hamiltonian function into three integrable parts is possible, and this approach proved to be the best choice for the model's integration with respect to both speed and accuracy. It is worth noting here that the numerical integration of the variational equations by SIs is done by the so-called Tangent Map Method ^[58,59,60].

The intention of this work is to present several efficient (symplectic and non-symplectic) integration techniques for the time evolution of both the equations of motion and the variational equations of multidimensional Hamiltonian systems. In particular, in our study we use these techniques to investigate the chaotic dynamics of the one-dimensional (1D) FPUT system, as well as one- and two-dimensional (2D) DDNLS models. We carry out this task by considering particular initial conditions for both the systems and for the deviation vectors, which we evolve in time requiring a predefined level of energy accuracy. Then, we record the needed CPU time to perform these numerical simulations and highlight those methods that ensure the smallest CPU times for obtaining accurate results. A similar analysis for the 1D and 2D DKG systems was performed in ^[56].

The paper is organized as follows. In Section 2 we introduce the models considered in our study. In Section 3 we introduce in detail the symplectic and non-symplectic schemes implemented in our investigations. In Section 4 we present our numerical results and report on the efficiency of each studied scheme. In Section 5 we review our findings and discuss the outcomes of our study. Finally, in Appendix B we provide the explicit forms of the used tangent map operators. Moreover, we present there the explicit expressions of the tangent map operators for some additional important many-body systems, the so-called $\beta$-FPUT chain, the KG system and the classical XY model (a Josephson junction chain-JJC) in order to facilitate the potential application of SIs for these systems from the interested reader, although we do not present numerical results for these models.

2.   Models and Hamiltonian functions

In this work we focus our analysis on the 1D FPUT system and the 1D and 2D DDNLS models. In what follows we briefly present the Hamiltonian functions of these systems.

2.1.   The $\alpha$-Fermi-Pasta-Ulam-Tsingou model

The 1D FPUT model ^[1,2,3] was first introduced in 1955 to study the road toward equilibrium of a chain of harmonic oscillators in the presence of weak anharmonic interactions. Since then, this system has been widely used as a toy model for investigating energy equipartition and chaos in nonlinear lattices. In the literature there exist two types of the FPUT model the so-called $\alpha$- and $\beta$-FPUT systems. In our study we consider the $\alpha$-FPUT model whose Hamiltonian function reads

In Eq.(2.1), $q_i$ and $p_i$ are respectively the generalized coordinates and momenta of the $i$ lattice site and $\alpha$ is a real positive parameter. In our study we consider fixed boundary conditions $q_0 = p_0 = p_{N+1} = q_{N+1} = 0$. We also note that this model conserves the value of the total energy $H_{\text{1F}}$.

In contrast to the $\alpha$-FPUT system, the $\beta$-FPUT model is characterized by a quartic nonlinear term [see Eq. (A.26) in Appendix B.4]. The fact that the value of the Hamiltonian function of the $\beta$-FPUT model is bounded from below leads to differences between the two models. For example, phenomena of chopping time which occur in the $\alpha$- model do not appear in the $\beta$-FPUT system ^[61]. Further discussion of the differences between the $\alpha$- and the $\beta$- models can be found in ^[62,63,64].

2.2.   The 1D and 2D disordered discrete nonlinear Schrödinger equations

The DDNLS model describes anharmonic interactions between oscillators in disordered media and has been used to study the propagation of light in optical media or Bose-Einstein condensates through the famous Gross-Pitaevskii equation ^[52], as well as investigate, at a first approximation, several physical processes (e.g., electron tight binding systems ^[39]). The Hamiltonian function of the 1D DDNLS model reads

where $q_i$ and $p_i$ are respectively the generalized coordinates and momenta of the $i$ lattice site and the onsite energy terms $\epsilon _i$ are uncorrelated numbers uniformly distributed in the interval $\left[-\frac{W}{2}, \frac{W}{2} \right]$. The real, positive numbers $W$ and $\beta$ denote the disorder and the nonlinearity strength respectively. We also consider here fixed boundary conditions i.e., $q_0 = p_0 = p_{N+1} = q_{N+1} = 0$.

The two-dimensional version of the DDNLS model was considered in ^[30,65]. Its Hamiltonian function is

where $q_{i, j}$ and $p_{i, j}$ are respectively the generalized positions and momenta at site $(i, j)$ and $\epsilon _{i, j}$ are the disorder parameters uniformly chosen in the interval $\left[-\frac{W}{2}, \frac{W}{2}\right]$. We again consider fixed boundary conditions i.e., $q_{0, j} = p_{0, j} = q_{N+1, j} = p_{N+1, j} = 0$ for $1\leq j\leq M$ and $q_{i, 0} = p_{i, 0} = q_{i, M+1} = p_{i, M+1} = 0$ for $1\leq i\leq N$.

Additionally to the energies $H_{1D}$ [Eq. (2.2)] and $H_{2D}$ [Eq. (2.3)] both systems conserve their respective norms $S_{\text{1D}}$ and $S_{\text{2D}}$:

3.   Numerical integration schemes

The Hamilton equations of motion

of the $N$ degree of freedom (dof) Hamiltonian $H = H(\boldsymbol{ q}, \boldsymbol{ p})$, with $\boldsymbol{ q} = (q_1, q_2, \ldots, q_N)$ and $\boldsymbol{ p} = (p_1, p_2, \ldots, p_N)$ being respectively the system's generalized positions and momenta, can be expressed in the general setting of first order ODEs as

where $\boldsymbol{ x} = \left(\boldsymbol{ q}, \boldsymbol{ p} \right) = (x_1, x_2, \ldots, x_{N}, x_{N+1}, \ldots, x_{2N}) = (q_1, q_2, \ldots, q_N, p_1, p_2, \ldots, p_N)$ is a vector representing the position of the system in its phase space and $(\, \dot{}\, )$ denotes differentiation with respect to time $t$. In Eq. (3.2)

is the symplectic matrix with $\boldsymbol{I}_N$ and $\boldsymbol{O}_N$ being the $N\times N$ identity and the null matrices respectively, and

with $(^T)$ denoting the transpose matrix.

The variational equations (see for example ^[37,58]) govern the time evolution of a small perturbation $\boldsymbol{ w}(t)$ to the trajectory $\boldsymbol{ x}(t)$ with $\boldsymbol{w}(t) = (\delta \boldsymbol{q}(t), \delta \boldsymbol{p}(t)) = \left(\delta q_1 (t), \delta q_2 (t), \ldots, \delta q_N (t), \delta p_1(t), \delta p_2(t), \ldots, \delta p_N (t) \right)$ (which can also be written as $\boldsymbol{w}(t) = \delta \boldsymbol{x}(t) = \left(\delta x_1 (t), \ldots, \delta x_N (t), \delta x_{N+1}(t), \ldots, \delta x_{2N} (t) \right)$) and have the following form

where

are the $2N\times 2N$ elements of the Hessian matrix $\boldsymbol{ D}_{ H}^2 (\boldsymbol{ x}(t))$ of the Hamiltonian function $H$ computed on the phase space trajectory $\boldsymbol{ x}(t)$. Eq. (3.5) is linear in $\boldsymbol{w}(t)$, with coefficients depending on the system's trajectory $\boldsymbol{ x}(t)$. Therefore, one has to integrate the variational equations (3.5) along with the equations of motion (3.2), which means to evolve in time the general vector $\boldsymbol{ X}(t) = (\boldsymbol{ x}(t), \delta \boldsymbol{ x}(t))$ by solving the system

In what follows we will briefly describe several numerical schemes for integrating the set of equations (3.7).

3.1.   Non-symplectic integration schemes

Here we present the non-symplectic schemes we will use in this work: the Taylor series method and the Runge-Kutta discretization scheme. These methods are referred to be non-symplectic because they do not preserve the geometrical properties of Hamiltonian equations of motion (e.g., their symplectic structure). The immediate consequence of that is that they do not conserve constants of motion (e.g., the system's energy).

3.1.1.   Taylor series method - $TIDES$

The Taylor series method consists in expanding the solution $\boldsymbol{ X}(t)$ at time $t_0 + \tau$ in a Taylor series of $\boldsymbol{X}(t)$ at $t = t_0$

Once the solution $\boldsymbol{X}$ at time $t_0 + \tau$ is approximated, $\boldsymbol{ X} (t_0 + \tau)$ is considered as the new initial condition and the procedure is repeated to approximate $\boldsymbol{X} (t_0 + 2 \tau)$ and so on so forth. Further information on this integrator can be found in [66,Sec. Ⅰ.8]. If we consider in Eq. (3.8) only the first $n + 1$ terms we then account the resulting scheme to be of order $n$. In addition, in order to explicitly express this numerical scheme, one has to perform $n-1$ successive differentiations of the field vector $\boldsymbol{ f}$. This task becomes very elaborate for complex structured vector fields. One can therefore rely on successive differentiation to automate the whole process (see e.g., ^[67]). For the simulations reported in this paper we used the software package called $TIDES$ ^[67,68] which is freely available ^[69]. We particularly focused on the implementation of the $TIDES$ package as it has been already used in studies of lattice dynamics ^[60]. The $TIDES$ package comes as a Mathematica notebook in which the user provides the Hamiltonian function, the potential energy or the set of ODEs themselves. It then produces FORTRAN (or C) codes which can be compiled directly by any available compiler producing the appropriate executable programs. In addition, the $TIDES$ package allows us to choose both the integration time step $\tau$ and the desired 'one-step' precision of the integrator $\delta$. In practice, the integration time step is accepted if the estimated local truncation error is smaller than $\delta$.

It is worth noting that there exists an equivalent numerical scheme to the Taylor series method, derived from Lie series ^[70] (for more details see e.g., the appendix of ^[71]). Indeed, Eq. (3.7) can be expressed as

where $L_{HZ}$ is the Lie operator ^[72] defined as

The formal solution of Eq. (3.9) reads

and can be expanded as

This corresponds to a Lie series integrator of order $n$. Similarly to the Taylor series method, one has to find the analytical expression of the successive action of the operator $L_{HZ}$ onto the vector $\boldsymbol{X} (t_0)$. For further details concerning Lie series one can refer to ^[72,73]. The equivalence between the Lie and Taylor series approaches can be seen in the following way: for each element $x_i$ and $\delta x_i$ of the phase space vector $\boldsymbol{X}$ we compute

Therefore $L_{HZ}^0 \boldsymbol{X} = \boldsymbol{X}$, $L_{HZ}^1 \boldsymbol{X} = \frac{d \boldsymbol{X}}{dt}$, $L_{HZ}^2\boldsymbol{X} = \frac{d^2 \boldsymbol{X}}{dt^2}$ etc.

3.1.2.   A Runge-Kutta family scheme $-$ $DOP853$

A huge hurdle concerning the applicability of the Taylor and Lie series methods can be the explicit derivation of the differential operators (see [73] and references therein). Over the years other methods have been developed in order to overcome such issues and to efficiently and accurately approximate Eq. (3.8) up to a certain order in $n$. One way to perform this task was through the use of the well-known 'Runge-Kutta family' of algorithms (see e.g., [40,73]) which nowadays is one of the most popular general purpose schemes for numerical computations. This is a sufficient motivation for us to introduce a $s-$stage Runge-Kutta method of the form

where the real coefficients $a_{i, j}$, $b_i$ with $i, j = 1, \ldots, s$ are appropriately chosen in order to obtain the desired accuracy (see e.g., [66,Sec. Ⅱ.1]). Eq. (3.13) can be understood in the following way: in order to propagate the phase space vector $\boldsymbol{X}$ from time $t = t_0$ to $t = t_0+\tau$, we compute $s$ intermediate values $\boldsymbol{X}_1, \boldsymbol{X}_2, \ldots, \boldsymbol{X}_s$ and $s$ intermediate derivatives $\boldsymbol{k}_1, \boldsymbol{k}_2, \ldots, \boldsymbol{k}_s$, such that $\boldsymbol{k}_i = \boldsymbol{f} (\boldsymbol{\boldsymbol{X}_i})$, at times $t_i = t_0 + \tau \sum _{j = 1}^{i - 1} a_{ij}$. Then each $\boldsymbol{X}_i$ is found through a linear combination of the known $\boldsymbol{k}_i$ values, which are added to the initial condition $\boldsymbol{X}(t_0)$.

In this work we use a $12$-stage explicit Runge-Kutta algorithm of order $8$, called $DOP853$ [66,Sec. Ⅱ.5]. This method is the most precise scheme among the Runge-Kutta algorithms presented in ^[66] (see Fig. $4.2$ of ^[66]). A free access Fortran77 implementation of $DOP853$ is available in ^[74] (see also the appendix of ^[66]). As in the case of the $TIDES$ package, apart from the integration time step $\tau$, $DOP853$ also admits a 'one-step' precision $\delta$ based on embedded formulas of orders $5$ and $3$ (see [66,Sec. Ⅱ.10], or ^[60] and references therein for more details).

3.2.   Symplectic integration schemes

Hamiltonian systems are characterized by a symplectic structure (see e.g., [40,Chap. VI] and references therein) and they may also possess integrals of motion, like for example the energy, the angular momentum, etc. SIs are appropriate for the numerical integration of Hamiltonian systems as they keep the computed energy (i.e., the value of the Hamiltonian) of the system almost constant over the integration in time. Let us remark that, in general, SIs do not preserve any additional conserved quantities of the system (for a specific exception see ^[75]). SIs have already been extensively used in fields such as celestial mechanics ^[76], molecular dynamics ^[77], accelerator physics ^[43], condensed matter physics ^[22,29,33], etc.

Several methods can be used to build SIs ^[40]. For instance, it has been proved that under certain circumstances the Runge-Kutta algorithm can preserve the symplectic structure (see e.g., ^[78,79]). However, the most common way to construct SIs is through the split of the Hamiltonian function into the sum of two or more integrable parts. For example, many Hamiltonian functions are expressed as the sum of the kinetic and potential energy terms, with each one of them corresponding to an integrable system (see Appendices B.1 and B.4). Let us remark that, in general, even if each component of the total Hamiltonian is integrable the corresponding analytical solution might be unknown. Nevertheless, we will not consider such cases in this work. In our study we also want to track the evolution of the deviation vector $\boldsymbol{ w}(t)$ by solving Eq. (3.7). Indeed this can be done by using SIs since upon splitting the Hamiltonian into integrable parts we know analytically for each part the exact mapping $\boldsymbol{ x} (t) \rightarrow \boldsymbol{ x}(t+\tau)$, along with the mapping $\delta \boldsymbol{ x} (t) \rightarrow \delta \boldsymbol{ x}(t+\tau)$.

Let us outline the whole process by considering a general autonomous Hamiltonian function $H\left(\boldsymbol{ q}, \boldsymbol{ p}\right)$, which can be written as sum of $I$ integrable intermediate Hamiltonian functions $A_i$, i.e., $H = \sum_{i = 1}^I A_i$. This decomposition implies that the operator $e^{ \tau L_{A_iZ}}$ in the formal solution in Eq. (3.11) of each intermediate Hamiltonian function $A_i$ is known. A symplectic integration scheme to integrate the Hamilton equations of motion from $t_0$ to $t_0+\tau$ consists in approximating the action of $e^{\tau L_{HZ}}$ in Eq. (3.11) by a product of the operators $e^{\gamma_i \tau L_{A_iZ}}$ for a set of properly chosen coefficients $\gamma_i$. In our analysis we will call as number of steps of a particular SI the total number of successive application of each individual operator $e^{ \tau L_{A_iZ}}$. Further details about this class of integrators can be found in ^[80,81,82] and references therein. In what follows, we consider the most common cases where the Hamiltonian function can be split into two ($I = 2$) or three ($I = 3$) integrable parts. Let us remark here that, in general, the split $H = \sum_{i = 1}^I A_i$ is not necessarily unique (see Appendix B.1). Studying the efficiency and the stability of different SIs upon different choices of splitting the Hamiltonian is an interesting topic by itself, which, nevertheless, is beyond the scope of our work.

3.2.1.   Two part split

Let us consider that the Hamiltonian $H\left(\boldsymbol{ q}, \boldsymbol{ p}\right)$ can be separated into two integrable parts, namely $H = A(\boldsymbol{ q}, \boldsymbol{ p}) + B(\boldsymbol{ q}, \boldsymbol{ p})$. Then we can approximate the action of the operator $e^{\tau L_{HZ}} = e^{\tau(L_{AZ} + L_{BZ})}$ by the successive actions of products of the operators $e^{\tau L_{AZ}}$ and $e^{\tau L_{BZ}}$ ^{[80,81,82,83]}

for appropriate choices of the real coefficients $c_j, d_j$ with $j = 1, \dots, p$. Different choices of $p$ and coefficients $c_j$, $d_j$ lead to schemes of different accuracy. In Eq. (3.14) the integer $n$ is called the order of a symplectic integrator.

The Hamiltonian function of the 1D $\alpha$-FPUT system [Eq. (2.1)] can be split into two integrable parts $H_{\text{1F}} = A\left(\boldsymbol{ p}\right) + B\left(\boldsymbol{ q}\right)$, with each part possessing $N$ cyclic coordinates. The kinetic part

depends only on the generalized momenta, whilst the potential part

depends only on the generalized positions. This type of split is the most commonly used in the literature, therefore a large number of SIs have been developed for such splits. Below we briefly review the SIs used in our analysis, based mainly on results presented in ^[56].

Symplectic integrators of order two. These integrators constitute the most basic schemes we can develop from Eq. (3.14)

$LF$: The simplest example of Eq. (3.14) is the so-called Störmer-Verlet or leap-frog scheme (e.g., see [40,Sect. I.3.1] and [83]) having $3$ individual steps

where $a_1 = \frac{1}{2}$ and $b_1 = 1$.

$SABA_2/SBAB_2$: We consider the $SABA_2$ and the $SBAB_2$ SIs with $5$ individual steps

where $a_1 = \frac{1}{2} - \frac{1}{2\sqrt{3}}$, $a_2 = \frac{1}{\sqrt{3}}$ and $b_1 = \frac{1}{2}$, and

with $a_1 = \frac{1}{2}$, $b_1 = \frac{1}{6}$ and $b_2 = \frac{2}{3}$. These schemes were presented in ^[84], where they were named the (4, 2) methods, and also used in ^[55]. We note that the $SABA_2$ and $SBAB_2$ SIs (as well as other two part split SI schemes) have been introduced for Hamiltonian systems of the form $H = A + \varepsilon B$, with $\varepsilon$ being a small parameter. Both the $SABA_2$ and $SBAB_2$ integrators have only positive time steps and are characterized by an accuracy of order $\mathcal{O}(\tau^4\varepsilon + \tau^2\varepsilon^2)$ ^[55]. Although these integrators are particularly efficient for small perturbations ($\varepsilon \ll1$), they have also shown a very good performance in cases of $\varepsilon = 1$ (see e.g., ^[24]).

$ABA82$: In addition, we use in our analysis the SI

with $9$ individual steps ^[84,85], where the constants $a_i, b_i$ with $i = 1, 2, 3$ can be found in Table 2 of ^[85]. We note that the $ABA82$ method is called $SABA_4$ in ^[55].

Symplectic integrators of order four. The order of symmetric SIs can be increased by using a composition technique presented in ^[80]. According to that approach starting from a symmetric SI $S_{2n} (\tau)$ of order $2n\; (n\geq 1)$, we can construct a SI $S_{2n+2} (\tau)$ of order $2n+2$ as ^*

^*We adopt the notation of Iserles and Quispel ^[86].

$FR4$: Using the composition given in Eq. (3.21) for the $LF$ SI of Eq. (3.17) we construct a SI which we name $FR4$ ^[80,87] having $7$ individual steps

with coefficients $a_1 = \frac{1}{2(2 - 2^{1/3})}$, $a_2 = \frac{1 - 2^{1/3}}{2(2 - 2^{1/3})}$, $b_1 = \frac{1}{2 - 2^{1/3}}$ and $b_2 = -\frac{2^{1/3}}{2 - 2^{1/3}}$.

$SABA_2Y4/SBAB_2Y4$: Applying the composition given in Eq. (3.21) to the $SABA_2$ [Eq. (3.18)] and the $SBAB_2$ [Eq. (3.19)] integrators we obtain the fourth order SIs $SABA_2Y4$ and $SBAB_2Y4$ having $13$ individual steps. In particular, we get

with coefficients $d_0 = - \frac{2^{1/3}}{2 - 2^{1/3} }$, $d_1 = \frac{1}{ 2 - 2^{1/3} }$, $a_1 = \frac{1}{2} - \frac{1}{2\sqrt{3}}$, $a_2 = \frac{1}{\sqrt{3}}$, $b_1 = \frac{1}{2}$, and $a_0 = d_1a_1 + d_0 a_1$, and

with coefficients $d_0 = -\frac{2^{1/3}}{ 2 - 2^{1/3} }$, $d_1 = \frac{1}{ 2 - 2^{1/3} }$, $a_1 = \frac{1}{2}$, $b_1 = \frac{1}{6}$, $b_2 = \frac{2}{3}$ and $b_0 = d_1 b_1 + d_0 b_1$.

$ABA82Y4$: Using the composition given in Eq. (3.21) for the second order scheme $ABA82$ of Eq. (3.20) we obtain an integrator with $25$ individual steps having the form

where $d_0 = -\frac{2^{1/3}}{ 2 - 2^{1/3} }$, $d_1 = \frac{1}{ 2 - 2^{1/3} }$, while $a_i, b_i$ with $i = 1, 2, 3$ can be found in Table 2 of ^[85]. Here $a_0 = d_1a_1 + d_0 a_1$.

$SABA_2K/SBAB_2K$: As was explained in ^[55] the accuracy of the $SABA_n$ (and the $SBAB_n$) class of SIs can be improved by a corrector term $K = \{ B, \{B, A\}\}$, defined by two successive applications of Poisson brackets ($\{\cdot, \cdot\}$), if $K$ corresponds to a solvable Hamiltonian function. In that case, the second order integration schemes can be improved by the addition of two extra operators with negative time steps in the following way

with analogous result holding for the $SBAB_n$ scheme. By following this approach for the $SABA_2$ and $SBAB_2$ SIs [which are of the order $\mathcal{O}(\tau^4\varepsilon + \tau^4\varepsilon^2)$] we produce the fourth order SIs $SABA_2K$ and $SBAB_2K$, with $g = (2 - \sqrt{3})/24$ and $g = 1/72$ respectively. These new integration schemes are of the order $\mathcal{O}(\tau^4\varepsilon + \tau^4\varepsilon^2)$ ^[55].

$ABA864/ABAH864$: The fourth order SIs $ABA864$ and $ABAH864$ were proposed in ^[57,85]. They have respectively $15$ and $17$ individual steps and have the form

with coefficients $a_i, b_i\; i = 1, 2, 3, 4$ taken from Table 3 of ^[57], and

with coefficients $a_i, b _i, \; i = 1, 2, \dots, 5$ found in Table 4 of ^[57]. We note that both schemes were designed for near-integrable systems of the form $H = A + \varepsilon B$, with $\varepsilon$ being a small parameter, but the construction of $ABAH864$ was based on the assumption that the integration of the $B$ part cannot be done explicitly, but can be approximated by the action of some second order SI, since $B$ is expressed as the sum of two explicitly integrable parts, i.e., $B = B_1+B_2$. The $ABA864$ and $ABAH864$ SIs are of order four, but their construction satisfy several other conditions at higher orders, improving in this way their performance.

Symplectic integrators of order six. Applying the composition technique of Eq. (3.21) to the fourth order SIs $FR4$ [Eq. (3.22)], $SABA_2Y4$ [Eq. (3.23)], $SBAB_2Y4$ [Eq. (3.24)], $ABA82Y4$ [Eq. (3.25)], $SABA_2K$ and $ABA864$ [Eq. (3.27)], we construct the sixth order SIs $FR4Y6$, $SABA_2Y4Y6$, $SBAB_2Y4Y6$, $ABA82Y4Y6$, $SABA_2KY6$ and $ABA864Y6$ with $19$, $37$, $37$, $73$, 19 and $43$ individual steps.

In ^[80] a composition technique using fewer individual steps than the one obtained by the repeated application of Eq. (3.21) to SIs of order two was proposed, having the form

whose coefficients $w_i, \; i = 0, 1, 2, 3$ are given in Table 1 in ^[80] for the case of the so-called 'solution A' of that table. Here $S_2$ and $S_6$ respectively represent a second and a sixth order symmetric SI. Note that Eq. (3.29) corresponds to the composition scheme $s6odr6$ of ^[88]. Applying the composition given in Eq. (3.29) to the $SABA_2$ [Eq. (3.18)], the $SBAB_2$ [Eq. (3.19)] and the $ABA82$ [Eq. (3.20)] SIs we generate the order six schemes $SABA_2Y6$, $SBAB_2Y6$ and $ABA82Y6$ having $29$, $29$ and $57$ individual steps respectively.

We also consider in our study the composition scheme $s9odr6b$ of ^[88] which is based on $9$ successive applications of $S_2$

The values of $\delta _i, \; i = 1, 2, \dots, 5$ in Eq. (3.30) can be found in the Appendix of ^[88]. Furthermore, we also implement the composition method

of ^[89], which involves 11 applications of a second order SI $S_2$, whose coefficients $\gamma_i$, $i = 1, 2, \dots, 6$ are reported in Section 4.2 of ^[89]. Using the $SABA_2$ of Eq. (3.18) as $S_2$ in Eqs. (3.30) and (3.31), we respectively build two SIs of order six, namely the $s9SABA_26$ and the $s11SABA_26$ SIs with $37$ and $45$ individual steps. In addition, using the $ABA82$ integrator of Eq. (3.20) as $S_2$ in Eqs. (3.30) and (3.31) we construct two other order six SIs namely the $s9ABA82\_6$ and $s11ABA82\_6$ schemes with $73$ and $89$ individual steps respectively.

Runge-Kutta-Nyström methods: In addition, we consider in our analysis two SIs of order six, belonging in the category of the so-called Runge-Kutta-Nyström (RKN) methods (see e.g., ^{[40,90,44,71]} and references therein), which respectively have $21$ and $29$ individual steps

and

The values of the coefficients appearing in Eqs. (3.32) and (3.33) can be found in Table 3 of ^[90]. This class of integrators has, for example, been successfully implemented in a recent investigation of the chaotic behavior of the DNA molecule ^[91].

Symplectic integrators of order eight. Following ^[80] we can construct an eighth order SI $S_8$, starting from a second order one $S_2$, by using the composition

In our study we consider two sets of coefficients $w_i$, $i = 1, \ldots, 7$, and in particular the ones corresponding to the so-called 'solution A' and 'solution D' in Table 2 of ^[80]. Using in Eq. (3.34) the $SABA_2$ [Eq. (3.18)] SI as $S_2$ we construct the eighth order SIs $SABA_2Y8\_A$ (corresponding to 'solution A') and $SABA_2Y8\_D$ (corresponding to 'solution D') with $61$ individual steps each. In a similar way the use of $ABA82$ [Eq. (3.20)] in Eq. (3.34) generates the SIs $ABA82Y8\_A$ and $ABA82Y8\_D$ with $121$ individual steps each.

In addition, considering the composition scheme $s15odr8$ of ^[88], having 15 applications of $S_2$, we construct the eighth order SI $s15SABA_28$ [$s15ABA82\_8$] having $61$ [$121$] individual steps when $SABA_2$ of Eq. (3.18) [$ABA82$ of Eq. (3.20)] is used in the place of $S_2$.

Furthermore, implementing the composition technique $s19odr8$ presented in ^[89], which uses 19 applications of a second order SI $S_2$, we construct the SI $s19SABA_28$ [$s19ABA82\_8$] with $77$ [$153$] individual steps, when $SABA_2$ of Eq. (3.18) [$ABA82$ of Eq. (3.20)] is used in the place of $S_2$.

The SIs of this section will be implemented to numerically integrate the $\alpha$-FPUT model [Eq. (2.1)] since this Hamiltonian system can be split into two integrable parts $A(\boldsymbol{ p})$ and $B(\boldsymbol{ q})$. In Section B.1 of the Appendix the explicit forms of the operators $e^{\tau L_{AZ} }$ and $e^{\tau L_{BZ} }$ are given, along with the operator $e^{\tau L_{KZ} }$ of the corrector term used by the $SABA_2K$ and $SBAB_2K$ SIs. In addition, in Section B.4 of the Appendix the explicit forms of the tangent map method operators of some commonly used lattice systems, whose Hamiltonians can be split in two integrable parts, are also reported.

3.2.2.   Three part split

Let us now consider the case of a Hamiltonian function $H\left(\boldsymbol{ q}, \boldsymbol{p}\right)$ which can be separated into three integrable parts, namely $H\left(\boldsymbol{ q}, \boldsymbol{p}\right) = \mathcal{A}\left(\boldsymbol{ q}, \boldsymbol{p}\right) + \mathcal{B}\left(\boldsymbol{ q}, \boldsymbol{ p}\right) + \mathcal{C}\left(\boldsymbol{ q}, \boldsymbol{ p}\right)$. This for example could happen because the Hamiltonian function may not be split into two integrable parts, or to simplify the solution of one of the two components, $A$ or $B$, of the two part split schemes discussed in Section 3.2.1. In such cases we approximate the action of the operator $e^{\tau L_{HZ}}$ of Eq. (3.11) by the successive application of operators $e^{\tau L_{\mathcal{A}Z}}$, $e^{\tau L_{\mathcal{B}Z}}$ and $e^{\tau L_{\mathcal{C}Z}}$ i.e.,

for appropriate choices of the real coefficients $c_j$, $d_j$ and $e_j$ with $j = 1, \dots, p$. As in Eq. (3.14), in Eq. (3.35) the integer $n$ is the order of a symplectic integrator.

As examples of Hamiltonians which can be split in three integrable parts we mention the Hamiltonian function of a free rigid body ^[92] and the Hamiltonian functions of the 1D [Eq. (2.2)] and 2D [Eq. (2.3)] DDNLS models we consider in this work. For example the 1D DDNLS Hamiltonian of Eq. (2.2) can be split in the following three integrable Hamiltonians: a system of $N$ independent oscillators

where $J_i = (q_i^2 + p_i^2)/2$, $i = 1, \ldots, N$ are $N$ constants of motion, and the Hamiltonian functions of the $\boldsymbol{q}-$ and $\boldsymbol{p}-$hoppings

with each one of them having $N$ cyclic coordinates. The three part split of the 2D DDNLS of Eq. (2.3) can be found in Section B.3 of the Appendix [Eq. (A.17)].

We note that a rather thorough survey on three part split SIs can be found in ^[49,50]. We decided to include in our study a smaller number of schemes than the one presented in these works, focusing on the more efficient SIs. We briefly present these integrators below.

Symplectic integrators of order two. We first present the basic three part split scheme obtained by the application of the Störmer-Verlet/leap-frog method to three-part separable Hamiltonians. This scheme has $5$ individual steps and we call it $\mathcal{ABC}2$ ^[93]

where $a_1 = \frac{1}{2}$, $b_1 = \frac{1}{2}$ and $c_1 = 1$.

Symplectic integrators of order four. In order to built higher order three part split SIs we apply some composition techniques on the basic $\mathcal{ABC}2$ SI of Eq. (3.38).

$\mathcal{ABC}Y4$: Using the composition given in Eq. (3.21) for $n = 1$, we construct

with $d_0 = \frac{- 2^{1/3}}{2 - 2^{1/3}}$ and $d_1 = \frac{1}{2 - 2^{1/3}}$. This integrator has $13$ individual steps, it has been explicitly introduced in ^[93] and implemented in ^[50], where it was called $\mathcal{ABC}^{4}_{[Y]}$.

$\mathcal{ABC}S4$: Implementing a composition scheme which was introduced in ^[81] and studied in ^[88] (where it was named $s5odr4$) we obtain the SI

where $p_2 = \frac{1}{4 -4^{1/3}}$ and $1 - 4p_2 = \frac{-4^{1/3}}{4 - 4^{1/3}}$, having $21$ individual steps. This integrator was denoted as $\mathcal{ABC}^{4}_{[S]}$ in ^[50].

$SS864S$: Using the $ABAH864$ integrator of Eq. (3.28) where $B$ is considered to be the sum of functions $\mathcal{B}_1$ and $\mathcal{C}_1$ of Eq. (3.37), i.e. $B = \mathcal{B}_1+\mathcal{C}_1$, and its solution is approximated by the second order $SABA_2$ SI of Eq. (3.18), we construct a SI with 49 steps, which we call $SS864S$. This integrator has been implemented for the integration of the equations of motion of the 1D DDNLS system [Eq. (2.2)] in ^[49], where it was called $SS_{864}^4$.

Symplectic integrators of order six.

$\mathcal{ABC}Y4Y6/\mathcal{ABC}S4Y6$: Applying the composition technique of Eq. (3.21) to the fourth order SIs $\mathcal{ABC}Y4$ [Eq. (3.39)] and $\mathcal{ABC}S4$ [Eq. (3.40)], we respectively construct the schemes $\mathcal{ABC}Y4Y6$ and $\mathcal{ABC}S4Y6$ with $37$ and $49$ individual steps.

$\mathcal{ABC}Y6\_A$: Using the composition given in Eq. (3.29) we build a sixth order SI with $29$ individual steps, considering the integrator $\mathcal{ABC}2$ in the place of $S_2$

In particular, we consider in this construction the coefficients $w_i$, $i = 0, 1, 2, 3$, corresponding to the 'solution A' of Table 1 in ^[80]. Note that this SI has already been implemented in ^[49,50], where it was denoted as $\mathcal{ABC}^{6}_{[Y]}$.

$s9\mathcal{ABC}6$: Implementing the composition given in Eq. (3.30), with $\mathcal{ABC}2$ in the place of $S_2$, we get

which was referred to as $\mathcal{ABC}^{6}_{[KL]}$ in ^[49,50].

$s11\mathcal{ABC}6$: From the composition given in Eq. (3.31) we create the SI scheme

which has 45 individual steps. This integrator was referred to as $\mathcal{ABC}^{6}_{[SS]}$ in ^[49,50].

Symplectic integrators of order eight.

$\mathcal{ABC}Y8\_A/\mathcal{ABC}Y8\_D$: Based on the composition given in Eq. (3.34) we construct the SI

setting $\mathcal{ABC}2$ in the place of $S_2$. This SI has $61$ individual steps. Considering the 'solution A' of Table 2 in ^[80] for the coefficients $w_i$, $0\leq i \leq 7$, we obtain the $\mathcal{ABC}Y8\_A$ SI, while the use of 'solution D' of the same table leads to the construction of the SI $\mathcal{ABC}Y8\_D$.

$s17\mathcal{ABC}8$: Consider the composition method $s17odr8b$ of ^[88] we build the SI (referred to as $\mathcal{ABC}^{8}_{[KL]}$ in ^[49,50]) $s17\mathcal{ABC}8$ having $69$ individual steps

$s19\mathcal{ABC}8$: Finally, we also implement the composition $s19odr8b$ reported in [89,Eq. (13)] and construct the SI

which has $72$ individual steps. We note that this scheme corresponds to the SI $\mathcal{ABC}^{8}_{[SS]}$ considered in ^[49,50].

The explicit forms of the operators related to the three part split of the DDNLS Hamiltonians [Eqs. (2.2) and (2.3)] are given in Sections. B.2 and B.3 of the Appendix.

4.   Numerical results

We test the efficiency of the integrators presented in Section 3 by using them to follow the dynamical evolution of the $\alpha$-FPUT model [Eq. (2.1)], the 1D DDNLS system [Eq. (2.2)] and the 2D DDNLS Hamiltonian [Eq. (2.3)]. For each model, given an initial condition $\boldsymbol{X} (t_0) = \left(\boldsymbol{ x} (t_0), \delta\boldsymbol{ x} (t_0)\right)$ we compute the trajectory $\{\boldsymbol{x}(t_n)\}_{n\in \mathbb{N}}$ with $\boldsymbol{ x} (t) = \left(q_1 (t), q_2 (t), \ldots, q_N (t), p_1(t), p_2(t), \ldots, p_N (t) \right)$ and check the integrators' efficiency through their ability to correctly reproduce certain observables of the dynamics. We also follow the evolution of a small initial perturbation to that trajectory $\boldsymbol{ w}(t_0) = \delta \boldsymbol{x} (t_0) = \left(\delta q_1 (t_0), \delta q_2 (t_0), \ldots, \delta q_N (t_0), \delta p_1(t_0), \delta p_2(t_0), \ldots, \delta p_N (t_0) \right)$ and use it to compute the time evolution of the finite time mLE ^[35,36,37]

in order to characterize the regular or chaotic nature of the trajectory through the estimation of the most commonly used chaos indicator, the mLE $\chi$, which is defined as $\chi = \lim_{t\rightarrow +\infty} X_1(t)$. In Eq. (4.1) $\|\cdot \|$ is the usual Euclidian norm, while $\boldsymbol{ w}(t_0)$ and $\boldsymbol{ w}(t_0 + \tau)$ are respectively the deviation vectors at $t = t_0$ and $t_0 + \tau > t_0$. In the case of regular trajectories $X_1(t)$ tends to zero following the power law ^[37,94]

whilst it takes positive values for chaotic ones.

4.1.   The $\alpha$-Fermi-Pasta-Ulam-Tsingou model

We present here results on the computational efficiency of the symplectic and non-symplectic schemes of Section 3 for the case of the $\alpha$-FPUT chain [Eq. (2.1)]. As this system can be split into two integrable parts we will use for its study the two part split SIs of Section 3.2.1. In our investigation we consider a lattice of $N = 2^{10}$ sites with $\alpha = 0.25$ and integrate up to the final time $t_f = 10^6$ two sets of initial conditions:

● Case $\text{I}_{\text{F}}$: We excite all lattice sites by attributing to their position and momentum coordinates a randomly chosen value from a uniform distribution in the interval $[-1, 1]$. These values are rescaled to achieve a particular energy density, namely $H_{1F}/N = 0.1$.

● Case $\text{II}_{\text{F}}$: Same as in case $\text{I}_{\text{F}}$, but for $H_{1F}/N = 0.05$.

We consider these two initial conditions in an attempt to investigate the potential dependence of the performance of the tested integrators on initial conditions [73,Sec. 8.3]. Since we have chosen non-localized initial conditions, we also use an initial normalized deviation vector $\boldsymbol{ w} (t)$ whose components are randomly selected from a uniform distribution in the interval $[-1, 1]$.

To evaluate the performance of each integrator we investigate how accurately it follows the considered trajectories by checking the numerical constancy of the energy integral of motion, i.e., the value of $H_\text{1F}$ [Eq. (2.1)]. This is done by registering the time evolution of the relative energy error

at each time step. In our analysis we consider two energy error thresholds $E_r \approx 10^{-5}$ and $E_{r} \approx 10^{-9}$. The former, $E_r \approx 10^{-5}$, is typically considered to be a good accuracy in many studies in the field of lattice dynamics, like for example in investigations of the DKG and the DDNLS models, as well as in systems of coupled rotors (see for example ^{[24,25,26,28,29,30]}). In some cases, e.g., for very small values of conserved quantities, one may desire more accurate computations. Then $E_r \approx 10^{-9}$ is a more appropriate accuracy level. In addition, in order to check whether the variational equations are properly evolved, we compute the finite time mLE $X_1(t)$ [Eq. (4.1)].

In Figure 1 we show the time evolution of the relative energy error $E_r(t)$ [panels (a) and (d)], the finite time mLE $X_1(t)$ [panels (b) and (e)], and the required CPU time $T_C$ [panels (c) and (f)], for cases $\text{I}_{\text{F}}$ and $\text{II}_{\text{F}}$ respectively, when the following four integrators were used: the fourth order SI $ABA864$ (blue curves), the sixth order SI $SABA_2Y6$ (red curves), the $DOP853$ scheme (green curves) and the $TIDES$ package (brown curves). These results are indicative of our analysis as in our study we considered in total 37 different integrators (see Tables 1 and 2). In Figure 1 the integration time steps $\tau$ of the SIs (reported in Tables 1 and 2) were appropriately chosen in order to achieve $E_r \approx 10^{-9}$, while for the $DOP853$ algorithm and the $TIDES$ package $E_r(t)$ eventually grows in time as a power law [Figure 1(a),(d)]. Nevertheless, all schemes succeed in capturing correctly the chaotic nature of the dynamics as they do not present any noticeable difference in the computation of the finite time mLE $X_1$ in Figure 1(b), (e). For both sets of initial conditions $X_1$ eventually saturates to a constant positive value indicating that both trajectories are chaotic. The CPU time $T_C$ needed for the integration of the equations of motion and the variational equations are reported in Figure 1(c), (f). From these plots we see that the SIs need less computational time to perform the simulations than the $DOP853$ and $TIDES$ schemes.

In Table 1 (Table 2) we present information on the performance of all considered integration schemes for the initial condition of case $\text{I}_{\text{F}}$ (case $\text{II}_{\text{F}}$). From the results of these tables we see that the performance and ranking (according to $T_C$) of the integrators do not practically depend on the considered initial condition. It is worth noting that although the non-symplectic schemes manage to achieve better accuracies than the symplectic ones, as their $E_r$ values are smaller [Figure 1(a),(d)], their implementation is not recommended for the long time evolution of the Hamiltonian system, because they require more CPU time and eventually their $E_r$ values will increase above the bounded $E_r$ values obtained by the symplectic schemes.

From the results of Tables 1 and 2 we see that the best performing integrators are the fourth order SIs $ABA864$ and $ABAH864$ for $E_r\approx 10^{-5}$, and the sixth order SIs $SRKN_{14}^a$ and $SRKN_{11}^b$ for $E_r \approx 10^{-9}$. We note that the best SI for $E_r\approx 10^{-5}$, the $ABA864$ scheme, shows a quite good behavior also for $E_r \approx 10^{-9}$, making this integrator a valuable numerical tool for dynamical studies of multidimensional Hamiltonian systems. We remark that the eighth order SIs we implemented to achieve the moderate accuracy level $E_r\approx 10^{-5}$ exhibited an unstable behavior failing to keep their $E_r$ values bounded. A similar behavior was also observed for the two RKN schemes $SRKN_{14}^a$, $SRKN_{11}^b$ when they were used to obtain $E_r\approx 10^{-5}$. Thus, the higher order SIs are best suited for more accurate computations. It is also worth mentioning here that the ranking presented in Tables 1 and 2 is indicative of the performance of the various SIs in the sense that small changes in the implementation (e.g., a change in the last digit of the used $\tau$ value) of integrators with similar behaviors (i.e., similar $T_C$ values) could interchange their ranking positions without any noticeable difference in the produced results.

4.2.   The 1D disordered discrete nonlinear Schrödinger equation system

We investigate the performance of various integrators of Section 3 for the 1D DDNLS system [{Eq. (2.2)] by considering a lattice of $N = 2^{10}$ sites and integrating two sets of initial conditions (for the same reason we did that for the $\alpha$-FPUT system) up to the final time $t_f = 10^6$. We note that, as was already mentioned in Section 3.2.2, this model can be split into three integrable parts, so we will implement the SIs presented in that section. In particular, we consider the following two cases of initial conditions:

● Case $\text{I}_{\text{1D}}$: We initially excite $21$ central sites by attributing to each one of them the same constant norm $s_j = (q_j^2 + p_j^2)/2 = 1$, $1\leq i \leq N$, for $W = 3.5$ and $\beta = 0.62$. This choice sets the total norm $S_{\text{1D}} = 21$. The random disorder parameters $\epsilon _i$, $1\leq i \leq N$, are chosen so that the total energy is $H_{1D} \approx 0.0212$.

● Case $\text{II}_{\text{1D}}$: Similar set of initial conditions as in case $\text{I}_{\text{1D}}$ but for $W = 3$, $\beta = 0.03$. The random disorder parameters $\epsilon _i$, $1\leq i \leq N$, are chosen such that $H_{\text{1D}}\approx 3.4444$.

We note that cases $\text{I}_{\text{1D}}$ and $\text{II}_{\text{1D}}$ have been studied in ^[33] and respectively correspond to the so-called 'strong chaos' and 'weak chaos' dynamical regimes of this model. As initial normalized deviation vector we use a vector having non-zero coordinates only at the central site of the lattice, while its remaining elements are set to zero.

To evaluate the performance of each implemented integrator we check if the obtained trajectory correctly captures the statistical behavior of the normalized norm density distribution $\zeta _j = s_j/S_{\text{1D}}$, $1\leq j\leq N$, by computing the distribution's second moment

where $\overline{j} = \sum_{j = 1}^N j \zeta _j$ is the position of the center of the distribution ^{[22,24,26,29,33,49,50]}. We also check how accurately the values of the system's two conserved quantities, i.e., its total energy $H_{\text{1D}}$ [Eq. (2.2)] and norm $S_{\text{1D}}$ [Eq. (2.4)], are kept constant throughout the integration by evaluating the relative energy error $E_r (t)$ [similarly to Eq. (4.3)] and the relative norm error

In addition, we compute the finite time mLE $X_1(t)$ [Eq. (4.1)] in order to characterize the system's chaoticity and check the proper integration of the variational equations.

We consider several three part split SIs, which we divide into two groups: (i) those integrators with order $n \leq 6$, which we implement in order to achieve an accuracy of $E_r \approx 10^{-5}$, and (ii) SIs of order eight used for $E_r \approx 10^{-9}$. In addition, the two best performing integrators of the first group are also included in the second group. We do not use higher order SIs for obtaining the accuracy level of $E_r \approx 10^{-5}$ because, as was shown in ^[49] and also discussed in Section 4.1, usually this task requires large integration time steps, which typically make the integrators unstable. Moreover, increasing the order $n$ of SIs beyond eight does not improve significantly the performance of the symplectic schemes for high precision ($E_r\approx 10^{-9}$) simulations ^[49]. Therefore we do not consider such integrators in our study.

In Figure 2 we show the time evolution of the relative energy error $E_r(t)$ [panel (a)], the relative norm error $S_r(t)$ [panel (b)], the second moment $m_2(t)$ [panel (d)], as well as the norm density distribution $\zeta_j$ at time $t_f\approx 10^6$ [panel (c)] for case $\text{I}_{\text{1D}}$ (we note that analogous results were also obtained for case $\text{II}_{\text{1D}}$, although we do not report them here). These results are obtained by the implementation of the the second order SI $\mathcal{ABC}2$ (red curves), the fourth order SI $\mathcal{ABC}Y4$ (blue curves), and the non-symplectic schemes $DOP853$ (green curves) and $TIDES$ (brown curves). The results of Figure 2 are indicative of the results obtained by the integrators listed in Table 3. The integration time step $\tau$ of the SIs was selected so that the relative energy error is kept at $E_r \approx 10^{-5}$ [Figure 2(a)]. From the results of Figure 2(b) we see that the SIs do not keep $S_r$ constant. Nevertheless, our results show that we lose no more than two orders of precision (in the worst case of the $\mathcal{ABC}2$ scheme) during the whole integration. On the other hand, both the relative energy [$E_r(t)$] and norm [$S_r(t)$] errors of the $TIDES$ and $DOP853$ integrators increase in time, with the $TIDES$ scheme behaving better than the $DOP853$ one. Figure 2(c), (d) show that all integrators correctly reproduce the dynamics of the system, as all of them practically produce the same norm density distribution at $t_f = 10^6$ [Figure 2(c)] and the same evolution of the $m_2(t)$ [Figure 2(d)]. We note that $m_2$ increases by following a power law $m_2 \propto t^\alpha$ with $\alpha = 1/2$, as was expected for the strong chaos dynamical regime (see for example ^[33] and references therein).

In Figure 3 we show the evolution of the finite time mLE $X_1$ [panels (a) and (c)] and the required CPU time $T_C$ [panels (b) and (d)] for the integration of the Hamilton equations of motion and the variational equations for cases $\text{I}_{\text{1D}}$ [panels (a) and (b)] and $\text{II}_{\text{1D}}$ [panels (c) and (d)] obtained by using the same integrators of Figure 2. Again the results obtained by these integrators are practically the same for both sets of initial conditions, reproducing the tendency of the finite time mLE to asymptotically decrease according to the power law $X_1(t) \propto t^{\alpha_L}$ with $\alpha _L \approx -0.3$ (case $\text{I}_{1D}$) and $\alpha _L \approx -0.25$ (case $\text{II}_{1D}$), in accordance to the results of ^[32,33].

We now check the efficiency of the used symplectic and non-symplectic methods by comparing the CPU time $T_C$ they require to carry out the simulations. These results are reported in Table 3 for case $\text{I}_{\text{1D}}$ and in Table 4 for case $\text{II}_{\text{1D}}$. These tables show that the comparative performance of the integrators does not depend on the chosen initial condition, as the ranking of the schemes is practically the same in both tables. As in the case of the $\alpha$-FPUT model, the $DOP853$ and $TIDES$ integrators required, in general, more CPU time than the SIs, although they produced more accurate results (smaller $E_r$ and $S_r$ values) at least up to $t_f = 10^{6}$, with $TIDES$ being more precise. The integrators exhibiting the best performance for $E_r \approx 10^{-5}$ are the sixth order SIs $s11\mathcal{ABC}6$ and $s9\mathcal{ABC}6$, while for $E_r \approx 10^{-9}$ we have the eighth order SIs $s19\mathcal{ABC}8$ and $s17\mathcal{ABC}8$, with the $s11\mathcal{ABC}6$ scheme performing quite well also in this case.

4.3.   The 2D disordered discrete nonlinear Schrödinger equation system

We now investigate the performance of the integrators used in Section 4.2 for the computationally much more difficult case of the 2D DDNLS lattice of Eq. (2.3), as its Hamiltonian function can also be split into three integrable parts. In order to test the performance of the various schemes we consider a lattice with $N\times M = 200\times 200$ sites, resulting to a system of $4\times 40\, 000 = 160\, 000$ ODEs (equations of motion and variational equations). The numerical integration of this huge number of ODEs is a very demanding computational task. For this reason we integrate this model only up to a final time $t_f = 10^5$, instead of the $t_f = 10^6$ used for the $\alpha$-FPUT and the 1D DDNLS systems. It is worth noting that due to the computational difficulty of the problem very few numerical results for the 2D DDNLS system exist in the literature (e.g., ^[39,65]). We consider again two sets of initial conditions:

● Case $\text{I}_{\text{2D}}$: We initially excite $7\times 7$ central sites attributing to each one of them the same norm $s_{i, j} = (q_{i, j} ^2 + p_{i, j}^2)/2 = 1/6$ so that the total norm is $S_{\text{2D}} = 49/6$, for $W = 15$ and $\beta = 6$. The disorder parameters $\epsilon_{i, j}$, $1\leq i \leq N$, $1\leq j \leq M$, are chosen so that the initial total energy is $H_{\text{2D}} \approx 1.96$.

● Case $\text{II}_{\text{2D}}$: We initially excite a single central site of the lattice with a total norm $S_{\text{2D}} = 1$, i.e., $\zeta _{100,100} = 1$, for $W = 16$, $\beta = 1.25$ and $H_{\text{2D}} = 0.625$.

The initial normalized deviation vector considered in our simulations has random non-zero values only at the $7\times 7$ initially excited sites for case $\text{I}_{\text{2D}}$, and only at site $i = 100$, $j = 100$ for case $\text{II}_{\text{2D}}$. In both cases, all others elements of the vectors are initially set to zero. Both considered cases belong to a Gibbsian regime where the thermalization processes are well defined by Gibbs ensembles ^[52,97]. Therefore, we expect a subdiffusive spreading of the initial excitations to take place for both cases, although their initial conditions are significantly different.

As was done in the case of the 1D DDNLS system (Section 4.2), in order to evaluate the performance of the used integrators we follow the time evolution of the normalized norm density distribution $\zeta _{i, j} = s_{i, j} /S_{\text{2D}}$, $1\leq i \leq N$, $1\leq j \leq M$ and compute the related second moment $m_2$ and participation number $P$

where $(\overline{i}, \overline{j})^T = \sum_{i, j} (i, j)^T \zeta _{i, j}$ is the mean position of the norm density distribution. We also evaluate the relative energy [$E_r (t)$] and norm [$S_r (t)$] errors and compute the finite time mLE $X_1(t)$.

In Figure 4 we present results obtained for case $\text{I}_{\text{2D}}$ by the four best performing SIs (see Table 5), the $s11\mathcal{ABC}6$ (red curves), $s9\mathcal{ABC}6$ (blue curves), $\mathcal{ABC}Y6\_A$ (green curves) and $\mathcal{ABC}Y4$ (brown curves) schemes, along with the $DOP853$ integrator (grey curves). The integration time steps $\tau$ of the SIs were adjusted in order to obtain an accuracy of $E_r\approx 10^{-5}$ [Figure 4(a)], while results for the conservation of the second integral of motion, i.e., the system's total norm, are shown in Figure 4(b). We see that for all SIs the $S_r$ values increase slowly, remaining always below $S_r \approx 10^{-4}$, which indicates a good conservation of the system's norm. As in the case of the 1D DDNLS system, the $E_r$ and $S_r$ values obtained by the $DOP853$ integrator increase, although the choice of $\delta = 10^{-16}$ again ensures high precision computations.

The norm density distributions at the final integration time $t_f = 10^5$ along the axis $i = 100$ [Figure 4(c)] and $j = 100$ [Figure 4(d)] obtained by the various integrators practically overlap indicating the ability of all numerical schemes to correctly capture the system's dynamics, as well as the fact that the initial excitations expand along all directions of the 2D lattice. From Figure 4(e) [Figure 4(f)] we see that $m_2 (t)$ [$P(t)$] is increasing according to the power law $m_2 = t^{1/3}$ [$P = t^{1/3}$] as expected from the analysis presented in ^[30], indicating that the 2D lattice is being thermalized. The results of Figure 4(e), (f) provide additional numerical evidences that all numerical methods reproduce correctly the dynamics. This is also seen by the similar behavior of the finite time mLE curves in Figure 4(g). From the results of this figure we see that $X_1$ exhibits a tendency to decrease following a completely different decay from the $X_1\propto t^{-1}$ power law observed for regular motion. This behavior was also observed for the 2D DKG model ^[56], as well as for the 1D DKG and DDNLS systems in ^[32,33] where a power law $X_1(t) \propto t^{\alpha_L}$ with $\alpha _L \approx -0.25$ and $\alpha _L \approx -0.3$ for, respectively, the weak and strong chaos dynamical regimes was established. Further investigations of the behavior of the finite mLE in 2D disordered systems are required in order to determine a potentially global behavior of $X_1$, since here and in ^[56] only some isolated cases were discussed. Such studies will require the statistical analysis of results obtained for many different disorder realizations, parameter sets and initial conditions. Thus, the utilization of efficient and accurate numerical integrators, like the ones presented in this study, will be of utmost importance for the realization of this goal.

From Tables 5 and 6, where the CPU times $T_C$ required by the tested integrators are reported, we see that, as in the case of the 1D DDNLS model, the SIs $s11\mathcal{ABC}6$ and $s9\mathcal{ABC}6$ have the best performance for $E_r\approx 10^{-5}$ and the SIs $s19\mathcal{ABC}8$ and $s17\mathcal{ABC}8$ for $E_r\approx 10^{-9}$.

5.   Conclusions

In this work we carried out a methodical and detailed analysis of the performance of several symplectic and non-symplectic integrators, which were used to integrate the equations of motion and the variational equations of some important many-body classical Hamiltonian systems in one and two spatial dimensions: the $\alpha$-FPUT chain, as well as the 1D and 2D DDNLS models. In the case of the $\alpha$-FPUT system we used two part split SIs, while for the integration of the DDNLS models we implemented several three part split SIs. In order to evaluate the efficiency of all these integrators we evolved in time different sets of initial conditions and evaluated quantities related to (a) the dynamical evolution of the studied systems (e.g., the second moment of norm density distributions for the DDNLS models), (b) the quantification of the systems' chaotic behavior (i.e., the finite time mLE), and (c) the accurate computation of the systems' integrals of motion (relative energy and norm errors), along with the CPU times needed to perform the simulations.

For the $\alpha$-FPUT system several two part split SIs showed very good performances, among which we mention the $ABA864$ and $ABAH864$ SIs of order four to be the best schemes for moderate energy accuracies ($E_r \approx 10^{-5}$), while the $SRKN^a_{14}$ and $SRKN^b_{11}$ SIs of order six were the best integration schemes for higher accuracies ($E_r \approx 10^{-9}$). In particular, the $ABA864$ scheme appears to be an efficient, general choice as it showed a quite good behavior also for $E_r \approx 10^{-9}$. Concerning the 1D and the 2D DDNLS models our simulations showed that the SIs $s9\mathcal{ABC}6$ and $s11\mathcal{ABC}6$ (order six), along with the SIs $s17\mathcal{ABC}8$ and $s19\mathcal{ABC}8$ (order eight) are the best integrators for moderate ($E_r \approx 10^{-5}$) and high ($E_r \approx 10^{-9}$) accuracy levels respectively.

The $DOP853$ and $TIDES$ non-symplectic integrators required, in general, much longer CPU times to carry out the simulations, although they produced more accurate results (i.e., smaller $E_r$ and $S_r$ values) than the symplectic schemes. Apart from the drawback of the high CPU times, the fact that $E_r$ (and $S_r$) values exhibit a constant increase in time signifies that such schemes should not be preferred over SIs when very long time simulations are needed.

It is worth noting that two part split SIs of order six and higher often do now not produce reliable results for relative low energy accuracies like $E_r \approx 10^{-5}$ for the $\alpha$-FPUT system (similar behaviors were reported in ^[56] for the DKG model). This happens because the required integration time step $\tau$ needed to keep the relative energy error at $E_r \approx 10^{-5}$ is typically large, resulting to an unstable behavior of the integrator i.e., the produced $E_r$ values do not remain bounded. Thus, SIs of order $n \geq 6$ are more suitable for calculations that require higher accuracies (e.g., $E_r \approx 10^{-9}$ or lower).

We note that we presented here a detailed comparison of the performance of several two and three part split SIs for the integration of the variational equations through the tangent map method and consequently for the computation of a chaos indicator (the mLE), generalizing, and completing in some sense, some sporadic previous investigation of the subject ^{[56,58,59,60]}, which were only focused on two part split SIs.

We hope that the clear description of the construction of several two and three part split SIs in Section 3, along with the explicit presentation in the Appendix of the related differential operators for many commonly used classical, many-body Hamiltonians will be useful to researchers working on lattice dynamics. The numerical techniques presented here can be used for the computation of several chaos indicators, apart from the mLE (e.g., the SALI and the GALI methods ^[98]) and for the dynamical study of various lattice models, like for example of arrays of Josephson junctions in regimes of weak non-integrability ^[53], granular chains ^[99] and DNA models ^[91], to name a few.

Appendix

A.   Explicit forms of tangent map method operators

We present here the exact expressions of the operators needed by the various SIs we implemented in our study to simultaneously solve the Hamilton equations of motion and the variational equations, or in order words to solve the system of Eq. (3.7)

A.1.   The $\alpha$-Fermi-Pasta-Ulam-Tsingou model

The Hamiltonian of the $\alpha$-FPUT chain [Eq. (2.1)] can be split into two integrable parts as

As we have already stated, the split into two integrable parts is not necessarily unique. In this particular case another possible choice of integrable splits for the $\alpha$-FPUT chain is to group together the quadratic terms of the Hamiltonian [i.e., $A({\bf p}, {\bf q}) = \sum_{ i = 0}^N\ \frac{p_i^2}{2} + \frac{1}{2}(q_{i+1} - q_{ i})^2$] and keep separately the nonlinear terms [i.e., $B({\bf q}) = \sum_{ i = 0}^N\ \frac{\alpha}{3}(q_{i+1} - q_{ i})^3$]. The set of equations of motion and variational equations for the Hamiltonian function $A\left(\boldsymbol{ p}\right)$ is

and the corresponding operator $e^{\tau L_{AZ}}$, which propagates the values of $q_i$, $p_i$, $\delta q_i$ and $\delta p_i$ for $\tau$ time units in the future, obtaining $q_i'$, $p_i'$, $\delta q_i'$ and $\delta p_i'$, takes the form

In a similar way for the $B\left(\boldsymbol{ q}\right)$ Hamiltonian of Eq. (A.2) we get

and

According to Eq. (3.26) the accuracy of the $SABA_n$ and $SBAB_n$ integrators can be improved by using a corrector Hamiltonian $K$ ^[55]. In the case of a separable Hamiltonian $H(\boldsymbol{q}, \boldsymbol{p}) = A(\boldsymbol{p})+B(\boldsymbol{q})$ with $A\left(\boldsymbol{ p}\right) = \sum_{i = 1}^N p_i^2/2$, the corrector $K$ becomes

For the $\alpha$-FPUT chain, the corrector Hamiltonian $K$ is

As the equations of motion and variational equations associated to the corrector Hamiltonian $K$ are cumbersome, we report here only the form of the operator $e^{ \tau L_{KZ}}$

where

We did not specify the range of index $i$ in Eqs. (A.5), (A.6) and (A.9) intentionally, because it depends on the type of the used boundary conditions. In particular, the expression of Eqs. (A.5), (A.6) and (A.9) are accurate for the case of periodic boundary conditions, i.e., $q_0 = q_N$, $p_0 = p_N$, $\delta q_0 = \delta q_N$, $\delta p_0 = \delta p_N$, $q_{N+1} = q_1$, $p_{N+1} = p_1$, $\delta q_{N+1} = \delta q_1$, $\delta p_{N+1} = \delta p_1$. In the case of fixed boundary conditions we considered in our numerical simulations, some adjustments have to be done for the $i = 1$ and $i = N$ equations, like the ones presented in the Appendix of ^[56] where the operators for the 1D and 2D DKG models were reported (with the exception of the corrector term).

For completeness sake in Section B.4 we provide the explicit expression of the $e^{ \tau L_{BZ}}$ and $e^{ \tau L_{KZ}}$ operators for some commonly used Hamiltonians, which can be split into two integrable parts, one of which is the usual kinetic energy $A\left(\boldsymbol{ p}\right) = \sum_{i = 1}^N p_i^2/2$.

A.2.   The 1D disordered discrete nonlinear Schrödinger equation

Here we focus on the 1D DDNLS system, whose Hamiltonian [Eq. (2.2)] can be split into three integrable parts as

The set of equations of motion and variational equations associated with the Hamiltonian function $\mathcal{A}_1$ is

with $\theta_{i} = \epsilon _{i} + \beta (q_{i} ^ 2 + p_{i} ^2)/2$ for $i = 1, 2, \ldots, N$ being constants of the motion. The corresponding operator $e^{\tau L_{\mathcal{A}_1Z}}$ takes the form

with $J_i \neq 0$ and

We note that in the special case of $J_i = 0$ we have $q_i = p_i = 0$. Then the system of Eq. (A.12) takes the simple form $\dot{q_i} = 0$, $\dot{p_i} = 0$, $\dot{\delta q}_{i} = \epsilon_i \delta p_i$, $\dot{\delta p}_{i} = - \epsilon_i \delta q_i$, leading to $q_i' = q_i$, $p_i' = p_i$, $\delta q_i' = \delta q_i \cos (\epsilon_i \tau) + \delta p_i\sin (\epsilon_i \tau)$, $\delta p_i' = \delta p_i \cos (\epsilon_i \tau) - \delta q_i\sin (\epsilon_i \tau)$.

The set of equations of motion and variational equations associated to the intermediate Hamiltonian functions $\mathcal{B}_1$ and $\mathcal{C}_1$ are respectively

These yield to the operators $e^{L_{\mathcal{B}_1Z}}$ and $e^{L_{\mathcal{C}_1Z}}$ given by

As in the case of the $\alpha$-FPUT model, Eqs. (A.15) and (A.16) correspond to the case of periodic boundary conditions and adjustments similar to the ones presented in the Appendix of ^[56] for the 1D DKG model, should be implemented when fixed boundary conditions are imposed.

A.3.   The 2D disordered discrete Nonlinear Schrödinger equation

The Hamiltonian $H_{\text{2D}}$ of Eq. (2.3) can be split into three integrable parts $\mathcal{A}_2$, $\mathcal{B}_2$ and $\mathcal{C}_2$ as

The equations of motion and the variational equations associated with the $\mathcal{A}_2$ Hamiltonian are

for $1 \leq i \leq N$, $1 \leq j \leq M$, with $\theta_{i, j} = \epsilon _{i, j} + \beta (q_{i, j} ^ 2 + p_{i, j} ^2)/2$ being constants of motion for the Hamiltonian $\mathcal{A}_2$. Then, the operator $e^{\tau L_{\mathcal{A}_2Z}}$ is

with $J_{i, j}\neq 0$ and

Again, in the special case of $J_{i, j} = 0$, the system of Eq. (A.18) takes the form $\dot{q_{i, j}} = 0$, $\dot{p_{i, j}} = 0$, $\dot{\delta q}_{i, j} = \epsilon_{i, j} \delta p_{i, j}$, $\dot{\delta p}_{i, j} = - \epsilon_{i, j} \delta q_{i, j}$, leading to $q_{i, j}' = q_{i, j}$, $p_{i, j}' = p_{i, j}$, $\delta q_{i, j}' = \delta q_{i, j} \cos (\epsilon_{i, j}\tau) + \delta p_{i, j}\sin (\epsilon_{i, j}\tau)$, $\delta p_{i, j}' = \delta p_{i, j}\cos (\epsilon_{i, j}\tau) - \delta q_{i, j}\sin (\epsilon_{i, j}\tau)$.

The equations of motion and the variational equations for Hamiltonians $\mathcal{B}_2$ and $\mathcal{C}_2$ are

and

while the corresponding operators $e^{L_{\mathcal{B}_2Z}}$ and $e^{L_{\mathcal{C}_2Z}}$ are respectively

and

Here again Eqs. (A.21)–(A.24) correspond to the case of periodic boundary conditions. For fixed boundary conditions adjustments similar to the ones report in the Appendix of ^[56] for the 2D DKG lattice should be performed.

A.4.   Other Hamiltonian models which can be split into two integrable parts

We present here the exact expressions of the tangent map operators needed in symplectic integration schemes which can be used to numerically integrate some important models in the field of classical many-body systems: the $\beta$-FPUT chain, the KG model, and the classical XY model (a JJC system) ^{[53,100,101,102]}. Similarly to the $\alpha$-FPUT chain of Eq. (2.1), the Hamiltonians $H(\boldsymbol{q}, \boldsymbol{p})$ of each of these systems can be split as

with appropriately defined potential terms $B\left(\boldsymbol{ q}\right)$:

where $\beta$, $k$ and $E_J$ are real coefficients. Obviously for all these systems the operator $e^{\tau L_{AZ}}$ of the kinetic energy part is the same as for the $\alpha$-FPUT chain in Eq. (A.4). Thus, we report below only the expressions of the operators $e^{\tau L_{BZ}}$ and $e^{\tau L_{KZ}}$ when periodic boundary conditions are imposed.

1. The $\beta$-Fermi-Pasta-Ulam-Tsingou chain

The operator $e^{\tau L_{BZ}}$ of the $\beta$-FPUT chain of Eq. (A.26) is

The corrector Hamiltonian $K$ of Eq. (A.7) becomes

while the corresponding operator is

with

2. The 1D Klein-Gordon chain model

The operator $e^{\tau L_{BZ}}$ of the Klein-Gordon chain [Eq. (A.27)] is

The corresponding corrector Hamiltonian $K$ [Eq.(A.7)] is written as

while $e^{\tau L_{KZ}}$ takes the form

with

3. The XY model of a Josephson junctions array

The operator $e^{\tau L_{BZ}}$ for the potential of Eq. (A.28) is

In this case the corrector Hamiltonian $K$ of Eq.(A.7) becomes

and the operator $e^{\tau L_{KZ}}$ is given by the following set of equations

with

Acknowledgments

We thank S. Flach for useful discussions. C.D. and M.T. acknowledge financial support from the IBS (Project Code No. IBS-R024-D1). Ch.S. and B.M.M. were supported by the National Research Foundation of South Africa (Incentive Funding for Rated Researchers, IFFR and Competitive Programme for Rated Researchers, CPRR). B.M.M. would also like to thank the High Performance Computing facility of the University of Cape Town (http://hpc.uct.ac.za) and the Center for High Performance Computing (https://www.chpc.ac.za) for the provided computational resources needed for the study of the two DDNLS models, as well as their user-support teams for their help on many practical issues.

Conflict of interest

All authors declare no conflict of interest in this paper.