Symplectic time-average propagators for the Schrödinger equation with a time-dependent Hamiltonian

Several symplectic splitting methods of orders four and six are presented for the step-by-step time numerical integration of the Schr¨odinger equation when the Hamiltonian is a general explicitly time-dependent real operator. They involve linear combinations of the Hamiltonian evaluated at some intermediate points. We provide the algorithm and the coefﬁcients of the methods, as well as some numerical examples showing their superior performance with respect to other available schemes. Published by AIP Publishing. [http://dx.doi.org/10.1063/1.4978410]


I. INTRODUCTION
Many quantum mechanical problems of physical interest involve time-dependent potentials.One has to deal then with the time-dependent Schrödinger equation ( = 1) where the Hamiltonian is given by where µ is the reduced mass, ∆ is the Laplacian operator, and V (t) depends explicitly on time.If the potential changes slowly with time, a usual technique for getting approximate numerical solutions consists of subdividing the time integration interval in a number of sufficiently short subintervals of length τ so that V is approximately constant in each subinterval and then solve (1) and (2) with the time-varying potential V (t) replaced by with time-independent potentials V k .In this way, one gets ψ k ≈ ψ(•, t k ) at times t k = kτ, k = 1, 2, 3, . . .and eventually reaches the desired final time.The time-independent potential V k can be chosen in several ways: as the average on the subinterval, V k = (1/τ) ∫ t k t k−1 V (t)dt, as the value at the midpoint, V k = V (t k−1 + τ/2), etc. [1][2][3] The procedure typically requires discretizing first the Schrödinger equation (1) in space, for which several techniques are widely used: finite differences, spectral methods based on collocation with trigonometric polynomials, Galerkin method with a Hermite basis, etc, both in one or several dimensions (see Ref. 4 and references therein).Once this process has been carried out, one is led to a linear ordinary differential a) Electronic mail: serblaza@imm.upv.es.b) Electronic mail: Fernando.Casas@uji.es.c) Electronic mail: Ander.Murua@ehu.es.

equation (ODE) of the form
where H(t) is a Hermitian piecewise constant N × N matrix, and each H k = T + V k is a discretized counterpart of the operator H k ≡ T + V k .Its solution u(t) ∈ C N represents a fully discretized version of the wave function ψ(x, t) at N space grid points, with N usually a large number.The values u k = u(kτ) of the solution u(t) of ( 4) evaluated at the time grid points t = k τ (k = 1, 2, . ..) can be computed as but instead of evaluating exactly the matrix exponentials in (6) for determining u k , two alternatives are widely used: (i) the exponential is approximated by the so-called unitary split operator algorithms, i.e., compositions of the form where {a i , b i } are appropriately chosen real coefficients, [5][6][7][8] and (ii) the action of the exponential on the vector u k 1 is approximated as where P m (y) is a polynomial of degree m in y that approximates the exponential e −i y .Standard choices for such P m (y) are truncated Taylor or Chebyshev series expansion of e −i y for an appropriate real interval of y or a Lanczos approximation. 9hen H k is an arbitrary real symmetric matrix, the use of symplectic splitting methods has allowed the present authors to develop an algorithm that approximately computes e −i τ H k for given τ ∈ R and ∈ C N within a prescribed error tolerance with a minimum number of matrix-vector products H k .The algorithm is more efficient than Chebyshev methods under the same conditions for all tolerances and time integration intervals. 10f the Hamiltonian does not evolve so slowly with time, the previous procedure is not particularly appropriate, however.Notice that either averaging Ĥ on each subinterval or taking each value at the midpoint produces only a second order of approximation in τ, whereas both the Chebyshev methods and the algorithm based on symplectic splitting methods have been designed to be used with large time steps.In consequence, the performance of such methods for approximately computing e −i τ H k degrades with a reduced τ, whereas large values of τ produce results with low accuracy due to a poor approximation of the original time-dependent potential V (t) by the piecewise constant potential (3).In that case, it would be more appropriate to use the numerical integrators designed to solve a linear system of the form (4) with H(t), a continuously time-dependent matrix.
Many algorithms specifically designed for Equation (4) exist indeed in the literature.We mention in particular splitoperator methods, 11 Runge-Kutta 12 and symplectic partitioned Runge-Kutta methods, 13 and a combination of a 4th-order Magnus method with the Lanczos algorithm 2 and the so-called (t, t ) method. 14In the paper, 2 an exhaustive comparison is carried out of these schemes.Polynomial approximations to the propagator based on Taylor 15 and Chebyshev polynomials [16][17][18][19] have also been proposed to deal with this problem.The review 19 provides a thorough presentation of one such class of polynomial approximations which are in principle also valid for nonlinear versions of the Schrödinger equation.
Motivated by the excellent performance exhibited by symplectic splitting methods for the computation of e −i τ H for real symmetric matrices H, our purpose here is to adapt that methodology to the more general case of a time-dependent matrix H(t) by using some techniques introduced in Refs.20 and 21.The resulting algorithms only involve matrix-vector products of the form H(t) at different values of t.Here we mainly focus on the application to linear ordinary differential equation (4) obtained by a space discretization of a Hamiltonian operator Ĥ(t) in ( 1) and ( 2) with a time-dependent potential, although the integration schemes presented in this paper are valid for arbitrary linear ODE systems of the form (4) with a time-dependent real matrix H(t), provided that H(t) is reasonably smooth in t.

