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

We investigate the computational performance of various numerical methods for the integration of the equations of motion and the variational equations for some typical classical many-body models of condensed matter physics: the Fermi-Pasta-Ulam-Tsingou (FPUT) chain and the one- and two-dimensional disordered, discrete nonlinear Schr\"odinger equations (DDNLS). In our analysis we consider methods based on Taylor series expansion, Runge-Kutta discretization and symplectic transformations. The latter have the ability to preserve particular properties of Hamiltonian models such as the constancy of the system's total energy. We perform extensive numerical simulations for several initial conditions of the studied models and compare the numerical efficiency of the used integrators by testing their ability to accurately reproduce characteristics of the systems' dynamics and quantify their chaoticity through the computation of the maximum Lyapunov exponent. We also report the expressions of the implemented symplectic schemes and provide the explicit forms of the used differential operators. Among the tested numerical schemes the symplectic integrators $ABA864$ and $s11SABA_26$ exhibit the best performance, respectively for moderate and high accuracy levels in the case of the FPUT chain, while for the DDNLS models $s9\mathcal{ABC}6$ and $s11\mathcal{ABC}6$ (moderate accuracy), along with $s17\mathcal{ABC}8$ and $s19\mathcal{ABC}8$ (high accuracy) proved to be the most efficient schemes.

maximum Lyapunov exponent (mLE) [33][34][35]. 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 [36,37]. Another category of integrators is the so-called symplectic integrators (SIs), which were particularly designed for Hamiltonian systems (see e.g. [38][39][40][41][42] and references therein). The main characteristic of SIs is that they conserve the value of the Hamiltonian function (which is usually referred to as the system's energy) by keeping its error bounded as time evolves. This 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 FPUT models, the DKG, the DDNLS and systems of coupled rotors (see e.g. [20-27, 31, 32, 43-51]). In these studies SIs belonging to the so-called SABA family [52] were mostly used, since the FPUT, the KG 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 [32,53] the ABA864 [54] 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 [26]. In [46,47] it was shown that for the DDNLS model, splitting 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 [55][56][57].
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 [53].
The paper is organized as follows. In Sec. II we introduce the models considered in our study. In Sec. III we introduce in detail the symplectic and non-symplectic schemes implemented in our investigations. In Sec. IV we present our numerical results and report on the efficiency of each studied scheme. In Sec. V we review our findings and discuss the outcomes of our study. Finally, in Appendix A 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 β-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.

II. 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.

A. The α-Fermi-Pasta-Ulam-Tsingou model
The one-dimensional 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 non-linear lattices. In the literature there exist two types of the FPUT model the so-called α− and β−FPUT systems. In our study we consider the α−FPUT model whose Hamiltonian function reads In this expression q i and p i are respectively the generalized coordinates and momenta of the i lattice site and α 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 1F .
B. 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 [49], as well as investigate, at a first approximation, several physical processes (e.g. electron tight binding systems [37]). 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 i are uncorrelated numbers uniformly distributed in the interval − W 2 , W 2 . The real, positive numbers W and β 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 [27,58]. Its Hamiltonian function is where q i,j and p i,j are respectively the generalized positions and momenta at site (i, j) and i,j are the disorder parameters uniformly chosen in the interval − W 2 , W 2 . 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 ≤ j ≤ M and q i,0 = p i,0 = q i,M +1 = p i,M +1 = 0 for 1 ≤ i ≤ N .

III. NUMERICAL INTEGRATION SCHEMES
The Hamilton equations of motion dq dt = ∂H ∂p , dp dt = − ∂H ∂q , (III. 1) of the N degree of freedom (dof) Hamiltonian H = H(q, p), with q = (q 1 , q 2 , . . . , q N ) and p = (p 1 , p 2 , . . . , p N ) being respectively the system's generalized positions and momenta, can be expressed in the general setting of first order ODEs as is a vector representing the position of the system in its phase space and (˙) denotes differentiation with respect to time t. In Eq. (III.2) is the symplectic matrix with I N and O N being the N × N identity and the null matrices respectively, and with ( T ) denoting the transpose matrix. The variational equations (see for example [35,55]) govern the time evolution of a small perturbation w(t) to the trajectory x(t) with w(t) = (δq(t), δp(t)) = (δq 1 (t), δq 2 (t), . . . , δq N (t), δp 1 (t), δp 2 (t), . . . , δp N (t)) (which can also be written as w(t) = δx(t) = (δx 1 (t), . . . , δx N (t), δx N +1 (t), . . . , δx 2N (t))) and have the following forṁ are the 2N ×2N elements of the Hessian matrix D 2 H (x(t)) of the Hamiltonian function H computed on the phase space trajectory x(t). Eq. (III.5) is linear in w(t), with coefficients depending on the system's trajectory x(t). Therefore, one has to integrate the variational equations (III.5) along with the equations of motion (III.2), which means to evolve in time the general vector X(t) = (x(t), δx(t)) by solving the systeṁ .
(III. 7) In what follows we will briefly describe several numerical schemes for integrating the set of equations (III.7).

A. Non-symplectic integration schemes
Here we present the non-symplectic schemes we are going to 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).

