Efficient geometric integrators for nonadiabatic quantum dynamics. I. The adiabatic representation

TL;DR

This paper introduces implicit geometric integrators for nonadiabatic quantum dynamics that are energy-conserving, unitary, and more accurate than traditional methods, demonstrated on a two-surface NaI model.

Contribution

The paper develops and analyzes new implicit geometric integrators applicable to non-separable Hamiltonians in the adiabatic representation, enhancing accuracy and efficiency.

Findings

Integrators conserve energy exactly and are symplectic and unitary.

Numerical tests show significant speedup over traditional methods.

Higher-order integrators achieve arbitrary accuracy with controlled computational cost.

Abstract

Geometric integrators of the Schr\"{o}dinger equation conserve exactly many invariants of the exact solution. Among these integrators, the split-operator algorithm is explicit and easy to implement, but, unfortunately, is restricted to systems whose Hamiltonian is separable into a kinetic and potential terms. Here, we describe several implicit geometric integrators applicable to both separable and non-separable Hamiltonians, and, in particular, to the nonadiabatic molecular Hamiltonian in the adiabatic representation. These integrators combine the dynamic Fourier method with recursive symmetric composition of the trapezoidal rule or implicit midpoint method, which results in an arbitrary order of accuracy in the time step. Moreover, these integrators are exactly unitary, symplectic, symmetric, time-reversible, and stable, and, in contrast to the split-operator algorithm, conserve energy exactly, regardless of the accuracy of the solution. The order of convergence and conservation of geometric properties are proven analytically and demonstrated numerically on a two-surface NaI model in the adiabatic representation. Although each step of the higher order integrators is more costly, these algorithms become the most efficient ones if higher accuracy is desired; a thousand-fold speedup compared to the second-order trapezoidal rule (the Crank-Nicolson method) was observed for wavefunction convergence error of $10^{-10}$. In a companion paper [J. Roulet, S. Choi, and J. Van\'{\i}\v{c}ek (2019)], we discuss analogous, arbitrary-order compositions of the split-operator algorithm and apply both types of geometric integrators to a higher-dimensional system in the diabatic representation.

Figures

Tables

Table 1. Table 1: Geometric properties and computational cost of various integrators. Cost is measured by the number of Fourier transforms required per time step (see Sec.  II.6 ). I 𝐼 I is the number of iterations for the implicit step and n = 0 , 1 , 2 , … 𝑛 0 1 2 … n=0,1,2,\ldots is the number of recursive compositions. C 𝐶 C is the total number of composition steps per time step ( C = 3 n 𝐶 superscript 3 𝑛 C=3^{n} for the triple jump Yoshida ( 1990 ); Suzuki ( 1990 ) , C = 5 n 𝐶 superscript 5 𝑛 C=5^{n} for Suzuki’s fractals Suzuki ( 1990 ) ). + + denotes that the geometric property of the exact evolution operator is preserved and − - denotes that it is not. SOD stands for the second-order differencing and SO for the split-operator algorithm.

	Order	Unitary	Sympl-	Commutes	Energy	Symm-	Time-	Stable	Cost
			ectic	with $\hat{H}$	cons.	etric	reversible
Elementary methods
$1^{\mathrm{st}}$ order SO	1	$+$	$+$	$-$	$-$	$-$	$-$	$+$	$2$
Expl. Euler	1	$-$	$-$	$+$	$-$	$-$	$-$	$-$	$4D$
Impl. Euler	1	$-$	$-$	$+$	$-$	$-$	$-$	$+$	$4D\left(2+I\right)$
SOD	2	$-$	$-$	$+$	$-$	$+$	$+$	$+$11footnotemark: 1	$4D$
Recursively composable methods
$2^{\mathrm{nd}}$ order SO	2$\left(n+1\right)$	$+$	$+$	$-$	$-$	$+$	$+$	$+$	$2C$
Midpoint	2$\left(n+1\right)$	$+$	$+$	$+$	$+$	$+$	$+$	$+$	$4D\left(3+I\right)C$
Trapezoidal	2$\left(n+1\right)$	$+$	$+$	$+$	$+$	$+$	$+$	$+$	$4D\left(3+I\right)C$

Equations

$$i\hslash \frac{d\psi (t)}{dt}=\hat{H}\psi (t)$$

$$\hat{U}(t)=e^{-i\hat{H}t/\hslash }$$

$$⟨\hat{U}\psi |\hat{U}\varphi ⟩\equiv ⟨\psi |{\hat{U}}^{†}\hat{U}\varphi ⟩=⟨\psi |\varphi ⟩,$$

$${\hat{U}}^{-1}={\hat{U}}^{†}.$$

$$\omega (\psi ,\varphi ):=-2\hslash \mathrm{Im}⟨\psi |\varphi ⟩;$$

$$E(t)$$

$$=⟨\psi (0)|\hat{U}(t{)}^{†}\hat{H}\hat{U}(t)|\psi (0)⟩$$

$$=⟨\psi (0)|\hat{U}(t{)}^{†}\hat{U}(t)\hat{H}|\psi (0)⟩$$

$$=⟨\psi (0)|\hat{H}|\psi (0)⟩=E(0).$$

$$\hat{U}(t{)}^{\ast }:=\hat{U}(-t{)}^{-1}.$$

$$\hat{U}(t{)}^{\ast }=\hat{U}(t).$$

$$\hat{U}(-t)\hat{U}(t)\psi (0)=\psi (0).$$

$$\Vert \psi (0)-\varphi (0)\Vert <\delta \text{ implies }\Vert \psi (t)-\varphi (t)\Vert <ϵ\text{ for all t;}$$

$$\Vert \psi (0)-\varphi (0)\Vert <\delta \text{ implies }\Vert \psi (t)-\varphi (t)\Vert \to 0\text{ as }t\to \infty ;$$

$$\Vert \psi (t)-\varphi (t)\Vert =\Vert \psi (0)-\varphi (0)\Vert .$$

$${\hat{U}}_{\mathrm{expl}}(\Delta t):=1-\frac{i}{\hslash }\Delta t\hat{H}$$

$${\hat{U}}_{\mathrm{impl}}(\Delta t):={\left(1+\frac{i}{\hslash }\Delta t\hat{H}\right)}^{-1}$$

$$\Vert {\psi }_{\mathrm{expl}}(t)-{\varphi }_{\mathrm{expl}}(t)\Vert \to \infty \text{ as }t\to \infty ,$$

$${\psi }_{\mathrm{sod}}(t+\Delta t)-{\psi }_{\mathrm{sod}}(t-\Delta t)=-2\frac{i}{\hslash }\Delta t\hat{H}{\psi }_{\mathrm{sod}}(t).$$

$$\frac{d\psi (t)}{dt}\approx \frac{\psi (t+\Delta t)-\psi (t-\Delta t)}{2\Delta t}.$$

$${\hat{U}}_{\mathrm{trap}}(\Delta t):={\hat{U}}_{\mathrm{impl}}(\Delta t/2){\hat{U}}_{\mathrm{expl}}(\Delta t/2),$$

$${\hat{U}}_{\mathrm{midp}}(\Delta t):={\hat{U}}_{\mathrm{expl}}(\Delta t/2){\hat{U}}_{\mathrm{impl}}(\Delta t/2).$$

$${\hat{U}}_{p+2}(\Delta t):={\hat{U}}_p({\gamma }_M\Delta t)\cdots {\hat{U}}_p({\gamma }_1\Delta t)$$

