Error analysis of splitting methods for the time dependent Schrödinger equation
Abstract
A typical procedure to integrate numerically the time dependent Schrödinger equation involves two stages. In the first one carries out a space discretization of the continuous problem. This results in the linear system of differential equations , where is a real symmetric matrix, whose solution with initial value is given by . Usually, this exponential matrix is expensive to evaluate, so that time stepping methods to construct approximations to from time to are considered in the second phase of the procedure. Among them, schemes involving multiplications of the matrix with vectors, such as Lanczos and Chebyshev methods, are particularly efficient.
In this work we consider a particular class of splitting methods which also involves only products . We carry out an error analysis of these integrators and propose a strategy which allows us to construct different splitting symplectic methods of different order (even of order zero) possessing a large stability interval that can be adapted to different space regularity conditions and different accuracy ranges of the spatial discretization. The validity of the procedure and the performance of the resulting schemes are illustrated on several numerical examples.

Instituto de Matemática Multidisciplinar, Universidad Politécnica de Valencia, E46022 Valencia, Spain.

Institut de Matemàtiques i Aplicacions de Castelló and Departament de Matemàtiques, Universitat Jaume I, E12071 Castellón, Spain.

Konputazio Zientziak eta A.A. saila, Informatika Fakultatea, EHU/UPV, Donostia/San Sebastián, Spain.
1 Introduction
To describe and understand the dynamics and evolution of many basic atomic and molecular phenomena, their time dependent quantum mechanical treatment is essential. Thus, for instance, in molecular dynamics, the construction of models and simulations of molecular encounters can benefit a good deal from time dependent computations. The same applies to scattering processes such as atomdiatom collisions and triatomic photodissociation and, in general, to quantum mechanical phenomena where there is an initial state that under the influence of a given potential evolves through time to achieve a final asymptotic state (e.g., chemical reactions, unimolecular breakdown, desorption, etc.). This requires, of course, to solve the time dependent Schrödinger equation ()
(1) 
where is the Hamiltonian operator, is the wave function representing the state of the system and the initial state is . Usually
(2) 
and the operators , are defined by their actions on as
The solution of (1) provides all dynamical information on the physical system at any time. It can be expressed as
(3) 
where represents the evolution operator, which is linear and satisfies the equation with . Since the Hamiltonian is explicitly time independent, the evolution operator is given formally by
(4) 
In practice, however, the Schrödinger equation has to be solved numerically, and the procedure involves basically two steps. The first one consists in considering a faithful discrete spatial representation of the initial wave function and the operator on an appropriately constructed grid. Once this spatial discretization is built, the initial wave function is propagated in time until the end of the dynamical event. It is on the second stage of this process where we will concentrate our analysis.
As a result of the discretization of eq. (1) in space, one is left with a linear equation , with a Hermitian matrix of large dimension and large norm. Since evaluating exactly the exponential is computationally expensive, approximation methods requiring only matrixvector products with are particularly appropriate [22]. Among them, the class of splitting symplectic methods has received considerable attention in the literature [13, 23, 27, 2, 3]. In this case is approximated by a composition of symplectic matrices. While it has been shown that stable high order methods belonging to this family do exist, such high degree of accuracy may be disproportionate in comparison with the error involved in the spatial discretization, and also inappropriate particularly when the problem at hand involves nonsmooth solutions, as high order methods make small phase errors in the low frequencies but much larger errors in the high frequencies. An error analysis of this family of integrators, in particular, could help one to design different efficient time integrators adapted to different accuracy requirements and spacial regularity situations.
The analysis carried out in the present paper could be considered a step forward in this direction. We present a strategy which allows us to construct different splitting symplectic methods of different order and large stability interval (with a large number of stages) that can be adapted to different space regularity conditions and different accuracy ranges of the spatial discretization. When this regularity degree is low, sometimes the best option is provided by methods of order zero.
Since the splitting methods we analyze here only involve products of the matrix with vectors, they belong to the same class of integrators as the Chebyshev and Lanczos methods, in the sense that all of them approximate by linear combinations of terms of the form ().
The plan of the paper is as follows. In section 2 we review first the Fourier collocation approach carrying out the spatial discretization of the Schrödinger equation, and then we turn our attention to the time discretization errors of symplectic splitting methods. The bulk of the paper is contained in section 3. There we carry out a theoretical analysis of symplectic splitting methods and obtain some estimates on the time discretization error. These estimates in turn allows us to build different classes of splitting schemes in section 4, which are then illustrated in section 5 on several numerical examples exhibiting different degrees of regularity. Here we also include, for comparison, results achieved by the Lanczos and Chebyshev methods. Finally, section 6 contains some conclusions.
2 Space and time discretization
Among many possible ways to discretize the Schrödinger equation in space, collocation spectral methods possess several attractive features: they allow a relatively small grid size for representing the wave function, are simple to implement and provide an extremely high order of accuracy if the solution of the problem is sufficiently smooth [11, 12]. In fact, spectral methods are superior to local methods (such as finite difference schemes) not only when very high spatial resolution is required, but also when long time integration is carried out, since the resulting spatial discretization does not cause a deterioration of the phase error as the integration in time goes on [16].
To simplify the treatment, we will limit ourselves to the onedimensional case and assume that the wave function is negligible outside an interval . In such a situation one may reformulate the problem on the finite interval with periodic boundary conditions. After rescaling, one may assume without loss of generality that the space interval is , and therefore
(5) 
with for all . In the Fouriercollocation (or pseudospectral) approach, one intends to construct approximations based on the equidistant interpolation grid
where is even (although the formalism can also be adapted to an odd number of points). Then one seeks a solution of the form [22]
(6) 
where the coefficients are related to the grid values through a discrete Fourier transform of length , [26]. Its computation can be accomplished by the Fast Fourier Transform (FFT) algorithm with floating point operations.
In the collocation approach, the grid values are determined by requiring that the approximation (6) satisfies the Schrödinger equation precisely at the grid points [22]. This yields a system of ordinary differential equations to determine the point values :
(7) 
where
(8) 
for and . Observe that the matrices on the righthand side of (7) are Hermitian.
An important qualitative feature of this space discretization procedure is that it replaces the original Hilbert space defined by the quantum mechanical problem by a discrete one in which the action of operators are approximated by (Hermitian) matrices obeying the same quantum mechanical commutation relations [18]. From a quantitative point of view, if the function is sufficiently smooth and periodic, then the coefficients exhibit a rapid decay (in some cases, faster than algebraically in , uniformly in ), so that typically the value of in the expansion (6) needs not to be very large to represent accurately the solution. Specifically, in [22] the following result is proved.
Theorem 1
When the problem is not periodic, the use of a truncated Fourier series introduces errors in the computation. In that case several techniques have been proposed to minimize its effects (see [1, 5] and references therein).
The previous treatment can be generalized to several spatial dimensions, still exploiting all the onedimensional features, by taking tensor products of onedimensional expansions. The resulting functions are then defined on the Cartesian product of intervals [6, 22].
We can then conclude that after the previous space discretization has been applied to eq. (5), one ends up with a linear system of ODEs of the form
(9) 
where is a real symmetric matrix. This is the starting point for carrying out an integration in time. Although a collocation approach has been applied here, in fact any space discretization scheme leading to an equation of the form (9) fits in our subsequent analysis.
The spatial discretization chosen has of course a direct consequence on the time propagation of the (discrete) wave function , since the matrix representing the Hamiltonian has a discrete spectrum which depends on the scheme. This discrete representation, in addition, restricts the energy range of the problem and therefore imposes an upper bound to the high frequency components represented in the propagation [19].
The exact solution of eq. (9) is given by
(10) 
but to compute the exponential of the complex and full matrix (typically also of large norm) by diagonalizing the matrix can be prohibitively expensive for large values of . In practice, thus, one turns to time stepping methods advancing the approximate solution from time to , so that the aim is to construct an approximation as a map .
Among them, exponential splitting schemes have been widely used when the Hamiltonian has the form given by (2) [9, 19, 22]. In that case, equation (9) reads
(11) 
where is a diagonal matrix associated with and is related to the kinetic energy . It turns out that the solutions and of equations and , respectively, can be easily determined [22], so that one may consider compositions of the form
(12) 
where . In (12) the number of exponentials (and therefore the number of coefficients ) has to be sufficiently large to solve all the equations required to achieve order (the so called order conditions).
Splitting methods of this class have several structurepreserving properties. They are unitary, so that the norm of is preserved along the integration, and timereversible when the composition (12) is symmetric. Error estimates of such methods applied to the Schrodinger equation [17, 24, 25] seem to suggest that, while they are indeed very efficient for high spatial regularity, they may not be very appropriate under conditions of limited regularity.
Here we will concentrate on another class of splitting methods that have been considered in the literature [13, 23, 27, 2, 3]. Notice that the corresponding in eq. (7) is a real symmetric matrix, and thus is not only unitary, but also symplectic with canonical coordinates and momenta . In consequence, equation (9) is equivalent to [13, 14]
(13) 
Alternatively, one may write
(14) 
with the matrices and given by
The solution operator corresponding to (14) can be written in terms of the rotation matrix
(15) 
as , which is an orthogonal and symplectic matrix. Computing exactly (by diagonalizing the matrix ) is, as mentioned before for its complex representation , computationally very expensive, so that one typically splits the whole time interval into subintervals of length and then approximate acting on the initial condition at each step. Since
it makes sense to apply splitting methods of the form
(16) 
Observe that the evaluation of the exponentials of and requires only computing the products and , and this can be done very efficiently with the FFT algorithm.
Several methods with different orders have been constructed along these lines [13, 21, 27]. In particular, the schemes presented in [13] use only exponentials and to achieve order for and . Furthermore, when the idea of processing is taken into account, it is possible to design families of symplectic splitting methods with large stability intervals and a high degree of accuracy [2, 3]. They have the general structure , where (the kernel) is built as a composition (16) and (the processor) is taken as a polynomial.
Although these methods are neither unitary nor unconditionally stable, they are symplectic and conjugate to unitary schemes. In consequence, neither the average error in energy nor the norm of the solution increase with time. In other words, the quantities and are both approximately preserved along the evolution, since the committed error is (as shown in Subsection 3.1 below) only local and does not propagate with time. The mechanism that takes place here is analogous to the propagation of the error in energy for symplectic integrators in classical mechanics [15]. In addition, the families of splitting methods considered here are designed to have large stability intervals and can be applied when no particular structure is required for the Hamiltonian matrix . Furthermore, they can also be used in more general problems of the form , , resulting, in particular, from the space discretization of Maxwell equations [3].
3 Analysis of symplectic splitting methods for time discretization
In this section we proceed to characterize the family of splitting symplectic methods (16), paying special attention to their stability properties. By interpreting the numerical solution as the exact solution corresponding to a modified differential equation, it is possible to prove that the norm and energy of the original system are approximately preserved along evolution. We also provide rigorous estimates of the time discretization error that are uniformly valid as both the space and time discretizations get finer and finer. The analysis allows us to construct new methods with large stability domains such that the error introduced is comparable to the error coming from the space discretization.
3.1 Theoretical analysis
It is clear that the problem of finding appropriate compositions of the form (16) for equation (14) is equivalent to getting coefficients , in the matrix
(17) 
such that approximates the solution , where denotes the rotation matrix (15). The matrix that propagates the numerical solution of the splitting method (17) can be written as
(18) 
where the entries and (respectively, and ) are even (repect., odd) polynomials in , and . It is worth stressing here that by diagonalizing the matrix with an appropriate linear change of variables, one may transform the system into uncoupled harmonic oscillators with frequencies . Although in practice one wants to avoid diagonalizing , numerically solving system (13) by a splitting method is mathematically equivalent to applying the splitting method to each of such onedimensional harmonic oscillators (and then rewriting the result in the original variables). Clearly, the numerical solution of each individual harmonic oscillator is propagated by the matrix with polynomial entries () for . We will refer to in the sequel as the propagation matrix, although other denominations have also been used [3, 23].
Moreover, for a given with polynomial entries, an algorithm has been proposed to factorize as (17) and determine uniquely the coefficients , of the splitting method [3, Proposition 2.3]. Thus, any splitting method is uniquely determined by its propagation matrix . For this reason, in the analysis that follows we will be only concerned with such matrices .
When applying splitting methods to the system (13) with time step size , the numerical solution is propagated by as an approximation to , which is bounded (with norm equal to 1) independently of . It then makes sense requiring that be also bounded independently of . This clearly holds if for each eigenvalue of , the corresponding matrix with is stable, i.e., if for some constant .
In our analysis, use will be made of the stability polynomial, defined for a given by
(19) 
The following proposition, whose proof can be found in [3], provides a characterization of the stability of .
Proposition 2
Let be a matrix with , and , with . Then, the following statements are equivalent:

