Splitting Methods for Semi-Classical Hamiltonian Dynamics of Charge Transfer in Nonlinear Lattices

Citation: Baj¯ ars, J.; Archilla, J.F.R. Splitting Methods for Semi-Classical Hamiltonian Dynamics of Charge Transfer in Nonlinear Lattices. Mathematics 2022,10, 3460. https:// doi.org/10.3390/math10193460 Academic Editors: Vicente Martínez and Pablo Gregori Received: 2 August 2022 Accepted: 18 September 2022 Published: 22 September 2022 Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations. Copyright: © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/). mathematics Article Splitting Methods for Semi-Classical Hamiltonian Dynamics of Charge Transfer in Nonlinear Lattices J¯ anis Baj¯ ars 1,† and Juan F. R. Archilla 2,*,† 1 Faculty of Physics, Mathematics and Optometry, University of Latvia, Jelgavas Street 3, LV-1004 Riga, Latvia 2Group of Nonlinear Physics, Universidad de Sevilla, ETSII, Avda Reina Mercedes s/n, 41012 Sevilla, Spain *Correspondence: ar[email protected] † These authors contributed equally to this work. Abstract: We propose two classes of symplecticity-preserving symmetric splitting methods for semi-classical Hamiltonian dynamics of charge transfer by intrinsic localized modes in nonlinear crystal lattice models. We consider, without loss of generality, one-dimensional crystal lattice models described by classical Hamiltonian dynamics, whereas the charge (electron or hole) is modeled as a quantum particle within the tight-binding approximation. Canonical Hamiltonian equations for coupled lattice-charge dynamics are derived, and a linear analysis of linearized equations with the derivation of the dispersion relations is performed. Structure-preserving splitting methods are constructed by splitting the total Hamiltonian into the sum of Hamiltonians, for which the individual dynamics can be solved exactly. Symmetric methods are obtained with the Strang splitting of exact, symplectic flow maps leading to explicit second-order numerical integrators. Splitting methods that are symplectic and conserve exactly the charge probability are also proposed. Conveniently, they require only one solution of a linear system of equations per time step. The developed methods are computationally efficient and preserve the structure; therefore, they provide new means for qualitative numerical analysis and long-time simulations for charge transfer by nonlinear lattice excitations. The properties of the developed methods are explored and demonstrated numerically considering charge transport by mobile discrete breathers in an example model previously proposed for a layered crystal. Keywords: semi-classical Hamiltonian dynamics; splitting methods; symplectic integrators; lattice models; charge transfer; intrinsic localized modes; discrete breathers MSC: 65P10; 37M05 1. Introduction In this introduction, we briefly recall the phenomenon of the coupling of electric charge and lattice vibrations, the concept of intrinsic localized modes in nonlinear lattices, first without charge and then with it. Then, we describe the numerical challenges and introduce the methods this paper proposes to overcome them. The introduction finishes with the outline of the paper. 1.1. Coupling of Electric Charge and Lattice Vibrations The coupling of electric charge with lattice vibrations has been of interest since long ago, perhaps starting with the pioneering work by Landau, where the self-trapping of an electron in a crystal due to the lattice deformation by the electron was described [ 1 ]. It was followed by Pekar, who proposed the term polaron [ 2 ], and subsequently in [ 3 ]. These works describe a localized deformation coupled to a localized charge and apply, for example, to ionic crystals, polar semiconductors, molecular crystals, and polymers. The development of the subject can be read in the review book by Alexandrov [ 4 ]. The concept of a polaron Mathematics 2022,10, 3460. https://doi.org/10.3390/math10193460 https://www.mdpi.com/journal/mathematics Mathematics 2022,10, 3460 2 of 43 can be extended to other lattice systems in which an electron or hole can be described by an atomic wave function and therefore bound to a specific particle, which can be an atom, ion, or molecule. This approach is named the tight-binding approximation (TBA) [ 5 ]. Often, the electron or hole is described by quantum mechanics, but the particle is treated classically, and these models are therefore called semi-classical models. The resulting Hamiltonian usually has three components, a quantum or electronic Hamiltonian, a classical or lattice Hamiltonian, and an interaction term. A well-known model was proposed by Holstein [ 6 , 7 ]. In these systems, localized solutions of the electron or hole probability appear. 1.2. Intrinsic Localized Modes without Charge Localized excitations without charge have also been studied in discrete nonlinear lattices. Important steps were the Frenkel–Kontorova model in the 1930s [ 8 , 9 ] and Davidov’s solitons in the 1970s [ 10 ]. If these excitations have a topological charge, they are called kinks and antikinks. In a lattice, they correspond to moving vacancies and interstitials, respectively, the latter also named crowdions. If they do not have a topological charge, they do not transport mass and are called solitons. If they also have a vibration, they are called breathers or intrinsic localized modes (ILMs) [ 11 ]. However, the term ILM, being unspecific, is also used as a synonym for nonlinear localized excitation. Kinks and breathers are sometimes associated with a finite amplitude extended wave, a wing, and are called pterokinks or nanopterons and pterobreathers, respectively [12,13]. 1.3. Nonlinear Lattice Excitations Transporting Charge In a tight-binding, semi-classical model in a nonlinear lattice, nonlinear lattice excitations can bind to the wave function and in this way transport charge, bringing about objects that have received different names from different authors, such as polarobreathers [14–17] , electron-vibron breather [ 17 ], and solectrons or soliton-elecrons [ 18 – 22 ]. For example, in Figure 1 , we illustrate the evolution of charge probability carried by a mobile discrete breather in the model example of Section 2.7, where we plot the amplified lattice particle displacements together with the charge probability indicated by the color plot. Note that a kink or crowdion in an ionic crystal transports charge by itself but could also be bound to an additional one [23,24]. Figure 1. Charge transfer by a mobile discrete breather in the model example of Section 2.7. Solid lines indicate the lattice particle displacements amplified by the factor 1.5, while the color plot illustrates the probability of finding the charge at a lattice site at a given time. Numerical results are obtained with the explicit, symplecticity-preserving, symmetric numerical method PQDABADQP, see Section 4.3 . See also Section 5for additional numerical results of charge transfer by discrete breathers. Mathematics 2022,10, 3460 3 of 43 Recently, charge transport in absence of an electric field was observed experimentally in layered silicates, a phenomenon called hyperconductivity [ 25 – 28 ]. This is a characteristic of moving localized excitations coupled to a charge, called quodons by the authors. The energy is given to the lattice by incident alpha particles, and quodons propagate long distances within the lattice without a potential bias. Semi-classical tight-binding models can, in principle, model this behavior as well. 1.4. Existing and Proposed Numerical Methods ILMs have been extensively studied from analytic and numerical points of view in crystal lattice models [ 12 , 21 , 29 – 35 ], to mention but a few references. Traditionally, such lattice models at constant energy are described by classical Hamiltonian dynamics with empirical particle interaction potentials and by thermostated Hamiltonian dynamics at a given temperature [ 36 ]. There exist very successful numerical methods, such as, for example, the Verlet method, that preserve the underlying structure of the Hamiltonian equations [36,37]. The transport of charge (electrons or holes) by ILMs poses new numerical simulation challenges due to the different oscillation time scales of the charge and the lattice dynamics because an electron mass is ≃ 10 4 times smaller than most lattice atoms or ions and the small value of the Planck constant. In addition, the dynamics is conservative. Thus, we may be required to use very small time steps to resolve high-frequency oscillations in charge and limit the amount of dissipation in (standard) explicit numerical integration. This motivates the consideration of structure-preserving splitting and composition methods [ 38 ] for the numerical integration of coupled lattice-charge differential equations. We demonstrate that lattice-charge dynamics can be stated into classical canonical Hamiltonian form, and while the total Hamiltonian is not separable in all variables, it is possible to construct explicit numerical methods that are symplectic, preserve the time-reversibility, and have good approximate energy conservation in long-time numerical simulations. This property can be attributed to the existence of the modified Hamiltonian of a symplectic method applied to a Hamiltonian system [ 37 ]. We show that different splitting strategies allow computations with larger time steps, e.g., compared to the standard fourth-order RungeKutta method. Computational analysis is performed in linear and nonlinear regimes, where, in the latter case, we consider charge transfer by mobile discrete breathers. The performed analysis demonstrates that explicit splitting methods not only approximately conserve the Hamiltonian up to the second order but also the total charge probability in long-time numerical simulations. Apart from computationally efficient explicit methods, we also develop semi-implicit splitting methods that are symplectic and conserve the charge probability exactly. Moreover, they require only one solution of a linear system per time step. They provide charge probability amplitude solutions that are unconditionally stable. In addition, semi-implicit methods exactly preserve the rotational invariance of the charge amplitude variables, which is not exactly preserved by the explicit symplectic splitting methods. We demonstrate that semi-implicit splitting methods have Hamiltonian conservation errors that are smaller or on par with the errors of the explicit methods. Such a splitting method approach will further allow the development of multiple timestepping schemes, such as the impulse method [ 36 , 37 ], to further improve numerical integration of the multiscale dynamics and construct higher-order composition methods [39,40] (see also [ 36 – 38 ]) to improve the numerical accuracy of charge transfer simulations. Importantly, the proposed methods can be easily incorporated into splitting methods for the thermostated dynamics to achieve efficient and accurate sampling [ 36 ], i.e., to study charge transport in thermalized crystal lattice models at different temperatures. 1.5. Outline The paper is organized in the following way. In Section 2we derive the general model of lattice-charge dynamics described in canonical Hamiltonian form. At the end Mathematics 2022,10, 3460 4 of 43 of Section 2 , we present an example model for numerical simulations in Section 5. Linear analysis of the linearized Hamiltonian system of equations is demonstrated in Section 3 with the derived lattice and charge dispersion relations. Two classes of splitting methods are proposed in Section 4, i.e., explicit and semi-implicit, both preserving symplecticity and time-reversibility. All methods in explicit representation forms are listed in Appendix A. In Section 5, we demonstrate numerical results and perform a detailed numerical study of the proposed splitting methods. Discussion and conclusions are provided in Sections 6and 7, respectively. 2. Semi-Classical Hamiltonian Dynamics In this section, we describe the classical Hamiltonian dynamics of a crystal lattice coupled with a quantum Hamiltonian to model a charge (electron or hole) transfer by nonlinear lattice excitations. After variable and parameter rescaling, we obtain the dimensionless total semi-classical Hamiltonian for which we derive canonical Hamiltonian equations. In addition, we present a model example for the numerical analysis and simulations of Sections 4and 5. 2.1. The Lattice Hamiltonian Classical one-dimensional nonlinear lattice dynamics of N particles are modeled by the lattice Hamiltonian in the following form [12,31,33]: Hlat =KE+UE+VE= N ∑ n=1 1 2Mn p2 n+U(qn) + 1 2 N ∑ n0=1 n06=n V(|qn−qn0|)!, (1) where KE is the kinetic energy, UE is the on-site potential energy, and VE is the radial interparticle potential energy. In our considerations, the crystal lattice model (1) is a onedimensional model of a three-dimensional layered crystal, where the on-site potential UE models the periodic energy of the particles in a close-packed line due to the rest of the crystal. The interparticle potential VE models the interaction between the particles in the close-packed line. This interaction is generically repulsive at short distances and attractive at larger ones. qn∈Ω⊂R is the n th particle’s position, where Ω is a bounded computational domain (segment) containing the crystal lattice imposed by the periodic boundary conditions on qn . Alternatively, some form of confining potentials of q1 and qN could be added to (1) to impose boundary conditions [ 36 ]. pn∈R and Mn> 0 are the n th particle’s momentum and mass, respectively. From (1) , we derive the system of Hamiltonian equations: ˙ qn=1 Mn pn, (2) ˙ pn=−U0(qn)− N ∑ n0=1 n06=n ∂qnV(|qn−qn0|), (3) for all n= 1, . . . , N , subject to initial conditions at time t= 0, where the dot indicates time derivative with respect to t≥0, i.e., qnand pnare functions of time variable t. A special case of (1) includes the so-called fixed close neighbor interaction model with the Hamiltonian: Hf ixed lat = N ∑ n=1 1 2Mn p2 n+U(qn) + 1 2∑ n0∈Nn V(|qn−qn0|)!, (4) where N n is a set of n0 indices specifying a particle’s qn fixed neighbors, i.e., if n0∈ N n , then n∈ N n0 . Commonly, in one-dimensional models N n={n− 1, n+ 1 } . Such model (4) provides a computationally efficient approximation of (1) when only close-range Mathematics 2022,10, 3460 5 of 43 interaction potentials V are considered. Alternatively, a smooth cut-off of the potentials V can be considered, with N n being adjusted in time along the solution of the Hamiltonian dynamics (2) and (3) . Note that for generality in the formulation (1) , we have included both types of empirical potentials, but their actual use and form will highly depend on the problem in question. Below, we list the most general properties of the on-site and interaction potentials commonly used in practice. We assume that U and V are smooth functions of qn and satisfy the conditions for any nthat we present below. 2.1.1. The On-Site Potential Conditions on Uare: U(q0 n) = 0, U0(q0 n) = 0, U00(q0 n)>0, (5) where q0 n=σ(n− 1 ) are the lattice equilibrium particle positions at which all forces add to zero, i.e., for all n there is an index pair (n0 , n00) , n06=n00 6=n , where also n0 , n00 ∈ N n in (4) , such that ∂qnV(|q0 n−q0 n0|) = −∂qnV(|q0 n−q0 n00|), (6) where σ> 0 is the lattice constant and defines the length scale of the lattice particle interaction. It is often convenient to define the on-site potential function U independently of the lattice site index n, i.e., U(qn) = Uqn−q0 n σ,∀n=1, . . . , N. (7) The definition (7) and conditions (5) imply that we have a harmonic approximation of the on-site potential Uat the particle equilibrium: Uh(qn) = 1 2ω2 0(qn−q0 n)2,ω0=qU00(q0 n), (8) where ω0 defines the oscillation frequency of isolated oscillators of unit mass for the linearized equations with V=0. 2.1.2. The Interaction Potential Modeling the particle interactions with a radial potential V(r), where r=rn,n0=|qn−qn0|>0, n06=n, we differentiate between two types of potentials, i.e., potentials with energy minimum, if both forces of repulsion and attraction are considered, and potentials without energy minimum, e.g., when only particle repulsion is considered. In the latter case, we achieve an equilibrium state, e.g., by considering periodic boundary conditions. In the former case, we impose that the functions V(r) and V00(r) are monotonically decreasing in the interval ( 0, σ) , such as V0(r)< 0, and the function V(r) is monotonically increasing for r>σ, such as V0(r)>0, and V(σ) = −V0,V0(σ) = 0, V00(σ)>0 as well as lim r→∞V(r) = 0−, lim r→∞V0(r) = 0+, (9) where V0> 0 describes the relative strength of the potential and σ> 0 is an equilibrium interaction distance between two particles. The conditions (9) (see also (11) ) imply that the interaction strength between particles diminishes as the distance between particles Mathematics 2022,10, 3460 6 of 43 increases. Examples of such radial interparticle potentials include the Lennard–Jones and Morse potentials: VLJ(r) = V0σ r12 −2σ r6(10) and VM(r) = V0exp−2br σ−1−2 exp−br σ−1,b>0, respectively, where the dimensionless parameter b controls the width of the Morse potential well. In the latter case, when only repulsion forces are considered in the models of equally charged particles, e.g., ions, the functions V(r)> 0 and V00(r)> 0 are assumed to be monotonically decreasing for all r> 0, while the function V0(r)< 0 is assumed to be monotonically increasing as well as V(σ) = V0and the following conditions hold: lim r→∞V(r) = 0+, lim r→∞V0(r) = 0−. (11) Examples of such radial potentials include the Coulomb and Pauli repulsion potentials: VC(r) = V0 σ r and VP(r) = V0 σ rexp−br σ−1,b>0, respectively, where the dimensionless parameter b models the rate of the potential decay, and σ>0, as already stated above, defines the length scale of the interaction. In either case, the interaction potentials can be expanded in the Taylor series at r=σ , i.e., V(r) = V(σ) + V0(σ)r σ−1+1 2V00(σ)r σ−12+. . . In the former case, since V0(σ) = 0, we obtain the harmonic approximation of the potential Vin the following form: Vh(r) = −V0+1 2V00(σ)r σ−12. To ensure that the Hamiltonians (1) and (4) at the lattice equilibrium (q0 n , p0 n) = (σ(n− 1 ) , 0 ) for all n= 1, . . . , N are equal to zero. Without loss of generality, we redefine the interaction potential Vin (1) and (4) as: V=V(|qn−qn0|)−V(|q0 n−q0 n0|), (12) which does not alter the Hamiltonian Equations (2) and (3). 2.2. The Charge Hamiltonian While the lattice dynamics is modeled by the classical Hamiltonian dynamics (1) , an extra charge is treated as a quantum particle. The coupled motion of the charge with the lattice (1) is modeled within the tight-binding approximation [ 5 ]. With index n , we indicate the internal quantum number of the charge state corresponding to the lattice (1) particle at the position qn . In addition, we assume that there is only one quantum state per particle that the charge can occupy. Mathematics 2022,10, 3460 7 of 43 2.2.1. The Wave Function and the Charge Probability The extra charge dynamics are described by the wave function, also called state function [41]: |φ(t)i= N ∑ n=1 cn(t)|ni, (13) where the ket |ni represents the element of the orthonormal basis describing the state in which the charge is bound to the site n and cn(t)∈C are time-dependent dimensionless complex coefficients known as probability amplitudes. The adjoint operator of |φ(t)iis: hφ(t)|= N ∑ n=1 c∗ n(t)hn|, (14) where c∗ n is the complex conjugate of cn . The ket and bra operators, |ni and hn| , satisfy the property hn|n0i=δn,n0 , where δn,n0 is the Kronecker delta function. Accordingly, we can define a space-discrete, time-dependent wave function ψn(t) = hn|φ(t)i= N ∑ n0=1 cn0(t)hn|n0i= N ∑ n0=1 cn0(t)δn,n0=cn(t), to obtain a probability of an electron or hole being at a state nat time t, i.e., |ψn(t)|2=ψ∗ n(t)ψn(t) = c∗ n(t)cn(t) = |cn(t)|2=:Pn(t), (15) which defines the probability function Pn . Thus, we have the total probability conservation in the following form: N ∑ n=1 Pn(t) = 1. (16) To simplify the notation, in what follows, we omit explicit dependence on t from the variables cn, i.e., cn=cn(t). 2.2.2. The Evolution of the Wave Function: The Schrödinger Equation The state function (13) is the solution to the discrete Schrödinger equation: i¯h|˙ φ(t)i=Hc|φ(t)i, (17) where ¯his the Planck constant, and Hcis the charge Hamiltonian operator: Hc= N ∑ n=1 En(q1, . . . , qN)|nihn|− N ∑ n0=1 n06=n J(qn,qn0)|nihn0|!, (18) where indices n and n0 indicate internal quantum numbers of the charge states corresponding to the lattice (1) particles at positions qn and qn0 , respectively. Notice that we have not added factor 1 / 2 in front of the second sum in the Hamiltonian operator (18) as in (1) , since |nihn0| 6=|n0ihn| for arbitrary different n and n0 . Similarly to the lattice Hamiltonian (4) , we can consider the charge Hamiltonian operator with fixed close neighbor interactions: Hf ixed c= N ∑ n=1 En(qn0|n0∈Nn∪{n})|nihn|− ∑ n0∈Nn J(qn,qn0)|nihn0|!(19) for which the following derivations hold as well. Mathematics 2022,10, 3460 8 of 43 2.2.3. The Local Charge Energy Smooth multivariable function En(q1 , . . . , qN)∈R describes the local charge energy at site n and, in general, will depend on the lattice particle positions, or can be modeled as a constant. In addition, we impose that at the lattice equilibrium ∂qnEn(q0 1, . . . , q0 N) = 0, ∀n=1, . . . , N, (20) and for any n0there exists n00 such that ∂qnEn0(q0 1, . . . , q0 N) = −∂qnEn00(q0 1, . . . , q0 N). (21) In practice, the charge energy function En can be modeled to consist of a sum of charge on-site and interaction potentials, i.e., En(q1, . . . , qN) = QUc(qn) + Q N ∑ n0=1 n06=n Vc(|qn−qn0|), (22) respectively, where −Uc (note the minus sign) and Vc are smooth functions and share the same properties of the on-site and interaction potentials U and V of Section 2.1, constant Q= 1 for a positive charge, and Q=− 1 for a negative charge. Note that, depending on the model, charge energy functions E1 and EN must be adjusted if non-periodic boundary conditions are used. 2.2.4. The Charge Transition Matrix The symmetric transition matrix elements Jn,n0=J(qn , qn0)≥ 0, i.e., Jn,n0=Jn0,n , of charge transfer between states n and n0 will depend on the particle interaction distance rn,n0and is often modeled with exponential decay: [18,19,21,22]: J(qn,qn0) = J(rn,n0) = J0exp−αrn,n0 σ, lim rn,n0→∞J(rn,n0) = 0, (23) where the dimensionless parameter α> 0 specifies the rate of exponential decay, while σ> 0 defines lattice particle interaction length scale, see Section 2.1. Constant J0≥ 0 models the relative strength of the charge transfer from one site to another and is model dependent. In a general setting, we assume that J is a smooth monotonically decreasing function of particle interaction distance r>0 or taken to be constant. 2.2.5. Dynamical Equations for the Charge Probability Amplitude Notice that the right-hand side of the Schrödinger Equation (17) Hc|φ(t)i= N ∑ m=1 cmHc|mi = N ∑ m=1 cm N ∑ n=1 En(q1, . . . , qN)|nihn|mi− N ∑ n0=1 n06=n J(qn,qn0)|nihn0|mi! = N ∑ n=1 En(q1, . . . , qN)N ∑ m=1 cm|nihn|mi− N ∑ n0=1 n06=n J(qn,qn0)N ∑ m=1 cm|nihn0|mi! = N ∑ n=1 En(q1, . . . , qN)cn|ni− N ∑ n0=1 n06=n J(qn,qn0)cn0|ni!, (24) Mathematics 2022,10, 3460 9 of 43 and the projection of the Equation (17) onto the state n: i¯hhn|˙ φ(t)i=hn|Hc|φ(t)i, i¯h N ∑ m=1 ˙ cmhn|mi= N ∑ m=1 Em(q1, . . . , qN)cmhn|mi− N ∑ n0=1 n06=m J(qm,qn0)cn0hn|mi!, leads to dynamical equations for the probability amplitude cn, i.e., i¯h˙ cn=En(q1, . . . , qN)cn− N ∑ n0=1 n06=n J(qn,qn0)cn0, (25) where n=1, . . . , N, with the charge classical Hamiltonian Hc=hφ(t)|Hc|φ(t)i = N ∑ m=1 N ∑ n=1 En(q1, . . . , qN)c∗ mcnhm|ni− N ∑ n0=1 n06=n J(qn,qn0)c∗ mcn0hm|ni! = N ∑ n=1 En(q1, . . . , qN)c∗ ncn− N ∑ n0=1 n06=n J(qn,qn0)c∗ ncn0!, (26) where we have applied the adjoint operator (14) to (24) and used the Kronecker delta property of the ket and bra operators. 2.3. The Coupled Hamiltonian In the previous sections, we derived the lattice and charge Hamiltonians (1) and (26) , respectively, with their dynamical Equations (2) , (3) and (25) . In this section, we formulate the coupled Hamiltonian system with the total Hamiltonian Htotal =Hlat +Hc and the associated lattice-charge coupled system of equations: ˙ qn=1 Mn pn, (27) ˙ pn=−U0(qn)− N ∑ n0=1 n06=n ∂qnV(|qn−qn0|)− N ∑ n0=1 ∂qnEn0(q1, . . . , qN)c∗ n0cn0+ N ∑ n0=1 n06=n ∂qnJ(qn,qn0)(c∗ ncn0+c∗ n0cn), (28) i¯h˙ cn=En(q1, . . . , qN)cn− N ∑ n0=1 n06=n J(qn,qn0)cn0, (29) for all n= 1, . . . , N , where we have obtained the additional charge-lattice force term entering Equation (28). Recall that c∗ ncn0+c∗ n0cn=2Re(c∗ ncn0). With the energy function Engiven by (22), we have that N ∑ n0=1 ∂qnEn0(q1, . . . , qN)c∗ n0cn0=QU0 c(qn)c∗ ncn+Q N ∑ n0=1 N ∑ m=1 m6=n0 ∂qnVc(|qn0−qmk)c∗ n0cn0, Mathematics 2022,10, 3460 16 of 43 all lattice particle masses are equal, i.e., M=Mn for all n= 1, . . . , N . Although, they are highly restrictive assumptions, we are still able to obtain very valuable information. 3.1. Linearization of Lattice Forces The Hamiltonian dynamics (4) admits the following linearized equations. The force from the interaction potential Vof a particle qn0acting on the particle qnin (3) is given by: S(qn,qn0) = −∂qnV(|qn−qn0|) = −1 rn,n0V0(rn,n0)(qn−qn0). (62) The linearized force of (62) at the lattice equilibrium is Slin(qn,qn0) = ∂qnS(q0 n,q0 n0)(qn−q0 n) + ∂qn0S(q0 n,q0 n0)(qn0−q0 n0), where the terms S(q0 n , q0 n0) have been omitted, even if S(q0 n , q0 n0)6= 0, since (see (6) ) there exists a particle qn00 acting on the particle qn with the force S(q0 n , q0 n00) = −S(q0 n , q0 n0) . Note also that ∂qnS(q0 n,q0 n0) = −∂qn0S(q0 n,q0 n0),∂qnS(q0 n,q0 n0) = −V00(r0 n,n0), where r0 n,n0=|q0 n−q0 n0|. With introduction of the particle displacement function un=qn−q0 n , we derive linear coupled oscillator equations: Mn¨ un=−ω2 0un−∑ n0∈Nn V00(r0 n,n0)(un−un0),∀n=1, . . . , N, (63) where we have added the linearized force from the on-site potential (8) . Equations (63) are a linearized lattice Hamiltonian system of (4) in Newton’s equation form with the Hamiltonian: Hlin lat = N ∑ n=1 1 2Mn ˙ u2 n+1 2ω2 0u2 n+1 2∑ n0∈Nn 1 2V00(r0 n,n0)(un−un0)2!. (64) Since N n={n− 1, n+ 1 } and r0 n,n+1=r0 n,n−1=σ , the linear Equation (63) simplifies to: Mn¨ un=−ω2 0un−V00(σ)(2un−un+1−un−1), (65) which is also the linearized lattice equation of the model example of Section 2.7. 3.2. The Lattice Dispersion Relation Valuable spectral information of the linear lattice wave solutions of (65) can be obtained from the dispersion relation when all masses are assumed to be equal. Alternatively, normal mode solutions, i.e., single frequency linear wave solutions, can be obtained for periodic primitive lattice cells with particles of different masses. To derive the dispersion relation for the linearized lattice Equations (65) we make an ansatz of harmonic wave solutions: un=Aexp(i(−ωt+kn)),A∈R6=0, (66) where ω=ˆ ω/√M∈R is the lattice frequency, and k∈R is its wavenumber. Inserting (66) into (65), we obtain the dispersion relation: −ˆ ω2A=−ω2 0A−V00(σ)A(2−exp(ki)−exp(−ki)), ˆ ω2=ω2 0+V00(σ)(2−2 cos(k)), ˆ ω2=ω2 0+4V00(σ)sin2k 2, Mathematics 2022,10, 3460 17 of 43 Mω2=ω2 0+4V00(σ)sin2k 2. (67) From our assumptions of the interaction potential V in Section 2.1, i.e., V00(σ)> 0, we can obtain upper and lower bounds on the lattice frequency ω: ω2 0≤Mω2≤ω2 0+4V00(σ), (68) where the case ω0= 0 leads to acoustic-like dispersion relations, while models with ω0> 0 lead to optic-like dispersion relations [ 31 ]. With respect to the physical scales set in Section 2.5 , we have that M[ M ] , ω2[ T −2] , ω2 0[ MT −2] , and V00[ E σ−2] . Thus, in the dimensionless form, the dispersion relation (67) reads: ω2=ω2 0+4V00(1)sin2k 2,ω2 0≤ω2≤ω2 0+4V00(1). For the model example of Section 2.7, we have that ω2=4π2U0+288V0sin2k 2. 3.3. The Charge Dispersion Relation The charge dynamical Equations (25) or in dimensionless form (32) and (33) – (34) are already linear with respect to the variables cn(t) or an(t) and bn(t) , respectively. To derive the charge dispersion relation, we consider the lattice at its equilibrium, i.e., Equation (25) in the following form: i¯h˙ cn=E0cn− N ∑ n0=1 n06=n J0 n,n0cn0, (69) where J0 n,n0=J(q0 n,q0 n0) = J0exp −αr0 n,n0 σ! and we have assumed that E0=En(q0 1 , . . . , q0 N)∈R for all n= 1, . . . , N . Notice that, for generality, we have not restricted Equation (69) to only the fixed neighbors model. As already stated in Section 2.1, at the lattice equilibrium q0 n=σ(n− 1 ) for any n06=n there is different n00 6=n such that r0 n,n0=r0 n,n00 . Thus, the Equation (69) can be more conveniently rewritten in the following form: i¯h˙ cn=E0cn−∑ (n0,n00) J0 n,n0(cn0+cn00), (70) where the sum is over all such index pairs (n0 , n00) . Then, we proceed with an ansatz of harmonic wave solutions, see also Section 3.2: cn=Bexp(i(−ωct+kcn)),B=±1 √N, (71) where the amplitude B is chosen such that the probability conservation (16) holds, ωc∈R is the charge frequency and kc∈R is its wavenumber. For any fixed n and associated index pair (n0,n00), there exists an index jn0∈Z6=0such that cn0=Bexp(i(−ωct+kc(n+jn0))) (72) and cn00 =Bexp(i(−ωct+kc(n−jn0))). (73) Mathematics 2022,10, 3460 18 of 43 Inserting expressions (71)–(73) into (70), we obtain the dispersion relation: ¯hωcB=E0B−∑ (n0,n00) J0 n,n0B(exp(kcjn0i) + exp(−kcjn0i)), ¯hωc=E0−2∑ (n0,n00) J0 n,n0cos(kcjn0), ¯hωc=E0−2 K ∑ m=1 J(mσ)cos(mkc), (74) which is real and the frequencies ωc can have positive as well as negative values. With K=|{(n0 , n00)}| , we identify the total number of index pairs (n0 , n00) , which is independent of the index n . For the general case (69) and a finite periodic lattice of N particles, we have that K=N− 1. From the assumptions of the function J together with definition (23) , we find that J0exp(−α) = J(σ)>J(2σ)>··· >J(Kσ)>0, J06=0, and can estimate the charge frequency ωcfrom the dispersion relation (74), i.e., |ωc| ≤ ¯h−1 |E0|+2 K ∑ m=1 J(mσ)!=¯h−1 |E0|+2J0 K ∑ m=1 exp(−α)m!=¯h−1 |E0|+2J0exp(−α)1−exp(−Kα) 1−exp(−α)!. For example, in the one-dimensional fixed close neighbor model of Section 2.7, when N n={n− 1, n+ 1 } for all n , there is only one index pair (n0 , n00)=(n− 1, n+ 1 ) , such as K=|{(n0,n00)}|=1, the dispersion relation (74) reads: ωc=¯h−1(E0−2J0exp(−α)cos(kc))(75) and |ωc| ≤ ¯h−1(|E0|+2J0exp(−α)). In addition, in the rescaled dimensionless form of Section 2.5, the dispersion relation (75) and estimate transform to: ωc=τ−1(E0−2J0exp(−α)cos(kc)),|ωc| ≤ τ−1(|E0|+2J0exp(−α)). Comparing the charge dispersion relation (75) to the lattice dispersion relation (67) multiple time scales become evident, especially, when |E0|¯h−1 1 or J0 1. This observation is important for designing efficient numerical integrators for solving coupled charge-lattice Equations (43)–(46). 3.4. The Charge Equilibrium States After coupling both dynamics and since the Hamiltonian (26) (and (36) ) also depends on the lattice particle positions qn , we have additional force term −∂qHa,b(q , a , b) entering Equation (45) . The force term is explicitly expressed in (28) (and (40) ), where from the definition (23), it follows that: Z(qn,qn0) = −∂qnJ(qn,qn0) = −1 rn,n0J0(rn,n0)(qn−qn0), (76) and, in addition, we find that: W(qn,qn0):=−∂qnVc(|qn−qn0|) = −1 rn,n0V0 c(rn,n0)(qn−qn0). Mathematics 2022,10, 3460 19 of 43 Recall that at the lattice equilibrium (q0 , p0) , where q0 n=σ(n− 1 ) and p0=0 , all lattice Hamiltonian forces in (45) add to zero, i.e., −∂qHq lat(q0) = 0, and for any n0there exists an index n00 such that ∂qnJ(q0 n,q0 n0) = −∂qnJ(q0 n,q0 n00), (77) and also for the potential Vc , the relation (6) holds. In addition, with assumptions (20) and (21) and with any eigenmode charge solution given by (71) , conveniently written in the vector form: (a0(t) , b0(t)) = ±(¯a0(t) , ¯ b0(t))/√N , for any wavenumber kc and frequency ωcsatisfying the dispersion relation (74), we obtain that (in (28)) −∂qHa,b(q0,a0(t),b0(t)) = 0. To observe that, notice (77) and the relations (37) applied to (71) and taking into account (72) and (73), i.e., c∗ n(t)cn(t) = 1 N, Re(c∗ n(t)cn0(t))=Re(c∗ n(t)cn00(t))=1 Ncos(kcjn0), where the expressions are time-independent but depend on the charge wavenumber kc . This observation motivates to explore linear properties at the lattice equilibrium of the subsystem: ˙q=M−1p, (78) ˙p=−∂qHq lat(q)−e∂qHa,b(q, ¯a0(t),¯ b0(t)),e=1 N, (79) which we discuss in the following section. Notice that the additional force term in (79) is of order e , which tends to zero when N→∞ . Thus, for a large lattice of particles the additional force term is negligible and the contributions to the lattice dispersion relations of Section 3.2 will also be of order e . Importantly, accordingly scaled (to satisfy probability conservation (16) ) linear combinations of different charge eigenmode solutions (71) are not charge equilibrium states of the Hamiltonian dynamics (43) – (46) at the lattice equilibrium, and the force term in (79) will be not of order e but will contain time-periodic forcing coefficients of different frequencies. Such induced forcing by linear combinations of charge eigenmode solutions may lead to resonance in lattice vibrations and localization effects due to nonlinearity. The study of this phenomenon is left for future research. 3.5. Semi-Coupled Lattice–Charge Dispersion Relation In this section, we discuss linearization and spectral properties of the fixed close neighbor model of the subsystem (78) and (79) . Considering only the e order term the equations can be written in the following form: Me−1¨ qn=−QU0 c(qn) + 2Q∑ n0∈Nn W(qn,qn0)−2∑ n0∈Nn Z(qn,qn0)cos(kc·jn0), (80) for all n= 1, . . . , N , where all particle masses are assumed to be constant, i.e., M=Mn for all n , such that the dispersion relation for harmonic wave solution can be derived. In addition, we will also assume that the on-site potential −Uc has a harmonic approximation in the form (8) with the isolated oscillator frequency ωc0per unit mass when Vc=0. Linearization of the force terms Z(qn , qn0) and W(qn , qn0) follows the linearization of the lattice force (62) of Section 3.1, i.e., ∂qnZ(q0 n,q0 n0) = −J00(r0 n,n0),∂qnW(q0 n,q0 n0) = −V00 c(r0 n,n0), Mathematics 2022,10, 3460 20 of 43 respectively. Then the linearized equations of (80) for the particle displacements unread: Me−1¨ un=Qωc2 0un−2Q∑ n0∈Nn V00 c(r0 n,n0)(un−un0) + 2∑ n0∈Nn J00(r0 n,n0)(un−un0)cos(kcjn0) =Qωc2 0un−2QV00 c(σ)−J00(σ)cos(kcjn0)(2un−un+1−un−1), since N n={n− 1, n+ 1 } and r0 n,n+1=r0 n,n−1=σ . Thus, together with the lattice dispersion relation of Section 3.2, we derive the semi-coupled charge–lattice dispersion relation in the following form: Mω2=ω2 0−eQωc2 0+4V00(σ) + 2eQV00 c(σ)−J00(σ)cos(kc)sin2k 2, or in the dimensionless form: ω2=ω2 0−eQωc2 0+4V00(1) + 2eQV00 c(1)−J00(1)cos(kc)sin2k 2, (81) where the contribution to the lattice dispersion relation (67) is of order e . For the model example of Section 2.7, the dispersion relation (81) simplifies even further: ω2=ω2 0+4V00(1)sin2k 2−eQωc2 0−8eJ00(1)cos(kc)sin2k 2 =4π2U0+288V0sin2k 2−eQU0 c−8eα2J0exp(−α)cos(kc)sin2k 2 ≤4π2U0+288V0+e(U0 c+8α2J0exp(−α)). With the analysis above, we can conclude that approximate linear wave dynamics of the Hamiltonian system (43) – (46) when N 1 is predominantly characterized by decoupled lattice and charge dispersion relations of Sections 3.2 and 3.3, which has to be taken into account for the absolute linear stability considerations of splitting methods discussed in the following section. 4. Structure-Preserving Splitting Methods In this section, we propose symplecticity-preserving symmetric splitting methods of semi-classical canonical Hamiltonian dynamics (43) – (46) , which we restate in the following form: ˙q=M−1p, (82) ˙a=D(q)b+L(q)b, (83) ˙p=F(q) + G(q,a,b), (84) ˙ b=−D(q)a−L(q)a, (85) where Π(q) = D(q) + L(q) , D(q) = diag(Π(q)) and L(q)T=L(q) , which follows from Π(q)T=Π(q). The system of differential Equations (82) – (85) is highly nonlinear and, thus, it is very desirable to obtain an explicit numerical integration scheme. In addition, from the linear analysis of Section 3, we observed that linear dynamics at the lattice equilibrium is (predominately) described by two decoupled lattice and charge dispersion relations (67) and (74) , respectively, which demonstrate the different time scales of lattice and charge dynamics. This naturally motivates to split the lattice dynamics from the charge dynamics in numerical integration, while at the same time attempting to preserve as many as possible structural properties of the Hamiltonian system (82) – (85) stated in Section 2.6. Such splitting method approach may further allow us to investigate and consider multiple time-stepping schemes, such as the impulse method [36,37], to further improve numerical integration of Mathematics 2022,10, 3460 21 of 43 the multi-scale dynamics, and higher-order methods [ 39 , 40 ], to improve numerical accuracy of nonlinear wave simulations, which is left for future work. In Section 3.3, we also observed an effect of the energy density values En divided by τ on the charge linear wave frequencies, see the charge dispersion relation (74) . This motivates us to split the symmetric matrix Π(q) into the sum of two matrices D(q) and L(q) consisting of diagonal and off-diagonal entries, respectively, where the matrix D(q) contains exactly En/τ values on the diagonal, i.e., Dnn(q) = En(q1 , . . . , qN)/τ for all n= 1, . . . , N . Such splitting naturally occurs by splitting the kinetic from the potential energy in splitting methods for solving continuous Schödinger equations. Importantly, as illustrated below, for given values of q , the Hamiltonian split lattice–charge dynamics associated with D(q)terms can be efficiently solved analytically. 4.1. Splitting of Lattice–Charge Dynamics With the stated motivation above, we consider the following pieces of the right-hand side vector field of the dynamics (82)–(85): ˙     q a p b    =    M−1p 0 0 0     | {z } Q +    0 0 F(q) 0     | {z } P | {z } L +    0 D(q)b GD(q,a,b) −D(q)a     | {z } D +    0 L(q)b GA(q,b) 0     | {z } A +    0 0 GB(q,a) −L(q)a     | {z } B | {z } W | {z } C , (86) where we have assigned letters Q, P, D, A, and B to identify each piece. We have split the right hand sides of (82) – (85) in such a way that each piece can be solved exactly with associated analytic flow maps φQ t,φP t,φD t,φA t, and φB t, for which explicit representations are stated in Table A1 of Appendix A.1. For the analytic solution of the piece D see the explanations and Equations (87)–(89) below. Each piece Q, P, D, A, and B represents Hamiltonian dynamics with the Hamiltonian obtained from splitting the total Hamiltonian (42) into the following terms, i.e., H(q,a,p,b) = HQ(p) + HP(q) + HD(q,a,b) + HA(q,b) + HB(q,a), where HQ(p) = Hp lat(p),M−1p=∂pHQ(p), HP(q) = Hq lat(q),F(q) = −∂qHP(q), HD(q,a,b) = 1 2aTD(q)a+1 2bTD(q)b,GD(q,a,b) = −∂qHD(q,a,b), HA(q,b) = 1 2bTL(q)b,GA(q,b) = −∂qHA(q,b), HB(q,a) = 1 2aTL(q)a,GB(q,a) = −∂qHB(q,a). In addition, with letters L and C, we have identified pieces of (semi-decoupled) lattice and charge dynamics with associated lattice and charge Hamiltonians Hlat and Ha,b , respectively. The letter W represents the combination of pieces A and B with the Hamiltonian: HW(q,a,b) = HA(q,b) + HB(q,a) and the force term GW(q,a,b) = GA(q,b) + GB(q,a). Mathematics 2022,10, 3460 22 of 43 Thus, all analytic flow maps associated with pieces Q, P, L, D A, B, W, and C are Hamiltonian and symplectic (49) , which implies phase volume preservation as well. In addition, they are also time-reversibility (50) preserving. Apart from the flow maps φA t and φB t , all remaining flow maps also preserve rotational invariance (51) and conserve the probability (16). Notice that in the pieces C, D, and W, equations for a and b can also be written down in complex form: i˙c=Π(q)c, i˙c=D(q)c, i˙c=L(q)c, respectively, where c=a+b i. Since the matrix D(q) is diagonal, for any given q the second equation above can be easily solved analytically with the (complex) flow map: φE(q) t=exp(−D(q)ti),φE(q) tn=exp(−Dnn(q)ti),∀n=1, . . . , N, (87) whereas computing analytic solutions of the remaining two equations poses a challenge for large systems of equations when N 1, but efficient linear algebra solvers can be used if, e.g., an implicit charge-preserving integrator is applied to solve C or W dynamics, since the system of equations for c is linear and both matrices Π(q) and L(q) are symmetric and potentially sparse as well. An important observation is that with the application of the flow map φE(q) t(87) , we find that GD(q, Re(φE(q) t(a+bi)), Im(φE(q) t(a+bi))) = GD(q,a,b)(88) for all values of q , a , b , and t , e.g., see Equation (28) or (40) . Thus, the piece D can be efficiently solved exactly with the flow map φD texpressed in the following explicit form: Q=q, A=Re(φE(q) t(a+bi)), P=p+tGD(q,a,b), B=Im(φE(q) t(a+bi)), (89) where, for a given state (q , a , p , b)T , we have found a new system’s state (Q , A , P , B)T at time t. In what follows, with φh and ψh , we identify exact and numerical flow maps of autonomous differential equations advancing a given state (q , a , p , b)T in time to a new state (Q , A , P , B)T with the time step h> 0. The numerical flow map ψh is symplectic if it satisfies (49) for all h> 0. Recall that composition of symplectic maps is also symplectic [ 37 ]; thus, we can obtain a symplectic numerical method just by the composition of the exact symplectic flow maps above. It is well known that symmetric numerical methods preserve time-reversibility (50) [ 37 ]. Thus, the construction of symmetric numerical methods, which are ρ-reversible, i.e., ρ◦ψh=ψ−1 h◦ρ,∀h>0, for numerical integration of (82) – (85) is highly desirable. The numerical method is called symmetric if ψ† h=ψh, where † indicates the adjoint method of ψhdefined as: ψ† h=ψ−1 −h, with respective properties [37]: ψ† h†=ψh,(ψh◦¯ ψh)†=¯ ψ† h◦ψ† h,ψh◦ψ† h†=ψh◦ψ† h, Mathematics 2022,10, 3460 23 of 43 where ψh and ¯ ψh are two different numerical flow maps, and the last relation tells us that the numerical method composed with its adjoint method, or vice versa, is symmetric. Notice that the property ψ† h=ψh holds for exact flows maps φt , see (48) . The symmetry of the method implies that exchanging (q,a,p,b)↔(Q,A,P,B)and h↔ −h leaves the method unaltered. Similarly to the composition of symplectic methods, the composition of symmetric flow maps is also symmetric. In addition, all symmetric methods are of even order [ 37 ]. Note that the preservation of rotational invariance (51) and charge probability (16) will depend on the choice of the exact flow maps above for the construction of a numerical method based on the flow map composition. In the following section, we describe semi-implicit splitting methods for the system (82) – (85) that are symplectic, symmetric, second-order, preserve time-reversibility, and conserve exactly the charge probability. Moreover, they require only one force F(q) and matrix Π(q) elements evaluation and at most two force term G(q , a , b) calculations per time step. 4.2. Semi-Implicit Methods with Exact Charge Probability Conservation Modeling charge particle transport by nonlinear lattice excitations is very important to preserve the conservation of the charge probability (16) , which is a quadratic invariant. We propose four semi-implicit splitting methods that are symplectic, symmetric, and also conserve the charge probability exactly. They are based on the symplectic and symmetric implicit midpoint method’s [ 37 ] solution of the pieces C and W. All methods require only one force F(q) and matrix Π(q) element evaluations, as well as one solution of a linear system of the charge Equations (83) and (85) per time step. The first method reads: ψPQCQP h=φP h/2 ◦φQ h/2 ◦ψC h◦φQ h/2 ◦φP h/2, where the symmetry follows from ψPQCQP† h=ψPQCQP h with ψC† h=ψC h and symplecticity follows from the composition of symplectic flow maps. The charge probability conservation (16) follows from the application of the implicit midpoint rule, i.e., the map ψC h: Q=q, C=I+ih 2Π(q)−1I−ih 2Π(q)c,c=a+bi, A=Re(C),B=Im(C), P=p+hGq,a+A 2,b+B 2, (90) which requires only one force Gq,a+A 2,b+B 2 calculation per time step. Notice that in the method PQCQP during the application of the lattice flow maps φQ h/2 and φP h/2 the charge variable cis kept constant. Thus, since η◦ψC h=ψC h◦η,∀h>0, then also η◦ψPQCQP h=ψPQCQP h◦η,∀h>0. Composition φP h◦φQ h is nothing more than the symplectic Euler method [ 37 ] applied to the decoupled lattice equations, i.e., Q=q+hM−1p, Mathematics 2022,10, 3460 24 of 43 P=p+hF(Q), and, disregarding charge equations, the numerical method based on the Strang splitting with the flow map ψLv h=φP h/2 ◦φQ h/2◦φP h/2 ◦φQ h/2†=φP h/2 ◦φQ h◦φP h/2, i.e., ¯p=p+h 2F(q), Q=q+hM−1¯p, P=p+h 2F(Q), is the symplectic second-order velocity Verlet method [ 36 ], while the numerical method with the flow map ψLp h=φP h/2 ◦φQ h/2†◦φP h/2 ◦φQ h/2=φQ h/2 ◦φP h◦φQ h/2 is the so-called symplectic second-order position Verlet method, which are staple methods in classical molecular dynamics with good energy conservation properties in long-time simulations due to the existence of the modified differential equation, which is also Hamiltonian. Interested readers in theoretical numerical analysis aspects of this subject are referred to [ 37 ]. Thus, we will refer to the method PQCQP as the velocity semi-implicit method and to the method defined by the flow map ψQPCPQ h=φQ h/2 ◦φP h/2 ◦ψC h◦φP h/2 ◦φQ h/2 as the position semi-implicit method. Note that we do not consider methods CPQPC and CQPQC, since those methods would require two solutions of the linear system of equations for the charge equations, see (90) , per time step. Recall that the time step restriction for the absolute linear stability of the symplectic Euler and Verlet methods is [36] |hω| ≤ 2, (91) where, for the lattice dynamics, frequency ω is given by the lattice dispersion relations (67) and (81). It is well known that the implicit midpoint method is an unconditionally stable method for stable linear dynamical systems. Despite this, large values of En/τ , as indicated by the charge dispersion relation (74) , may require to use small time steps to acquire the necessary accuracy, i.e., to resolve high-frequency oscillations of the charge particle, which is demonstrated in numerical results Section 5. Thus, we propose splitting the charge dynamics piece C into two pieces, D and W, where the piece D can be solved exactly as already indicated in (89) . Note that the analytic charge solution given by (87) conserves charge probability (16) as well as the properties (50) and (51). The remaining two proposed semi-implicit splitting methods are of the following form: ψPQDWDQP h=φP h/2 ◦φQ h/2 ◦φD h/2 ◦ψW h◦φD h/2 ◦φQ h/2 ◦φP h/2, ψQPDWDPQ h=φQ h/2 ◦φP h/2 ◦φD h/2 ◦ψW h◦φD h/2 ◦φP h/2 ◦φQ h/2, Mathematics 2022,10, 3460 25 of 43 where the map ψW hof the implicit midpoint rule is Q=q, C=I+ih 2L(q)−1I−ih 2L(q)c,c=a+bi, A=Re(C),B=Im(C), P=p+hGWq,a+A 2,b+B 2. (92) As for the methods PQCQP and QPCPQ, we have that ψPQDWDQP† h=ψPQDWDQP h,η◦ψPQDWDQP h=ψPQDWDQP h◦η, ψQPDWDPQ† h=ψQPDWDPQ h,η◦ψQPDWDPQ h=ψQPDWDPQ h◦η, for all h>0. Comparing methods PQDWDQP and QPDWDPQ to methods PQCQP and QPCPQ, we still require only one force term GW(q , a , b) evaluation, but two force term GD(q , a , b) calculations with the same value of q per time step, unless the force GD(q , a , b) is equal to zero when all Enare constant, see Appendix A.2. In our implementation of the methods PQDWDQP and QPDWDPQ in MATLAB with double-precision floating-point arithmetic, we observed a gradual linear increase in relative error in charge probability conservation (16) due to the accumulation of round-off errors. To address this problem, we propose the normalization of the charge variable c after each solution of the piece D in (87), i.e., ¯c=exp(−tD(q)i)c,c=a+bi, ¯c=r2τ S¯c,S=aTa+bTb, (93) thus, we obtain ¯aT¯a+¯ bT¯ b= 2 τ to machine double-precision. We identify these methods with the flow maps ψPQ ¯ DW ¯ DQP h and ψQP ¯ DW ¯ DPQ h . We do not add this rescaling of c values after the implicit midpoint steps (90) and (92) , where we did not observe a linear increase in relative error. In Appendix A.2, we state all semi-implicit methods in explicit forms, see Tables A2 and A3. In the numerical results in Section 5, we compare the proposed semi-implicit splitting methods above, i.e., PQCQP, QPCPQ, PQDWDQP, QPDWDPQ, PQ ¯ D W ¯ D QP and QP ¯ D W ¯ D PQ, with the standard fourth-order explicit Runge-Kutta method (RK4) [ 37 ], which is neither a symplectic, symmetric, nor a charge-conserving method, although RK4 is of order four. We identify each method’s strengths and weaknesses, in addition to in comparison to the explicit symplectic splitting methods derived in the following section. 4.3. Explicit Splitting Methods There are multiple ways we can construct fully explicit symplecticity-preserving symmetric splitting methods for Hamiltonian dynamics (82) – (85) considering pieces in (86) . We propose four symmetric Strang splitting methods derived from the composition of two first-order, one-step splitting methods: ψPQABD h=φP h◦φQ h◦φA h◦φB h◦φD h,ψPQDAB h=φP h◦φQ h◦φD h◦φA h◦φB h, (94) which are formed from the composition of exact symplecticity-preserving flow maps. Note that φA h◦φB h6=φB h◦φA h . Since, structurally, pieces A and B are equivalent (see (86) ), in both methods above, we only consider the composition φA h◦φB h . In addition, in numerical simulations (not shown), we did not observe qualitative differences in numerical solutions Mathematics 2022,10, 3460 32 of 43 them accordingly such that (16) holds. With these initial conditions, the charge probability is spread across the whole computational domain, and the lattice solutions are dominated by linear nonlocalized phonon waves. In the nonlinear regime, we generate eleven different initial initial conditions using patterns (100) and (101), where γj=0.4 +jhγ,hγ=0.05, j=0, . . . , 10, i.e., with γ∈[ 0.4, 0.9 ] . With these initial conditions, we can generate nonlinear and localized solutions of charge transfer by mobile discrete breathers, e.g., see Figure 1. In addition, in both Figures 5and 6with the reduction in the time step by a factor of 10, we have also reduced the y -axis scale by a factor of 100 to illustrate the second-order approximation error. We have also removed the probability errors for the semi-implicit splitting methods PQ ¯ D W ¯ D QP and QP ¯ D W ¯ D PQ from the probability error plots (bottom rows), since the total probability is conserved up to the machine precision, see Figure 3c,f. In Figure 5a–c, it can clearly be seen that the Hamiltonian is conserved up to the second-order and that the errors of the semi-implicit methods PQ ¯ D W ¯ D QP and QP ¯ D W ¯ D PQ coincide very well with the errors of the explicit splitting methods DBAQPQABD and PQDABADQP for all considered ratio E0/τ values. Interestingly, as the value of E0/τ increases, we observe significant differences in the Hamiltonian errors between the explicit methods. Such a large discrepancy cannot be seen in Figure 4b,e, since the results there were demonstrated for E0/τ< 10 4 . In Figure 4c,f, we observed that explicit methods DBAQPQABD and PQDABADQP can have much smaller probability conservation errors compared to the methods PQABDBAQP and BADQPQDAB. This observation is more evident in Figure 5d–f, where, in fact, the PQABDBAQP and BADQPQDAB methods’ Hamiltonians and probability errors are indistinguishable. Notice that the probability errors of the methods DBAQPQABD and PQDABADQP precisely obey the second-order reduction in error only for E0/τ values satisfying the condition (102) . Despite that, the DBAQPQABD and PQDABADQP methods’ probability errors are significantly smaller compared to the errors of the methods PQABDBAQP and BADQPQDAB. In addition, the PQDABADQP method’s produced errors are slightly smaller compared to errors of the method DBAQPQABD. The linear regime results, Figure 5, do not directly reflect the nonlinear regime results, see Figure 6. First of all, in the nonlinear regime, we observe larger Hamiltonian and probability errors comparing y -axis scales in Figures 5and 6. In all plots of Figure 6, we can see that two explicit methods DBAQPQABD and PQDABADQP, overall, outperform other two explicit methods PQABDBAQP and BADQPQDAB, especially for larger values of E0/τ and smaller time steps, with method PQDABADQP having slightly smaller errors compared to the method DBAQPQABD. We can also observe slightly smaller errors for the semi-implicit velocity method PQ ¯ D W ¯ D QP compared to the position method QP ¯ D W ¯ D PQ, see Figure 6a–c, which can also be seen in Figure 3b,e, but not in the linear regime results of Figure 5a–c. Mathematics 2022,10, 3460 33 of 43 (a) (b) (c) (d) (e) (f) Figure 5. Hamiltonian and total probability conservation errors of the splitting methods: PQABDBAQP, BADQPQDAB, DBAQPQABD, PQDABADQP, PQ ¯ D W ¯ D QP, and QP ¯ D W ¯ D PQ. Top ( a – c ): maximal relative error of the Hamiltonian (42) for t∈[ 0, 100 ] . Bottom ( d – f ): maximal relative error of the total charge probability (16) for t∈[ 0, 100 ] , where the results for semi-implicit methods are excluded. Left ( a , d ): time step h= 0.01. Middle ( b , e ): time step h= 0.001. Right ( c , f ): time step h= 0.0001. The Hamiltonian and probability errors are averaged over eleven numerical simulations with randomly generated initial conditions from the uniform distribution U(−0.01, 0.01). To summarize, after performing the computational analysis of the proposed splitting methods’ convergence and the conserved quantities’ conservation properties, we advocate the velocity semi-implicit, exactly charge-conserving method PQ ¯ D W ¯ D QP over the position method QP ¯ D W ¯ D PQ due to the smaller errors in the Hamiltonian conservation. From all the explicit splitting methods, we advocate the method PQDABADQP, which showed overall smaller errors and conserves the total probability to a high degree in long-time simulations of charge transfer by discrete breathers as demonstrated in the following section. Mathematics 2022,10, 3460 34 of 43 (a) (b) (c) (d) (e) (f) Figure 6. Hamiltonian and total probability conservation errors of the splitting methods: PQABDBAQP, BADQPQDAB, DBAQPQABD, PQDABADQP, PQ ¯ D W ¯ D QP, and QP ¯ D W ¯ D PQ. Top ( a – c ): maximal relative error of the Hamiltonian (42) for t∈[ 0, 100 ] . Bottom ( d – f ): maximal relative error of the total charge probability (16) for t∈[ 0, 100 ] , where results for semi-implicit methods are excluded. Left ( a , d ): time step h= 0.01. Middle ( b , e ): time step h= 0.001. Right ( c , f ): time step h= 0.0001. The Hamiltonian and probability errors are averaged over eleven numerical simulations with initial conditions (100) and (101) and different values of γ∈[0.4, 0.9]. 5.3. Charge Transfer by Mobile Discrete Breathers In the example in Figure 1, we demonstrate charge transfer by a mobile discrete breather in time until Tend = 50 initiated by the patterns (100) and (101) with γ= 0.6 and computed with the explicit method PQDABADQP and time step h= 0.01. In this section, we demonstrate the long-time charge transfer until Tend = 10 4 with the three proposed splitting methods of Section 4, i.e., with three lattice–charge–lattice methods PQ ¯ D W ¯ D QP, PQDABADQP, and PQABDBAQP. Numerical simulations are performed in a lattice with N= 254 particles, time step h= 0.01, and E0/τ= 1000. We excite the charge transfer with the patterns (100) and (101) set at n∗=64 with γ=0.6. In Figure 7a–c, we plot contours of charge probability function Pn(t) = |cn(t)|2(15) . All three methods illustrate localized charge transport by a mobile discrete breather, where short-time solutions appear indistinguishable. As time progresses, numerical solutions start to differ while producing qualitatively acceptable solutions. Importantly, such results are not possible to obtain with the methods RK4, PQCQP and QPCPQ when h= 0.01, see Figure 2d–f. All three methods produce Hamiltonian (42) conservation errors to an equivalent degree over the whole computational time segment [ 0, Tend] , see Figure 7d. The same is not true for the total probability conservation (16) , as can be seen in Figure 7e. Semiimplicit method PQ ¯ D W ¯ D QP exactly conserves the total probability, while the explicit method PQDABADQP has smaller errors compared to the method PQABDBAQP, which is consistent with the results of Figure 6d. Mathematics 2022,10, 3460 35 of 43 (a) (b) (c) (d)(e)(f) Figure 7. Simulation of charge transfer by a mobile discrete breather with three splitting methods: PQ ¯ D W ¯ D QP, PQDABADQP, and PQABDBAQP. The solutions are initiated with patterns (100) and (101) , where γ= 0.6. Other parameter values: h= 0.01, N= 254 and Tend = 10 4 . Top ( a – c ): contours of the charge probability Pn . Bottom: ( d ) relative error of the Hamiltonian (42) at each time step; ( e ) relative error of the total charge probability (16) at each time step; ( f ) participation ratio value at each time step. To gain more information on the charge probability localization in charge transport by discrete breathers, in Figure 7f, we plot the normalized participation ratio (adopted from [19]) Pr(t) = 1 N−1 N N ∑ n=1 Pn(t)2−1!/ N ∑ n=1 Pn(t)2∈[0, 1], as a function of time for all three methods. The participation ratio characterizes the localization of the charge probability, i.e., the larger the value of Pr , the more localized the charge probability is. For example, if the charge probability is completely localized at one site, then Pr= 1, but in the situation when the charge probability is equally spread across the whole lattice, then Pr= 0. With our charge probability’s initial excitation pattern (101) , we have that Pr( 0 ) = 1. From Figure 7f, we can see that as the time progresses, the charge probability becomes more and more delocalized and eventually disperses within the lattice. We have verified this (not shown) by performing simulations with smaller time steps, and the phenomenon persists, which could be attributed to the properties of our example model of Section 2.7. Results of Figure 7f do not imply that charge transfer by mobile discrete breathers may not persist for longer times and distances in different lattice models. In addition, the excitation patterns (100) – (101) do not provide a means to produce exact time-periodic and space-translated traveling wave solutions of charge transfer by discrete breathers. It is still an open question if such solutions even exist. Overall, Figure 7 demonstrates that even the explicit, symplecticity-preserving, symmetric methods, which do not exactly preserve the charge probability (16) , may produce qualitatively good long- Mathematics 2022,10, 3460 36 of 43 time results with good approximate probability conservation properties compared to the exactly charge-probability conserving semi-implicit methods. 6. Discussion In this work, we have constructed structure-preserving numerical methods for an important class of Hamiltonian equations to model charge transfer by intrinsic localized modes in nonlinear crystal lattice models. The importance of the study of charge transfer phenomena by lattice excitations cannot be overstated and traces back to original works by Landau L.D. and Pekar S.I. at the beginning of the last century. We demonstrate that such Hamiltonian equations can be written into a canonical form and address the question of solving the system of equations by direct numerical integration. The coupling of classically modeled lattice equations with the charge dynamics described by quantum mechanics presents serious challenges due to the different time scales of the charge and lattice dynamics, thus also inducing great important difficulties for the construction of stable, accurate, and efficient numerical integration schemes. In this work, we were able to address the problem by effective splitting strategies considering the structural properties of the equations splitting the piece C of charge dynamical equations into pieces D and W or into pieces D, A, and B, as explained in Section 4.1. In the construction of the splitting methods presented in Section 4, we have considered geometric structural properties of the Hamiltonian dynamics, i.e., symplecticity, which also implies phase volume preservation, time-reversibility, rotational invariance, and conserved quantities. We have derived two classes of computationally efficient splitting methods, i.e., semi-implicit and fully explicit methods. All the methods are symplectic, symmetric, second-order, and require a single calculation of the lattice force per time step and a single evaluation per time step of the lattice variable dependent terms in the charge equations. We have split the dynamical equations into pieces that correspond to explicit flow maps, which are Hamiltonian, symplectic, and time-reversible. The composition of flow maps will automatically inherit the same properties. However, to preserve the total charge probability, which is a quadratic invariant, and the rotational invariance of the probability amplitude, we require implicit methods for the pieces corresponding to the charge equations. One convenient method is, for example, the implicit midpoint rule. Since charge equations are linear with respect to the charge variables, we require only one solution of a linear system of equations per time step to exactly preserve the total charge probability and the rotational invariance of charge variables. Performing numerical computational analysis for all the proposed splitting methods, we found that the Hamiltonian is conserved up to second-order in long-time simulations, which we would anticipate from symplectic integrators. The total charge probability is also preserved with high accuracy even with symplecticity-preserving symmetric explicit methods, which do not exactly conserve the total charge probability by design. With both classes of methods, we were able to obtain qualitatively good numerical results of charge transfer by mobile discrete breathers, but it remains to be seen if the approximate charge probability conservation by explicit methods is sufficient to reproduce the physical properties of the phenomenon. Incoherencies may appear when post-processing the numerical results. To address this question we are developing multiple time-stepping methods based on the proposed symplectic explicit methods. In addition, explicit splitting methods do not preserve rotational invariance compared to the proposed semi-implicit methods. At this stage, it is unclear how important it is for the numerical method to preserve rotational invariance in connection to the study of charge transfer, but we conjecture that reproducing the physical properties of the equations is extremely beneficial. We will explore it in future work. It is also unclear at this stage if simulations with large time steps and second-order accuracy are sufficient to provide meaningful answers regarding charge transfer properties, despite the ability to obtain qualitatively good numerical results in long-time simulations. These considerations are of particular importance for spectral Mathematics 2022,10, 3460 37 of 43 representation and theory, including simulations to find numerically exact periodic wave solutions with charge transfer and their stability, see [12,13]. When explicit splitting methods are compared, we find a peculiar sudden increase in accuracy in the total probability conservation at a certain time step, which is dependent on the charge oscillation frequency value E0/τ . We plan to address this phenomenon in a future work from an analytic point of view, in addition to multiple time-stepping or the development of the impulse method to improve probability conservation by explicit methods. Since all the proposed methods are only second-order, we plan to explore higher-order composition methods. It is well known that the number of equations to evaluate significantly increases with higher-order composition methods. An interesting prospect to potentially overcome this is to explore the so-called processing and effective order methods [37,38]. While the computational analysis was performed on an idealized example model, we have ensured sufficient generality by considering the charge energy function En(q1 , . . . , qN) to be dependent on the positions of the lattice particles. The developed methods can be directly applied to higher dimensional lattice models and with long-range interaction potentials. The flexibility of the methods stems from the fact that they can be directly incorporated into splitting methods of thermostated Hamiltonian dynamics and provide a means to construct other methods with multiple time-stepping and higher-order. 7. Conclusions In this article, we have addressed the problem of numerical integration of conservative semi-classical systems, that is, a lattice of atoms, ions, or molecules that are described classically and a charge carrier that is described as a quantum particle within the tightbinding approximation. This approximation implies that the charge carrier states can be expressed in a basis of localized states at each atom or ion. The much faster dynamics of the charge carrier due to the small value of the mass of an electron and the very small value of the Planck constant provide serious challenges to numerical methods that attempt to resolve mathematical models of realistic systems. Importantly, numerical methods that do not conserve energy and the probability one of finding the charge carrier in the whole system can be viewed as unreliable. We have solved the problem expressing the whole system in a canonical form and designing integration schemes based on the splitting methods approach. We have proposed different algorithms with the objective of conserving the underlying properties of the original mathematical model, i.e., symplecticity, time-reversibility, conservation of the total charge probability and the Hamiltonian, rotational invariance of the charge variables, and allowing computations with large time steps and small numerical errors. We have tested these methods producing polarobreathers, i.e., breathers transporting charge probability, in a phenomenological model for silicates. Polarobreathers are good candidates to explain experimental results on charge transport without an electric field. We have demonstrated that explicit splitting methods conserve the Hamiltonian up to the second order and approximately conserve the charge probability in long-time simulations. However, semi-implicit methods conserve the charge probability and the rotational invariance of the charge-probability amplitude exactly. They also have smaller or on par numerical errors compared to the explicit methods. Moreover, the proposed semi-implicit methods allow larger time steps with smaller errors than other conventional methods. Therefore, at least for small models, they are the method of choice for integration of semiclassical systems. Further developments of the explicit methods are highly motivated by the numerical integration of charge transfer in twoand three-dimensional crystal lattice models, where solutions of a large linear system of equations for charge variables per time can become very cumbersome. The methods are not limited to the systems tested as they rely on the generic properties of tight-binding semi-classical models, which are a consequence of the first principles and the laws of physics, particularly classical Hamiltonian dynamics and the Schrödinger Mathematics 2022,10, 3460 38 of 43 equations. We think that the proposed methods will provide efficient computational means to study the phenomenon of charge transfer by nonlinear lattice excitations in physical and biological lattice models with realistic potentials and applications. Author Contributions: Conceptualization, J.B. and J.F.R.A.; methodology, J.B. and J.F.R.A.; software, J.B.; formal analysis, J.B. and J.F.R.A.; investigation, J.B. and J.F.R.A.; writing—original draft preparation, J.B.; writing—review and editing, J.B. and J.F.R.A.; visualization, J.B.; project administration, J.B.; funding acquisition, J.B. and J.F.R.A. All authors have read and agreed to the published version of the manuscript. Funding: J. Baj ¯ ars acknowledges support from the PostDoc Latvia grant No.1.1.1.2/VIAA/4/20/617 funded by the Latvian Council of Science. J.F.R. Archilla acknowledges support from projects MICINN PID2019-109175GB-C22 and Junta de Andalucía US-1380977 and a travel grant from VIIPPITUS 2022. Institutional Review Board Statement: Not applicable. Informed Consent Statement: Not applicable. Data Availability Statement: Not applicable. Conflicts of Interest: The authors declare no conflict of interest. Appendix A In this appendix, we list all the exact flow maps and numerical methods stated in this manuscript in the explicit vector component forms, where a given system’s state (q , a , p , b)T is advanced in time to a new state (Q , A , P , B)T at time t or with the time step h, respectively. Appendix A.1. Exact Flow Maps All exact flow maps φQ t , φP t , φD t , φA t , and φB t of pieces Q, P, D, A, and B, respectively (see Equation (86) ), are listed in Table A1, where the complex flow map φE(q) t is given in (87) . Notice that we have applied the property (88) in the application of the flow map φD t. Table A1. Explicit representations of exact flow maps φQ t,φP t,φD t,φA t, and φB t. φQ tφP tφD tφA tφB t Q=q+tM−1pQ=q Q =q Q =q Q =q A=a A =aA=Re(φE(q) t(a+bi)) A=a+tL(q)b A =a P=p P =p+tF(q)P=p+tGD(q,a,b)P=p+tGA(q,b)P=p+tGB(q,a) B=b B =bB=Im(φE(q) t(a+bi)) B=b B =b−tL(q)a Appendix A.2. Semi-Implicit Numerical Methods In Tables A2 and A3, we list all symplecticity-preserving, symmetric, semi-implicit, exactly charge-conserving, numerical methods PQCQP, QPCPQ, PQDWDQP, QPDWDPQ, PQ ¯ D W ¯ D QP, and QP ¯ D W ¯ D PQ. Explicit representations of the numerical flow maps ψC h and ψW hare stated in (90) and (92), respectively. Mathematics 2022,10, 3460 39 of 43 Table A2. Explicit representation of semi-implicit numerical flow maps ψPQCQP hand ψQPCPQ h. ψPQCQP hψQPCPQ h ¯p=p+h 2F(q)¯q=q+h 2M−1p ¯q=q+h 2M−1¯p¯p=p+h 2F(¯q) c=a+bic=a+bi C=I+ih 2Π(¯q)−1I−ih 2Π(¯q)c C =I+ih 2Π(¯q)−1I−ih 2Π(¯q)c A=Re(C)A=Re(C) B=Im(C)B=Im(C) ˜p=¯p+hG¯q,a+A 2,b+B 2˜p=¯p+hG¯q,a+A 2,b+B 2 Q=¯q+h 2M−1˜p P =˜p+h 2F(¯q) P=˜p+h 2F(Q)Q=¯q+h 2M−1P Table A3. Explicit representations of semi-implicit numerical flow maps ψPQDWDQP h , ψPQ ¯ DW ¯ DQP h , ψQPDWDPQ h, and ψQP ¯ DW ¯ DPQ h. ψPQDWDQP hψPQ ¯ DW ¯ DQP h ¯p=p+h 2F(q)¯p=p+h 2F(q) ¯q=q+h 2M−1¯p¯q=q+h 2M−1¯p ˜p=¯p+h 2GD(¯q,a,b)˜p=¯p+h 2GD(¯q,a,b) ¯a=Re(φE(¯ q) h/2 (a+bi)) ¯a=Re(φE(¯ q) h/2 (a+bi)) ¯ b=Im(φE(¯ q) h/2 (a+bi)) ¯ b=Im(φE(¯ q) h/2 (a+bi)) ¯c=¯a+¯ bi ¯c=¯a+¯ bi ˜c=I+ih 2L(¯q)−1I−ih 2L(¯q)¯c¯c=q2τ S¯c,S=aTa+bTb ˜a=Re(˜c)˜c=I+ih 2L(¯q)−1I−ih 2L(¯q)¯c ˜ b=Im(˜c)˜a=Re(˜c) ˆp=˜p+hGW¯q,¯ a+˜ a 2,¯ b+˜ b 2˜ b=Im(˜c) ˘p=ˆp+h 2GD(¯q, ˜a,˜ b)ˆp=˜p+hGW¯q,¯ a+˜ a 2,¯ b+˜ b 2 A=Re(φE(¯ q) h/2 (˜a+˜ bi)) ˘p=ˆp+h 2GD(¯q, ˜a,˜ b) B=Im(φE(¯ q) h/2 (˜a+˜ bi)) A=Re(φE(¯ q) h/2 (˜a+˜ bi)) Q=¯q+h 2M−1˘pB=Im(φE(¯ q) h/2 (˜a+˜ bi)) P=˘p+h 2F(Q)C=A+Bi C=q2τ SC,S=˜aT˜a+˜ bT˜ b A=Re(C) B=Im(C) Q=¯q+h 2M−1˘p P=˘p+h 2F(Q) Mathematics 2022,10, 3460 40 of 43 Table A3. Cont. ψQPDWDPQ hψQP ¯ DW ¯ DPQ h ¯q=q+h 2M−1p¯q=q+h 2M−1p ¯p=p+h 2F(¯q)¯p=p+h 2F(¯q) ˜p=¯p+h 2GD(¯q,a,b)˜p=¯p+h 2GD(¯q,a,b) ¯a=Re(φE(¯ q) h/2 (a+bi)) ¯a=Re(φE(¯ q) h/2 (a+bi)) ¯ b=Im(φE(¯ q) h/2 (a+bi)) ¯ b=Im(φE(¯ q) h/2 (a+bi)) ¯c=¯a+¯ bi ¯c=¯a+¯ bi ˜c=I+ih 2L(¯q)−1I−ih 2L(¯q)¯c¯c=q2τ S¯c,S=aTa+bTb ˜a=Re(˜c)˜c=I+ih 2L(¯q)−1I−ih 2L(¯q)¯c ˜ b=Im(˜c)˜a=Re(˜c) ˆp=˜p+hGW¯q,¯ a+˜ a 2,¯ b+˜ b 2˜ b=Im(˜c) ˘p=ˆp+h 2GD(¯q, ˜a,˜ b)ˆp=˜p+hGW¯q,¯ a+˜ a 2,¯ b+˜ b 2 A=Re(φE(¯ q) h/2 (˜a+˜ bi)) ˘p=ˆp+h 2GD(¯q, ˜a,˜ b) B=Im(φE(¯ q) h/2 (˜a+˜ bi)) A=Re(φE(¯ q) h/2 (˜a+˜ bi)) P=˘p+h 2F(¯q)B=Im(φE(¯ q) h/2 (˜a+˜ bi)) Q=¯q+h 2M−1PC=A+Bi C=q2τ SC,S=˜aT˜a+˜ bT˜ b A=Re(C) B=Im(C) P=˘p+h 2F(¯q) Q=¯q+h 2M−1P Appendix A.3. Explicit Numerical Methods In Table A4, we list both symplecticity-preserving, symmetric, explicit, lattice–charge– lattice, numerical methods PQABDBAQP and PQDABADQP, while in Table A5, we list both symplecticity-preserving, symmetric, explicit, charge–lattice–charge, numerical methods DBAQPQABD and BADQPQDAB. Mathematics 2022,10, 3460 41 of 43 Table A4. Explicit representations of explicit numerical flow maps ψPQABDBAQP h and ψPQDABADQP h . ψPQABDBAQP hψPQDABADQP h ¯p=p+h 2F(q)¯p=p+h 2F(q) ¯q=q+h 2M−1¯p¯q=q+h 2M−1¯p ¯a=a+h 2L(¯q)b˜p=¯p+h 2GD(¯q,a,b) ˜p=¯p+h 2(GA(¯q,b) + GB(¯q, ¯a))¯a=Re(φE(¯ q) h/2 (a+bi)) ¯ b=b−h 2L(¯q)¯a¯ b=Im(φE(¯ q) h/2 (a+bi)) ˆp=˜p+hGD(¯q, ¯a,¯ b)˜a=¯a+h 2L(¯q)¯ b ˜a=Re(φE(¯ q) h(¯a+¯ bi)) ˆp=˜p+h 2GA(¯q,¯ b) ˜ b=Im(φE(¯ q) h(¯a+¯ bi)) ˜ b=¯ b−hL(¯q)˜a B=˜ b−h 2L(¯q)˜a˘p=ˆp+hGB(¯q, ˜a) ˘p=ˆp+h 2(GA(¯q,B) + GB(¯q, ˜a))ˆa=˜a+h 2L(¯q)˜ b A=˜a+h 2L(¯q)B p0=˘p+h 2GA(¯q,˜ b) Q=¯q+h 2M−1˘p p00 =p0+h 2GD(¯q, ˆa,˜ b) P=˘p+h 2F(Q)A=Re(φE(¯ q) h/2 (ˆa+˜ bi)) B=Im(φE(¯ q) h/2 (ˆa+˜ bi)) Q=¯q+h 2M−1p00 P=p00 +h 2F(Q) Table A5. Explicit representations of explicit numerical flow maps ψDBAQPQABD h and ψBADQPQDAB h . ψDBAQPQABD hψBADQPQDAB h ¯p=p+h 2GD(q,a,b)¯ b=b−h 2L(q)a ¯a=Re(φE(q) h/2 (a+bi)) ¯p=p+h 2GA(q,¯ b) + GB(q,a) ¯ b=Im(φE(q) h/2 (a+bi)) ¯a=a+h 2L(q)¯ b ˜ b=¯ b−h 2L(q)¯a˜p=¯p+h 2GD(q, ¯a,¯ b) ˜p=¯p+h 2GA(q,˜ b) + GB(q, ¯a)˜a=Re(φE(q) h/2 (¯a+¯ bi)) ˜a=¯a+h 2L(q)˜ b˜ b=Im(φE(q) h/2 (¯a+¯ bi)) ¯q=q+h 2M−1˜p¯q=q+h 2M−1˜p ˆp=˜p+hF(¯q)ˆp=˜p+hF(¯q) Q=¯q+h 2M−1ˆp Q =¯q+h 2M−1ˆp ˆa=˜a+h 2L(Q)˜ b˘p=ˆp+h 2GD(Q, ˜a,˜ b) ˘p=ˆp+h 2GA(Q,˜ b) + GB(Q, ˆa)ˆa=Re(φE(Q) h/2 (˜a+˜ bi)) ˆ b=˜ b−h 2L(Q)ˆaˆ b=Im(φE(Q) h/2 (˜a+˜ bi)) P=˘p+h 2GD(Q, ˆa,ˆ b)A=ˆa+h 2L(Q)ˆ b A=Re(φE(Q) h/2 (ˆa+ˆ bi)) P=˘p+h 2GA(Q,ˆ b) + GB(Q,A) B=Im(φE(Q) h/2 (ˆa+ˆ bi)) B=ˆ b−h 2L(Q)A