Geometric numerical integration of semi-classical Hamiltonain lattice dynamics

Geometric numerical integration of semi-classical Hamiltonain lattice dynamics J¯ anis Baj¯ ars†and Juan F. R. Archilla‡ †Faculty of Physics, Mathematics and Optometry, University of Latvia Jelgavas Street 3, Riga, LV-1004, Latvia ‡Group of Nonlinear Physics, Universidad de Sevilla ETSII, Avda Reina Mercedes s/n, 41012-Sevilla, Spain Email: [email protected], [email protected] Abstract—In this work we provide a brief overview of recently proposed symplecticity-preserving symmetric splitting methods for semi-classical Hamiltonian dynamics of charge transfer by intrinsic localized modes in nonlinear crystal lattice models [1]. Without loss of generality, we consider one-dimensional crystal lattice models described by classical Hamiltonian dynamics, whereas charge particle is modeled as quantum particle within the tight-binding approximation. Canonical Hamiltonian equations for the coupled lattice-charge dynamics are derived. Structurepreserving splitting methods are constructed by splitting the total Hamiltonian into the sum of Hamiltonians which individual dynamics can be solved exactly. Exactly charge conserving symplectic splitting methods are also proposed which require only one solution of a linear system of equations per time step. Developed computationally efficient non-dissipative methods provide new means for long-time simulations of charge transfer by nonlinear lattice excitations. 1. Introduction Study of charge transfer by nonlinear lattice excitations is of particular interest in solid state physics, where transport of charge in silicates by moving nonlinear localized exictations is experimentally confirmed, the phenomenon called hyperconductivity [2]. Intrinsic localized modes (ILMs), such as discrete breathers,kinks and solitons, have been extensively studied from analytic and numerical points of view in crystal lattice models [3, 4]. Traditionally, such lattice models at zero temperature are described by classical Hamiltonian dynamics with empirical particle interaction potentials, and by thermostated Hamiltonian dynamics at a given temperature [5]. The transport of charge (electrons or holes) by ILMs [4, 6, 7, 8] poses new numerical simulation challenges due to different oscillation time scales of charge and lattice particles which motivates to consider splitting methods for numerical integration of coupled lattice-charge dynamical equations. We demonstrate that lattice-charge dy- ORCID iDs First Author: 0000-0001-7601-8694, Second Author: 0000-0001-6583-6114 namics can be stated into classical canonical Hamiltonian form and, while the Hamiltonian is not separable in all variables, explicit symplectic and time-reversibility preserving numerical integrators can be constructed with good approximate energy and charge conservation in long-time numerical simulations. We also develop semi-implicit exactly charge conserving symplectic splitting methods. In addition, such splitting method approach may further allow us to develop multiple time stepping schemes, such as the impulse method [5, 9], to further improve numerical integration of the multiscale dynamics, and high order methods [5, 9], to improve numerical accuracy of nonlinear localized wave simulations. 2. Mathematical Problem In this section we describe mathematical model in dimensionless form of coupled lattice-charge dynamics modeled by the total (classical) Hamiltonian: H=Hlat +Hc,(1) i.e., the sum of lattice and charge Hamiltonians, respectively, where Hlat = N X n=1        1 2Mn p2 n+U(qn)+1 2 N X n0=1 n0,n V(|qn−qn0|)      ,(2) Hc= N X n=1        1 2τEn(q1,...,qN)a2 n+b2 n −1 2τ N X n0=1 n0,n J(qn,qn0)(anan0+bnbn0)      .(3) Nis the number of lattice particles, qn∈R,pn∈Rand Mn>0 are the nth particle’s position, momentum and rescaled mass, respectively. In the lattice Hamiltonian (2) Uis the on-site potential and Vis the radial interparticle potential energy. In our considerations lattice model (2) is thought to be a one-dimensional crystal lattice model of a three-dimensional layered crystal where the on-site potential Umodels forces from the upper and lower layers of the – 576 – 2022 International Symposium on Nonlinear Theory and Its Applications, NOLTA2022, Virtual, December 12-15, 2022 This work is licensed under a Creative Commons Attribution NonCommercial, No Derivatives 4.0 License. crystal, whereas the interaction potential Vmodels repulsive and, potentially, as well as attractive forces of lattice particles. To derive Hamiltonian (1) equations in the dimensionless form we consider characteristic mass M >0, length scale σ > 0, time scale T >0 and energy scale E >0 values such that the relation E =Mσ2T−2holds. In the charge Hamiltonian (3) τ=~(ET)−1>0 is the dimensionless Planck’s 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: E =~2(τMσ)−2and T =τMσ2~−1. In the charge Hamiltonian (3) variables anand bnare the real and imaginary parts of √2τcn, where cnis the probability amplitude and Pn:=|cn|2is the probability of finding an electron or hole at site n. The total probability N X n=1 Pn=1 (4) is conserved as it is the sum PN n=1a2 n+b2 n=2τ. The introduction of this normalization is convenient from the problem formulation point of view as the system becomes canonical. Symmetric charge hopping function J(qn,qn0)≥0, i.e., J(qn,qn0)=J(qn0,qn), describes charge transfer between states nand n0. Naturally, the transfer function Jdepends on the particle interaction distance rn,n0:=|qn−qn0|and is often modeled with exponential decay [4, 6]: J(qn,qn0) := J0exp −αrn,n0, where the dimensionless parameter α > 0 specifies the rate of exponential decay, while the constant J0≥0 models relative strength of charge transfer from one site to another and is model dependent. Smooth multivariable function En(q1,...,qN)∈Rin (3) describes charge energy at site nand, in general, will depend on the lattice particle positions, or can be modeled as constant. For completeness we state rescaled canonical Hamiltonian equations derived from (1) with (2)–(3) in the component form: ˙qn=1 Mn pn,(5) ˙an=1 τEn(q1,...,qN)bn−1 τ N X n0=1 n0,n J(qn,qn0)bn0,(6) ˙pn=−U0(qn)− N X n0=1 n0,n ∂qnV(|qn−qn0|) −1 2τ N X n0=1 ∂qnEn0(q1,...,qN)a2 n+b2 n +1 τ N X n0=1 n0,n ∂qnJ(qn,qn0)(anan0+bnbn0),(7) ˙ bn=−1 τEn(q1,...,qN)an+1 τ N X n0=1 n0,n J(qn,qn0)an0,(8) for all n=1,...,N. Canonical Hamiltonian equations are well known to be symplectic, and, thus, they also conserve phase volume. Apart from the conservation of the total Hamiltonian (1), charge probability (4) is also conserved. Dynamics (5)– (8) is time-reversible, i.e., equations are invariant under the transformation: ρ(q,a,−p,−b,−t)=(q,a,p,b,t),(9) where q:=(q1,q2,...,qN)T,a:=(a1,a2,...,aN)T,p:= (p1,p2,...,pN)Tand b:=(b1,b2,...,bN)T. In addition, Hamiltonian equations (5)–(8) are invariant under a constant rotation of charge variables anand bn, i.e., for any given angle θ∈Requations (5)–(8) are invariant under the following transformation: η(q,¯a,p,¯ b,t)=(q,a,p,b,t),(10) where ¯a=cos(θ)a−sin(θ)band ¯ b=sin(θ)a+cos(θ)b. The system of differential equations (5)–(8) is highly nonlinear and, thus, it is very desirable to obtain an explicit, or at least linearly implicit, numerical integration scheme, while at the same time attempting to preserve as many as possible structural properties of the Hamiltonian system (5)–(8) stated above, which we discuss in the following section. 3. Splitting Methods In this section we review recently proposed symplecticity-preserving symmetric splitting methods for the Hamiltonian dynamics (5)–(8) [1]. We consider two classes of methods, i.e., semi-implicit exactly charge probability (4) and rotational invariance (10) conserving and fully explicit methods. All proposed methods are symplectic, symmetric and of second order. From the computational efficiency point of view all methods require only one lattice force: first two force terms in (7), and function En(q1,...,qN) and J(qn,qn0) calculations in (6) and (8) per time step, and at most two evaluations of the last two force terms in (7) with the same values of qn, i.e., only values of anand bndiffer. Compared to the semi-implicit methods, which require one solution of a linear system of equations per time step, the fully explicit methods do not exactly conserve charge probability (4) and rotational invariance (10). Recall that symmetric methods preserve time-reversibility property (9) [9]. To construct symplectic splitting methods we consider the total Hamiltonian (1) as the sum of the following Hamiltonians: H=HQ+HP+HD+HA+HB, where HQ= N X n=1 1 2Mn p2 n,(11) – 577 – HP= N X n=1       U(qn)+1 2 N X n0=1 n0,n V(|qn−qn0|)      ,(12) HD=1 2τ N X n=1 En(q1,...,qN)a2 n+b2 n,(13) HA=−1 2τ N X n=1 N X n0=1 n0,n J(qn,qn0)bnbn0,(14) HB=−1 2τ N X n=1 N X n0=1 n0,n J(qn,qn0)anan0.(15) In addition, for construction of charge probability conserving splitting methods we also consider the Hamiltonians: HC=HD+HA+HB,(16) HW=HA+HB.(17) Importantly, corresponding Hamiltonian systems associated to each Hamiltonian (11)–(15) can be solved exactly. Each such solution we identify with analytic symplectic flow maps: φQ t,φP t,φD t,φA tand φB t, respectively. What follows, with q0 n,a0 n,p0 nand b0 nwe identify initial conditions at time t=0, where t≥0. Solution of the Hamiltonian system associated to split dynamics (11) reads: an(t)=a0 n,pn(t)=p0 n,bn(t)=b0 n, qn(t)=q0 n+1 Mn tp0 n, for all t>0 and n. Similarly, solution of the Hamiltonian system associated to split dynamics (12) reads: qn(t)=q0 n,an(t)=a0 n,bn(t)=b0 n, pn(t)=p0 n−t      U0(q0 n)+ N X n0=1 n0,n ∂qnV(|q0 n−q0 n0|)      . Exact solutions of the Hamiltonian systems with Hamiltonians (14) and (15) are, i.e., the flow maps φA tand φB t: qn(t)=q0 n,bn(t)=b0 n, an(t)=a0 n−1 τt N X n0=1 n0,n J(q0 n,q0 n0)b0 n0, pn(t)=p0 n+1 τt N X n0=1 n0,n ∂qnJ(q0 n,q0 n0)b0 nb0 n0, and qn(t)=q0 n,an(t)=a0 n, pn(t)=p0 n+1 τt N X n0=1 n0,n ∂qnJ(q0 n,q0 n0)a0 na0 n0, bn(t)=b0 n+1 τt N X n0=1 n0,n J(q0 n,q0 n0)a0 n0, respectively, for all t>0 and n=1,...,N. To find exact solution of the Hamiltonian system associated to the Hamiltonian (13) we first write down the system of equations: ˙qn=0,(18) ˙an=1 τEn(q1,...,qN)bn,(19) ˙pn=−1 2τ N X n0=1 ∂qnEn0(q1,...,qN)a2 n+b2 n,(20) ˙ bn=−1 τEn(q1,...,qN)an,(21) and make an important observation that with constant values of qncoupled equations (19) and (21) for variables an and bncan be solved exactly for all n, i.e., an(t)=a0 ncos(ωnt)+b0 nsin(ωnt),(22) bn(t)=−a0 nsin(ωnt)+b0 ncos(ωnt),(23) where ωn:=τ−1En(q0 1,...,q0 N), and the equation (20) is invariant under the solution (22)–(23), i.e., an(t)2+bn(t)2= a0 n 2+b0 n 2for all nand t. Thus, the exact (explicit) solution of the Hamiltonian system (18)–(21) is in the following form: qn(t)=q0 n, an(t)=a0 ncos(ωnt)+b0 nsin(ωnt), pn(t)=p0 n−1 2τt N X n0=1 ∂qnEn0(q0 1,...,q0 N)a0 n 2+b0 n 2, bn(t)=−a0 nsin(ωnt)+b0 ncos(ωnt). In conclusion, all analytic symplectic maps φQ t,φP t,φD t,φA t and φB tare explicit maps which will provide means to construct explicit symplecticity-preserving numerical methods. Note that flow maps φA tand φB tdo not preserve charge probability (4) and rotational invariance (10). To preserve (4) and (10) we could solve split Hamiltonian systems with Hamiltonians (16) or (17) with symplectic implicit midpoint-rule which preserves quadratic invariants [9]. Notice that all three Hamiltonians (13)–(15) are quadratic with respect to charge variables anand bn. In [1] we have demonstrated that numerical stability and accuracy can be greatly improved by splitting Hamiltonian (16) into two Hamiltonians HDand HW. Such splitting naturally occurs by splitting the kinetic from the potential energy in splitting methods for solving continuous Schr¨ odinger equations. With the analysis above we list the splitting methods proposed in [1]. What follows, with ψhwe identify a numerical flow map which advances a given system’s state (q,a,p,b)Tinto a new state (Q,A,P,B)Tin time after the time step h>0. All four proposed semi-implicit methods are listed in Table 1, while all four fully explicit splitting methods are shown in Table 2. In Table 1 with ψC hand ψW h – 578 – Numerical flow maps ψh PQCQP φP h/2◦φQ h/2◦ψC h◦φQ h/2◦φP h/2 QPCPQ φQ h/2◦φP h/2◦ψC h◦φP h/2◦φQ h/2 PQDWDQP φP h/2◦φQ h/2◦φD h/2◦ψW h◦φD h/2◦φQ h/2◦φP h/2 QPDWDPQ φQ h/2◦φP h/2◦φD h/2◦ψW h◦φD h/2◦φP h/2◦φQ h/2 Table 1: Semi-implicit symplecticity-preserving symmetric and charge probability conserving splitting methods. Numerical flow maps ψh PQABDBAQP φP h/2◦φQ h/2◦φA h/2◦φB h/2◦φD h◦φB h/2◦φA h/2◦φQ h/2◦φP h/2 QPABDBAPQ φQ h/2◦φP h/2◦φA h/2◦φB h/2◦φD h◦φB h/2◦φA h/2◦φP h/2◦φQ h/2 PQDABADQP φP h/2◦φQ h/2◦φD h/2◦φA h/2◦φB h◦φA h/2◦φD h/2◦φQ h/2◦φP h/2 QPDABADPQ φQ h/2◦φP h/2◦φD h/2◦φA h/2◦φB h◦φA h/2◦φD h/2◦φP h/2◦φQ h/2 Table 2: Fully explicit symplecticity-preserving symmetric splitting methods. we have identified numerical flow maps of implicit midpoint rule applied to respective Hamiltonian systems. For example, the implicit midpoint rule ψW hreads: Qn=qn,(24) An=an+h τEn(q)Bn+bn 2−h τ N X n0=1 n0,n J(qn,qn0)Bn0+bn0 2,(25) Pn=pn+h τ N X n0=1 n0,n ∂qnJ(qn,qn0)ζn,n0,(26) Bn=bn−h τEn(q)An+an 2+h τ N X n0=1 n0,n J(qn,qn0)An0+an0 2,(27) where ζn,n0=1 4[(An+an)(An0+an0)+(Bn+bn)(Bn0+bn0)], En(q)=En(q1,...,qN) and momentum value Pnin (26) is found explicitly after the linear equation system (equations (25) and (27) for all n=1,...,N) is solved for the charge variables Anand Bn. 4. Summary and Conclusions In this work we have reviewed recently proposed geometric structure-preserving splitting methods for the Hamiltonian dynamics (5)–(8) [1]. In [1] we have performed extensive numerical study (not shown) to investigate and compare all proposed methods in Tables 1 and 2. Numerical results, which will be presented at the conference, show that the best results from all semi-implicit methods are obtained with the method PQDWDQP, while from all explicit methods in Table 2 the best results were obtained with the method PQDABADQP, measured by numerical errors in the Hamiltonian (1) and charge probability (4) conservation in numerical simulations of charge transfer by discrete breathers. Acknowledgments 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. JFR Archilla thanks projects MICINN PID2019-109175GB-C22 and Junta de Andalucia US- 1380977, and a travel grant from VIIPPITUS-2022. References [1] J. Baj¯ ars and J.F.R. Archilla, “Splitting methods for semi-classical Hamiltonian dynamics of charge transfer in nonlinear lattices,” submitted, 2022. [2] F. Russell, J.F.R. Archilla, and S. Medina-Carrasco, “Localized waves in silicates. What do we know from experiments?,” in: C. Skiadas, Y. Dimotikalis (Eds.), 13th Chaotic Modeling and Simulation International Conference, Springer International Publishing, Cham, pp. 721–734, 2021. [3] S. Flach and A. Gorbach, “Discrete breathers — Advances in theory and applications,” Physics Reports, vol.467, pp.1–116, 2008. [4] J.F.R. Archilla et al, Eds, Quodons in Mica, Springer, Cham, 2015. [5] B. Leimkuhler and C. Matthews, Molecular Dynamics: With Deterministic and Stochastic Numerical Methods, Springer, Cham, 2015. [6] O.G. Cantu Ros, L. Cruzeiro, M.G. Velarde, and W. Ebeling, “On the possibility of electric transport mediated by long living intrinsic localized solectron modes,” Eur. Phys. J. B, vol.80, pp. 545–554, 2011. [7] L. A. Cisneros-Ake, “Stability of multi-hump localized solutions in the Holstein model for linear acoustic and soft nonlinear optical interactions,” Physica D, vol.431, 133138, 2022. [8] L. Cruzeiro-Hansson, J.C. Eilbeck, J.L. Mar´ın, and F.M. Russell, “Breathers in systems with intrinsic and extrinsic nonlinearities,” Physica D, vol.142, pp. 101– 112, 2000. [9] E. Hairer, C. Lubich, and G. Wanner, Geometric Numerical Integration: Structure-Preserving Algorithms for Ordinary Differential Equations, Springer Science &Business Media, 2006. – 579 –