$${\gamma }_1=\frac{1}{2-2^{1/(p+1)}},{\gamma }_2=-\frac{2^{1/(p+1)}}{2-2^{1/(p+1)}},$$

$${\gamma }_1={\gamma }_2=\frac{1}{4-4^{1/(p+1)}},{\gamma }_3=-\frac{4^{1/(p+1)}}{4-4^{1/(p+1)}}.$$

$$({\hat{U}}_2{\hat{U}}_1{)}^{†}={\hat{U}}_1^{†}{\hat{U}}_2^{†}={\hat{U}}_1^{-1}{\hat{U}}_2^{-1}=({\hat{U}}_2{\hat{U}}_1{)}^{-1}.$$

$$\omega ({\hat{U}}_2{\hat{U}}_1\psi ,{\hat{U}}_2{\hat{U}}_1\varphi )=\omega ({\hat{U}}_1\psi ,{\hat{U}}_1\varphi )=\omega (\psi ,\varphi ).$$

$$⟨\hat{H}{⟩}_{{\hat{U}}_2{\hat{U}}_1\psi }=⟨\hat{H}{⟩}_{{\hat{U}}_1\psi }=⟨\hat{H}{⟩}_{\psi }.$$

$$({\hat{U}}_2{\hat{U}}_1{)}^{\ast }={\hat{U}}_1^{\ast }{\hat{U}}_2^{\ast }={\hat{U}}_1{\hat{U}}_2\ne {\hat{U}}_2{\hat{U}}_1.$$

$$({\hat{U}}_1{\hat{U}}_2{\hat{U}}_1{)}^{\ast }={\hat{U}}_1^{\ast }{\hat{U}}_2^{\ast }{\hat{U}}_1^{\ast }={\hat{U}}_1{\hat{U}}_2{\hat{U}}_1.$$

$${\hat{U}}_2(-\Delta t_2){\hat{U}}_1(-\Delta t_1){\hat{U}}_2(\Delta t_2){\hat{U}}_1(\Delta t_1)={\hat{U}}_2(\Delta t_2{)}^{-1}{\hat{U}}_1(\Delta t_1{)}^{-1}{\hat{U}}_2(\Delta t_2){\hat{U}}_1(\Delta t_1)\ne 1.$$

$${\hat{U}}_1(-\Delta t_1){\hat{U}}_2(-\Delta t_2){\hat{U}}_1(-\Delta t_1){\hat{U}}_1(\Delta t_1){\hat{U}}_2(\Delta t_2){\hat{U}}_1(\Delta t_1)$$

Full text

I Introduction

Separating electronic from nuclear degrees of freedom leads to the Born–Oppenheimer approximationBorn and Oppenheimer (1927); Heller (2018) and the intuitive picture of electronic potential energy surfaces. However, many chemical, physical, and biological processes can only be described by taking into account the correlation between the nuclear and electronic motions,Mokhtari et al. (1990) which is reflected in the nonadiabatic couplings between different Born–Oppenheimer surfaces.Baer (2006); Domcke and Yarkony (2012); Nakamura (2012); Takatsuka et al. (2015); Bircher et al. (2017) To address such processes, one can forget the Born–Oppenheimer picture and treat electrons and nuclei on the same footing,Shin and Metiu (1995); Albert, Kaiser, and Engel (2016) use an exact factorizationAbedi, Maitra, and Gross (2010); Cederbaum (2008) of the molecular wavefunction, or, most commonly, determine which Born–Oppenheimer states are significantly coupledZimmermann and Vaníček (2010, 2012) and then solve the time-dependent Schrödinger equation with a molecular Hamiltonian that contains the nonadiabatic couplings. Below, we will only consider the third, yet the most traditional way to treat quantum nonadiabatic dynamics.

An approach particularly suited to study the nonadiabatic population dynamics of large chemical systems is the ab initio multiple spawningBen-Nun, Quenneville, and Martínez (2000); Curchod and Martínez (2018) and related methods, all of which represent the wavefunction by a superposition of time-dependent Gaussian basis functions moving along classicalShalashilin and Child (2001); Makhov et al. (2017) or variationalWorth, Robb, and Burghardt (2004); Richings et al. (2015) trajectories. If high accuracy is required and especially if the Hamiltonian can be expressed as a sum of products of one-dimensional operators, a nonadiabatic algorithm of choice is the multiconfigurational time-dependent Hartree (MCTDH) methodMeyer, Manthe, and Cederbaum (1990); Worth et al. (2008) or its multilayer extension,Wang and Thoss (2003) which expand the state using orthogonal time-dependent basis functions. The power of the MCTDH method relies on the fact that only a small fraction of the tensor-product Hilbert space is typically accessible during the time of interest; sparse-grid methodsAvila and Carrington Jr (2017); Lubich (2008) also take advantage of this phenomenon. However, there are systems, in which the full Hilbert space is accessible, and then full grid or time-independent basis sets are preferable.Kosloff (1988); Lubich (2008)

There also exist situations, where, in addition to prescribed accuracy, it pays off to conserve certain invariants of the exact solution exactly, regardless of the accuracy of the wavefunction. Because the above-mentioned methods typically conserve none or only some of these invariants, other methods, called geometric integrators,Hairer, Lubich, and Wanner (2006) are needed in this setting. The geometric integrators acknowledge that the Schrödinger equation is special, and not just another general differential equation. Using these integrators can be likened to realizing that the Earth is not flat but round, and even approximate models of its surface should take this curvature into account. Geometric integrators are highly exploited in classical molecular dynamics, where the deceptively simple Verlet algorithm,Verlet (1967); Frenkel and Smit (2002) despite its only second-order accuracy, results in exact conservation of $D$ invariants in a $D$-dimensional system, where $D$ can easily reach thousands or millions in state-of-the art simulations of proteins.

Time propagation schemes based on geometric integrators have also been applied to the time-dependent Schrödinger equation.McCullough and Wyatt (1971); Feit, Fleck, and Steiger (1982); Lubich (2008); Tannor (2007) Symmetric compositions of the first-order split-operator algorithms,Lubich (2008); Tannor (2007) including the standard second-order splitting,Feit, Fleck, and Steiger (1982) are unitary, symplectic, stable, symmetric, and time-reversible regardless of the size of the time step. Moreover, the symmetric split-operator algorithms can be recursively composed to obtain efficient methods of arbitrary order in the time step.Yoshida (1990); McLachlan (1995); Suzuki (1990); Hairer, Lubich, and Wanner (2006); Wehrle, Šulc, and Vaníček (2011); Roulet, Choi, and Vaníček In a companion paperRoulet, Choi, and Vaníček (below referred to as Paper II), we implement such higher-order compositions for the nonadiabatic quantum molecular dynamics in the diabatic representation.

Although the split-operator algorithms preserve numerous geometric properties of interest of the exact evolution operator, their use is limited to systems with Hamiltonians separable into a sum $\hat{H}=T(\hat{p})+V(\hat{q})$ of two terms, the first depending only on the momentum operator and the second only on the position operator. One must use a different time propagation scheme for systems with a non-separable Hamiltonian; for example, the nonadiabatic dynamics in the adiabatic representation or particles in crossed electric and magnetic fields.

The explicit Euler method is the simplest integrator applicable to non-separable Hamiltonians; it is, however, unstable.Leimkuhler and Reich (2004); Hairer, Lubich, and Wanner (2006) The implicit Euler method is stable regardless of the size of the time step but requires solving a large, although sparse, system of linear equations at every time step; furthermore, the method fails to preserve the unitarity, time reversibility, energy conservation, and other geometric properties of the exact evolution operator. The second-order differencing methodAskar and Cakmak (1978); Kosloff and Kosloff (1983a); Leforestier et al. (1991) introduces symmetry by combining the forward and backward step of the explicit Euler method. It is explicit and stable for small enough time steps, but does not conserve the norm or energy exactly.

