A time-reversible integrator for the time-dependent Schr\"{o}dinger equation on an adaptive grid

One of the most accurate methods for solving the time-dependent Schr\"{o}dinger equation uses a combination of the dynamic Fourier method with the split-operator algorithm on a tensor-product grid. To reduce the number of required grid points, we let the grid move together with the wavepacket, but find that the na\"ive algorithm based on an alternate evolution of the wavefunction and grid destroys the time reversibility of the exact evolution. Yet, we show that the time reversibility is recovered if the wavefunction and grid are evolved simultaneously during each kinetic or potential step; this is achieved by using the Ehrenfest theorem together with the splitting method. The proposed algorithm is conditionally stable, symmetric, time-reversible, and conserves the norm of the wavefunction. The preservation of these geometric properties is shown analytically and demonstrated numerically on a three-dimensional harmonic model and collinear model of He-H$_{2}$ scattering. We also show that the proposed algorithm can be symmetrically composed to obtain time-reversible integrators of an arbitrary even order. We observed $10000$-fold speedup by using the tenth- instead of the second- order method to obtain a solution with a time discretization error below $10^{-9}$. Moreover, using the adaptive grid instead of the fixed grid resulted in a 64-fold reduction in the required number of grid points in the harmonic system and made it possible to simulate the He-H$_{2}$ scattering for six times longer, while maintaining reasonable accuracy. Applicability of the algorithm to high-dimensional quantum dynamics is demonstrated using the strongly anharmonic eight-dimensional H\'{e}non--Heiles model.


I. INTRODUCTION
Understanding many dynamical phenomena in chemical physics requires the solution of the time-dependent Schrödinger equation. [1][2][3][4][5][6][7] This equation can be often solved both accurately and efficiently by employing a combination of the dynamic Fourier method with a highorder split-operator algorithm 8-12 on a tensor-product grid. 8,9,13 However, for simulations that access only a small portion of the tensor-product Hilbert space, more suitable methods exist. These methods focus the available computational resources on the important portions of the full Hilbert space.
The multiconfigurational time-dependent Hartree (MCTDH) method [14][15][16] and its multilayer extension 17 reduce the required number of basis functions by employing an optimized time-dependent basis set. The recently developed pruned, collocation-based MCTDH approach 18 mitigates the exponential scaling of the MCTDH method with the number of dimensions by employing a pruned basis 19 and simultaneously extends the applicability of the MCTDH method to general potential energy surfaces by using a collocation grid. 20 Another way of reducing the required computational resources is by appropriately truncating a lattice of Gaussian basis functions [21][22][23][24] or a set of grid points. 25,26 In particular, sparse-grid methods 10,[27][28][29] reduce the number of required grid points, e.g., by employing the Smolyak quadrature. 30 Making the grid adaptive 31 is another coma) Electronic mail: jiri.vanicek@epfl.ch mon approach to reduce the required number of grid points. For example, the adaptive moving grid has been used to improve the quantum trajectory method 32,33 near wavefunction nodes, [34][35][36] to treat the interaction of molecules with intense time-dependent electromagnetic fields, 37 and to reduce the size 38,39 of grids employed in discrete variable representation (DVR), 40,41 which has been widely employed to compute vibrational spectra. 42,43 The parallelized 44,45 time-dependent DVR 39,46 has enabled quantum dynamics simulations in various multi-state and multi-dimensional systems. 47 To reduce the number of grid points and memory required for quantum simulations of systems that occupy only a small part of the accessible phase space at any given time, in this paper, we use an adaptive tensorproduct grid that moves according to the wavepacket expectation values of position and momentum. The most naïve approach is to evolve the grid only after each time step of the wavefunction propagation. However, this naïve grid adaptation breaks the symmetry and, therefore, also the time reversibility of the time propagation scheme.
We find that the time reversibility is recovered when the grid is evolved simultaneously with the wavepacket. The resulting algorithm is not only symmetric and timereversible, but also norm-preserving and conditionally stable (i.e., stable for small enough time steps). In addition, because of its symmetry, this algorithm can be composed by various symmetric composition schemes to obtain higher-order integrators. 11,12,[48][49][50][51] As well as having favorable geometric properties, the proposed algorithm is also very simple to implement because its implementation does not depend on the form of the potential energy arXiv:1909.06412v2 [quant-ph] 8 Nov 2019 function and because the grid adaptation requires no adjustable parameters.
The remainder of this paper is organized as follows: In Sec. II, we give a brief overview of the splitoperator algorithm and dynamic Fourier method, including their discretized implementation on a tensorproduct grid. Then, we demonstrate the breakdown of the time reversibility by the naïve grid adaptation, and the recovery of the time reversibility by employing a combination of the Ehrenfest theorem 52 and splitting method. 10,11 In Sec. III, we numerically confirm the geometric and convergence properties of the proposed algorithm using a three-dimensional harmonic model of electronic excitation and a two-dimensional modified Secrest-Johnson 53,54 model of He-H 2 scattering. Using the highly nonlinear eight-dimensional Hénon-Heiles system, we demonstrate that the proposed algorithm can be also used for high-dimensional quantum dynamics, which is beyond the reach of the conventional split-operator algorithm on a fixed tensor-product grid. Section IV concludes the paper.