Taylor series method −T IDES
The Taylor series method consists in expanding the solution X(t) at time t 0 + τ in a Taylor series of X(t) at t = t 0 Once the solution X at time t 0 + τ is approximated, X(t 0 + τ ) is considered as the new initial condition and the procedure is repeated to approximate X(t 0 + 2τ ) and so on so forth. Further information on this integrator can be found in [59,Sec. I.8]. If we consider in Eq. (III.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 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. [60]). For the simulations reported in this paper we used the software package called T IDES [60,61] which is freely available [62]. We particularly focused on the implementation of the T IDES package as it has been already used in studies of lattice dynamics [57]. The T IDES algorithm 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 T IDES algorithm allows us to choose both the integration time step τ and the desired 'one-step' precision of the integrator δ. In practice, the integration time step is accepted if the estimated local truncation error is smaller than δ.
It is worth noting that there exists an equivalent numerical scheme to the Taylor series method, derived from Lie series [63]. Indeed, Eq. (III.7) can be expressed as where L HZ is the Lie operator [64] defined as The formal solution of Eq. (III.9) reads X(t 0 + τ ) = e τ L HZ X(t 0 ), (III.11) 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 X(t 0 ). For further details concerning Lie series one can refer to [64,65]. The equivalence between the Lie and Taylor series approaches can be seen in the following way: for each element x i and δx i of the phase space vector X we compute

A Runge-Kutta family scheme − DOP 853
In the 1960s during the 'golden age' of the Taylor and Lie series method, the main difficulty was the derivation of the explicit expressions from the expansion of the differential operators as the forms of these numerical schemes depend strongly on the system considered (see [65] and references therein). Therefore, issues related to the portability and reusability of the generated numerical scheme was of serious concern. Fortunately, other methods have been developed to obtain Eq. (III.8). One way to perform this task was through the use of the well-known 'Runge-Kutta family' of algorithms (see e.g. [38,65]) which nowadays consist one of the most popular general purpose schemes for numerical computations. This is a sufficient motivation for us to introduce a p−stage Runge-Kutta method of the form where the real coefficients a ij , b i with i, j = 1, . . . , p are appropriately chosen in order to obtain the desired accuracy. In our study a 12-stage explicit Runge-Kutta algorithm called DOP 853 is used [66]. This numerical scheme is based on the Dormand and Prince's method [59]. As in the case of the T IDES algorithm, apart from the integration time step τ , DOP 853 also admits a 'one-step' precision δ (see [57] and references therein for more details).

