Error analysis of splitting methods for the time dependent Schrodinger equation

A typical procedure to integrate numerically the time dependent Schr\"o\-din\-ger 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 $i du/dt = H u$, where $H$ is a real symmetric matrix, whose solution with initial value $u(0) = u_0 \in \mathbb{C}^N$ is given by $u(t) = \e^{-i t H} u_0$. Usually, this exponential matrix is expensive to evaluate, so that time stepping methods to construct approximations to $u$ from time $t_n$ to $t_{n+1}$ are considered in the second phase of the procedure. Among them, schemes involving multiplications of the matrix $H$ 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 $Hu$. 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.


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 atom-diatom collisions and triatomic photo-dissociation 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, ψ : R d × R −→ C is the wave function representing the state of the system and the initial state is ψ(x, 0) = ψ 0 (x). UsuallyĤ =T (P ) +V (X) ≡ 1 2µP 2 +V (X) (2) and the operatorsX,P are defined by their actions on ψ(x, t) aŝ Xψ(x, t) = x ψ(x, t),P ψ(x, t) = −i ∇ψ(x, t).
The solution of (1) provides all dynamical information on the physical system at any time. It can be expressed as whereÛ represents the evolution operator, which is linear and satisfies the equation i dÛ (t)/dt =ĤÛ(t) withÛ (0) = I. Since the Hamiltonian is explicitly time independent, the evolution operator is given formally bŷ 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 ψ 0 (x) 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 idu/dt = Hu, with a Hermitian matrix H of large dimension and large norm. Since evaluating exactly the exponential exp(−itH) is computationally expensive, approximation methods requiring only matrix-vector products with H 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 exp(−itH) 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 H with vectors, they belong to the same class of integrators as the Chebyshev and Lanczos methods, in the sense that all of them approximate exp(−itH)u 0 by linear combinations of terms of the form H j u 0 (j ≥ 1).
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.

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 one-dimensional 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 [0, 2π], and therefore with ψ(0, t) = ψ(2π, t) for all t. In the Fourier-collocation (or pseudospectral) approach, one intends to construct approximations based on the equidistant interpolation grid where N is even (although the formalism can also be adapted to an odd number of points). Then one seeks a solution of the form [22] where the coefficients c n (t) are related to the grid values ψ N (x j , t) through a discrete Fourier transform of length N , F N [26]. Its computation can be accomplished by the Fast Fourier Transform (FFT) algorithm with O(N log N ) floating point operations.
In the collocation approach, the grid values ψ N (x j , t) are determined by requiring that the approximation (6) satisfies the Schrödinger equation precisely at the grid points x j [22]. This yields a system of N ordinary differential equations to determine the N point values ψ N (x j , t): where for n = −N/2, . . . , N/2 − 1 and j = 0, . . . , N − 1. Observe that the matrices on the right-hand side of (7) are Hermitian. An important qualitative feature of this space discretization procedure is that it replaces the original Hilbert space L 2 (0, 2π) defined by the quantum mechanical problem by a discrete one in which the action of operators are approximated by N × N (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 c n exhibit a rapid decay (in some cases, faster than algebraically in n −1 , uniformly in N ), so that typically the value of N 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 Suppose that the exact solution ψ(x, t) of (5) is such that, for some s ≥ 1, ∂ s+2 x ψ(·, t) ∈ L 2 (0, 2π) for every t ≥ 0. Then the error due to the approximation ψ N (x, t) defined by (6) in the collocation approach is bounded by where C depends only on s.
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 one-dimensional features, by taking tensor products of one-dimensional 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 where H 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 u(t), since the matrix H 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 but to compute the exponential of the N × N complex and full matrix −itH (typically also of large norm) by diagonalizing the matrix H can be prohibitively expensive for large values of N . In practice, thus, one turns to time stepping methods advancing the approximate solution from time t n to t n+1 = t n + ∆t, so that the aim is to construct an approximation u n+1 ≈ u(t n+1 ) = e −i∆t H u(t n ) as a map u n+1 = φ ∆t u n . 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 where V is a diagonal matrix associated withV and T is related to the kinetic energyT . It turns out that the solutions e −itT u 0 and e −itV u 0 of equations iu = T u and iu = V u, respectively, can be easily determined [22], so that one may consider compositions of the form where τ ≡ ∆t. In (12) the number of exponentials m (and therefore the number of coefficients {a i , b i } m i=1 ) has to be sufficiently large to solve all the equations required to achieve order r (the so called order conditions).
Splitting methods of this class have several structure-preserving properties. They are unitary, so that the norm of u is preserved along the integration, and time-reversible 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 H in eq. (7) is a real symmetric matrix, and thus e −itH is not only unitary, but also symplectic with canonical coordinates q = Re(u) and momenta p = Im(u). In consequence, equation (9) is equivalent to [13,14] q = Hp,ṗ = −Hq.
Alternatively, one may write d dt with the 2N × 2N matrices A and B given by The solution operator corresponding to (14) can be written in terms of the rotation matrix O(y) = cos(y) sin(y) − sin(y) cos(y) (15) as O(t H), which is an orthogonal and symplectic 2N × 2N matrix. Computing O(t H) exactly (by diagonalizing the matrix H) is, as mentioned before for its complex representation e −i t H , computationally very expensive, so that one typically splits the whole time interval into subintervals of length τ ≡ ∆t and then approximate O(τ H) acting on the initial condition at each step. Since it makes sense to apply splitting methods of the form u n+1 = e τ bmB e τ amA · · · e τ b 1 B e τ a 1 A u n .
Observe that the evaluation of the exponentials of A and B requires only computing the products Hp and Hq, 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 m = r exponentials e τ a i A and e τ b i B to achieve order r for r = 4, 6, 8, 10 and 12. 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 P (τ H)K(τ H)P −1 (τ H), where K (the kernel) is built as a composition (16) and P (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 u 2 = u Tū /N and u T Hū/(2N ) 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 H. Furthermore, they can also be used in more general problems of the formq = M 1 p,ṗ = −M 2 q, resulting, in particular, from the space discretization of Maxwell equations [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.

Theoretical analysis
It is clear that the problem of finding appropriate compositions of the form (16) for equation (14) is equivalent to getting coefficients a i , b i in the matrix (15). The matrix K(τ H) that propagates the numerical solution of the splitting method (17) can be written as where the entries K 1 (y) and K 4 (y) (respectively, K 2 (y) and K 3 (y)) are even (repect., odd) polynomials in y ∈ R, and det K(y) = K 1 (y)K 4 (y)−K 2 (y)K 3 (y) ≡ 1. It is worth stressing here that by diagonalizing the matrix H with an appropriate linear change of variables, one may transform the system into N uncoupled harmonic oscillators with frequencies ω 1 , . . . , ω N . Although in practice one wants to avoid diagonalizing H, numerically solving system (13) by a splitting method is mathematically equivalent to applying the splitting method to each of such one-dimensional harmonic oscillators (and then rewriting the result in the original variables). Clearly, the numerical solution of each individual harmonic oscillator is propagated by the 2 × 2 matrix K(y) with polynomial entries K j (y) (j = 1, 2, 3, 4) for y = τ ω j . We will refer to K(y) in the sequel as the propagation matrix, although other denominations have also been used [3,23]. Moreover, for a given K(y) with polynomial entries, an algorithm has been proposed to factorize K(y) as (17) and determine uniquely the coefficients a i , b i of the splitting method [3, Proposition 2.3]. Thus, any splitting method is uniquely determined by its propagation matrix K(y). For this reason, in the analysis that follows we will be only concerned with such matrices K(y).
When applying splitting methods to the system (13) with time step size τ , the numerical solution is propagated by (K(τ H)) n as an approximation to O(τ H) n = O(nτ H), which is bounded (with L 2 norm equal to 1) independently of n. It then makes sense requiring that (K(τ H)) n be also bounded independently of n ≥ 1. This clearly holds if for each eigenvalue ω j of H, the corresponding 2×2 matrix K(y) with y = τ ω j is stable, i.e., if (K(τ ω j )) n ≤ C for some constant C > 0.
In our analysis, use will be made of the stability polynomial, defined for a given K(y) by The following proposition, whose proof can be found in [3], provides a characterization of the stability of K(y).
Proposition 2 Let K(y) be a 2 × 2 matrix with det K(y) = 1, and p(y) = 1 2 tr K(y), with y ∈ R. Then, the following statements are equivalent: (a) The matrix K(y) is stable.
(b) The matrix K(y) is diagonalizable with eigenvalues of modulus one.
We define the stability threshold y * as the largest non negative real number such that K(y) is stable for all y ∈ (−y * , y * ). Thus, (K(τ H)) n will be bounded independently of n ≥ 1 if all the eigenvalues of τ H lie on the stability interval (−y * , y * ), that is, if τ ρ(H) < y * , where ρ(H) is the spectral radius of the matrix H. For instance, if a Fourier-collocation approach based on N nodes is applied to discretize (5) in space, the spectral radius is of size ρ(H) = O(N 2 ) (cf. eqs. (7)-(8)), which shows that τ must decrease proportionally to N −2 as the number of nodes N increases.
The stability threshold y * depends on the coefficients {a i , b i } of the method (16) and verifies y * ≤ 2m, since 2m is the optimal value for the stability threshold achieved by the concatenation of m steps of length τ /m of the leapfrog scheme [7].
The stability of the matrix K(y) of a splitting method for a given y ∈ R can alternatively be characterized as follows.

Proposition 3
The matrix K(y) is stable for a given y ∈ R if and only if there exist real quantities φ(y), ǫ(y), γ(y), with γ(y) = 0, such that Proof. If K(y) is of the form (21) then, obviously, trK(y) = 2 cos(φ(y)) and thus cos(φ(y)) = p(y). Moreover, it is straightforward to check that (20) holds with so that K(y) is stable in that case. Let us assume now that K(y) is stable, so that from the third characterization given in Proposition 2, φ(y) = arccos(p(y)) ∈ R, where p(y) is the stability polynomial. We now consider two cases: • p(y) = 1 (resp. p(y) = −1), so that K(y) (resp. −K(y)) is similar to the identity matrix, which implies that K(y) (resp. −K(y)) is also the identity matrix. In that case, (21) holds with ǫ(y) = 0 and γ(y) = 1.
Since det(K(y)) = 1, one has which implies that γ(y) = 0 and Notice that, for a given splitting method with a non-empty stability interval (−y * , y * ), Proposition 3 determines two odd functions φ(y) and ǫ(y) and an even function γ(y) defined for y ∈ (−y * , y * ) which characterize the accuracy of the method when applied with step size τ to a harmonic oscillator of frequency ω, with y = τ ω. An accurate approximation will be obtained if |φ(y)−y|, |γ(y)−1|, and |ǫ(y)| are all small quantities. In particular, if the splitting method is of order r, then as y → 0. For instance, for the simple first order splitting e τ A e τ B one has and one can easily check that It is worth stressing that (21) implies that (20) holds with Q(y) 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, If τ ρ(H) < y * , then it holds that In other words,ũ n is the exact solution at t n = nτ of the initial value problem With this backward error analysis interpretation at hand, it readily follows the preservation of both the discrete L 2 norm of and the energy corresponding to (23). This implies that the discrete L 2 norm of u = q + i p and the energy of the original system will be approximately preserved (that is, their variation will be uniformly bounded for all times t n ).

Error estimates
Our goal now is to obtain meaningful estimates of the time discretization error that are uniformly valid as N → ∞ and τ → 0, 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 τ < y * /ρ(H), where the stability threshold must verify y * ≤ 2m for an m-stage splitting method. Since the Hermitian matrix H comes from the space discretization of an unbounded selfadjoint operator, the spectral radius ρ(H) will tend to infinity as N → ∞, and thus inevitably τ → 0. It seems then reasonable to introduce the parameter and analyze the time-integration error corresponding to a fixed value of θ < y * . The fact that, for each y ∈ (−y * , y * ), (20) holds with Q(y) given by (22) implies that, for each n ≥ 1, This will allow us to obtain rigorous estimates for the error of approximating e −itH u 0 by applying n steps of a splitting method with time-step τ = t/n. Specifically, we have where O(y) denotes the rotation matrix (15) and the 2 × 2 matrix E(y) is given by so that the following theorem can be stated.
Theorem 4 Given u 0 = q 0 + ip 0 , let u n = q n + ip n be the approximation to u(nτ ) = e −i nτ H u 0 obtained by applying n steps of length (24) of a splitting method with stability threshold y * . Then one has Proof. From (25), we can write where, for clarity, φ ≡ φ(τ H). For the first contribution we have As for the second contribution, one has and thus the proof is complete. Notice that the error estimate in the previous theorem does not guarantee that, for a given t, the error in approximating e −i t H u 0 by applying n steps of the method is bounded as ρ(H) → ∞. As a matter of fact, since τ = t/n must satisfy the stability restriction θ = τ ρ(H) < y * , so that n > t ρ(H)/y * , one has that n (and hence the error bound above) goes to infinity as ρ(H) → ∞. This can be avoided by estimating the error in terms of Hu 0 in addition to u 0 .
The assumption that Hu 0 can be bounded uniformly as the space discretization parameter N → ∞, implies that the initial state ψ(x, 0) of the continuous time dependent Schrödinger equation is such that ∂ 2 x ψ(x, 0) is square-integrable. The converse will also be true for reasonable space semi-discretizations and a sufficiently smooth potential V (x).
More generally, the assumption that ψ(x, 0) has sufficiently high spatial regularity (together with suitable conditions on the potential V (x)) is related to the existence of bounds of the form H k u 0 ≤ C k that hold uniformly as ρ(H) → ∞. In this sense, it is useful to introduce the following notation: • Given k ≥ 0, we denote for each u ∈ C N u k := H k u .
• For a m-stage splitting method with stability threshold y * , given k ≥ 0 and θ ∈ [0, y * ) we denote Clearly, µ k (θ) and ν k (θ) are bounded if and only if the method is of order r ≥ k.
We are now ready to state the main result of this section.
Theorem 5 Given u 0 = q 0 + ip 0 and t ∈ R + , let n be such that τ = t/n = θ/ρ(H) (with θ < y * ), and let u n = q n + ip n be the approximation to u(t) = e −i t H u 0 obtained by applying n steps of length τ of a r-th order splitting method with stability threshold y * . Then, for each k ∈ [0, r], Proof. We proceed as in the proof of Theorem 4. First we bound Then the second term in (25) verifies from which (29) is readily obtained.
Some remarks are in order at this point: 1. 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 N of collocation points goes to infinity. Our estimate (29) shows in turn the behavior of the time discretization error as N → ∞, provided that τ = θ/ρ(H) with a fixed θ < y * . In that case, it can be shown that ρ(H) −1 ≤ L N −2 uniformly for all N , and thus the error of the full discretization admits the estimate Notice the similarity of both the space and time discretization errors ( u 0 k+1 = H k+1 u 0 is a discrete version of a continuous norm ||ψ(·, 0)|| k+1 which is equivalent to the Sobolev norm ∂ 2k+2 x ψ(·, 0) ).
2. Given a splitting method with stability threshold y * of order r for the harmonic oscillator, consider µ r (θ) and ν r (θ) in (27)-(28) for a fixed θ < y * . If instead of analyzing the behavior as N increases of the time discretization error committed by the splitting method when applied with τ = θ/ρ(H), one is interested in analyzing the error with fixed H and decreasing τ ≤ θ/ρ(H), 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 ρ(H) 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 semi-discretization 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 µ k (θ).

On the construction of new symplectic splitting methods
Observe that when comparing the error estimates in Theorem 5 for a given k ≥ 0 corresponding to two methods with different number of stages m and m ′ respectively, one should consider time steps τ and τ ′ that are proportional to m and m ′ respectively. In this way, the same computational effort is needed for both methods to obtain a numerical approximation of u(t) for a given t > 0. It makes sense, then, to consider a scaled time step of the application with time step τ of a m-stage splitting method to the system (13). This can be defined as so that the relevant error coefficients associated to the error estimates in Theorem 5 are µ k (θ ′ m) and ν k (θ ′ m).
The task of constructing a splitting method in this family can be thus precisely formulated as follows.
Problem. Given a fixed number m of stages in (17), and for prescribed values of k ≥ 0 and scaled time step θ ′ ∈ (0, 2), design some splitting method having order r ≥ k and stability threshold y * > θ ′ m, which tries to optimize the main error coefficient µ k (θ ′ m) while keeping ν k (θ ′ m) reasonably small.
• Given an even polynomial p(y) and an odd polynomial q(y), there exist a finite number of propagation matrices K(y) such that (19) and (31) hold. Indeed, the entries of such stability matrices are of the form K 1 (y) = p(y) + d(y), K 2 (y) = q(y) + e(y), where d(y) and e(y) are respectively even and odd polynomials satisfying It is not difficult to see that there is a finite number of choices for such polynomials d(y) and e(y). The ill-conditioning mentioned before seems to come mainly from the ill-conditioning of the problem of determining d(y) and e(y) from prescribed polynomials p(y) and q(y). Obviously, a necessary condition for the existence of such polynomials d(y) and e(y) with real coefficients is that p(y) 2 + q(y) 2 − 1 ≥ 0 for all y. It is also straightforward to see that, for a rth order method, d(y) = O(y r+1 ) and e(y) = O(y r+1 ) (as y → 0), and thus In addition, if the method has stability threshold y * > 0, then there exists 0 < y 1 < · · · < y l < y * such that φ(y j ) = jπ, and thus p(y j ) = (−1) j , p ′ (y j ) = 0, q(y j ) = 0, for j = 1, . . . , l.
3. Apply the algorithm given in [3] to each of the matrices K(y) obtained in the previous step. In this way we will get the vector of coefficients (a i , b i ) of all the splitting methods having a progagation matrix K(y) satisfying (19) and (31) for the pair of polynomials (p(y), q(y)) determined in the first step. Since Theorem 5 gives exactly the same error estimate (29) for all the splitting schemes obtained in that way, we choose (with the aim of reducing the effect of round-off errors) one that minimizes m j=1 (|a j | + |b j |).
This procedure has been applied to construct several splitting methods of different orders r, number of stages m and scaled time steps θ ′ . We collect in Table 1 the relevant parameters of some of them, whereas the actual coefficients a j , b j can be found at www.gicas.uji.es/Research/splitting1.html. As a matter of fact, all the methods have m + 1 pairs of coefficients a i , b i , but b m+1 = 0. In consequence, the last stage at a given step can be concatenated with the first one at the next step, so that the overall number of stages is m. This property is called FSAL (first-same-as-last) in the numerical analysis m r  literature. According to the previous comments, the new schemes are aimed at integrating equation (13) under very different conditions of regularity. The first three methods in Table 1 are designed to be applied with the same scaled time-step θ ′ = τ ρ(H)/m = 1, and thus the three of them have the same computational cost. The order r of the methods is increased by adding more stages, while keeping reasonably small error coefficients µ r (θ) and ν r (θ). This will be advantageous, according to Theorem 5, for sufficiently regular initial states. Alternatively, one may want to use the additional number of stages to reduce the computational cost while keeping the same order r = 6. This can be illustrated with the fourth method in Table 1, which has been optimized for scaled time-step θ ′ = 1.4, and thus is substantially cheaper than the first method (optimized for θ ′ = 1), while having smaller error coefficients µ 6 (θ) and ν 6 (θ).
We now turn our attention to the methods of order zero in Table 1, which according to Theorem 5, are the methods of choice for very low regularity conditions. Although they have comparatively smaller error coefficients than the methods of order r > 1, one should bear in mind that the error estimates (29) for r ≥ k > 1 decrease with ρ(H) −k as the spectral radius ρ(H) increases, while for methods of order r = 0 the same estimate holds independently of the size of ρ(H). Comparing the first three methods of order r = 0 and m = 30, we see that, not surprisingly, the accuracy of the methods can be improved by increasing the computational cost (by considering methods optimized for lower values of θ ′ ). This is analogous to increasing the accuracy of the application of a Chebyshev polynomial of degree m = 30 by decreasing the time-step size. Now, by comparing the first 30-stage method of order 0 with the method with m = 40 and order 0, we see that the accuracy can be increased also by increasing the number of stages from m = 30 to m = 40 while keeping the same computational cost (with θ ′ = 1). This is similar to increasing the accuracy of Chebyshev approximations, while keeping the same computational cost, by increasing the degree of the polynomial from m = 30 to m = 40.
Recall that, if ω 1 , . . . , ω N are the eigenvalues of H, numerically integrating (13) by a splitting method is mathematically equivalent to applying the splitting method to N uncoupled harmonic oscillators with frequencies ω j . Particularizing the proof of Theorem 5 to this case, it is quite straightforward to conclude that, when integrating the system with scaled time step θ ′ (that is, with τ = mθ ′ /ρ(H), where m is the number of stages of the scheme,) the relative error made in each oscillator can be bounded by and thus |y j | ≤ θ ′ . When such a system of harmonic oscillators originates from a continuous problem possessing a high degree of regularity, the highest frequency oscillators have much smaller amplitude and thus can be approximated less accurately than the lower frequency oscillators without compromising the overall precision. For lower regularity conditions, the overall precision will be more affected by the accuracy of the approximations corresponding to higher frequency oscillators.
With the aim of illustrating the relative error made in each harmonic oscillator, in Figure 1 we represent (in double logarithmic scale) |φ(my)/(my)−1| and E(my) (which are even functions of y) for y ∈ [0, θ ′ ] for some of the methods collected in Table 1, identified by appropriate labels indicating their respective number of stages, order and scaled time step (m, r, θ ′ ). Observe that both functions of y exhibit a similar behavior for each splitting method, although the values taken by the second one are several orders of magnitude smaller, since the methods are designed to minimize mainly the phase error coefficient µ k (θ). The order r of each of the methods is reflected in the slope of the curves as y approaches 0.
We also include in Figure 1 the 12th order 12-stage scheme presented in [13], which is perhaps the most efficient when applied to harmonic oscillators among those (non-processed) splitting methods currently found in the literature. We denote it by GM12. It has a relative stability threshold y * /12 = 0.2618, so that, strictly speaking, it should be used with θ ′ = τ ρ(H)/12 < 0.2618 to guarantee stability. However, it seems in practice that the method can be safely used with θ ′ = 0.932 (for a larger value of the scaled time step θ ′ , the method becomes very unstable), because ||p(y)| − 1| < 10 −6 provided that |y|/12 < 0.932183. The theoretical instability for θ ′ ∈ (0.2618, 0.932183) is only relevant after a very large number of steps, and reveals itself as resonance peaks (which are clearly visible in the graph of E(my) in Figure 1 for k = 2, 3) near the values θ ′ = kπ/12, k = 1, 2, 3.
We can see in Figure 1 that GM12 is less accurate than the 30-stage 24th order method for the whole frequency range, and thus the former will show a poorer performance than the later for any regularity conditions. If the amplitudes at higher frequencies decrease fast enough, the 24th order method will   Table 1 and the 12th-order scheme GM12. Each new splitting method is identified by the triad (m, r, θ ′ ), indicating its number of stages m, order r and scaled time step θ ′ , as in Table 1. The resonances associated with the instability of GM12 at y ≈ kπ (k = 2, 3) are visible, especially in the second graph.
give very accurate approximations of u(t) = e −itH u 0 at a relatively low cost (since θ ′ = 1). The 6th order method with m = 30 stages gives correct approximations for all harmonic oscillators within the range y ∈ [−1.4, 1.4], and hence it is expected to give excellent approximations under mild regularity with a comparatively lower cost than the 24th order method (θ ′ = 1.4 compared with θ ′ = 1). Clearly, the methods of order 0 will be the right choices for low regularity conditions, since the corresponding phase errors |φ(my)/(my) − 1| are uniformly bounded for all y ∈ [−1, 1]. Among them, the method with m = 30 and θ ′ = 1 can be accurate enough in many practical computations. If more precision is required, one can either consider the method with θ ′ = 0.75 (with result in a 25% increase of the computational cost), or use the method with m = 40 and θ ′ = 1, without any increase of the computational cost (at the expense of having less frequent output).

Numerical examples
The purpose of this section is twofold. On the one hand, since the symplectic splitting methods we have presented here to approximate e −itH u 0 involve only products of the matrix H with vectors, it makes sense to compare them with other well established schemes of this kind, such as the Chebyshev and Lanczos methods. Although a thorough comparison with the family of splitting methods proposed in this work will be the subject of a forthcoming paper [4], we include here some results which show that the new schemes are indeed competitive for evaluating exp(−itH)u 0 , at least in the example considered.
On the other hand, since Theorem 5 provides a rigorous a priori estimate on the error committed when using a splitting method of the form (17) in the time integration of equation (9), it is interesting to check how this theoretical error estimate behave in practice for some of the methods constructed here.

A preliminary comparison with Chebyshev and Lanczos methods
As is well known, Chebyshev and Lanczos methods provide high order polynomial approximations to e −itH u 0 requiring only matrix-vector products. The former is neither unitary nor symplectic, whereas the later is unitary, but symplectic only in the Krylov subspace (which changes from one time step to the next).
To carry out this comparison we choose the very simple example previously considered in [22]. The problem consists in approximating y = e −iA v, where v is a random vector of unit norm and A is the tridiagonal matrix A = ω 2 tridiag(−1, 2, −1) of dimension 10000. The eigenvalues of A are contained in the interval [0, 2ω].
We have implemented the Chebyshev and Lanczos algorithms in the usual way (see [22]) with the particularity that, since the range of values for the eigenvalues is known, both the Chebyshev and the new splitting methods are used with a shift to the midpoint of the spectrum. In other words, y = e −iωI e −i(A−ωI) v, with I the identity matrix. This shift allows us to take ρ(A − ωI) ≃ ω. Figure 2 shows the error, y −y ap for different degrees m of the polynomials used and for ω = 15, 20, 30, 40. Here y is computed numerically to high accuracy and y ap corresponds to the approximate solution obtained by each scheme. Each particular value of m in the Lanczos and Chebyshev methods corresponds to a different mth-order polynomial approximation (denoted by small crosses and circles, respectively). We clearly observe that, for this irregular problem, the Lanczos method converges to the optimal Chebyshev method, the main difference between both schemes being the number of vectors to be kept in memory.
To apply the new splitting methods, we notice that for this problem the time step τ = 1, and the corresponding spectral radius ρ ≃ ω. Therefore we shall consider splitting methods whose value of θ ′ given by (30) satisfies In other words, for each ω the method (m, r, θ ′ ) is such that m θ ′ ≥ ω.
From the graphs of Figure 2, it is clear that, for each value of ω, we can always select one particular splitting method (big dots) which outperforms both Chebyshev and Lanczos. High-order splitting methods show a worst performance than schemes of order zero for this problem, since they require typically a higher degree of regularity.

The Pöschl-Teller potential
We next illustrate the error estimate provided by Theorem 5 for the class of splitting methods proposed here. We also compare the error in the time integration with the error coming from the space discretization for different values of the mesh size N . For that purpose we choose a well known anharmonic quantum potential leading to analytical solutions and consider, for clarity, only the 30-stage splitting methods of order six and order zero collected in Table 1. Specifically, we consider the Pöschl-Teller potential with λ > 1. It has been frequently used in polyatomic molecular simulation and is also of interest in supersymmetry, group symmetry, the study of solitons, etc. [8,10,20]. The parameter λ gives the depth of the well, whereas α is related to the range of the potential. The energies are We take the following values for the parameters (in atomic units, a.u.): the reduced mass µ = 1745 a.u., α = 2, λ = 24.5 (leading to 24 bounded states), x ∈ [−5, 5], and assume the system is periodic. The periodic potential is continuous and very close to differentiable. For N = 128 we have ρ(H) ≃ 0.635, whereas for N = 256 one gets ρ(H) ≃ 1.85.
We take as initial condition the Gaussian function, ψ(x, 0) = σ e −b 2 x 2 , where σ is a normalizing constant. With b = 3 the function and all its derivatives of practical interest vanish up to round off accuracy at the boundaries. The initial conditions contain part of the continuous spectrum, but this fact does not cause any trouble due to the smoothness of the periodic potential and wave function. As an illustration, some of the corresponding values of u 0 k for N = 128 are: u 0 1 = 0.629909, u 0 6 = 0.0722513, u 0 7 = 0.0478388. They decrease only moderately for the first values of k (before they increase again due to the contributions coming from higher energies). The corresponding values for N = 256 are quite similar to the previous ones. For this problem, both large and very small spatial errors are expected from spectral methods, depending on the mesh employed. It is then useful to have different methods with large values of θ ′ when low accuracy is desired and smaller values of θ ′ for high accuracy.
We integrate for t ∈ [0, 128 T ] with T = 333 and measure the 2-norm error in the discrete wave function, u ex (2 i T ) − u ap (2 i T ) , for i = 0, 1, . . . , 7. The values u ex are computed using the same spatial discretization and an accurate time integration (using a very small time step), whereas u ap stand for the numerical approximations obtained with splitting methods. For a given spatial discretization, we choose the time step for each of the new methods such that τ ≤ m θ ′ /ρ(H). In consequence, a period T has to be divided into M steps such that M = T /τ ≥ T ρ(H)/(mθ ′ ). In particular, for the 6th-order method (30, 6, 1.4) and N = 256, since ρ(H) ≃ 1.85, we take M = 15 ≥ (333 × 1.85)/(30 × 1.4) (each period T requires 15 steps, 450 stages or 1800 FFTs calls).
The results obtained are shown in Figure 3. Solid lines represent the error with respect to the exact solution for the same spatial discretization, whereas dashed lines correspond to the total error with respect to the exact solution obtained with a finer mesh (it is the sum of the spatial error and the error from the time integration). Dotted lines are obtained with the estimate (29).
We observe that the spatial error decreases exponentially with N due to the smoothness and periodicity of the problem. To estimate this spatial error we take the results obtained with N = 512 and an accurate time integration as the exact solution, and compare with the solution computed up to a high accuracy for N = 128 and N = 256. For N = 128 the spatial error dominates the total error, so that the most convenient time integration scheme is one able to provide such accuracy with a large time step. These requirements are fulfilled by the (30, 6, 1.4) method, especially designed to be used with θ ′ = 1.4, whereas scheme (30, 0, 1) gives us higher accuracy than necessary and with more computational cost. Method (30, 6, 1.4) can be used with a time step τ about a 40% larger than method (30, 0, 1), and thus its computational cost is reduced approximately by this factor .
However, for N = 256 the spatial error reaches nearly round off accuracy, and it could be convenient to employ methods able to provide this accuracy with the minimal computational cost. Notice that in this case the error committed by the 30-stage method with θ ′ = 1, (30, 0, 1), is still larger than the spatial error. In consequence, it makes sense integrating in time with a method specially designed to be used with a smaller time step. Thus, in particular, we reach round off accuracy with the 30-stage method (30, 0, 0.75) (θ ′ = 0.75) which is nearly twice more expensive than the method with θ ′ = 1.4.
We have also performed here the time integration with the 12th-order scheme GM12, which in the case of N = 128 requires a scaled time step of θ ′ = 0.49 to give a precision similar to that obtained by (30, 6, 1.4) with θ ′ = 1.4. With N = 256, GM12 must be applied with θ ′ = 0.19 to achieve the precision obtained by (30, 0, 0.75) with θ ′ = 0.75, thus requiring approximately four times more FFT calls.

Concluding remarks
The time integration of the Schrödinger equation previously discretized in space has been extensively studied in the literature. This is essentially equivalent to approximate u(t) = e −itH u(0), where H is a real symmetric matrix and u(0) represents the discrete wave function. In this work we propose using symplectic splitting integration methods to get this approximation. The main difference with standard polynomial approximations is that in the products Hv = H Re(v) + i H Im(v), the real and imaginary parts are computed sequentially instead of simultaneously (i.e., the computation of the real part is used in the computation of the imaginary part and vice versa in consecutive stages). These schemes are conjugate to unitary methods, so that the errors in norm and energy do not grow secularly [3].
To carry out the integration, one divides the whole time interval into n steps of length τ = t/n and applies an m-stage method at each time step. The total computational cost of the method is measured by the product n m instead of m. The analysis carried out in this paper allows us, in particular, to construct a particular symplectic splitting scheme of the form (16) which minimizes the total cost n m, given a a prescribed tolerance, the spectral radius ρ(H) of the corresponding Hamiltonian matrix H and the norm of its action on the initial condition, H k u(0) . We have observed that the optimal methods in this sense have relatively large values of m. We can choose the most appropriate method for each problem, i.e. the method, (m, r, θ ′ ), with the largest value of θ ′ which provides the desired accuracy for a given problem.
The error analysis of splitting methods provided here allows one to get a priori bounds on the propagating error when numerically integrating with a given time step which are comparable to similar estimates for the space discretization error. Moreover, it permits to construct new classes of schemes with a large stability interval specifically designed to be used with a certain (large) time step in such a way that the accuracy is similar to the spatial discretization error for a given space regularity. The main ingredients in the process are again the values of ρ(H), H k u 0 , and the estimate provided by Theorem 5. The numerical examples considered illustrate the validity of our approach. In particular, they show that there are methods in this family which are competitive with other standard procedures, such as Chebyshev and Lanczos methods.
By following this procedure it is indeed possible to generate a list of integration schemes specifically designed to be used under different regularity conditions on the initial state and the Hamiltonian matrix which involve in each case an error comparable to that coming from the spatial discretization. It is our purpose in a forthcoming paper [4] to elaborate an algorithm in such a way that, given a prescribed tolerance, an initial state u 0 and a Hamiltonian matrix H, automatically selects the most efficient time integration method in this family fulfilling the requirements supplied by the user. Moreover, we will also carry out a detailed numerical study of this family of splitting methods and the proposed automatic algorithm in comparison with the Chebyshev polynomial expansion scheme and the Lanczos iteration method.