Another issue with the second-order differencing is that a much too small time step is required to obtain an accurate solution.Lubich (2002) This problem has been addressed by using the ChebyshevTal-Ezer and Kosloff (1984) and short iterative Lanczos algorithms;Lanczos (1950); Park and Light (1986); Leforestier et al. (1991) both methods increase remarkably the efficiency of numerical integration by effectively approximating the exact evolution operator. However, these two methods are neither time-reversible nor symplectic, and the Chebyshev propagation scheme does not even conserve the norm.

To address either the low accuracy or nonconservation of geometric properties by various nonadiabatic integrators, we propose time propagation schemes based on symmetric compositions of the trapezoidal rule (also known as the Crank-Nicolson methodCrank and Nicolson (1947); McCullough and Wyatt (1971)) or implicit midpoint method. As we show below, because these elementary methods are unitary, symplectic, energy conserving, stable, symmetric, and time-reversible, so are their symmetric compositions. Furthermore, like any other symmetric second-order algorithm, the trapezoidal rule and implicit midpoint method can be recursively composed to obtain integrators of arbitrary order of accuracy in the time step.Yoshida (1990); McLachlan (1995); Suzuki (1990); Hairer, Lubich, and Wanner (2006) Methods with higher orders of accuracy are useful for obtaining highly accurate solutions because, for that purpose, they are more efficient than the second-order algorithms. Although each time step of a higher-order method costs more, the solution with the same accuracy can be obtained using a larger time step and, hence, a smaller total number of steps in comparison to lower-order methods. The final benefit of the proposed geometric integrators is the simple, abstract, and general implementation of the compositions of the trapezoidal rule and implicit midpoint methods; indeed, even these “elementary” methods are, themselves, compositions of simpler explicit and implicit Euler methods.

In the adiabatic representation, the proposed integrators cannot be fully compared with the integrators based on the compositions of the split-operator algorithm, which are only applicable to separable Hamiltonians. Both types of integrators, however, can be used in the diabatic representation, which is the focus of Paper II.Roulet, Choi, and Vaníček We, therefore, compare the two integrators there, using a one-dimensional modelEngel and Metiu (1989) of NaI and a three-dimensional modelStock et al. (1995) of pyrazine.

The remainder of this paper is organized as follows: In Section II, after defining geometric properties of the exact evolution operator, we discuss their breakdown in elementary methods and recovery in the proposed symmetric compositions of the trapezoidal rule and implicit midpoint methods. Next, we present the dynamic Fourier method for its ease of implementation and the exponential convergence with the number of grid points. Yet, the proposed integrators can be combined with any other basis or grid representation. We conclude Section III by discussing the relationship between the molecular Hamiltonians in the adiabatic and diabatic representations. In Section III, the convergence properties and conservation of geometric invariants by various methods are analyzed numerically on a two-surface NaI modelEngel and Metiu (1989) in the adiabatic representation. This system has a non-separable Hamiltonian due to an avoided crossing between its potential energy surfaces and a corresponding region of large nonadiabatic momentum coupling. Section IV concludes the paper.

II Theory

For a time-independent Hamiltonian $\hat{H}$, the time-dependent Schrödinger equation

[TABLE]

has the formal solution $\psi(t)=\hat{U}(t)\psi(0)$, where $\psi(0)$ is the initial state and $\hat{U}(t)$ the so-called evolution operator. The exact evolution operator

[TABLE]

is linear (in particular, independent of the initial state), reversible, stable, and, moreover, conserves both the norm and energy of the quantum state. Let us define and discuss these and other geometric properties of the exact evolution operator because they are also desirable in approximate numerical evolution operators $\hat{U}_{\text{appr}}(t)$.

II.1 Geometric properties of the exact

evolution operator

An operator $\hat{U}$ on a Hilbert space is said to preserve the norm $\|\psi\|:=\langle\psi|\psi\rangle^{1/2}$ if $\|\hat{U}\psi\|=\|\psi\|$. For linear operators $\hat{U}$, preserving the norm is equivalent to preserving the inner product,

[TABLE]

where $\hat{U}^{{\dagger}}$ is the Hermitian adjoint of $\hat{U}$. The preservation of inner product is, therefore, equivalent to the condition that $\hat{U}^{\dagger}\hat{U}$ be the identity operator, or that

[TABLE]

Such an operator $\hat{U}$ is said to be unitary. The exact evolution operator is unitary since $\hat{U}(t)^{{\dagger}}=\exp(i\hat{H}t/\hbar)=\hat{U}(t)^{-1}$.

An operator $\hat{U}$ is said to be symplectic if it preserves the symplectic two-form $\omega(\psi,\phi)$, i.e., a nondegenerate skew-symmetric bilinear form on the Hilbert space, if $\omega(\hat{U}\psi,\hat{U}\phi)=\omega(\psi,\phi)$. In classical mechanics, conservation of the symplectic two-form has many far-reaching consequences, one of which is Liouville’s theorem—the conservation of phase space volume. In quantum mechanics, a symplectic two-form can be defined asLubich (2008)

[TABLE]

obviously, it is conserved if the inner product $\langle\psi|\phi\rangle$ itself is. The exact evolution operator is therefore symplectic.

The expectation value of energy is conserved if the evolution operator is unitary and commutes with the Hamiltonian:

[TABLE]

The exact evolution operator is unitary, and because $\hat{U}(t)=\exp(-i\hat{H}t/\hbar)$ can be Taylor expanded into a convergent series in powers of $\hat{H}$, $\hat{U}(t)$ also commutes with $\hat{H}$. As a result, the exact evolution conserves energy.

An adjoint $\hat{U}(t)^{\ast}$ of an evolution operator $\hat{U}(t)$ is defined as its inverse taken with a reversed time step:

[TABLE]

An evolution operator is said to be symmetric if it is equal to its own adjoint:Hairer, Lubich, and Wanner (2006)

[TABLE]

An evolution is time-reversible if a forward propagation for time $t$ is exactly cancelled by an immediately following backward propagation for the same time, i.e., ifHairer, Lubich, and Wanner (2006)

[TABLE]

Time reversibility in quantum dynamics is, therefore, a direct consequence of symmetry. The exact evolution operator is both symmetric and time-reversible because $\hat{U}(t)^{\ast}=\exp(-i\hat{H}t/\hbar).$

Finally, the time evolution is said to be:

(i) stableAuslander, Bhatia, and Seibert (1964); Leimkuhler and Reich (2004); Bhatia and Szegö (2006) if for every $\epsilon>0,$ there exists $\delta(\epsilon)>0$ such that

[TABLE]

(ii) attractingAuslander, Bhatia, and Seibert (1964); Bhatia and Szegö (2006) if there exists a $\delta>0$ such that

[TABLE]

(iii) *asymptotically stable *if it is both stable and attracting.

These conditions are visualized in Fig. 1. The exact evolution operator is stable but not asymptotically stable because, due to norm conservation,

[TABLE]

II.2 Loss of geometric properties by approximate

methods

In approximate propagation methods, the state $\psi(t+\Delta t)$ at time $t+\Delta t$, where $\Delta t$ is the numerical time step, is obtained from the state $\psi(t)$ at time $t$ by applying an approximate time evolution operator $\hat{U}_{\mathrm{appr}}(\Delta t)$. This operator is