II. THEORY
The time-dependent Schrödinger equation, whereĤ is a time-independent Hamiltonian and |ψ t is the quantum state at time t, has the solution |ψ t = U (t)|ψ 0 with the exact evolution operatorÛ (t) := e −itĤ/ . In general, this exact solution must be approximated by one of many possible time propagation schemes. Here, we will only discuss the split-operator algorithms [8][9][10] in detail because the splitting method 10,11 is crucial for the time-reversible grid adaptation that we derive in Sec. II G.
A. Split-operator algorithm and dynamic Fourier method The splitting method requires the Hamiltonian to be separable into a sum of kinetic and potential energy operators:Ĥ where p and q are D-dimensional momentum and position, respectively. For orthogonal coordinates, the kinetic energy T ( p ) has a simple form where m is a real symmetric (and often diagonal) D × D mass matrix. For Hamiltonians of form (2), all splitoperator algorithms can be expressed as a composition 11 of kinetic [ÛT (t)] and potential [ÛV (t)] evolution operators, whereÛÂ(t) := e −itÂ/ is the exact evolution operator for a HamiltonianĤ =Â. The simplest first-order split-operator algorithm 55 has the approximate evolution operatorÛ This "TV" algorithm can be composed with its adjoint, U VT (∆t) :=Û TV (−∆t) −1 , to obtain the second-order TVT algorithm, 56 which is symmetric, i.e., satisfiesÛ TVT (∆t) = U TVT (−∆t) −1 . Because it is symmetric, the TVT algorithm can be recursively composed by symmetric schemes 10,11,[48][49][50][51] to obtain symmetric algorithmŝ of arbitrary even orders of accuracy in the time step, where γ k is the kth composition coefficient (with Ncomp k=1 γ k = 1, γ Ncomp+1−k = γ k ) and N comp is the number of composition steps. [Exchanging V and T in algorithms (4)-(6) results in another set of split-operator algorithms.] The higher-order algorithms, despite their higher computational cost per time step, are often more efficient if high accuracy is desired. Moreover, all symmetric split-operator algorithms are examples of geometric integrators because they preserve many important geometric properties of the exact solution of Eq. (1), namely linearity, unitarity, symplecticity, stability, symmetry, and time reversibility. 10,11,57 For additional details about the properties and numerical implementation of the higher-order split-operator algorithms, we refer the reader to Ref. 12.
The action of compositions ofÛT (∆t) andÛV (∆t) on |ψ is simply evaluated using the dynamic Fourier method 8,9,13 in which the action of a function g(ˆ x) of an operatorˆ x on the wavepacket |ψ is evaluated as g( x)ψ( x) in the x-representation, where x is either the position q or momentum p. The wavefunction is transformed, if required, to the x-representation by either the Fourier or inverse Fourier transformation: where p · q := D l=1 p l q l .

