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

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


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.

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

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

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][15][16][17], electron-vibron breather [17], and solectrons or soliton-elecrons [18][19][20][21][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]. 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 5 for additional numerical results of charge transfer by discrete breathers.
Recently, charge transport in absence of an electric field was observed experimentally in layered silicates, a phenomenon called hyperconductivity [25][26][27][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.

Existing and Proposed Numerical Methods
ILMs have been extensively studied from analytic and numerical points of view in crystal lattice models [12,21,[29][30][31][32][33][34][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 Runge-Kutta 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][37][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.

Outline
The paper is organized in the following way. In Section 2 we derive the general model of lattice-charge dynamics described in canonical Hamiltonian form. At the end 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 6 and 7, respectively.

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 4 and 5.

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]: where K E is the kinetic energy, U E is the on-site potential energy, and V E 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 U E models the periodic energy of the particles in a close-packed line due to the rest of the crystal. The interparticle potential V E models the interaction between the particles in the close-packed line. This interaction is generically repulsive at short distances and attractive at larger ones. q n ∈ Ω ⊂ R is the nth particle's position, where Ω is a bounded computational domain (segment) containing the crystal lattice imposed by the periodic boundary conditions on q n . Alternatively, some form of confining potentials of q 1 and q N could be added to (1) to impose boundary conditions [36]. p n ∈ R and M n > 0 are the nth particle's momentum and mass, respectively. From (1), we derive the system of Hamiltonian equations:q 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., q n and p n are functions of time variable t. A special case of (1) includes the so-called fixed close neighbor interaction model with the Hamiltonian: where N n is a set of n ′ indices specifying a particle's q n fixed neighbors, i.e., if n ′ ∈ N n , then n ∈ N n ′ . 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 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 q n and satisfy the conditions for any n that we present below.

The On-Site Potential
Conditions on U are: where q 0 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 (n ′ , n ′′ ), n ′ ̸ = n ′′ ̸ = n, where also n ′ , n ′′ ∈ N n in (4), such that 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., The definition (7) and conditions (5) imply that we have a harmonic approximation of the on-site potential U at the particle equilibrium: where ω 0 defines the oscillation frequency of isolated oscillators of unit mass for the linearized equations with V = 0.

The Interaction Potential
Modeling the particle interactions with a radial potential V(r), where r = r n,n ′ = |q n − q n ′ | > 0, n ′ ̸ = 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 V ′′ (r) are monotonically decreasing in the interval (0, σ), such as V ′ (r) < 0, and the function V(r) is monotonically increasing for r > σ, such as V ′ (r) > 0, and where V 0 > 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 increases. Examples of such radial interparticle potentials include the Lennard-Jones and Morse potentials: and 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 V ′′ (r) > 0 are assumed to be monotonically decreasing for all r > 0, while the function V ′ (r) < 0 is assumed to be monotonically increasing as well as V(σ) = V 0 and the following conditions hold: Examples of such radial potentials include the Coulomb and Pauli repulsion potentials: 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., In the former case, since V ′ (σ) = 0, we obtain the harmonic approximation of the potential V in the following form: To ensure that the Hamiltonians (1) and (4) at the lattice equilibrium (q 0 n , p 0 n ) = (σ(n − 1), 0) for all n = 1, . . . , N are equal to zero. Without loss of generality, we redefine the interaction potential V in (1) and (4) as: which does not alter the Hamiltonian Equations (2) and (3).

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 q n . In addition, we assume that there is only one quantum state per particle that the charge can occupy.

The Wave Function and the Charge Probability
The extra charge dynamics are described by the wave function, also called state function [41]: where the ket |n⟩ represents the element of the orthonormal basis describing the state in which the charge is bound to the site n and c n (t) ∈ C are time-dependent dimensionless complex coefficients known as probability amplitudes. The adjoint operator of |ϕ(t)⟩ is: where c * n is the complex conjugate of c n . The ket and bra operators, |n⟩ and ⟨n|, satisfy the property ⟨n|n ′ ⟩ = δ n,n ′ , where δ n,n ′ is the Kronecker delta function. Accordingly, we can define a space-discrete, time-dependent wave function to obtain a probability of an electron or hole being at a state n at time t, i.e., which defines the probability function P n . Thus, we have the total probability conservation in the following form: To simplify the notation, in what follows, we omit explicit dependence on t from the variables c n , i.e., c n = c n (t).

The Evolution of the Wave Function: the Schrödinger Equation
The state function (13) is the solution to the discrete Schrödinger equation: whereh is the Planck constant, and H c is the charge Hamiltonian operator: where indices n and n ′ indicate internal quantum numbers of the charge states corresponding to the lattice (1) particles at positions q n and q n ′ , 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 |n⟩ ⟨n ′ | ̸ = |n ′ ⟩ ⟨n| for arbitrary different n and n ′ . Similarly to the lattice Hamiltonian (4), we can consider the charge Hamiltonian operator with fixed close neighbor interactions: for which the following derivations hold as well.

The Local Charge Energy
Smooth multivariable function E n (q 1 , . . . , q N ) ∈ 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 and for any n ′ there exists n ′′ such that In practice, the charge energy function E n can be modeled to consist of a sum of charge on-site and interaction potentials, i.e., respectively, where −U c (note the minus sign) and V c 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 E 1 and E N must be adjusted if non-periodic boundary conditions are used.

The Charge Transition Matrix
The symmetric transition matrix elements J n,n ′ = J(q n , q n ′ ) ≥ 0, i.e., J n,n ′ = J n ′ ,n , of charge transfer between states n and n ′ will depend on the particle interaction distance r n,n ′ and is often modeled with exponential decay: [18,19,21,22]: J(q n , q n ′ ) = J(r n,n ′ ) = J 0 exp −α r n,n ′ σ , lim r n,n ′ →∞ J(r n,n ′ ) = 0, where the dimensionless parameter α > 0 specifies the rate of exponential decay, while σ > 0 defines lattice particle interaction length scale, see Section 2.1. Constant J 0 ≥ 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.

Dynamical Equations for the Charge Probability Amplitude
Notice that the right-hand side of the Schrödinger Equation (17) and the projection of the Equation (17) onto the state n: leads to dynamical equations for the probability amplitude c n , i.e., where n = 1, . . . , N, with the charge classical Hamiltonian where we have applied the adjoint operator (14) to (24) and used the Kronecker delta property of the ket and bra operators.

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 H total = H lat + H c and the associated lattice-charge coupled system of equations: for all n = 1, . . . , N, where we have obtained the additional charge-lattice force term entering Equation (28). Recall that c * n c n ′ + c * n ′ c n = 2Re(c * n c n ′ ).
With the energy function E n given by (22), we have that where the double sum on the right-hand side is equal to and we obtain that Thus, it is also easy to see that we can obtain the model with fixed close nearest neighbor interaction with Equations (27)-(29) just by restricting the sums to the indices n ′ ∈ N n .

Rotational Invariance and Redefinition of the Charge Local-Energy Origin
Notice that under the transformation where E 0 and θ are arbitrary real constants, the charge probability conservation (16) and Equation (28) do not change, since |c n | = |d n | and c * n c n ′ = d * n d n ′ for any n and n ′ , but only the Equation (29) transforms into the following form: and is independent of the constant θ. Thus, if all E n in (22) are constant and equal to E 0 at the lattice equilibrium, then, without loss of generality, in (22), we can redefine the charge interaction potential V c = V c (|q n − q n ′ |) − V c (|q 0 n − q 0 n ′ |), equivalently to the lattice interaction potential in (12), such that at the lattice equilibrium E n (q 0 1 , . . . , q 0 N ) = 0, ∀ n = 1, . . . , N.

Scaled Hamiltonian Equations
To obtain the system of Equations (27)- (29) in dimensionless form we consider the characteristic mass M > 0, length scale σ > 0, time scale T > 0, and energy scale E > 0 values such that the relation E = Mσ 2 T −2 holds. Thus, we define dimensionless variables through the following rescaling: Substituting relations above into the system (27)- (29) and dropping the tildes from the variables, Equations (27)- (29) do not change apart from the Planck constant in the Equation (29), i.e., we obtain where τ =h(ET) −1 > 0 is the dimensionless Planck constant rescaled with respect to the energy and time scales. Thus, depending on the model we can define either energy scale E or time scale T to obtain the value of τ, or, on the other hand, for any value of τ, we can find associated energy and time scales given by the following expressions: Note that the rescaled coupled Equations (27), (28) and (32) are not in the canonical Hamiltonian form, which we derive in the following section.

Canonical Hamiltonian Equations
The charge Equation (32) is in the complex form for the variable c n and can be rewritten in the real form for c n = (a n + ib n )/ √ 2τ, i.e., where the scaling of c n by 1/ √ 2τ does not alter the Equation (32) but will add a factor 1/(2τ) into the Hamiltonian (26), see also (36), and in front of the last two force terms in (28). Note that the total probability conservation (16) in the new variables a n and b n now reads: Thus, we have obtained the charge Hamiltonian Equations (33)-(34) for the canonical variables a n and b n with the classical Hamiltonian derived from (26), i.e., where the factor 1/2 in front of the second sum is trivially justified, since a n a n ′ = a n ′ a n and b n b n ′ = b n ′ b n . To obtain (36) from (26) notice that (a n a n ′ + b n b n ′ ) (37) and the double sum in (26) since Im(c * n c n ′ ) = −Im(c * n ′ c n ), J(q n , q n ′ ) = J(q n ′ , q n ).

Hamiltonian Equations in Real Form
For completeness, before stating the equations in vector form, we state the rescaled canonical Hamiltonian equations in the component form: for all n = 1, . . . , N.
The canonical Hamiltonian Equations (38)-(41) are derived from the Hamiltonian: i.e., the sum of (1) and (36), where (36) is also separable and quadratic in variables a and b, i.e., and the canonical Hamiltonian Equations (38)- (41) in vector form read: where M ∈ R N×N is the nonsingular diagonal mass matrix of masses M n , F(q) + G(q, a, b) ∈ R N and describes the total forces acting on the lattice particles and both Equations (44) and (46) are system of Equations (39) and (41) written in the matrix-vector form with the symmetric system matrix Π(q) ∈ R N×N .

Conservation Laws and Properties of the Semi-Classical Hamiltonian Dynamics
Apart from the energy conservationḢ = 0 and the second invariant, i.e., the charge probability conservation law (35), which in the vector form is expressed as the canonical Hamiltonian dynamics (43)-(46) is also symplectic. Recall that symplecticity also implies phase volume preservation [37] by the flow map ϕ t of the dynamical system (43)-(46). For future reference, we recall the basic properties of analytic flow maps ϕ t of which not all are commonly shared by numerical flow maps, i.e., for all t and s the flow map exists and is invariable, where id is an identity map. With •, we identify the composition operator. We say that the flow map ϕ t of the canonical Hamiltonian dynamics (43)-(46) is symplectic if the following equation holds: for all states (q, a, p, b) T and t the flow map is defined, where I ∈ R 2N×2N is an identity matrix and ∂ϕ t (q, a, p, b)/∂(q, a, p, b) denotes the Jacobian matrix of the flow map ϕ t . In addition to symplecticity, the Hamiltonian dynamics (43)-(46) are also time-reversible under the transformation which leaves the equations unchanged. Accordingly, the flow map ϕ t of (43)-(46) is ρreversible and satisfies the relation and for all t, the map ϕ t and its inverse exist. We already noted the transformation (30), i.e., the rotation of the charge variable c n by a constant angle θ, which leaves the Hamiltonian Equations (43)-(46) unchanged under the following transformation The property can be restated as To summarize, we aim to construct efficient numerical integration methods for solving Hamiltonian Equations (43)-(46), while preserving as many as possible structural properties of the system mentioned above, which we discuss in Section 4.

Model Example: Frenkel-Kontorova with Lennard-Jones Interaction Coupled to a Charge
In this section, we describe a model example of the canonical Hamiltonian Equations (38)-(41) ((43)-(46)) for numerical analysis and simulations in Sections 4 and 5, see also Figure 1.

Lattice Classical Hamiltonian for the Model Example
Without loss of generality, we consider a one-dimensional periodic crystal lattice of N particles of equal masses, i.e., M = M n for all n = 1, . . . , N, coupled with close neighbor interactions. Note that in dimensionless form M = 1 and N n = {n − 1, n + 1}. The model is thought to be a reduced-order model of a row of particles within a three-dimensional layered crystal lattice, where the periodic lattice on-site potential models forces acting on the particle q n from the rest of the crystal. This system is known as the Frenkel-Kontorova model [8,9]. Particle equilibrium positions are q 0 n = n − 1 with σ = 1 in the dimensionless form, such that r 0 = r 0 n+1,n = |q 0 n+1 − q 0 n | = 1 for all n, where q N+1 = q 1 . We find that For the lattice particle interactions, we consider the Lennard-Jones potential (10) and find that Thus, the lattice Hamiltonian (4) is where we have subtracted V(r 0 ) from the potential V, see (12), such that H f ixed lat = 0 at the lattice equilibrium.

Coupling to an Electric Charge in the Model Example
We model the charge energy function E n (22) as: where the charge on-site potential is chosen to be a harmonic potential of the following form: The charge Hamiltonian for the fixed close neighbor example model reads: )(a n a n+1 + b n b n+1 ) .

Hamiltonian Equations for the Model Example
Thus, the canonical Hamiltonian Equations (38)- (41) for the one-dimensional model example above are in the following form: where n = 1, . . . , N. For numerical simulations, the obtained Hamiltonian Equations (54)-(57) can be simplified further by the introduction of the particle displacement function , and r n+1,n = q n+1 − q n = 1 + u n+1 − u n > 0 for all n, as long as the particles do not cross, which would give us a nonphysical solution. The reformulated system of Equations (54)-(57) with respect to the particle displacements u n reads: ,n )(a n a n+1 + b n b n+1 ) + 1 τ J ′ (r n,n−1 )(a n a n−1 + b n b n−1 ), The Hamiltonian system (58)-(61) contains seven parameters to be specified, i.e., U 0 , V 0 , Q, U 0 c , τ, J 0 and α, in addition to the energy values E 0 , which we specify in the numerical results Section 5.