II. SYMPLECTIC SPLITTING METHODS
Before presenting the proposed numerical integrators for systems of the form (4), we first recall the relatively simpler case of a constant real matrix H(t) ≡ H and the existing schemes within the same family designed to approximate u(t) = e −i t H u 0 .

A. Time-independent Hamiltonian
The class of symplectic splitting methods developed in Refs.22 and 23 as an alternative to Chebyshev polynomial approximations of e −i τ H is obtained by applying special purpose integrators to the equation expressed in terms of the real and imaginary part of the vector u ∈ C N , that is, so that Eq. ( 9) is equivalent to Application of a splitting scheme with appropriate coefficients leads (see Refs. 10 and 24-26) to an approximation of the solution operator of (10) of the form If this procedure is applied to get a numerical approximation of u(t) in the time interval t ∈ [0, t f ] using M steps with time step τ = t f /M, then the following algorithm results: The coefficients of method (12) usually verify a m+2j = a j and b m+1j = b j .In other words, composition ( 12) is left-right palindromic so that the scheme preserves time-symmetry.Although it is neither unitary nor unconditionally stable (as is the case of method ( 7)), it is symplectic and conjugate to a unitary scheme so that neither the average error in energy nor the norm of the approximate solution grows with time.Symplectic integrators preserve most qualitative properties of the exact solution and belong to a class of geometric numerical integrators [27][28][29] showing, in general, a more favorable error propagation than standard (e.g., Runge-Kutta, Taylor, multistep, etc.) integrators.
The error analysis carried out in Ref. 26, in particular the sharp error estimates obtained there, provides appropriate criteria to construct splitting methods of the form (12) possessing a large stability interval and especially adapted to different regularity conditions by optimizing the main error terms.This analysis has allowed the present authors to construct a set of symplectic splitting methods designed to be used under different precision requirements and time intervals, as well as a technique to select among them the most efficient scheme to approximate e −iτ H with the desired accuracy and a reduced computational cost.The resulting algorithm is presented in Ref. 10, where its performance is also compared with methods based on Chebyshev polynomials.