B. Dynamic Fourier method on a grid
The dynamic Fourier method on a grid follows the same approach as in Sec. II A except that g( x)ψ( x) and the integral transforms (7) and (8) are now discretized on a grid, consisting of points x I for all I ∈ I. Here, the multi-index I = (i 1 , . . . , i D ) is an ordered D-tuple of integers from the set I of all admissible multi-indices, where and N l is the number of grid points in the lth dimension. When iterating over all admissible multi-indices, we will simply write I ∈ I. The D coordinates of the grid point x I are given by x I l := x ctr,l + (i l − N l /2)∆x l for l ∈ {1, . . . , D}, (10) where x ctr is the x-grid center and ∆ x are the x-spacings of the grid. Note that ∆q l ∆p l = 2π /N l .
Application of the operator g(ˆ x) to ψ in xrepresentation is given by where grid = denotes "is represented on a grid as." The integral transforms (7) and (8) are discretized on a grid asψ with prefactors C x := D l=1 (∆x l / √ 2π ). To express Eqs. (12) and (13) in terms of the standard discrete Fourier transform (DFT), we scale the wavefunctions asψ K :=ψ( p K )/ C q , ψ J := ψ( q J )/ C p and use Eq. (10). As a result, we obtaiñ where N := D l=1 N l denotes the total number of grid points, [(p ctr,l − ∆p l N l /2)(q ctr,l − ∆q l N l /2) + p ctr,l j l ∆q l + q ctr,l k l ∆p l − π (j l + k l )], (16) and the multi-index inner product K, J := D l=1 k l j l /N l . The scaled wavefunctionψ K can be viewed as a standard DFT of e −it KJ / ψ J and ψ J as a standard inverse DFT of e it KJ / ψK . In practice, the DFT is implemented using the celebrated fast Fourier transform algorithm, which has been parallelized via the message passing interface (MPI) and the open multi-processing (openMP) interface. 58 C. Shifting the grid In Sec. II B, we assumed the grid centers to be fixed, but this assumption will now be dropped to allow for grid shifting. Moreover, all quantities defined in Sec. II B will be re-expressed in a more compact matrix form. Equations (14) and (15) thus becomẽ where the "vectors" ψ( q ctr ) andψ( p ctr ) are rank-D tensors with N components ψ J ( q ctr ) andψ K ( p ctr ), respectively, and the "matrices" representing the Fourier transforms are In Eqs. (17) and (18), a compact notation for the "matrix-vector" multiplication, defined by (aψ) K := J∈I a KJ ψ J , K ∈ I, is employed, and˜above a matrix denotes that it is applied to a wavefunction in the p-representation.
Equation (18) expresses that the momentum wavefunctionψ( p ctr ) represented on the p-grid centered at p ctr is transformed to the position wavefunction ψ( q ctr ) represented on the q-grid centered at q ctr by applying f ( q ctr , p ctr ). Similarly, ψ( q ctr ) is transformed toψ( p ctr ) by applying f ( q ctr , p ctr ) according to Eq. (17).

D. Split-operator algorithm on a grid
Kinetic and potential evolution operatorsÛT (∆t) and UV (∆t), which are composed to obtain any split-operator algorithm (see Sec. II A), are discretized on a grid as diagonal finite-dimensional tensors Therefore, the time-evolved wavefunctions are p |ÛT (∆t)|ψ grid = C qŨT (∆t, p ctr )ψ( p ctr ); (24) note that the scaled wavefunctions, ψ J andψ K , must be scaled back to ψ( q J ) andψ( p K ) at the end of the propagation with factors C p and C q , respectively (see Sec II B).