Analysis of the Linearized Equations
In this section, we discuss the linearized equations of the canonical Hamiltonian Equations (38)- (41) and derive the so-called dispersion relations, which relate harmonic wave solution wavenumbers to the wave oscillation frequencies, and obtain time step restrictions for the linear absolute stability of the numerical splitting methods of Section 4. Initially, lattice Equations (2) and (3) are treated separately from the charge Equation (32), thus deriving decoupled lattice and charge dispersion relations, while at the end of this section we discuss particular charge solutions at which coupled dispersion relations can be derived for coupled linearized equations. In what follows, linearization of the lattice equations is performed for the fixed close neighbor interaction model of the Hamiltonian (4) with N n = {n − 1, n + 1}. In addition, to derive the dispersion relations, we assume that all lattice particle masses are equal, i.e., M = M n for all n = 1, . . . , N. Although, they are highly restrictive assumptions, we are still able to obtain very valuable information.

Linearization of Lattice Forces
The Hamiltonian dynamics (4) admits the following linearized equations. The force from the interaction potential V of a particle q n ′ acting on the particle q n in (3) is given by: The linearized force of (62) at the lattice equilibrium is where the terms S(q 0 n , q 0 n ′ ) have been omitted, even if S(q 0 n , q 0 n ′ ) ̸ = 0, since (see (6)) there exists a particle q n ′′ acting on the particle q n with the force S(q 0 n , q 0 where r 0 n,n ′ = |q 0 n − q 0 n ′ |. With introduction of the particle displacement function u n = q n − q 0 n , we derive linear coupled oscillator equations: 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: Since N n = {n − 1, n + 1} and r 0 n,n+1 = r 0 n,n−1 = σ, the linear Equation (63) simplifies to: which is also the linearized lattice equation of the model example of Section 2.7.

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: where ω =ω/ √ M ∈ R is the lattice frequency, and k ∈ R is its wavenumber. Inserting (66) into (65), we obtain the dispersion relation: , From our assumptions of the interaction potential V in Section 2.1, i.e., V ′′ (σ) > 0, we can obtain upper and lower bounds on the lattice frequency ω: 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 . Thus, in the dimensionless form, the dispersion relation (67) reads: For the model example of Section 2.7, we have that

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 c n (t) or a n (t) and b n (t), respectively. To derive the charge dispersion relation, we consider the lattice at its equilibrium, i.e., Equation (25) in the following form: where and we have assumed that E 0 = E n (q 0 1 , . . . , q 0 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 q 0 n = σ(n − 1) for any n ′ ̸ = n there is different n ′′ ̸ = n such that r 0 n,n ′ = r 0 n,n ′′ . Thus, the Equation (69) can be more conveniently rewritten in the following form: where the sum is over all such index pairs (n ′ , n ′′ ). Then, we proceed with an ansatz of harmonic wave solutions, see also Section 3.2: where the amplitude B is chosen such that the probability conservation (16) holds, ω c ∈ R is the charge frequency and k c ∈ R is its wavenumber. For any fixed n and associated index pair (n ′ , n ′′ ), there exists an index j n ′ ∈ Z ̸ =0 such that and c n ′′ = B exp(i(−ω c t + k c (n − j n ′ ))).
Inserting expressions (71)-(73) into (70), we obtain the dispersion relation: which is real and the frequencies ω c can have positive as well as negative values. With K = |{(n ′ , n ′′ )}|, we identify the total number of index pairs (n ′ , n ′′ ), 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 and can estimate the charge frequency ω c from the dispersion relation (74), i.e., 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 (n ′ , n ′′ ) = (n − 1, n + 1), such as K = |{(n ′ , n ′′ )}| = 1, the dispersion relation (74) reads: and |ω c | ≤h −1 (|E 0 | + 2J 0 exp(−α)).
In addition, in the rescaled dimensionless form of Section 2.5, the dispersion relation (75) and estimate transform to: Comparing the charge dispersion relation (75) to the lattice dispersion relation (67) multiple time scales become evident, especially, when |E 0 |h −1 ≫ 1 or J 0 ≫ 1. This observation is important for designing efficient numerical integrators for solving coupled charge-lattice Equations (43)-(46).

The Charge Equilibrium States
After coupling both dynamics and since the Hamiltonian (26) (and (36)) also depends on the lattice particle positions q n , we have additional force term −∂ q H a,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: and, in addition, we find that: W(q n , q n ′ ) := −∂ q n V c (|q n − q n ′ |) = − 1 r n,n ′ V ′ c (r n,n ′ )(q n − q n ′ ).
Recall that at the lattice equilibrium (q 0 , p 0 ), where q 0 n = σ(n − 1) and p 0 = 0, all lattice Hamiltonian forces in (45) add to zero, i.e., −∂ q H q lat (q 0 ) = 0, and for any n ′ there exists an index n ′′ such that and also for the potential V c , 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: (a 0 (t), b 0 (t)) = ±(ā 0 (t),b 0 (t))/ √ N, for any wavenumber k c and frequency ω c satisfying the dispersion relation (74), we obtain that (in (28) To observe that, notice (77) and the relations (37) applied to (71) and taking into account (72) and (73), i.e., where the expressions are time-independent but depend on the charge wavenumber k c . This observation motivates to explore linear properties at the lattice equilibrium of the subsystem:q which we discuss in the following section. Notice that the additional force term in (79) is of order ϵ, 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 ϵ. 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 ϵ 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.

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 ϵ order term the equations can be written in the following form: for all n = 1, . . . , N, where all particle masses are assumed to be constant, i.e., M = M n 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 −U c has a harmonic approximation in the form (8) with the isolated oscillator frequency ω c0 per unit mass when V c = 0. Linearization of the force terms Z(q n , q n ′ ) and W(q n , q n ′ ) follows the linearization of the lattice force (62) of Section 3.1, i.e., ∂ q n Z(q 0 n , q 0 n ′ ) = −J ′′ (r 0 n,n ′ ), ∂ q n W(q 0 n , q 0 n ′ ) = −V ′′ c (r 0 n,n ′ ), respectively. Then the linearized equations of (80) for the particle displacements u n read: since N n = {n − 1, n + 1} and r 0 n,n+1 = r 0 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: or in the dimensionless form: where the contribution to the lattice dispersion relation (67) is of order ϵ. For the model example of Section 2.7, the dispersion relation (81) simplifies even further: 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.

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 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 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 E n 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 E n /τ values on the diagonal, i.e., D nn (q) = E n (q 1 , . . . , q N )/τ 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.

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): 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 q, a, b), In addition, with letters L and C, we have identified pieces of (semi-decoupled) lattice and charge dynamics with associated lattice and charge Hamiltonians H lat and H a,b , respectively. The letter W represents the combination of pieces A and B with the Hamiltonian: , a) and the force term q, a).
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ċ = Π(q)c, iċ = D(q)c, iċ = L(q)c, respectively, where c = a + bi. 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: 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 ϕ 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 t expressed in the following explicit form: (a + bi)), q, a, b), 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., for numerical integration of (82)-(85) is highly desirable. The numerical method is called symmetric if ψ † h = ψ h , where † indicates the adjoint method of ψ h defined as: with respective properties [37]: 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.

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: 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 : 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 c is kept constant. Thus, since 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 −1 p, P = p + hF(Q), and, disregarding charge equations, the numerical method based on the Strang splitting with the flow map is the symplectic second-order velocity Verlet method [36], while the numerical method with the flow map 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 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, 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 E n /τ, 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: where the map ψ W h of the implicit midpoint rule is Q = q, As for the methods PQCQP and QPCPQ, we have that for all h > 0.
Comparing methods PQDWDQP and QPDWDPQ to methods PQCQP and QPCPQ, we still require only one force term G W (q, a, b) evaluation, but two force term G D (q, a, b) calculations with the same value of q per time step, unless the force G D (q, a, b) is equal to zero when all E n are 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., thus, we obtainā Tā +b Tb = 2τ to machine double-precision. We identify these methods with the flow maps ψ PQDWDQP h and ψ QPDWDPQ 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, PQDWDQP and QPDWDPQ, 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.

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: which are formed from the composition of exact symplecticity-preserving flow maps. Note 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 between considering compositions ϕ A h • ϕ B h and ϕ B h • ϕ A h , but we find that it does matter when the map ϕ D h is applied, see Section 5. The adjoint methods of (94) read: Thus, compositions of methods (94) and (95) lead to four fully explicit symplecticitypreserving symmetric methods: where we have applied the property (48) of exact flow maps, i.e., The methods PQABDBAQP (96) and PQDABADQP (98) can be viewed as latticecharge-lattice-splitting methods, while the methods DBAQPQABD (97) and BADQPQDAB (99) can be viewed as charge-lattice-charge-splitting methods. Note that none of the methods (96)-(99) conserve the probability exactly (16) or preserve the rotational invariance (51). All four methods in explicit forms are stated in Tables A4 and A5 of Appendix A.3, respectively.
Notice that the composition ϕ A h • ϕ B h is the symplectic Euler numerical method of the piece W, i.e., the composition in the compact explicit form reads: (q, B) + G B (q, a)), Similarly, the composition ϕ B h • ϕ A h is the adjoint method of the symplectic Euler method above and reads: (98) is the Verlet method applied to the piece W. Thus, the absolute linear stability for the numerical integration of the piece W is provided by the condition (91), where the charge frequency ω c is given by the charge dispersion relation (74) with E 0 = 0, since we have split the piece C into two pieces D and W, where the piece D containing E n /τ values is solved exactly, see Equations (87)-(89).