B. Time-dependent Hamiltonian
In the following, we address the more general problem (4) with a generic time-dependent real matrix H(t) and show how the previous splitting scheme (12) can be generalized to this setting.Here we only provide a brief summary of the main steps involved and refer the reader to Refs.20 and 21 for a more detailed treatment.
We again express u ∈ C N in terms of the real and imaginary part so that Eq. ( 4) is equivalently written as Let K(t, τ) denote the solution operator so that any solution (q(t), p(t)) of ( 14) is formally expressed as for all t, τ ∈ R.
To determine the solution u(t) of the linear ODE system ( 14) at times t = t k (k = 1, 2, 3, . ..), we will obtain approximations (starting from q 0 = Re(u 0 ), p 0 = Im(u 0 )), where K(t, τ) is some approximation of the solution operator K(t, τ).More specifically, our proposal is to take K(t, τ) as where for some n and appropriately chosen coefficients a i, j , b i, j , c j ∈ R. Notice that the composition ( 15) is similar in form to (12), but now the matrix H in each factor in (15) is replaced by a different (weighted) time-average of H(t).The approximation ( 15) is said to have m stages.When schemes (15) and ( 16) with n = 3 are considered as the approximate solution operator, then the following algorithm results for determining the approximations u k ≈ u(t k ), k = 1, 2, 3, . . .

C. New schemes of order 4 and 6
We present next several time-symmetric schemes ( 15) and ( 16) of orders 4 and 6.In all of them, n = 3 and the coefficients c j , j = 1, 2, 3, correspond to nodes of the Gauss-Legendre (GL) quadrature rule of order 6, that is, The remaining coefficients {a ij , b ij } must be determined to guarantee that (19) holds for r = 6.As a matter of fact, condition (19) translates into a set of polynomial equations on the coefficients (the so-called order conditions), whose solutions provide actual integration schemes.These equations have been obtained in Ref. 20 and are collected, for completeness, in a more simplified but equivalent form in Subsection 1 of the Appendix.
The process of obtaining appropriate values for the coefficients a ij , b ij is made simpler by noticing that the schemes (15) and (16) reduce in fact to (12) when applied to a constant matrix H(t) = H as long as It is then clear that a necessary condition for schemes (15) and (16), to be of order r for arbitrary smooth matrices H(t), is that (12) with coefficients (21) be an approximation of order r to the solution operator of system (10).In the particular case of 6th-order time-symmetric methods, such conditions are the first six order conditions listed in Subsection 1 of the Appendix, which only depend on the coefficients a i , b i .Hence, a convenient starting point in the construction process of 6thorder integrators of the form ( 15) and ( 16) for Equation ( 14) is to design previously a symplectic splitting method of the form (12) for the corresponding autonomous problem (10).In other words, to determine appropriate values for the coefficients (11)  satisfying the first six order conditions given in Subsection 1 of the Appendix.Once this is done, one has to solve the remaining 18 order conditions and (21) for the coefficients a i,j , b i,j .We have explored compositions (12) ranging from 8 to 24 stages and different orders, and tried to obtain time-symmetric schemes (15)  and ( 16) based on them.In many cases, the remaining order conditions cannot be solved for real valued coefficients a i,j and b i,j satisfying (21).In particular, all the 6th-order schemes (12)  with m = 10 stages we have analyzed lead to complex coefficients a i,j and b i,j .The most efficient schemes we have found are time-symmetric compositions of order 6 involving m = 11 stages.We have chosen as a representative the following: • An 11-stage scheme, denoted in the sequel as SM [6] 11 .
• A second 11-stage scheme whose coefficients provide in the autonomous case an 8th-order approximation.It will be denoted as SM [8] 11 .We have constructed other compositions with more stages, but they turn out to be clearly less efficient than SM [6] 11 and SM [8] 11 on all the numerical tests (not reported here) we have carried out with them.For completeness, we have also considered 8-stage schemes ( 15) and ( 16) of order 4 that is of order 6 for constant H(t) = H, which we denote as SM [4] 8 .The coefficients of the three proposed integrators are collected in Subsection 2 of the Appendix.The same coefficients with more significant digits as well as their implementation in several examples can be found in Ref. 30.
Although only schemes with coefficients c j , j = 1, 2, 3 corresponding to the nodes (20) of the 6th-order Gauss-Legendre quadrature rule, are considered here, other choices are of course possible.In Subsection 3 of the Appendix, we provide the necessary transformations to adapt the present schemes to any other quadrature rule of order r ≥ 6.
It is worth remarking that the m-stage splitting schemes presented here are particularly advantageous with respect to other time integrators (such as those considered in Section III) when the computational cost of evaluating the product H(t 1 )u is essentially the same as the cost of the product (H(t 1 ) + H(t 2 ))u.This is certainly the case when the Hamiltonian is of the form with H i constant and s N, as occurs in many applications.Notice that the algorithm requires 2m products with real vectors, and this is equivalent to m products with complex vectors (and so roughly the same cost as approximating the action of the exponential by a polynomial of degree m).In addition, only one new real vector, , needs to be kept in memory.