B. Symplectic integration schemes
Hamiltonian flows possess integrals of motion (e.g. Hamiltonian energy, angular momentum, etc) and are recognized to exhibit a symplectic nature (see e.g. [38,Chap. VI] and references therein). These properties embedded in the Hamiltonian formalism play an important role in the explanation of the real world physical processes these models try to describe. Therefore, it is natural to look for numerical schemes that conserve these invariants to get as close as possible to real situations. The so-called SIs fulfill these requirements. In the case of autonomous Hamiltonian systems, SIs have the ability of conserving the Hamiltonian energy of the system. For this reason, they already have been extensively used in fields such as celestial mechanics [67], molecular dynamics [68], accelerator physics [41], condensed matter physics [20,26,53], etc.
In order to build SIs several methods can be used [38]. For instance, it has been proved that under certain circumstances the Runge-Kutta algorithm can preserve symplecticity [69,70]. However, many Hamiltonian functions are made up of the sum of the kinetic and potential energies, each one of which is described by an integrable expression. Thus, one commonly used way to approximate the system's exact solution is based on the construction of SIs through the split of the Hamiltonian function into integrable parts. In our study we also want to track the evolution of the deviation vector w(t) by solving Eq. (III.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 x(t) → x(t + τ ), along with the mapping δx(t) → δx(t + τ ).
Let us outline the whole process by considering a general autonomous Hamiltonian function H (q, p), which can be written as sum of I integrable intermediate Hamiltonian functions This decomposition implies that the resolvent operator e τ L A i Z in the formal solution in Eq. (III.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 + τ consists in approximating the action of e τ L HZ in Eq. (III.11) by a product of the operators e γiτ L A i Z for a set of properly chosen coefficients γ i . In our analysis we will call as number of steps of the particular SI the total number of successive application of each individual resolvent e τ L A i Z . Further details about this class of integrators can be found in [71][72][73] 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.

Two part split
Let us consider that the Hamiltonian H (q, p) can be separable into two integrable parts, namely H = A(q, p) + B(q, p). Then we can approximate the action of the operator e τ L HZ = e τ (L AZ +L BZ ) by the successive actions of products of the operators e τ L AZ and e τ L BZ [71][72][73][74]  depends only on the generalized momenta, whilst the potential part 16) 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 [53]. a. Symplectic integrators of order two. These integrators constitute the most basic schemes we can develop from Eq. (III.14) LF : The simplest example of Eq. (III.14) is the so-called leapfrog-Verlet scheme (LF ) (e.g. see [38,Sect. I.3.1] and [74]) having 3 individual steps LF (τ ) = e a1τ L AZ e b1τ L BZ e a1τ L AZ , (III.17) where a 1 = 1 2 and b 1 = 1.
ABA82: In addition, we use in our analysis the SI ABA82(τ ) = e a1τ L AZ e b1τ L BZ e a2τ L AZ e b2τ L BZ e a3τ L AZ e b2τ L BZ e a2τ L AZ e b1τ L BZ e a1τ L AZ , (III.20) with 9 individual steps [75,76], where the constants a i , b i with i = 1, 2, 3 can be found in Table 2 of [75]. We note that the ABA82 method is equivalent to the SABA 4 SI in [52]. b. Symplectic integrators of order four. The order of symmetric SIs can be increased by using a composition technique presented in [71]. According to that approach starting from a symmetric SI S 2n (τ ) of order 2n (n ≥ 1), we can construct a SI S 2n+2 (τ ) of order 2n + 2 as Here we adopt the notation of Ariech and Quispel [77].   19)] integrators we obtain the fourth order SIs SABA 2 Y 4 and SBAB 2 Y 4 having 13 individual steps. In particular, we get SABA 2 Y 4(τ ) = e d1a1τ L AZ e d1b1τ L BZ e d1a2τ L AZ e d1b1τ L BZ e a0τ L AZ e d0b1τ L BZ e d0a2τ L AZ e d0b1τ L BZ × e a0τ L AZ e d1b1τ L BZ e d1a2τ L AZ e d1b1τ L BZ e d1a1τ L AZ , (III.23) , and a 0 = d 1 a 1 + d 0 a 1 , and where d 0 = −2 1/3 / 2 − 2 1/3 , d 1 = 1/ 2 − 2 1/3 and a i , b i with i = 1, 2, 3 that can be found in Table 2 of [75]. Here SABA 2 K/SBAB 2 K: As was explained in [52] 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 ({·, ·}), 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 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 we produce the fourth order SIs SABA 2 K and SBAB 2 K, with g = (2 − √ 3)/24 and g = 1/72 respectively.
ABA864/ABAH864: The fourth order SIs ABA864 and ABAH864 were proposed in [54,75]. 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 [54], and In [71] a composition technique using fewer individual steps than the one obtained by the repeated application of Eq. (III.21) to SIs of order two was proposed, having the form (III.29) whose coefficients w i , i = 0, 1, 2, 3 are given in Table 1 in [71] for the case of the so-called 'solution A' of that table.
Here S 2 and S 6 respectively represent a second and sixth order symmetric SI. Note that Eq. (III.29) corresponds to the composition scheme s6odr6 of [79]. We also consider in our study the composition scheme s9odr6b of [79] which is based on 9 successive applications of S 2 The values of δ i , i = 1, 2, . . . , 5 in Eq. (III.30) can be found in the Appendix of [79]. Furthermore, we also implement the composition method of [80], which involves 11 applications of a second order SI S 2 , whose coefficients γ i , i = 1, 2, . . . , 6 are reported in Section 4.2 of [80]. d. Symplectic integrators of order eight. Following [71] 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, . . . , 7, and in particular the ones corresponding to the so-called 'solution A' and 'solution D' in Table 2  The SIs of this section will be implementing to numerically integrate the α-FPUT model [Eq. (II.1)] since this Hamiltonian system can be split into two integrable parts A(p) and B(q). In Sec. A 1 of the Appendix the explicit forms of the operators e τ L AZ and e τ L BZ are given, along with the operator e τ L KZ of the corrector term used by the SABA 2 K and SBAB 2 K SIs. In addition, in Sec. A 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.