and the composition ϕ
Similarly to the semi-implicit splitting methods of Section 4.2, the explicit splitting methods (96)-(99) also require only one lattice force F(q) calculation per time step as well as one evaluation of the matrix Π(q) elements, but at most two evaluations of the force terms G A (q, b), G B (q, a), and G D (q, a, b), essentially for the same value of q, see Tables A4  and A5. In the following section, we perform numerical simulations of the model example of Section 2.7 with the explicit splitting methods (96)-(99) and compare them to the semiimplicit integrators from Section 4.2.

Numerical Results
In this section, we perform numerical simulations of the example model of Section 2.7 and investigate the performance of the proposed splitting methods in Section 4. In what follows, for the models (58)-(61), we set the following parameter values: U 0 = 1, V 0 = 0.05, Q = 1, τ = 0.001, α = 15, J 0 /τ = exp(α)/2, and U 0 c = τ, where E 0 or the ratio E 0 /τ vary. With this choice of parameter values in our example model, we can observe the charge transfer by mobile discrete breathers. We were also able to obtain qualitatively similar results with Q = −1. All numerical simulations are performed with N = 64, unless stated otherwise, and we integrate system (58)-(61) in time until T end , which we vary as well as the time step value h.
To excite charge transfer by mobile discrete breathers, see Figure 1, we consider lattice particles in their dynamical equilibrium states while exciting three neighboring particle momenta with the pattern: for any chosen n * = 1, . . . , N (taking into account periodic boundary conditions), γ > 0 leads to a mobile breather solution traveling to the right, while with γ < 0, we obtain a traveling breather solution moving to the left. For example, the result in Figure 1 was obtained with γ = 0.6 and n * = 16. The localized charge solution at the initial time t = 0 is given by the pattern (a n * , b n * ) T = √ τ(1, −1) T (101) set at the middle particle of the pattern (100) and with the remaining values of a n and b n being equal to zero. Notice that the chosen initial pattern (101) satisfies charge probability conservation (16).