III. NUMERICAL ILLUSTRATIONS
We now illustrate the performance of the proposed methods on a pair of academic examples, namely, a high dimensional generalization of the Rosen-Zener model and the well-known Walker-Preston of a diatomic molecule.Although admittedly simple, they nevertheless represent adequately many typical applications and thus may serve as an indicative test bench of the methods proposed in this paper.
In our experiments, we compare with other well known schemes adapted to the linear problem (14).In particular, and guided by the conclusions of the analysis carried out in Ref. 2, we consider some highly efficient symplectic partitioned Runge-Kutta (splitting) methods.For our first example, we also include the efficiency diagrams obtained with a pair of explicit Runge-Kutta schemes as a standard reference.Specifically, the following integrators are taken for the numerical experiments: FIG. 1. Two-norm error in the evolution operator at the final time for the high dimensional Rosen-Zener model (23) in double logarithmic scale.FIG. 2. Two-norm error in the vector solution of the discretized onedimensional SE (24) versus the number of FFT calls in double logarithmic scale.
• RK [6] 7 : The 7-stage 6th-order explicit RK method that uses the Lobatto quadrature rule at internal stages to advance the time.
• S [4] 6 and S [6] 10 : The 6-stage 4th-order and the 10-stage 6th-order splitting methods, respectively, designed in Ref. 31 for general separable systems and adapted to the present problem by considering the time as two new dependent variables as proposed in Ref. 13.In this way, the linear problem ( 14) is transformed into a nonlinear autonomous system.One step of these schemes is evaluated as follows: (22)   with d i = i j=1 a j , e i = i−1 j=0 b j and b 0 = 0.These methods are independent of the Hamiltonian structure and have been tested in Ref. 2, showing a high performance.
In all cases, the computational cost of the algorithms is estimated as the number of products required to produce the final result.

A. A high dimensional Rosen-Zener model
For our first illustration, we consider a generalization of the well known Rosen-Zenner model for a quantum system of two levels 32 closely related to the model studied in Ref. 33.Specifically, the system is directly formulated as ( 4) with Here I k is the k × k identity matrix, D [k ] is a tridiagonal symmetric matrix with entries D [k] i,i = 0, D [k] i+1,i = D [k] i,i+1 = 1, i = 1, . . ., k − 1, D [k] k,k = 0, ω(t) is a time-dependent function, and σ i are the Pauli matrices.
We take N = 2k = 20, ω(t) = 0 + ε cos(δt), V (t) = V 0 / cosh(t/T 0 ), and 0 = 5, V 0 = 1/2 and integrate along the interval t ∈ [t 0 , t f ] with t 0 = 2 and t f = t 0 + 8π.We study the performance of the methods in the following cases: The results obtained are shown in Figure 1.We observe that the new methods show the best performance and their superiority increases when the explicitly time-dependent functions are smooth.In that case, the main error originates from terms associated with the splitting method for the autonomous case, and the new methods are optimized for this case.