Three part split
Let us now consider the case of a Hamiltonian function H (q, p) which can be separated into three integrable parts, namely H (q, p) = A (q, p) + B (q, p) + C (q, p). This for example could happen because the Hamiltonian function may require a split into more than two parts to be integrated, or to simplify the resolvent of one of the two components A or B of a two part split schemes discussed in Sec. III B 1. In such cases we approximate the action of the operator e τ L HZ of Eq. (III.11) by the successive application of operators e τ L AZ , e τ L BZ and e τ L CZ i.e.
. . , N are N constants of motion, and the Hamiltonian functions of the q− and p−hoppings and We note that a rather thorough survey on three part split SIs can be found in [46,47]. 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.
a. Symplectic integrators of order two. We present first the basic three part split scheme, which we call ABC2 . This integrator has 13 individual steps and has been already explicitly introduced in [82] and implemented in [47], where it was called ABC 4 [Y ] .
ABCS4: Implementing a composition scheme which was introduced in [72] and studied in [79] (where it was named s5odr4) we obtain the SI (III. 39) 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 [71]. Note that this SI has already been implemented in [46,47], where it was denoted as ABC 6 [Y ] .
s9ABC6: Implementing the composition of Eq. (III.30), with ABC2 in the place of S 2 , we get which was referred to as ABC 6 [KL] in [46,47].
s11ABC6: From the composition of Eq. (III.31) we create the SI scheme which has 45 individual steps. This integrator was referred to as ABC 6 [SS] in [46,47]. d. Symplectic integrators of order eight.
(III. 43) s19ABC8: Finally, we also implement the composition s19odr8b reported in [80,Eq. (13)] and construct the SI which has 72 individual steps. We note that this scheme corresponds to the SI ABC 8 [SS] considered in [46,47]. The explicit forms of the operators related to the three part split of the DDNLS Hamiltonians [Eqs. (II.2) and (II.3)] are given in Secs. A 2 and A 3 of the Appendix.