[TABLE]

in the explicit Euler method and

[TABLE]

in the implicit Euler method. Both Euler methods are of the first order in the time step $\Delta t$, and both are neither unitary nor symplectic. Due to their lack of unitarity, the methods do not conserve energy, even though their evolution operators commute with the Hamiltonian. Neither method is symmetric; in fact, they are adjoints of each other. Hence, neither method is time-reversible. The explicit Euler method is unstable with the distance between two wavefunctions diverging,

[TABLE]

whereas the implicit Euler method is asymptotically stable.

The second-order differencing methodKosloff and Kosloff (1983a); Leforestier et al. (1991); Askar and Cakmak (1978) recovers symmetry by combining a forward and backward steps of the explicit Euler method:

[TABLE]

The method can be also obtained directly from the time-dependent Schrödinger equation by using a finite-difference approximation

[TABLE]

While it is almost as simple as the explicit Euler method to implement, the second order differencing has a higher order of accuracy and, in contrast to the explicit Euler method, it is symmetric, time-reversible, and at least conditionally stable, meaning that it remains stable for sufficiently small time steps $\Delta t$. The second order differencing does not conserve the inner product, norm, energy, or symplectic two-form exactly. Yet, it conserves quantities analogous to the inner product [see Eq. (49)], norm,Kosloff and Kosloff (1983a); Leforestier et al. (1991) energyKosloff and Kosloff (1983a); Leforestier et al. (1991) [see Eq. (60)], and symplectic two-form [see Eq. (56)]. The corresponding exact quantities are conserved only up to the fourth order in $\Delta t$ (see Propositions 5 and 6 of Appendix A).

The properties of the different methods are summarized in Table 1; a more thorough justification of these properties is provided in Appendix A. Although the explicit and implicit Euler methods are not geometric, composing them in a specific way leads to arbitrary-order integrators that preserve many important geometric properties of the exact solution. Obviously, the compositions are applicable to systems with non-separable Hamiltonians just like the elementary methods themselves.

II.3 Recovery of geometric properties by

composed methods

Composing the explicit and implicit Euler methods, each for a time step $\Delta t/2,$ yields a symmetric second-order method (see Proposition 7 of Appendix A). Depending on the order of composition, one obtains either the trapezoidal rule

[TABLE]

or implicit midpoint method

[TABLE]

The trapezoidal rule is also known as the Crank-Nicolson method,Crank and Nicolson (1947) although the latter frequently implies a second-order finite-difference approximation to the spatial derivative in the kinetic energy operator, whereas we use the dynamic Fourier method (see Sec. II.5), which has exponential convergence with grid density (see Appendix B).

Both the trapezoidal rule and implicit midpoint methods are Cayley transformsGolub and Van Loan (1996) of $(i\Delta t/2\hbar)\hat{H}$ and, therefore, are unitary; in addition, both are second-order, symplectic, symmetric, time-reversible, and stable regardless of the size of the time step. Both methods also commute with the Hamiltonian, are energy conserving, and can be further recursively composed to obtain arbitrary-order methods (see Sec. II.4). The summary of the properties is given in Table 1 and a detailed justification provided in Appendix A.

It is necessary to stress that the geometric properties of the trapezoidal rule and implicit midpoint method are only preserved if the implicit step, which involves solving a set of linear equations, is executed exactly (or, in practice, to machine accuracy). We solved the system of equations using the generalized minimal residual method,Press et al. (2007); Saad and Schultz (1986); Saad (2003) an iterative method based on Arnoldi process.Arnoldi (1951); Saad (1980) It was an appropriate choice since the matrix being inverted was not positive-definite, symmetric, skew-symmetric, Hermitian, or skew-Hermitian, and therefore neither conjugate gradient nor minimal residual method was applicable.Saad (2003) The initial guess for the implicit step was approximated with the explicit Euler method since for small time steps, the solutions from the explicit and implicit Euler methods differ only by $(\Delta t/\hbar)^{2}\hat{H}^{2}|\psi(t)\rangle$.

II.4 Symmetric composition schemes for

symmetric methods

Recursively composing symmetric methods with appropriately chosen time steps leads to symmetric integrators of arbitrary orders.Hairer, Lubich, and Wanner (2006); Suzuki (1990); Yoshida (1990) More precisely, there exist a natural number $M$ and real numbers $\gamma_{n}$, $n=1,\ldots,M$, called composition coefficients, satisfying $\sum_{n=1}^{M}\gamma_{n}=1$ and such that if $\hat{U}_{p}(\Delta t)$ is any symmetric integrator of (necessarily even) order $p$, then

[TABLE]

is a symmetric integrator of order $p+2$. The most common composition schemes (see Fig. 2) are the triple jumpCreutz and Gocksch (1989); Forest and Ruth (1990); Suzuki (1990); Yoshida (1990) with $M=3$,

[TABLE]

and Suzuki’s fractalsSuzuki (1990) with $M=5$,

[TABLE]

The remaining coefficients are obtained from the relation $\gamma_{M+1-n}=\gamma_{n}$, which expresses that both of these are symmetric compositions.

Because each triple jump is formed of three steps while each Suzuki’s fractal is composed of five steps, the $p^{\mathrm{th}}$-order integrator obtained using Suzuki’s fractals has $(5/3)^{\frac{p}{2}-1}$ times more composition steps than the one obtained from the same symmetric second-order method using the triple jump. Therefore, the $p^{\mathrm{th}}$-order method obtained from Suzuki’s fractals takes $(5/3)^{\frac{p}{2}-1}$ times longer to execute per time step $\Delta t$ than does the method of the same order achieved through the triple jump. Yet, the leading order error coefficient of the $p^{\mathrm{th}}$-order integrator based on Suzuki’s fractal is smaller because the magnitude of each composition step is smaller in Suzuki’s fractal. Consequently, to achieve the same accuracy at a final time $t$, larger time steps can be typically used for calculations using Suzuki’s fractals compared to those based on the triple jump.

Non-recursive composition schemes, which require fewer composition steps and are also more efficient, have been obtained for various specific orders. We will refer to these as “optimal” schemes because they minimize the “magnitude” of composition steps. The magnitude of composition steps can be defined as either $\max_{n}|\gamma_{n}|$ or $\sum_{n=1}^{M}|\gamma_{n}|$. With either definition, Suzuki’s fractal is the optimal fourth-order scheme. The optimal sixth- and eighth-order schemes,Kahan and Li (1997) found by Kahan and Li by minimizing $\max_{n}|\gamma_{n}|$, have two more composition steps ($M=9$ and $M=17$, respectively) than the minimum number possible for the respective order; the optimal tenth-order scheme,Sofroniou and Spaletta (2005) obtained by Sofroniou and Spaletta by minimizing $\sum_{n=1}^{M}|\gamma_{n}|$, has four more ($M=35$).

Theorem. All compositions of the trapezoidal rule or implicit midpoint method are unitary, symplectic, stable, energy-conserving, and their evolution operators commute with the Hamiltonian; all symmetric compositions are symmetric and therefore time-reversible.

**Proof. **We prove the theorem in much greater generality. Indeed, a composition of any unitary operators $\hat{U}_{1}$ and $\hat{U}_{2}$ is unitary since

[TABLE]

A composition of any symplectic operators is symplectic since

[TABLE]

Proposition 3 of Appendix A shows that a composition of any operators commuting with the Hamiltonian commutes with the Hamiltonian. A composition of any energy-conserving operators conserves energy since