E. Loss of linearity by the grid adaptation
To be specific, we now assume that the initial wavefunction is provided in the q-representation and that the solution at time t is desired also in the q-representation. Let the adaptive x-grid be centered at the wavefunction's x-expectation value. The resulting equations of motion for the wavefunction and grid centers, is the Hamiltonian in the q-representation, containing appropriate Fourier transforms, and represented on the grid cen- ; the inner product between ψ and φ is defined as ).
(28) The grid adaptation leads to the loss of some geometric properties even if Eqs. (25)- (27) are solved exactly. In Eq. (25), the Hamiltonian is nonlinear due to its dependence on ψ t (via q t and p t ); the corresponding evolution operator is, therefore, also nonlinear, and does not preserve the inner product. 59 As a consequence, the symplectic two-form 10 ω(ψ, φ) := −2 Im ψ|φ is not preserved, either. In contrast, the exact solution of Eqs. (25)- (27) does preserve the norm and is both symmetric and timereversible.
[In practice, U(∆t; q 0 , p 0 )ψ 0 ( q 0 ) is approximated numerically using time propagation schemes, such as the splitoperator algorithm, short iterative Lanczos scheme, 60,61 or Crank-Nicolson 62,63 method.] The grid centers are then updated using the propagated wavefunction ψ ∆t ( q 0 ): Finally, the wavefunction ψ ∆t ( q 0 ) is represented on the updated grid: The overall evolution operator for the naïve adaptive grid is, therefore, where the dependence of U naïve on ψ 0 comes from the dependence of q ∆t and p ∆t on ψ 0 . The time propagation on the naïve adaptive grid preserves the norm ψ t ( q t ) := ψ t ( q t )|ψ t ( q t ) 1/2 because U naïve is a composition of three norm-preserving operators: That f andf preserve the norm follows from Eq. (20), and U(∆t, q 0 , p 0 ) preserves the norm because A symmetric operator is time-reversible [i.e., satisfies U(−∆t)U(∆t) = 1] and vice versa. Both of these properties are lost in the naïve adaptive grid approach because where Note that inequality (35) would still hold even in the unlikely situation that, by chance, q 0 = q 0 and p 0 = p 0 . As we shall see below, to preserve the symmetry and time reversibility, the grid must be evolved simultaneously with the wavefunction.

G. Recovery of time reversibility by a combination of the splitting method and Ehrenfest theorem
The Ehrenfest theorem 52 states that the time derivatives of the position and momentum expectation values satisfy˙ The system of differential and algebraic Eqs. (25), (26), (27) for ψ t , q t , p t is equivalent to and, hence, can be replaced with the system of differential Eqs. (25), (40), (41). These equations can be solved analytically ifĤ = V (ˆ q) orĤ = T (ˆ p). This is the essence of the splitting method 10,11 (see Sec. II A).
Time reversibility of U V,adpt (∆t; ψ 0 , q 0 ) follows be- Similarly, U T,adpt (∆t; ψ 0 , q 0 , p 0 ) is time-reversible because (56) where we used Eq. (20) and the identity U T (−∆t, p 0 )U T (∆t, p 0 ) = 1. A symmetric composition of time-reversible operators is time-reversible. 64 Therefore, all symmetric split-operator algorithms of form (6) that are composed from U V,adpt (∆t; ψ 0 , q 0 ) and U T,adpt (∆t; ψ 0 , q 0 , p 0 ) are time-reversible. The Verlet algorithm applied to the harmonic oscillator is stable for time steps that satisfy where T osc is the oscillation period. 57 In higherdimensional harmonic models, the restriction (57) on the time step must hold for the period T osc of the fastest normal mode. 57