Approximation Order and Convergence
In this section, we numerically demonstrate the second-order approximation error and convergence of the proposed splitting methods of Section 4, and the importance of splitting the piece C into two pieces D and W for the semi-implicit splitting methods of Section 4.2 and splitting the piece C into three pieces, D, A, and B, for the explicit splitting methods of Section 4.3. We integrate the model example (58)-(61) in time until T end = 1. For the reference "exact" solution, we consider the numerically computed solution of the explicit method PQDABADQP with the time step h = 10 −7 . In Figures 2-4 we illustrate numerical errors of the solution, the Hamiltonian (42), and the total charge probability (16) for different time step h values, i.e., left, middle, and right columns of the figures, respectively. To measure errors, we compute the maximal absolute error of the numerical solution at the final time T end , while for the Hamiltonian and total charge probability, we consider maximal relative error over the whole computational time segment [0, T end ]. We demonstrate results for two different ratio E 0 /τ values, i.e., when E 0 /τ = 100 (top row) and E 0 /τ = 1000 (bottom row).
In Figure 2, we plot errors of three semi-implicit splitting methods PQCQP, PQDWDQP and PQDWDQP of Section 4.2. Not shown but equivalent results were obtained with the method QPCPQ in place of PQCQP and with Q = −1. In addition to the proposed splitting methods, we have also added simulation results with the standard fourth-order Runge-Kutta method RK4, which requires four force evaluations per time step and it is neither symmetric nor symplectic. Figure 2a,d illustrate the second-order convergence of the splitting methods, while, as expected, RK4 is of order four. In order to better observe that, we have supplemented the figures with lines of slopes two and four, represented by the solid and dashed black lines, respectively. Figure 2a,d demonstrate that the solution errors of both methods PQDWDQP and PQDWDQP, which differ only by the rescaling (93), are indistinguishable. See also Figure 2b,e. They have much smaller solution errors compared to the method PQCQP, where we have not split the piece C, and to RK4, especially for larger time steps h and ratio value E 0 /τ. Recall that the value E 0 /τ appears in the charge dispersion relation (74) and produces high-frequency oscillations of the charge. From Figure 2, it is easy to see that to obtain stable numerical solutions with RK4 for large values of E 0 /τ, we require smaller time steps compared to the methods PQCQP, PQDWDQP, and PQDWDQP. When E 0 /τ → 0 (not shown), both methods PQCQP and PQDWDQP coincide and RK4 has smaller errors compared to the splitting methods, which can be attributed to fourth-order accuracy. Then, only in long-time simulations, the performance of splitting methods can be appreciated over RK4, since RK4 is a dissipative method. In addition, the case E 0 /τ → 0 may not be of particular interest in modeling charge transfer in realistic crystals. In Figure 2b,e, we demonstrate the maximal relative Hamiltonian errors and observe that with RK4, the Hamiltonian error converges at most with order four, and despite not splitting the piece C into two pieces D and W, the PQCQP method's Hamiltonian errors are almost indistinguishable from the Hamiltonian errors of the two methods PQDWDQP and PQDWDQP. Figure 2c,f show, as expected, that RK4 does not exactly conserve the total charge probability (16) compared to all semi-implicit splitting methods. Notice that the method PQDWDQP conserves the total probability to machine precision due to the rescaling (93), while the errors for the other two methods PQCQP and PQDWDQP fluctuate close to double machine precision and slightly increase with smaller time steps. It is important to mention that in numerical results we did not observe differences between methods PQDWDQP and PQDWDQP and between methods QPDWDPQ and QPDWDPQ. Thus, we postulate that the scaling of the charge variable c in (93) does not affect the numerical solution, at least in short-time simulations. In addition, the charge probability is conserved up to machine double-precision, as can be seen in Figure 2c,f, and does not increase with an increase in T end .
The results of Figure 2 illustrate the importance of splitting the piece C such that numerical simulations can be performed with larger time steps. Thus, we will disregard the methods RK4, PQCQP and QPCPQ, from further numerical study and will concentrate on the methods PQDWDQP and QPDWDPQ in place of the methods PQDWDQP and QPDWDPQ, respectively.
In Figure 3, we investigate the differences between the velocity and position semiimplicit splitting methods PQDWDQP and QPDWDPQ, respectively, for two different ratio E 0 /τ values. As in Figures 2 and 3, we plot the solution, Hamiltonian, and total probability errors for different time step h values. Figure 3a,d illustrate that the solution errors are very comparable with slightly smaller errors for the position method QPDWDPQ. On the contrary, the velocity splitting method PQDWDQP, see Figure 3b,e, gives smaller Hamiltonian conservation errors of greater magnitude compared to the differences in the solution errors. Figure 3c,f show that both methods exactly conserve the total charge probability (16). Thus, we advocate the method PQDWDQP over the position splitting method QPDWDPQ if Hamiltonian conservation is of higher importance. Recall that both methods are of the same computational complexity, see the discussion in Section 4.2 and Appendix A.2.
So far, in Figures 2 and 3, we have investigated semi-implicit splitting methods that conserve exactly the charge probability, as presented in Section 4.2. In Figure 4, we demonstrate the approximation and the convergence order of the explicit splitting methods proposed in Section 4.3, i.e., we compare methods PQABDBAQP, BADQPQDAB, DBAQPQABD, and PQDABADQP for the same two ratio E 0 /τ values. Figure 4 demonstrates that all explicit methods are second order. We do not observe significant differences between all four methods when the solution and Hamiltonian conservation errors are compared, see Figure 4a,b,d,e. Interesting differences appear in the total probability conservation, see is satisfied. In future work, we plan to explore this from the analytic point of view. This observation also motivates to explore in future work multiple time-stepping approaches, where the charge part, e.g., DABAD or ABA in method PQDABADQP, is applied multiple times with smaller time step. To even better indicate the differences between semi-implicit and all explicit splitting methods, we investigate in more detail the Hamiltonian and the total probability conservation errors for different ratio E 0 /τ values in the following section.