[TABLE]

However, a composition of two symmetric operators is, in general, not symmetric:

[TABLE]

It is symmetric if the two operators commute or if it is a symmetric composition, e.g.,

[TABLE]

Finally, a composition of time-reversible operators is not necessarily time-reversible since

[TABLE]

The composition is time-reversible if the two operators commute or if it is a symmetric composition, e.g.,

[TABLE]

II.5 Dynamic Fourier method

To propagate the wavepacket using the explicit or implicit Euler method, or one of their compositions (see Sec. II.2–II.4), only the action of the Hamiltonian operator $\hat{H}$ on $\psi(t)$ is required provided that the implicit steps are solved iteratively. The dynamic Fourier methodFeit, Fleck, and Steiger (1982); Kosloff and Kosloff (1983b, a); Tannor (2007) is an efficient approach to compute $f(\hat{x})\psi(t)$, where $f(\hat{x})$ is an arbitrary function of $\hat{x}$, which denotes either the nuclear position ($\hat{q}$) or momentum ($\hat{p}$) operator. Each action of $f(\hat{x})$ on $\psi(t)$ is evaluated in the $x$-representation (in which $\hat{x}$ is diagonal) by a simple multiplication, after Fourier-transforming $\psi(t)$ to change the representation if needed. On a grid of $N$ points, $f(\hat{x})\psi(t)$ is evaluated as $f(x_{i})\psi(x_{i},t)$, $1\leq i\leq N$, where $\psi(x,t)$ is the wavepacket in the $x$-representation and $x_{i}$ are either the position or momentum grid points.

II.6 Molecular Hamiltonian in the adiabatic basis

The molecular Hamiltonian in the adiabatic basis can be expressed as

[TABLE]

where $m$ is the diagonal $D\times D$ nuclear mass matrix, $D$ the number of nuclear degrees of freedom, $\mathbf{V}$ the diagonal $S\times S$ potential energy matrix, $S$ the number of considered electronic states, and $\mathbf{F}$ the nonadiabatic coupling vector (more precisely, a $D$-vector of $S\times S$ matrices). In Eq. (22), the dot $\cdot$ denotes the matrix product in nuclear $D$-dimensional vector space, the hat $\hat{}$ represents a nuclear operator, and the bold font indicates an electronic operator, i.e., an $S\times S$ matrix. Using the dynamic Fourier method, each evaluation of the action of $\hat{\mathbf{H}}$ on a molecular wavepacket $\bm{\psi}(t)$, which now becomes an $S$-component vector of nuclear wavepackets (one on each surface), involves $4D$ changes of the wavepacket’s representation.

In two-state models (i.e., for $S=2$), it is possible to obtain $\hat{\mathbf{H}}$ in the adiabatic representation analytically from the one in the diabatic representation,Sadygov and Yarkony (1998); Baer and Englman (1992); Hobey and McLachlan (1960)

[TABLE]

in which $\mathbf{W}(q)$ is the (real) diabatic potential energy matrix and in which the nonadiabatic vector couplings vanish. The adiabatic potential energy matrix $\mathbf{V}(q)$ is obtained by diagonalizing its diabatic analog $\mathbf{W}(q)$,

[TABLE]

and the molecular wavepacket in the adiabatic basis $\bm{\psi}(t)$ is obtained from its counterpart in the diabatic basis $\bm{\psi}_{\mathrm{diab}}(t)$ as

[TABLE]

using an orthogonal matrix

[TABLE]

with $\Delta(q)=V_{1}(q)-W_{11}(q)=-[V_{2}(q)-W_{22}(q)]$. The two adiabatic energies are given by

[TABLE]

where $\Delta W:=W_{22}-W_{11}$ and $\bar{W}:=(W_{11}+W_{22})/2$. Finally, the transformed momentum operator is

[TABLE]

By comparing with Eq. (22), we see that, in the adiabatic basis, the nonadiabatic coupling vector satisfies $\mathbf{F}(\hat{q})=\mathbf{O}(\hat{q})^{T}\mathbf{O}^{\prime}(\hat{q})$; in particular,

[TABLE]

III Numerical examples

To test the geometric and convergence properties of the integrators presented in Sections II.2–II.4, we used these integrators to simulate the nonadiabatic quantum dynamics in a two-surface modelEngel and Metiu (1989) of the NaI molecule. This one-dimensional model, motivated by the experiment by Mokhtari et al.,Mokhtari et al. (1990) has two electronic states, and therefore an analytical transformation between diabatic and adiabatic representations is available (see Sec. II.6). This allowed us to compare the proposed integrators, applied in the adiabatic basis, with the split-operator algorithm, which can only be used in the diabatic basis. Such a rigorous comparison would only be possible for a two-surface model potential because the split-operator algorithm requires the diabatization of the Hamiltonian formulated in the adiabatic representation and this cannot be done exactly for higher-dimensional ab initio potential energy surfaces with more electronic states.

Before the electronic excitation, the NaI molecule was assumed to be in the ground vibrational eigenstate of a harmonic fit to the ground-state potential energy surface at the equilibrium geometry. This vibrational wavepacket was then lifted to the excited-state surface, in order to obtain an initial Gaussian wavepacket ($q_{0}=4.9889$ a.u., $p_{0}=0$ a.u., $\sigma_{0}=0.110436$ a.u.) for the nonadiabatic dynamics. This use of the sudden approximation assumes an impulsive excitation, i.e., the simultaneous validity of the time-dependent perturbation theory and Condon and ultrashort pulse approximations during the excitation process. After that, the nonadiabatic dynamics was performed by solving the time-dependent Schrödinger equation using the dynamic Fourier method (see Sec. II.5) on a uniform grid with $2048$ points between $q=3.8$ a.u. and $q=47.0$ a.u.; Appendix B shows wavepacket represented on such a grid is converged for the duration of the dynamics. A long-enough propagation time, $t_{f}=10500$ a.u., was chosen so that the wavepacket traverses the avoided crossing and simultaneously witnesses the change of the nature of the excited adiabatic state from covalent to ionic. The top panel of Fig. 3 shows the two adiabatic potential energy surfaces as well as the initial wavepacket at $t=0$ and the final wavepacket at $t=t_{f}$. The population dynamics of NaI, displayed in the bottom panel of Fig. 3, shows that after crossing the region of highest nonadiabatic coupling, most of the wavepacket remains in the bound excited adiabatic state, while a small population transfer occurs to the dissociative ground state. The figure also confirms that the converged populations obtained with different integrators agree on the scale visible in the figure; in particular, the results obtained with integrators designed for the adiabatic basis agree with each other and also with the result of the split-operator algorithm in the diabatic basis.

To compare various integrators quantitatively, it is essential to “zoom in” and inspect the convergence error at the final time $t_{f}$; after all, the dynamic Fourier methodFeit, Fleck, and Steiger (1982); Kosloff and Kosloff (1983a); Tannor (2007) is expected to describe the wavepacket with a high degree of accuracy. In our setting, the convergence error at time $t$ is defined as the $L_{2}$-norm error $\left\|\psi_{\Delta t}(t)-\psi_{\Delta t/2}(t)\right\|$, where $\psi_{\Delta t}(t)$ denotes the wavepacket propagated with a time step $\Delta t$. We omit the split-operator algorithm, which served as a benchmark in Fig. 3, from the following analysis because this algorithm is not applicable to time propagation in the adiabatic representation. Note, however, that for separable Hamiltonians, such as the nonadiabatic Hamiltonian in the diabatic basis, the split-operator algorithms are more efficient than the present integrators of the corresponding order (see Table I and Paper II).