B. The Walker-Preston model
This constitutes a standard model for a diatomic molecule in a strong laser field. 34The system is described by the onedimensional Schrödinger equation (in units such that = 1) with ψ(x, 0) = ψ 0 (x).Here V (x) = D 1 − e −αx 2 is the Morse potential and f (t)x = A cos(ω(t))x accounts for the laser field.
As explained in the Introduction, the usual procedure in its numerical treatment consists of defining the wave function ψ in a certain domain x ∈ [x 0 , x N ] that is subdivided into N parts of length ∆x = (x N − x 0 )/N with x i = x 0 + i∆x, and then the vector u(t) ∈ C N with components u i = (∆x) 1/2 ψ(x i−1 , t), i = 1, . . ., N is formed.Notice that the 2-norm u(t) remains constant for all t and, in order to use the Fast Fourier Transform (FFT) algorithm, periodic boundary conditions ψ(x 0 , t) = ψ(x N , t) are assumed.
Here we take the interval x ∈ [−0.8, 4.32], subdivided into N = 64 and N = 128 parts of length ∆x = 0.08 and ∆x = 0.04, respectively, and the parameters are chosen as µ = 1745 a.u., D = 0.2251 a.u., and α = 1.1741 a.u.(corresponding to the HF molecule).For each choice of N, we take for time-dependent interaction A = A 0 = 0.011 025 a.u. and laser frequency ω = ω 0 = 0.017 87.The numerical experiments are then repeated with A = A 0 /2 and ω = ω 0 /2, corresponding to a reduction in the intensity and frequency of the laser.As initial conditions we take the ground state of the Morse potential, and σ is a normalizing constant, and integrate over the time interval t ∈ [0, 10 T 0 ] with T 0 = 2π/ω (the choice ω = ω 0 /2 requires a twice longer time integration).
To check accuracy, we measure the two-norm error in the solution at the final time, u(10 T 0 ).As usual, the exact solution is accurately approximated using a sufficiently small time step.
Figure 2 shows the efficiency plots for the methods.The largest time step (i.e., the smaller number of FFTs) corresponds to the stability limit of the method (an overflow or an exceedingly large error appears when the time step is slightly increased).We observe that the relative performance of the new methods increases both when taking a finer mesh or when the time-dependent interaction is weaker and smoother.We observe that in such limit, the 6th-order method that corresponds to an 8th-order one in the autonomous case shows the best performance when high accuracies are desired.The 4thorder method shows less accuracy but allows for longer time steps and more frequent outputs.

IV. CONCLUDING REMARKS
We have presented several symplectic splitting methods for the numerical integration in time of the Schrödinger equation previously discretized in space when the Hamiltonian is an explicitly time-dependent real operator.Our approach differs from others existing in the literature in that we do not need to impose beforehand a slow time-dependence in the potential so that it can be taken as approximately constant in each time subinterval along the integration.It is based on the fact that the semidiscretized equation can be formulated as a classical linear Hamiltonian system for which the solution operator can be approximated as a product of matrices involving linear combinations of the (semidiscretized) time-dependent Hamiltonian on quadrature nodes at each time step.Here we have considered methods of orders four and six based on the Gauss-Legendre quadrature rule of order six, but the extension to other quadrature rules (as shown in Subsection 3 of the Appendix) is rather straightforward.Although originally intended for the numerical integration of the Schrödinger equation with a time-dependent potential, the new methods are still valid for arbitrary linear ODE systems of the form u = −iH(t)u for an arbitrary real matrix H(t) provided its time dependence on t is reasonably smooth.Moreover, the proposed method could also be used to solve systems like x = M(t)y, y = N(t)x arising, e.g., in elastic waves analysis. 35

Order conditions
To derive the order conditions for the schemes (15) and (16) with n = 3 and the coefficients c j given by (20), it is advantageous to first write (15) in terms of a new set of variables in order to avoid redundancies and limit the presence of contributions at orders higher than six.Specifically, instead of considering the matrices H i ≡ H(t + c i τ), i = 1, 2, 3, so that τH i = O(τ), we introduce the variables As a matter of fact, the γ i are the (i 1)th derivatives, evaluated at the midpoint t + τ/2 and divided by (i 1)!, of a matrix that interpolates H(t) at the nodes.Then, schemes ( 15) and ( 16) can be expressed as with b m+1, j = b m+1 = b m+1 = b m+1 = 0, j = 1, 2, 3.In the autonomous case, we have that γ 1 = τ H, γ 2 = γ 3 = 0 and the composition ( 12) is recovered.If the composition satisfies the symmetry conditions (18), the scheme will be of order 6 when applied to a system of the form ( 14) whenever the following conditions are satisfied (see Ref. 20 for the derivation of equivalent conditions for symmetric schemes (12) of order 6):