IV. NUMERICAL RESULTS
We test the efficiency of the integrators presented in Sec. III by using them to follow the dynamical evolution of the α- 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 ) and use it to compute the time evolution of the finite time mLE [33][34][35] in order to characterize the regular or chaotic nature of the trajectory through the estimation of the most commonly used chaos indicator, the mLE χ, which is defined as χ = lim t→+∞ X 1 (t). In Eq. (IV.1) · is the usual Euclidian norm, while w(t 0 ) and w(t 0 + τ ) are respectively the deviation vectors at t = t 0 and t 0 + τ > t 0 . In the case of regular trajectories χ tends to zero following the power law [35,83] whilst it takes positive values for chaotic ones.
A. The α-Fermi-Pasta-Ulam-Tsingou model We present here results on the computational efficiency of the symplectic and non-symplectic schemes of Sec. III for the case of the α-FPUT chain [Eq. (II.1)]. As this system can be split into two integrable parts we will use for its study the two part split SIs of Sec. III B 1. In our investigation we consider a lattice of N = 2 10 sites with α = 0.25 and integrate up to the final time t f = 10 6 two sets of initial conditions: • Case I 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 II F : Same as in case I 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 [65,Sec. 8.3]. Since we have chosen non-localized initial conditions, we also use an initial normalized deviation vector 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 1F [Eq. (II.1)]. This is done by registering the time evolution of the absolute relative energy error at each time step. In our analysis we consider two energy error thresholds E r ≈ 10 −5 and E r ≈ 10 −9 . The former, E r ≈ 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  [22][23][24][25][26][27]). In some cases, e.g. for very small values of conserved quantities, one may desire more accurate computations. Then E r ≈ 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. (IV.1)].
In Fig. 1 we show the time evolution of the absolute relative energy error E r (t) [panels (a) and (d)], of the finite time mLE X 1 (t) [panels (b) and (e)], and the required CPU time T C [panels (c) and (f)], for the cases I F and II F respectively, when the following four integrators were used: the fourth order SI ABA864 (blue curves), the sixth order SI SABA 2 Y 6 (red curves), the DOP 853 scheme (green curves) and the T IDES algorithm (brown curves). These results are indicative of our analysis as in our study we considered in total 35 different integrators (see Tables I and  II). In Fig. 1 the integration time steps τ of the SIs (reported in Tables I and II) were appropriately chosen in order to achieve E r ≈ 10 −9 , while for the DOP 853 and the T IDES algorithms E r (t) eventually grows in time as a power law [ Figs. 1(a) and (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 Figs. 1(b) and (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 Figs. 1(c) and (f). From these plots we see that the SIs need less computational time to perform the simulations than the DOP 853 and T IDES schemes.   In Table I [II] we present information on the performance of all considered integration schemes for the initial condition of case I F [II 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 [ Figs. 1(a) and (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 I and II we see that the best performing integrators are the fourth order SIs ABA864 and ABAH864 for E r ≈ 10 −5 , and the sixth and eighth order SIs s11SABA 2 6 and s19SABA 2 8 for E r ≈ 10 −9 . We note that the best SI for E r ≈ 10 −5 , the ABA864 scheme, shows a very good behavior also for E r ≈ 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 ≈ 10 −5 exhibited an unstable behavior failing to keep their E r values bounded. Thus, the higher order SIs are best suited for more accurate computations. It is also worth mentioning here that the ranking presented in Tables I and II 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 τ 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.

B. The 1D disordered discrete nonlinear Schrödinger equation system
We investigate the performance of various integrators of Sec. III for the 1D DDNLS system [Eq. (II.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 α-FPUT system) up to the final time t f = 10 6 . We note that, as was already mentioned in Sec. III B 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 I 1D : We initially excite 21 central sites by attributing to each one of them the same constant norm ζ j = (q 2 j + p 2 j )/(2S 1D ) = 1, 1 ≤ i ≤ N , for W = 3.5 and β = 0.62. This choice sets the total norm S 1D = 21. The random disorder parameters i , 1 ≤ i ≤ N . are chosen so that the total energy is H 1D ≈ 0.0212.
• Case II 1D : Similar set of initial conditions as in the case I 1D but for W = 3, β = 0.03 and the random disorder parameters i , 1 ≤ i ≤ N are chosen such that H 1D ≈ 3.4444.
We note that cases I 1D and II 1D have been studied in [32] 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 norm distribution ζ j , 1 ≤ j ≤ N , by computing the distribution's second moment where j = N j=1 jζ j is the position of the center of the distribution [20,22,24,26,32,46,47]. We also check how accurately the values of the system's two conserved quantities, i.e. its total energy H 1D [Eq. (II. In addition, we compute the finite time mLE X 1 (t) [Eq. (IV.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 ≤ 6, which we implement in order to achieve an accuracy of E r ≈ 10 −5 , and (ii) SIs of order eight used for E r ≈ 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 ≈ 10 −5 because, as was shown in [46] and also discussed in Sec. IV A, 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 ≈ 10 −9 ) simulations [46]. Therefore we do not consider such integrators in our study.
In Fig. 2 Table III. The integration time step τ of the SIs was selected so that the relative energy error is kept at E r ≈ 10 −5 [Fig 2(a)]. From the results of Fig. 2(b) we see that the SIs do not keep S r constant as these integrators were designed to accurately conserve the value of the Hamiltonian. Nevertheless, our results show that we lose no more than two orders of precision (in the worst case of the ABC2 scheme) during the whole integration. On the other hand, both the absolute relative energy [E r (t)] and norm [S r (t)] error of the T IDES and DOP 853 integrators increase in time, with the T IDES scheme behaving better than the DOP 853 one. Figs. 2(c) and (d) show that all integrators correctly reproduce the dynamics of the system, as all of them practically produce the same norm distribution at t f = 10 6 [ Fig. 2(c)] and the same evolution of the m 2 (t) [Fig. 2(d)]. We note that m 2 increases by following a power law m 2 ∝ t α with α = 1/2, as was expected for the strong chaos dynamical regime (see for example [32] and references therein).
In Fig. 3 Fig. 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) ∝ t α L with α L ≈ −0.3 (case I 1D ) and α L ≈ −0.25 (case II 1D ), in accordance to the results of [31,32]. 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 III for case I 1D and in Table IV for case II 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 α-FPUT model, the DOP 853 and T IDES 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 T IDES being more precise. The integrators exhibiting the best performance for E r ≈ 10 −5 are the sixth order SIs s11ABC6 and s9ABC6, while for E r ≈ 10 −9 we have the eighth order SIs s19ABC8 and s17ABC8, with the s11ABC6 scheme performing quite well also in this case.

C. The 2D disordered discrete nonlinear Schrödinger equation system
We now investigate the performance of the integrators used in Sec. IV B for the computationally much more difficult case of the 2D DDNLS lattice of Eq. (II.3), as its Hamiltonian function can also be splitted into three integrable parts. In order to test the performance of the various schemes we consider a lattice with N × M = 200 × 200 sites, resulting to a system of 4 × 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 α-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. [37,58]). We consider again two sets of initial conditions: • Case I 2D : We initially excite 7 × 7 central sites attributing to each one of them the same norm density ζ i,j = (q 2 i,j + p 2 i,j )/(2S 2D ) = 1/6 so that the total norm is S 2D = 49/6, for W = 16 and β = 6. The disorder parameters i,j , 1 ≤ i ≤ N , 1 ≤ j ≤ M , are chosen so that the initial total energy is H 2D ≈ 1.96.
The initial normalized deviation vector considered in our simulations has random non-zero values only at the 7 × 7 initially excited sites for case I 2D , and only at site i = 100, j = 100 for case II 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 [49,84]. 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 (Sec. IV B), in order to evaluate the performance of the used integrators we follow the time evolution of the norm density ζ i,j , 1 ≤ i ≤ N , 1 ≤ j ≤ M and compute the related second moment m 2 and participation number P where (i, j) T = i,j (i, j) T ζ i,j is the mean position of the norm distribution. We also evaluate the absolute relative energy [E r (t)] and norm [S r (t)] errors and compute the finite time mLE X 1 (t).
In Fig. 4 we present results obtained for case I 2D by the four best performing SIs (see Table V), the s11ABC6 (red curves), s9ABC6 (blue curves), ABCY 6 A (green curves) and ABCY 4 (brown curves) schemes, along with the DOP 853 integrator (grey curves). The integration time time steps τ of the SIs were adjusted in order to obtain an accuracy of E r ≈ 10 −5 [ Fig. 4(a)], while results for the conservation of the second integral of motion, i.e. the system's total norm, are shown in Fig. 4(b). We see that for all SIs the S r values increase slowly, remaining always below S r ≈ 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 DOP 853 integrator increase, although the choice of δ = 10 −16 again ensures high precision computations.
The norm profiles at the final integration time t f = 10 5 along the axis i = 100 [ Fig. 4(c)] and j = 100 [ Fig. 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 Fig. 4(e) [Fig. 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 [27] indicating that the two-dimensional lattice is being thermalized. The results of Figs. 4(e) and (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 Fig. 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 ∝ t −1 power law observed for regular motion. This behavior was also observed for the 2D DKG model [53], as well as for the 1D DKG and DDNLS systems in [31,32] where a power law X 1 (t) ∝ t α L with α L ≈ −0.25 and α L ≈ −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 [53] 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 V and VI, 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 s11ABC6 and s9ABC6 have the best performance for E r ≈ 10 −5 and the SIs s19ABC8 and s17ABC8 for E r ≈ 10 −9 .

V. CONCLUSIONS
In this work we carried out a methodical and detailed analysis of the performance of several symplectic and nonsymplectic 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 α-FPUT chain, as well as the 1D and 2D DDNLS models. In the case of the α-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 interrogators we evolved in time different sets of initial conditions and evaluated quantities related to (a) the  ABA864 and ABAH864 SIs of order four to be the best schemes for moderate energy accuracies (E r ≈ 10 −5 ), while the s11SABA 2 6 and s19SABA 2 8 SIs of order six and eight respectively were the best integration schemes for higher accuracies (E r ≈ 10 −9 ). In particular, the ABA864 scheme appears to be an excellent, general choice as it showed a very good behavior also for E r ≈ 10 −9 . Concerning the 1D and the 2D DDNLS models our simulations showed that the SIs s9ABC6 and s11ABC6 (order six), along with the SIs s17ABC8 and s19ABC8 (order eight) are the best integrators for moderate (E r ≈ 10 −5 ) and high (E r ≈ 10 −9 ) accuracy levels respectively. The DOP 853 and T IDES 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 ≈ 10 −5 for the α-FPUT system (similar behaviors were reported in [53] for the DKG model). This happens because the required integration time step τ needed to keep the absolute relative energy error at E r ≈ 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 ≥ 6 are more suitable for calculations that require higher accuracies (e.g. E r ≈ 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 [53,[55][56][57], 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 Sec. III, along with the explicit presentation in the Appendix of the related differential operators for many commonly used classical, manybody 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 [85]) and for the dynamical study of various lattice models, like for example of arrays of Josephson junctions in regimes of weak non-integrability [50], granular chains [86] and DNA models [87], to name a few.

APPENDIX 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. (III.7)Ẋ H (x(t)) · δx(t) . (A.1)

The α-Fermi-Pasta-Ulam-Tsingou model
The Hamiltonian of the α-FPUT chain [Eq. (II.1)] can be split into two integrable parts as The set of equations of motion and variational equations for the Hamiltonian function A (p) is and the corresponding operator e τ L AZ , which propagates the values of q i , p i , δq i and δp i for τ time units in the future, obtaining q i , p i , δq i and δp i , takes the form e τ L AZ : In a similar way for the B (q) Hamiltonian of Eq. (A.2) we get dX dt = L BZ X : and e τ L BZ : According to Eq. (III.26) the accuracy of the SABA n K and SBAB n K integrators can be improved by using a corrector Hamiltonian K [52]. In the case of a separable Hamiltonian H(q, p) = A(p) + B(q) with A (p) = N i=1 p 2 i /2, the corrector K becomes (A.7) For the α-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 τ L KZ e τ L KZ : 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 , 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 reported in the Appendix of [53] where the operators for the 1D and 2D DKG models were reported (with the exception of the corrector term). For completeness sake in Sec. A 4 we provide the explicit expression of the e τ L BZ and e τ L KZ operators for some commonly used Hamiltonians, which can be split into two integrable parts, one of which is the usual kinetic energy

The one-dimensional disordered discrete nonlinear Schrödinger equation
Here we focus on the 1D DDNLS system, whose Hamiltonian [Eq. (II. 2)] can be split into three integrable parts as The set of equations of motion and variational equations associated with the Hamiltonian function A 1 is with θ i = i + β(q 2 i + p 2 i )/2 for i = 1, 2, . . . , N being constants of the motion. The corresponding operator e τ L A 1 Z takes the form e τ L A 1 Z : with J i = 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 The set of equations of motion and variational equations associated to the intermediate Hamiltonian functions B 1 and C 1 are respectively dX dt = L B1Z X : , and dX dt = L C1Z X : (A. 15) These yield to the operators e L B 1 Z and e L C 1 Z given by .
(A. 16) As in the case of the α-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 [53] for the 1D DKG model, should be implemented when fixed boundary conditions are imposed.

Other Hamiltonian models which can be split into two integrable parts
We present here the exact expressions of the resolvent 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 β-FPUT chain, the KG model, and the classical XY model (a JJC system) [50,[88][89][90]. Similarly to the α-FPUT chain of Eq. (II.1), the Hamiltonians H(q, p) of each of these systems can be split as with appropriately defined potential terms B (q): where β, k and E J are real coefficients. Obviously for all these systems the operator e τ L AZ of the kinetic energy part is the same as for the α-FPUT chain in Eq. (A.4). Thus, we report below only the expressions of the operators e τ L BZ and e τ L KZ when periodic boundary conditions are imposed.