Figure 4 plots the convergence error as a function of the time step and confirms, for each algorithm, the asymptotic order of convergence predicted in Sections II.2–II.4. Recall that the trapezoidal rule and implicit midpoint method are obtained by composing the explicit and implicit Euler methods, and that the higher order methods are obtained from the trapezoidal rule or implicit midpoint method using the triple jump, Suzuki’s fractal, or optimal composition. The top panel of Fig. 4 compares the convergence of all methods, while, for clarity, the bottom left-hand panel only compares the different orders of composition for the Suzuki’s fractal and the bottom right-hand panel compares different composition schemes with the sixth order of convergence. (In Fig. 4 and all following figures, we have omitted the results of the implicit midpoint method and of its compositions because they overlap almost perfectly with the corresponding results for the trapezoidal rule.) It is clear that, for a given order of convergence, the prefactor of the error is the largest for the triple jump composition, Yoshida (1990); Suzuki (1990) intermediate for the optimally composedKahan and Li (1997); Sofroniou and Spaletta (2005) method, and smallest for Suzuki’s fractalSuzuki (1990) composition. To guarantee the correct order of convergence of all composed methods, the composed elementary second-order method must be exactly symmetric, which requires that the systems of linear equations arising from implicit steps must be solved numerically exactly.

In Section II.2, we mentioned the instability of the explicit Euler method and the conditional stability of the second-order differencing. Both properties are reflected in the top panel of Fig. 4 in the divergence of the errors of the two methods for large time steps. The critical time step for the second-order differencing is $\Delta t\approx 0.5$ a.u., whereas the explicit Euler method is unstable regardless of $\Delta t$ but the effect of instability is more visible for larger time steps. In contrast, the trapezoidal rule, implicit midpoint method, and their compositions are stable, but implicit, and, therefore, require the solution of systems of linear equations. These methods could not be used beyond a certain time step ($\max_{n}|\gamma_{n}|\Delta t\approx 100$ a.u. for both the trapezoidal rule and implicit midpoint method) because the iterative generalized minimal residual algorithm did not converge for very large $\Delta t$.

Convergence of the wavepacket’s phase, which is very important, e.g., in the evaluation of spectra, is shown in Fig. 5. As a measure of the convergence error of the phase, we use $\lvert\varphi_{\Delta t}-\varphi_{\Delta t/2}\rvert$, where $\varphi_{\Delta t}:=\mathrm{arg}[\psi_{\Delta t}(q_{\mathrm{\max}},t_{f})]$ and $q_{\mathrm{\max}}:=\mathrm{arg\,max}_{q}(\lvert\psi_{\Delta t}(q,t_{f})\rvert)$, i.e., $\varphi_{\Delta t}$ is the phase of the wavefunction propagated with time step $\Delta t$ at the position $q_{\mathrm{max}}$, for which the amplitude of the wavefunction achieves its maximum. Note that the order of convergence is identical to that of the wavepacket.

The efficiency of an algorithm cannot be judged solely from the convergence error for a given time step $\Delta t$ because the number of composition steps depends on the composition scheme, and, indeed, grows exponentially for the triple-jump and Suzuki’s fractal compositions. Figure 6, therefore, displays the convergence error $\left\|\psi_{\Delta t}(t)-\psi_{0}(t)\right\|$, where $\psi_{0}(t)$ is the wavepacket propagated using the optimally composed $10^{\mathrm{th}}$-order trapezoidal rule with an infinitesimal time step (in practice, $\Delta t=0.01$ a.u.), as a function of the computational (CPU) time. Among the elementary first- and second-order methods, compared in the top panel of Fig. 6, the second-order differencing, which does not require the solution of a system of linear equations, is the most efficient. Comparison of the geometric integrators based on the trapezoidal rule in the middle and bottom panels of Fig. 6 shows that the fourth-order Suzuki composition takes less CPU time to achieve convergence error as high as $10^{-2}$ than does the elementary trapezoidal rule (i.e., Crank-Nicolson method). To reach errors below $10^{-2}$, it is already more efficient to use the higher-order integrators. For a more dramatic example, note that the CPU time required to reach an error of $10^{-10}$ is roughly $1000$ times longer for the original trapezoidal rule than for its optimal $8^{\mathrm{th}}$-order composition. (In Paper II, we confirm that this gain in efficiency holds in higher dimensions by applying the compositions of the trapezoidal rule to the nonadiabatic dynamics in a three-dimensional model of pyrazine in the diabatic representation.Stock et al. (1995)) The bottom panel of Fig. 6 confirms the prediction that the optimal compositions are the most efficient among composition methods of the same order. Finally, note that the dependence of CPU time on the error in Fig. 6 is not monotonous for the integrators with implicit steps because the convergence of the numerical solution to the system of linear equations required more iterations for larger time steps; as a result, both the error and CPU time increased for time steps larger than a certain critical value (see Fig. 6).

Besides increased efficiency, another benefit of the algorithms based on the composition of the trapezoidal rule or implicit midpoint method is the conservation of the geometric properties of the exact evolution operator. Conservation of the energy, norm, symplectic two-form, and time reversibility by the trapezoidal rule and their compositions is demonstrated in Fig. 7. Time reversibility is measured by the distance of an initial state $\psi(0)$ from $\hat{U}(-t)\hat{U}(t)\psi(0)$, i.e., a state propagated first forward in time for time $t$ and then backward in time for time $t$. The tiny residual errors ($<2\cdot 10^{-12}$ in all cases) of the invariants result from accumulated numerical errors of the fast Fourier transform and generalized minimal residual algorithm. In contrast, the second-order differencing conserves energy, norm, and symplectic two-form only approximately with much larger, $O(\Delta t^{4})$ errors (see Propositions 5 and 6 of Appendix A). Although the second-order differencing is time-reversible in theory, its practical implementation is not. [For the second-order differencing to be exactly time-reversible, the wavepackets at time $t=0$ and $t=-\Delta t$ would have to be known exactly before the start of the simulation. However, because only $\psi(0)$ is typically available, $\psi(-\Delta t)$ must be approximated with explicit methods such as the second-order Runge–Kutta scheme.Kosloff and Kosloff (1983a)] None of the four geometric properties or analogous quantities is conserved by the Euler methods. The explicit Euler method is unstable regardless of $\Delta t$, and will, for long enough times $t_{f}$, result in a norm divergent to infinity [see Fig. 7(b), top panel] even for very small $\Delta t$, implying that also the wavefunction will have an error increasing beyond any bound. As for the implicit Euler method, its error of the norm converges to $-1$ because $\|\psi_{\mathrm{impl}}(t)\|\rightarrow 0$ as $t\rightarrow\infty$ [see Fig. 7(b), top panel].

IV Conclusion

We have described geometric integrators for nonadiabatic quantum dynamics in the adiabatic representation, in which the popular split-operator algorithms cannot be used due to nonseparability of the Hamiltonian into a kinetic and potential terms. The proposed methods are based on the symmetric composition of the trapezoidal rule or implicit midpoint method, and as a result, are symmetric, stable, conserve the energy exactly and, in addition, are exactly unitary, symplectic, and time-reversible. We have shown that unlike the original trapezoidal rule or implicit midpoint method, which are only second-order, their recursive symmetric compositions can achieve accuracy of arbitrary even order in the time step.