Numerical Accuracy of Conserved Quantities
In the previous section, we numerically demonstrated the approximation and convergence orders of the proposed symplecticity-preserving symmetric splitting methods of Section 4 and showed the importance of splitting the piece C, which allows us to obtain more accurate and stable results with larger time steps. In this section, we continue with a comparison of the semi-implicit and explicit splitting methods by investigating the conservation of several quantities for different values of E 0 /τ. In Figures 5 and 6 and errors are averaged over eleven numerical simulations with different initial conditions, which are described below. In Figures 5 and 6, we illustrate the results for two regimes of solutions, i.e., (approximately) linear and nonlinear, respectively. In the linear regime simulations, we generate eleven different initial conditions with u n , a n , p n , and b n values for all n drawn from the uniform distribution U (−0.01, 0.01). After the generation of values a n and b n , we rescale 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 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 5 and 6 with 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 PQDWDQP and QPDWDPQ 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 PQDWDQP and QPDWDPQ coincide very well with the errors of the explicit splitting methods DBAQPQABD and PQDABADQP for all considered ratio E 0 /τ values. Interestingly, as the value of E 0 /τ 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 E 0 /τ < 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 E 0 /τ 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 5 and 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 E 0 /τ 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 PQDWDQP compared to the position method QPDWDPQ, see Figure 6a-c, which can also be seen in Figure 3b 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 PQDWDQP over the position method QPDWDPQ 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.