III. NUMERICAL EXAMPLES
A. Three-dimensional harmonic model To analyze the geometric and convergence properties of the algorithm proposed in Sec. II G, we devised a twosurface three-dimensional harmonic model of electronic excitation of a molecule. The initial vibrational state, determined using the ground-state potential energy surface, was propagated solely on the excited-state surface, following an impulsive electronic excitation. More precisely, the initial state for the propagation was the ground vibrational eigenstate, of the ground-state Hamiltonian, where q l is the lth ground-state normal mode coordinate, p l is its conjugate momentum, and ω l is the associated vibrational frequency. After the electronic excitation, ψ( q ) was propagated with the excited-state Hamiltonian where q 0 is the displacement of the excited-state potential energy surface and K is a symmetric positive definite matrix; K is not diagonal because the excited-state normal modes were chosen to be Duschinsky rotated 68 with respect to the ground-state normal modes. For the dynamics, natural units (n.u.) were used: = ω 2 = m H = 1, where m H is the mass of a hydrogen atom. The diagonal (K ll ) and off-diagonal (K lm ) elements of the K matrix, displacement q 0 , and ground-state vibrational frequency ω l in Eq. (60) are listed in Table I, which also contains the initial parameters of the adaptive grid and the total propagation time t f . −7 t f 50 q0,2 = q0, 3 7 To verify that grid adaptation does not decrease the accuracy of the solution, we compared the wavefunction ψ (∆t) t propagated using the adaptive grid with the time step ∆t to the corresponding "benchmark" wavefunction Ψ (∆t) t propagated using a fixed grid. Indeed, the errors ψ were minuscule (the errors were 5 × 10 −11 at t = 0, and 2 × 10 −10 at t = t f ). The wavefunctions were propagated with the optimally composed tenth-order TVT split-operator algorithm with ∆t = t f /2 9 . (See Ref. 64 and the references therein for a detailed discussion of composition schemes.) We used a high-order integrator with a small time step so that the error was dominated by grid adaptation, and not by time discretization. Both q-and p-ranges of the fixed grid were chosen to be twice larger than the ranges of  the adaptive grid because the amplitude of the adaptive grid's oscillation was approximately equal to its range. In order that the fixed and adaptive grids had the same density, the fixed grid was chosen to have 128 × 128 × 128 points (see Appendix A for the exponential convergence of the wavefunction with the number of grid points). Figure 1(a) shows that the expectation value of position is computed correctly with the adaptive grid, even when the wavefunction moves beyond the range of the initial grid. In fact, Fig. 1(b) shows that the error of the position expectation value is minuscule (of the order of 10 −11 ) for all times; the slow linear increase in the error is due to the accumulation of roundoff errors. Figure 2 demonstrates that the compositions 11,48-51 of the proposed algorithm from Sec. II G achieve the predicted higher orders of accuracy. The figure also demonstrates the divergence of the discretization errors ψ    Results of the second-order TVT algorithm were extrapolated using the line of best fit beyond CPU time = 7×10 4 s to highlight the higher efficiency of the higher-order methods. Top: all discussed methods; bottom: optimally composed methods (Suzuki's fractal is the optimal fourth-order composition scheme 50,51,64 ). Line labels are the same as in Fig. 2. second-order algorithm is much greater if a small error is desired. To reach an error of 10 −9 , 500-fold speedup is achieved by using the Suzuki fourth-order algorithm. Moreover, 2000-, 5000-, and 10000-fold speedups are achieved by using the optimal sixth-, eighth-, and tenthorder algorithms, respectively. Because a split-operator propagation is equivalent to an exact propagation with an effective, time-dependent Hamiltonian, the energy E is conserved only approximately. Figure 4(a) shows that the energy is conserved to the same order of accuracy [O(∆t m )] as the wavefunction. Figures 4(b) and (d) demonstrate that the compositions of the proposed algorithm are exactly norm-preserving and time-reversible as already justified analytically in Sec. II G. In Sec. II E, we showed that the grid adaptation leads to the nonconservation of the inner product. Figure 4(c) may, therefore, be misleading because the inner product appears to be conserved. However, this is not true in general, as shown later on the example of collinear He-H 2 scattering (see Sec. III B).
In all panels of Fig. 4 the slow increase in the error for decreasing time steps is due to the accumulation of roundoff errors; therefore, the (minuscule) errors are larger for methods with more composition steps per time step. 11 Panels (b), (c), and (d) show that, on the other hand, the errors diverge for large time steps ∆t because of the instability of the Verlet algorithm (see Sec. II H); larger errors result from methods with a larger maximum composition coefficient [max k |γ k |, see Eq. (57)]. Specifically, for large ∆t, the time reversibility of the proposed algorithm is lost because the centers of the position and momentum grids diverge due to the instability of the Verlet algorithm, used for propagating the grid centers. Similarly, beyond a certain time step size, the norm is no longer preserved because it is evaluated on a grid whose center has diverged to infinity. Figure 5 confirms that the naïve adaptive grid approach is not time-reversible while the algorithm proposed in Sec. II G is. Note that for very small time steps (∆t ≤ 10 −2 ), the solution is essentially exact, and even the naïve adaptive grid approach becomes effectively time-reversible. The bottom panel of Fig. 5 shows, however, that for a fixed time step ∆t, the time propagation on the naïve adaptive grid is not time-reversible already after a short propagation time t; the breaking of time reversibility is an inherent property of the naïve adaptive grid.

