An efficient algorithm based on splitting for the time integration of the Schrödinger equation

We present a practical algorithm based on symplectic splitting methods intended for the numerical integration in time of the Schrödinger equation when the Hamiltonian operator is either time-independent or changes slowly with time. In the later case, the evolution operator can be effectively approximated in a step-by-step manner: first divide the time integration interval in sufficiently short subintervals, and then successively solve a Schrödinger equation with a different time-independent Hamiltonian operator in each of these subintervals. When discretized in space, the Schrödinger equation with the time-independent Hamiltonian operator obtained for each time subinterval can be recast as a classical linear autonomous Hamiltonian system corresponding to a system of coupled harmonic oscillators. The particular structure of this linear system allows us to construct a set of highly efficient schemes optimized for different precision requirements and time intervals. Sharp local error bounds are obtained for the solution of the linear autonomous Hamiltonian system considered in each time subinterval. Our schemes can be considered, in this setting, as polynomial approximations to the matrix exponential in a similar way as methods based on Chebyshev and Taylor polynomials. The theoretical analysis, supported by numerical experiments performed for different time-independent Hamiltonians, indicates that the new methods are more efficient than schemes based on Chebyshev polynomials for all tolerances and time interval lengths. The algorithm we present automatically selects, for each time subinterval, the most efficient splitting scheme (among several new optimized splitting methods) for a prescribed error tolerance and given estimates of the upper and lower bounds of the eigenvalues of the discretized version of the Hamiltonian operator.