We have proven all these properties analytically as well as demonstrated them numerically on a two-surface model of NaI photodissociation. As expected, the higher-order integrators significantly sped up calculations when higher accuracy was required. For example, even to achieve a moderate wavefunction convergence error of $10^{-5}$, tenfold reduction in the computational time was observed by using higher-order methods compared to the elementary trapezoidal rule. It is probable that ChebyshevTal-Ezer and Kosloff (1984); Tal-Ezer (1989) and short iterative Lanczos schemesLanczos (1950); Park and Light (1986) would be more efficient in this and other typical systems, but these methods do not conserve exactly all the invariants conserved by the described geometric integrators.

Finally, we hope that the ability to run “geometric” quantum molecular dynamics in the adiabatic representation will be useful especially in conjunction with potential energy surfaces obtained from ab initio electronic structure calculations because this will avoid the tedious diabatization process necessary for the applicability of the split-operator algorithm.

The authors acknowledge the financial support from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement No. 683069 – MOLEQULE), and thank Nikolay Golubev and Rob Parrish for useful discussions.

Appendix A Geometric properties of various integrators

To shorten formulas, we set $\hbar=1$ and denote the increment $\Delta t$ with $\epsilon$ throughout the Appendix. The $\hbar$ can be reintroduced by replacing each occurrence of $t$ with $t/\hbar$ (and $\epsilon$ with $\epsilon/\hbar$). To analyze geometric properties of various integrators, we will use the following operator identities:

**Proposition 1. **Let $\hat{A}$ and $\hat{B}$ be invertible operators on a Hilbert space, and let $\hat{A}^{{\dagger}}$ and $\hat{B}^{{\dagger}}$ their Hermitian adjoints. Then both $\hat{A}^{{\dagger}}$ and $\hat{A}\hat{B}$ are invertible, and the following identities hold:

[TABLE]

The first property expresses the compatibility of the inverse and Hermitian adjoint operations, while the last three properties express that these two operations are involutive antiautomorphisms on the group of invertible operators. All four properties are easy to prove in finite-dimensional spaces;Halmos (1942) the proofs for infinite-dimensional spaces can be found in textbooks on advanced linear algebra or functional analysis.Halmos (1951)

**Proposition 2. **Let $\hat{A}$ and $\hat{B}$ be commuting operators on a vector space, i.e., $[\hat{A},\hat{B}]:=\hat{A}\hat{B}-\hat{B}\hat{A}=0$. If $\hat{A}$ is invertible, then $[\hat{A}^{-1},\hat{B}]=0$. If both $\hat{A}$ and $\hat{B}$ are invertible, then $[\hat{A}^{-1},\hat{B}^{-1}]=0$.

The first statement follows from the sequence of identities

[TABLE]

The second statement follows from the first by applying it twice, the second time for $\hat{A}:=\hat{B}$ and $\hat{B}:=\hat{A}^{-1}$, or directly from

[TABLE]

by using property (30).

Proposition 3. Let $\hat{A}$, $\hat{B}$, and $\hat{H}$ be operators on a vector space. If $[\hat{H},\hat{A}]=[\hat{H},\hat{B}]=0$, then $[\hat{H},\hat{A}\hat{B}]=0$.

This follows immediately from the identity $[\hat{H},\hat{A}\hat{B}]=\hat{A}[\hat{H},\hat{B}]+[\hat{H},\hat{A}]\hat{B}$.

A.1 Local error

The local error of an approximate evolution operator, defined as $\hat{U}_{\mathrm{appr}}(\epsilon)-\hat{U}(\epsilon)$, is obtained by comparing the Taylor expansion of $\hat{U}_{\mathrm{appr}}(\epsilon)$ with the Taylor expansion of the exact evolution operator:

[TABLE]

If the local error is $\mathcal{O(\epsilon}^{n+1})$, the method is said to be of order $n$ because the global error for a finite time $t=N\epsilon$ is $\mathcal{O}(\epsilon^{n})$.

For the explicit Euler method, the Taylor expansion is identical to the evolution operator (13) itself, and therefore the leading order local error is $(\epsilon\hat{H})^{2}/2$. The Taylor expansion of the implicit Euler method (14) is the Neumann seriesStewart (1998)

[TABLE]

the leading order local error is $-(\epsilon\hat{H})^{2}/2$.

The Taylor expansions of the trapezoidal rule and implicit midpoint method are obtained by composing Eqs. (13) and (34) with time steps $\epsilon/2$:

[TABLE]

the leading order local error of both methods is $i(\epsilon\hat{H})^{3}/12$.

Finally, the local error of the second-order differencing method is

[TABLE]

which is found by Taylor expanding $\psi_{\mathrm{sod}}(t-\epsilon)$, assumed to be exact, in Eq. (16) to obtain

[TABLE]

and

[TABLE]

Subtracting Eq. (33) from Eq. (38) gives the local error (36).

A.2 Unitarity

Neither Euler method is unitary because

[TABLE]

and

[TABLE]

In contrast, both the trapezoidal rule and implicit midpoint methods are unitary because

[TABLE]

(Proposition 1 was used in the first and Proposition 2 in the second line) and because

[TABLE]

(Proposition 1 was used in the first line).

The analysis of its geometric properties is simplified if the second-order differencing is represented by a $2\times 2$ propagation matrix

[TABLE]

acting on a vector of $\psi$ at two different times:Leforestier et al. (1991)

[TABLE]

Comparing the Hermitian conjugate $\hat{U}_{\mathrm{sod}}(\epsilon)^{{\dagger}}$ of $\hat{U}_{\mathrm{sod}}(\epsilon)$ with its inverse,

[TABLE]

found using $\det\hat{U}_{\mathrm{sod}}(\epsilon)=1,$ shows that the second-order differencing is not unitary.

When $\hat{U}(\epsilon)$ is not unitary, we can obtain the time dependence of the norm from

[TABLE]

For the explicit and implicit Euler methods, we find that

[TABLE]

where $\langle\hat{A}\rangle_{\psi}:=\langle\psi|\hat{A}|\psi\rangle$ denotes the expectation value of operator $\hat{A}$ in state $\psi$.

Although the second-order differencing is not unitary, a conserved quantity analogous to the inner product exists:

Proposition 4. The second-order differencing conserves the quantity

[TABLE]

The proof starts by projecting $\langle\phi_{\mathrm{sod}}(t)|$ on Eq. (16), which yields

[TABLE]

Adding the complex conjugate of Eq. (50) to the analogue of Eq. (50), in which $\psi$ and $\phi$ are exchanged, gives

[TABLE]

completing the proof. As an immediate corollary, obtained by taking $\phi=\psi$, the second-order differencing conserves the quantity $\mathrm{Re}\langle\psi_{\mathrm{sod}}(t)|\psi_{\mathrm{sod}}(t-\epsilon)\rangle$, which is an analogue of the norm.Kosloff and Kosloff (1983a)

Proposition 5. The global error of the inner product between two quantum states propagated by the second-order differencing is fourth-order in the time step, i.e., $\langle\psi_{\mathrm{sod}}(t_{f})|\phi_{\mathrm{sod}}(t_{f})\rangle-\langle\psi(0)|\phi(0)\rangle=\mathcal{O}(\epsilon^{4}).$

Assuming that the wavepackets at $t=-\epsilon$ are known exactly, Proposition 4 implies

[TABLE]

By Taylor expanding $\psi(-\epsilon)$ and $\phi(-\epsilon)$, and using Eq. (37), we obtain

[TABLE]

Rearranging the two sides gives

[TABLE]

The global error of the second-order differencing method is second-order in the time step and, therefore,