B. Collinear He-H2 scattering
As a more challenging test, we also applied the algorithm proposed in Sec. II G to a very anharmonic system. Following Ref. 38, we simulated the collinear He-H 2 scattering using a modified 54 Secrest-Johnson 53 potential energy surface, where β = 0.158 n.u., D = 20 n.u., and α = 0.3 n.u. The natural units (n.u.) are different from those defined in Sec. III A: = 1 as before, but 2Dβ 2 /m 1 = m 1 = 1 instead. In Eq. (61), q 1 is the vibrational coordinate of H 2 , and q 2 is the distance between the He atom and the center of mass of H 2 . 53 In this coordinate system, m 1 = 1 n.u. and m 2 = 2/3 n.u. 38,53 The Hamiltonian for this problem isĤ scat = T (ˆ p) + V SJ (ˆ q), where T (ˆ p) has the form (3). Identically to Ref. 38, the initial state is a product of two one-dimensional Gaussian wavepackets The Gaussian wavepacket ψ (2) (q) is sufficiently narrow and far from the interaction region so that there is no significant initial interaction between He and H 2 (σ 2 0 = 8 n.u. and q 0 = 24 n.u.). A negative initial momentum (p 0 = −3.56 n.u.) ensures a collision at a later time. Error ψt − Ψt of the wavefunction ψt propagated on either the adaptive or fixed grid (both with 128 × 128 points) in the He-H2 scattering from Sec. III B. Ψt is the "exact" reference wavefunction propagated on a fixed grid with 128 × 2048 points. The initial q1 range of all grids was (−14 n.u., 14 n.u.). The initial q2 range was (0 n.u., 48 n.u.) for both grids with 128 × 128 points and (−50 n.u., 400 n.u.) for the reference fixed grid. Figures 6-8 were produced using the optimal tenth-order composition of the VTV algorithm with ∆t = 0.1 n.u. Figure 6 shows the error of the wavefunction propagated using either the adaptive or fixed grid, both with the same number of grid points. The error of the wavefunction propagated using the adaptive grid remains reasonably small (< 10 −3 ) for six times longer than the error of the corresponding wavefunction on the fixed grid. The significance is that for a given number of grid points, determined, e.g., by the available memory, the time scale of a simulation can be extended by grid adaptation. Figure 7 displays the time dependence of the expectation value of q and of its error. Panel (a) shows that the collision between He and H 2 induces the vibration of H 2 , which was originally in its ground vibrational state. Panel (b) shows that the error of the position expectation value was reasonably small until t ≈ 20 n.u. on the adaptive grid. After t ≈ 20 n.u., however, the error starts to grow rapidly in the second dimension because the width of the wavepacket in this dimension increases approximately linearly from t ≈ 10 n.u. Thus, a significant portion of the wavepacket eventually escapes through the boundaries of the adaptive grid in the second dimension.
Finally, in Fig. 8, we show that the proposed algorithm preserves the geometric invariants even in the collinear scattering of He-H 2 , where the wavepacket is more delocalized than in the harmonic example from Sec. III A. As expected, the norm [ Fig. 8(b)] and time reversibility [ Fig. 8(d)] are preserved exactly (the slow linear increase of the invariants is again due to the accumulation of roundoff errors). On the other hand, the energy [ Fig. 8(a)] and inner product [ Fig. 8(c)] are not conserved. In particular, the apparent conservation of the inner product observed in Sec. III A is indeed not general.

C. Eight-dimensional Hénon-Heiles model
To demonstrate the applicability of the algorithm from Sec. II G to high-dimensional quantum dynamics, we applied it to the eight-dimensional Hénon-Heiles model, whose kinetic energy T (ˆ p) is of the form (3) and potential energy is given by with D = 8. Simulating the quantum dynamics of the Hénon-Heiles system is challenging because the potential V HH ( q) is anharmonic, unbound, and contains inter-mode couplings. Following Ref. 69, the parameters were chosen to be λ = 0.111803 n.u. and m l = 1 n.u. for l = 1, . . . , D, where the natural units (n.u.) were defined by fixing = m H = κ = 1; the initial state is ψ( q) = π −D/4 exp[−( q − q 0 ) 2 /2] with q 0,l = 2 n.u. for l = 1, . . . , D. Panels (a) and (b) of Fig. 9 show that the autocorrelation function ψ 0 |ψ t is obtained more accurately by the adaptive grid algorithm than by the standard splitoperator algorithm on a fixed grid when both the adaptive and fixed grids have the same number (8 8 ) of grid points. Especially, the detailed shape of the first recurrence of the autocorrelation function is correctly described only by the proposed algorithm [see insets of panels (a) and (b) of Fig. 9]. Consequently, only the spectrum calculated using the adaptive grid algorithm has the correct shape of the envelope [see Fig. 9(c)].