The matrix is stable.

The matrix is diagonalizable with eigenvalues of modulus one.
We define the stability threshold as the largest non negative real number such that is stable for all . Thus, will be bounded independently of if all the eigenvalues of lie on the stability interval , that is, if , where is the spectral radius of the matrix . For instance, if a Fouriercollocation approach based on nodes is applied to discretize (5) in space, the spectral radius is of size (cf. eqs. (7)(8)), which shows that must decrease proportionally to as the number of nodes increases.
The stability threshold depends on the coefficients of the method (16) and verifies , since is the optimal value for the stability threshold achieved by the concatenation of steps of length of the leapfrog scheme [7].
The stability of the matrix of a splitting method for a given can alternatively be characterized as follows.
Proposition 3
The matrix is stable for a given if and only if there exist real quantities , with , such that
(21) 
Proof. If is of the form (21) then, obviously, and thus . Moreover, it is straightforward to check that (20) holds with
(22) 
so that is stable in that case.
Let us assume now that is stable, so that from the third characterization given in Proposition 2, , where is the stability polynomial. We now consider two cases:

(resp. ), so that (resp. ) is similar to the identity matrix, which implies that (resp. ) is also the identity matrix. In that case, (21) holds with and .

If , then , and we set
Since , one has
which implies that and
Notice that, for a given splitting method with a nonempty stability interval , Proposition 3 determines two odd functions and and an even function defined for which characterize the accuracy of the method when applied with step size to a harmonic oscillator of frequency , with . An accurate approximation will be obtained if , , and are all small quantities. In particular, if the splitting method is of order , then
as . For instance, for the simple first order splitting one has
and one can easily check that
It is worth stressing that (21) implies that (20) holds with given by (22). This feature, in particular, allows us to interpret the numerical result obtained by a splitting method of the form (16) applied to (14) in terms of the exact solution corresponding to a modified differential equation. Specifically, assume that is obtained (for , ) as
If , then it holds that
trivially verifies . In other words, is the exact solution at of the initial value problem
(23) 
where . With this backward error analysis interpretation at hand, it readily follows the preservation of both the discrete norm of
and the energy corresponding to (23). This implies that the discrete norm of and the energy of the original system will be approximately preserved (that is, their variation will be uniformly bounded for all times ).
3.2 Error estimates
Our goal now is to obtain meaningful estimates of the time discretization error that are uniformly valid as and , that is, as both the space discretization and the time discretization get finer and finer. We know that, by stability requirements, the time step used in the time integration by a splitting method of system (13) must be chosen as , where the stability threshold must verify for an stage splitting method. Since the Hermitian matrix comes from the space discretization of an unbounded selfadjoint operator, the spectral radius will tend to infinity as , and thus inevitably . It seems then reasonable to introduce the parameter
(24) 
and analyze the timeintegration error corresponding to a fixed value of .
The fact that, for each , (20) holds with given by (22) implies that, for each ,
This will allow us to obtain rigorous estimates for the error of approximating by applying steps of a splitting method with timestep . Specifically, we have
(25)  
where denotes the rotation matrix (15) and the matrix is given by
(26) 
so that the following theorem can be stated.
Theorem 4
Given , let be the approximation to obtained by applying steps of length (24) of a splitting method with stability threshold . Then one has
(in the Euclidean norm), where
Proof. From (25), we can write
where, for clarity, . For the first contribution we have
since is Hermitian. Now
As for the second contribution, one has
whereas
and thus the proof is complete.
Notice that the error estimate in the previous theorem does not guarantee that, for a given , the error in approximating by applying steps of the method is bounded as . As a matter of fact, since must satisfy the stability restriction , so that , one has that (and hence the error bound above) goes to infinity as . This can be avoided by estimating the error in terms of in addition to . The assumption that can be bounded uniformly as the space discretization parameter , implies that the initial state of the continuous time dependent Schrödinger equation is such that is squareintegrable. The converse will also be true for reasonable space semidiscretizations and a sufficiently smooth potential .
More generally, the assumption that has sufficiently high spatial regularity (together with suitable conditions on the potential ) is related to the existence of bounds of the form that hold uniformly as . In this sense, it is useful to introduce the following notation:

Given , we denote for each

For a stage splitting method with stability threshold , given and we denote
(27) (28) Clearly, and are bounded if and only if the method is of order .
We are now ready to state the main result of this section.
Theorem 5
Given and , let be such that (with ), and let be the approximation to obtained by applying steps of length of a th order splitting method with stability threshold . Then, for each ,
(29) 
Proof. We proceed as in the proof of Theorem 4. First we bound
with
Then the second term in (25) verifies
and
from which (29) is readily obtained.
Some remarks are in order at this point:

Recall that the estimate in Theorem 1 shows the behavior of the space discretization error (of a spectral collocation method applied to the 1D Schrödinger equation) as the number of collocation points goes to infinity. Our estimate (29) shows in turn the behavior of the time discretization error as , provided that with a fixed . In that case, it can be shown that uniformly for all , and thus the error of the full discretization admits the estimate
Notice the similarity of both the space and time discretization errors ( is a discrete version of a continuous norm which is equivalent to the Sobolev norm ).

Given a splitting method with stability threshold of order for the harmonic oscillator, consider and in (27)(28) for a fixed . If instead of analyzing the behavior as increases of the time discretization error committed by the splitting method when applied with , one is interested in analyzing the error with fixed and decreasing , one proceeds as follows. Since by definition
then reasoning as in the proof of Theorem 5, one gets the estimate
From a practical point of view, (29) can be used to obtain a priori error estimates just by replacing the exact by an approximation (obtained for instance with some generalization of the power method), or by an estimation based on the knowledge of bounds of the potential and the eigenvalues of the discretized Laplacian.
The error estimates in Theorem 5 provide us appropriate criteria to construct splitting methods to be applied for the time integration of systems of the form (13) that result from the spatial semidiscretization of the time dependent Schrödinger equation. Such error estimates suggest in particular that different splitting methods should be used depending on the smoothness of initial state in the original equation. Also, Theorem 5 indicates that for sufficiently long time integrations, the actual error will be dominated by the phase errors, that is, the errors corresponding to .
4 On the construction of new symplectic splitting methods
Observe that when comparing the error estimates in Theorem 5 for a given corresponding to two methods with different number of stages and respectively, one should consider time steps and that are proportional to and respectively. In this way, the same computational effort is needed for both methods to obtain a numerical approximation of for a given . It makes sense, then, to consider a scaled time step of the application with time step of a stage splitting method to the system (13). This can be defined as
(30) 
so that the relevant error coefficients associated to the error estimates in Theorem 5 are and .
The task of constructing a splitting method in this family can be thus precisely formulated as follows.
Problem. Given a fixed number of stages in (17), and for prescribed values of and scaled time step , design some splitting method having order and stability threshold , which tries to optimize the main error coefficient while keeping reasonably small.
We have observed, however, that trying to construct such optimized methods in terms of the coefficients of the polynomial entries () of the propagation matrix leads us to very illconditioned systems of algebraic equations. That difficulty can be partly overcome by taking into account the following observations:

The functions and () determining the error estimates in Theorem 5 uniquely depend on two polynomials: the stability polynomial given in (19) and
(31) Indeed, from one hand, , so that uniquely depends on . On the other hand, according to Proposition 3,
and one can get by straigthforward algebra that is (in Euclidean norm)
(and thus as ).

Given an even polynomial and an odd polynomial , there exist a finite number of propagation matrices such that (19) and (31) hold. Indeed, the entries of such stability matrices are of the form
where and are respectively even and odd polynomials satisfying
(32) It is not difficult to see that there is a finite number of choices for such polynomials and