[TABLE]

Noting that under exact evolution, $\langle\psi(t_{f})|\hat{H}^{2}|\phi(t_{f})\rangle=\langle\psi(0)|\hat{H}^{2}|\phi(0)\rangle$, we obtain Proposition 5 by substituting Eq. (53) into Eq. (52).

A.3 Symplecticity

Using a shorthand notation $\omega_{\text{appr}}|_{t}:=\omega(\psi_{\text{appr}}(t),\phi_{\text{appr}}(t))$ and expressions $\hat{U}_{\text{appr}}(\epsilon)^{\dagger}\hat{U}_{\text{appr}}(\epsilon)$ from Appendix A.2 for the two Euler methods gives

[TABLE]

showing that neither first-order method is symplectic. In contrast, both the trapezoidal rule and implicit midpoint methods are symplectic because they are unitary.

Finally, the second-order differencing is strictly not symplectic, but Proposition 4 implies that the quantity

[TABLE]

analogous to the symplectic two-form, is conserved. In fact, Proposition 5 shows that the global error of the symplectic two-form is $\mathcal{O}(\epsilon^{4})$.

A.4 Commutation of the evolution operator with the

Hamiltonian

Evolution operators of both Euler methods commute with the Hamiltonian:

[TABLE]

where in Eq. (58), Proposition 2 was used. Applying Proposition 3 to $\hat{A}=\hat{U}_{\mathrm{expl}}(\epsilon/2)$ and $\hat{B}=\hat{U}_{\mathrm{impl}}(\epsilon/2)$ (or vice versa) then shows that the evolution operators of both the trapezoidal rule and implicit midpoint methods commute with the Hamiltonian. As for the second-order differencing, all entries in $\hat{U}_{\mathrm{sod}}$ are polynomials in $\hat{H}$ and hence commute with $\hat{H}$; as a result, $[\hat{H},\hat{U}_{\mathrm{sod}}]=0$ as well.

A.5 Energy conservation

Neither Euler method is unitary and hence neither conserves the energy. In contrast, both the trapezoidal rule and implicit midpoint methods conserve energy because their evolution operators are unitary and commute with the Hamiltonian.

The second-order differencing does not conserve energy exactly but a conserved quantity analogous to the energy has been defined:Leforestier et al. (1991) Applying $\langle\psi_{\mathrm{sod}}(t)|\hat{H}$ to Eq. (16) gives

[TABLE]

Because $\langle\hat{H}^{2}\rangle_{\psi_{\mathrm{sod}}(t)}$ is real, taking the real part of Eq. (59) shows that

[TABLE]

is conserved.

Proposition 6. The global error of the expectation value of energy of the quantum state propagated by the second-order differencing is fourth-order in the time step, i.e., $\langle\hat{H}\rangle_{\psi_{\mathrm{sod}}(t_{f})}-\langle\hat{H}\rangle_{\psi(0)}=\mathcal{O}(\epsilon^{4})$.

From Eq. (60),

[TABLE]

Assuming, as in the proof of Proposition 5, that $\psi(-\epsilon)$ is known exactly, by Taylor expanding $\psi(-\epsilon)$ and using Eq. (37), we obtain

[TABLE]

Invoking Eq. (53) and identity $\langle\hat{H}^{3}\rangle_{\psi(t_{f})}=\langle\hat{H}^{3}\rangle_{\psi(0)}$ completes the proof of Proposition 6.

A.6 Symmetry

Proposition 7. The adjoint of an evolution operator has the following properties:

[TABLE]

The first and second properties mean, respectively, that the adjoint operation ∗ is an involution and an antiautomorphism on the group of invertible operators, while the last property provides the simplest recipe for constructing a symmetric method—by composing a general method with its adjoint, with both composition coefficients of $1/2$. All three properties follow directly from the definition: the first because $(\hat{U}(\epsilon)^{\ast})^{\ast}=(\hat{U}(-\epsilon)^{\ast})^{-1}=(\hat{U}(\epsilon)^{-1})^{-1},$ the second because $(\hat{U}_{1}(\epsilon)\hat{U}_{2}(\epsilon))^{\ast}=(\hat{U}_{1}(-\epsilon)\hat{U}_{2}(-\epsilon))^{-1}=\hat{U}_{2}(-\epsilon)^{-1}\hat{U}_{1}(-\epsilon)^{-1}$, and the third by applying Eq. (64) to the product of $\hat{U}$ and $\hat{U}^{\ast}$, and using Eq. (63).

The explicit and implicit Euler methods are adjoints of each other, which follows from

[TABLE]

and Eq. (63). Therefore, neither method is symmetric. In contrast, the trapezoidal rule and implicit midpoint methods are both symmetric, which follows from Eq. (65) applied to the composition of the explicit and implicit Euler methods with composition coefficients $1/2$.

Taking the inverse of $\hat{U}_{\mathrm{sod}}(-\epsilon)$ gives

[TABLE]

implying that the second order differencing is symmetric if the sequence of wavefunctions is reversed when taking the inverse.

A.7 Time reversibility

For an elementary time step $\epsilon$, time reversibility is a direct consequence of the symmetry of the operator: if the operator is symmetric, i.e., if $\hat{U}(-\epsilon)^{-1}=\hat{U}(\epsilon)$, then a forward propagation is exactly cancelled by the immediately following backward propagation:

[TABLE]

This argument is easily extended, by induction, to a forward propagation for $N$ steps followed by a backward propagation for $N$ steps:

[TABLE]

As a result, the Euler methods are not time-reversible, whereas the trapezoidal rule, implicit midpoint, and second-order differencing methods are.

A.8 Stability

The explicit Euler method is unstable because, using Eq. (47),

[TABLE]

as long as $\hat{H}$ has no eigenvalue in the finite interval $(-E_{\text{min}},E_{\text{min}})$; composing the above inequality $N$ times shows that

[TABLE]

as $N\rightarrow\infty$ for $\psi(0)\neq\phi(0)$.

Asymptotic stability of the implicit Euler method follows, using Eq. (48), from an analogous inequality

[TABLE]

which implies

[TABLE]

as $N\rightarrow\infty$.

Both the trapezoidal rule and implicit midpoint methods are unitary, and therefore

[TABLE]

as a result, both methods are stable but not asymptotically stable.

Following Leforestier et al.Leforestier et al. (1991) and slightly abusing notation, the stability of the second-order differencing is analyzed by examining the eigenvalues

[TABLE]

of $\hat{U}_{\mathrm{sod}}(\epsilon)$. For the method to be stable, the eigenvalues must be complex units (i.e., $\left|\lambda_{1,2}\right|=1$), which is equivalent to requiring

[TABLE]

Otherwise, the magnitude of one of the eigenvalues is greater than unity and the method is unstable.Leforestier et al. (1991) For the stability criterion to be met, the condition (75) must be satisfied for all energy eigenstates and, therefore, the method is stable only for time stepsLeforestier et al. (1991)

[TABLE]

where $E_{\mathrm{max}}$ is the largest eigenvalue of the Hamiltonian operator approximated by a finite matrix.

Appendix B Exponential convergence with grid density

Figure 8 demonstrates the exponential convergence of the wavefunction with the increasing number of grid points. In order to have balanced position and momentum grids, the ranges as well as the densities of both the position and momentum grids were increased by a factor of $\sqrt{2}$ for every increase in the number of grid points by a factor of two. Comparison of wavepackets on grids with different densities was carried out by trigonometric interpolation of the wavepacket on the sparser grid.