Charge Transfer by Mobile Discrete Breathers
In the example in Figure 1, we demonstrate charge transfer by a mobile discrete breather in time until T end = 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 T end = 10 4 with the three proposed splitting methods of Section 4, i.e., with three lattice-charge-lattice methods PQDWDQP, PQDABADQP, and PQABDBAQP. Numerical simulations are performed in a lattice with N = 254 particles, time step h = 0.01, and E 0 /τ = 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 P n (t) = |c n (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, T end ], see Figure 7d. The same is not true for the total probability conservation (16), as can be seen in Figure 7e. Semiimplicit method PQDWDQP 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. 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]) 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 P r , the more localized the charge probability is. For example, if the charge probability is completely localized at one site, then P r = 1, but in the situation when the charge probability is equally spread across the whole lattice, then P r = 0. With our charge probability's initial excitation pattern (101), we have that P r (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-time results with good approximate probability conservation properties compared to the exactly charge-probability conserving semi-implicit methods.

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 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 E 0 /τ. 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 E n (q 1 , . . . , q N ) 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.

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 two-and 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 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. Q = q + tM −1 p Q = q Q = q Q = q Q = q A = a A = a A = Re(ϕ E(q) t (a + bi)) A = a + tL(q)b A = a P = p P = p + tF(q) P = p + tG D (q, a, b) P = p + tG A (q, b) P = p + tG B (q, a)