IV. CONCLUSION
We have described a split-operator algorithm combined with an adaptive phase space grid whose center moves according to the wavepacket's expectation values of position and momentum. By propagating the grid center exactly and simultaneously with the wavefunction, the symmetry and time reversibility were built into the proposed algorithm. Adapting the grid reduces the number of required grid points while maintaining high accuracy in situations where the wavepacket remains localized. Examples include harmonic systems or short-time dynamics in moderately anharmonic systems. On the example of He-H 2 scattering, we showed that the proposed algorithm is also suitable for longer-time dynamics if only a moderate accuracy of the wavepacket is required, i.e., when one can ignore small parts of the wavepacket escaping through the boundaries of the adaptive grid. The algorithm allowed us to compute accurately the medium-resolution spectrum of the eight-dimensional Hénon-Heiles system, which is not only high-dimensional but also highly nonlinear.
We showed both analytically and numerically that the time reversibility is lost by the naïve grid adaptation. Then, we introduced an amendment that recovered the time reversibility. The geometric properties of the resulting algorithm, namely norm preservation, conditional stability, symmetry, and time reversibility, were demonstrated analytically as well as numerically on two different model systems. Note that because neither the Chebyshev 70 nor the Lanczos 60,61 method is time with t damp = 30 n.u. as the damping function. The fixed grid was defined between q l = −3.5 n.u. and q l = 3.5 n.u., and the initial adaptive grid between q l = −3.0 n.u. and q l = 7.0 n.u.; both grids had 8 8 points. The benchmark was calculated using MCTDH. 14 reversible, 10 it is only reasonable to combine them with the naïve adaptive grid. Moreover, combining the Chebyshev method with adaptive grids would be inappropriate because the advantage of the Chebyshev algorithm as a global, long time propagator would be lost by changing the grid and, therefore, the Hamiltonian matrix regularly at small time intervals.
Because of its symmetry, the proposed algorithm can be composed to obtain higher-order integrators. We verified that these higher-order integrators are more efficient compared to the second-order integrator if high accuracy is desired. As an additional benefit, the proposed algorithm requires no adjustable parameters for the grid adaptation because the grid center follows the exact trajectory of the wavepacket's expectation values of position and momentum.
Finally, we hope that the proposed time-reversible integrator for the time-dependent Schrödinger equation on an adaptive grid could serve as a benchmark for more approximate methods, such as the thawed Gaussian approximation, 71-73 that rely on the wavepacket remaining localized for relevant time scales.
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 Tomislav Begušić and Lipeng Chen for useful discussions.
Appendix A: Exponential convergence with the number of grid points Figure 10 shows the exponential convergence of the wavefunction with the increasing number of grid points. As expected, the convergence is slower for the fixed grid, which must cover roughly four times larger phase space area in each dimension to account for the movement of the wavepacket. ψM (t f ) is evaluated on either the adaptive or fixed grid of N = M D points in the harmonic system from Sec. III A. For balanced q-and p-grids, the ranges as well as the densities of both q-and p-grids were increased by a factor of √ 2 for every doubling of the number of grid points. Wavefunctions on grids of different densities were compared through the trigonometric interpolation of the wavefunction on the sparser grid. Both the proposed adaptive grid algorithm from Sec. II G and the split-operator algorithm on a fixed grid were composed to the tenth order using the optimal scheme; ∆t = 0.1 n.u.
The fast Fourier transform, which is the bottleneck of the combination of the dynamic Fourier method with the split-operator algorithm on a tensor-product grid, scales as N log N . Because N = M D , where M is the geometric average number of grid points in each dimension, the method scales exponentially with the number of dimensions. Similarly to the MCTDH or time-dependent DVR methods, the proposed algorithm does not remove the exponential scaling, but slows down the exponential growth significantly compared to the conventional split-operator algorithm on a fixed grid (see, e.g., Fig. 10). This is, indeed, what allowed us to treat the 8-dimensional Hénon-Heiles system in Sec. III C.