Explicit schemes for time propagating many-body wavefunctions

Accurate theoretical data on many time-dependent processes in atomic and molecular physics and in chemistry require the direct numerical solution of the time-dependent Schr\"odinger equation, thereby motivating the development of very efficient time propagators. These usually involve the solution of very large systems of first order differential equations that are characterized by a high degree of stiffness. We analyze and compare the performance of the explicit one-step algorithms of Fatunla and Arnoldi. Both algorithms have exactly the same stability function, therefore sharing the same stability properties that turn out to be optimum. Their respective accuracy however differs significantly and depends on the physical situation involved. In order to test this accuracy, we use a predictor-corrector scheme in which the predictor is either Fatunla's or Arnoldi's algorithm and the corrector, a fully implicit four-stage Radau IIA method of order 7. We consider two physical processes. The first one is the ionization of an atomic system by a short and intense electromagnetic pulse; the atomic systems include a one-dimensional Gaussian model potential as well as atomic hydrogen and helium, both in full dimensionality. The second process is the decoherence of two-electron quantum states when a time independent perturbation is applied to a planar two-electron quantum dot where both electrons are confined in an anharmonic potential. Even though the Hamiltonian of this system is time independent the corresponding differential equation shows a striking stiffness. For the one-dimensional Gaussian potential we discuss in detail the possibility of monitoring the time step for both explicit algorithms. In the other physical situations that are much more demanding in term of computations, we show that the accuracy of both algorithms depends strongly on the degree of stiffness of the problem.

Accurate theoretical data on many time-dependent processes in atomic and molecular physics and in chemistry require the direct numerical ab initio solution of the time-dependent Schrödinger equation, thereby motivating the development of very efficient time propagators. These usually involve the solution of very large systems of first order differential equations that are characterized by a high degree of stiffness. In this contribution, we analyze and compare the performance of the explicit one-step algorithms of Fatunla and Arnoldi. Both algorithms have exactly the same stability function, therefore sharing the same stability properties that turn out to be optimum. Their respective accuracy however differs significantly and depends on the physical situation involved. In order to test this accuracy, we use a predictor-corrector scheme in which the predictor is either Fatunla's or Arnoldi's algorithm and the corrector, a fully implicit four-stage Radau IIA method of order 7. In this contribution, we consider two physical processes. The first one is the ionization of an atomic system by a short and intense electromagnetic pulse; the atomic systems include a one-dimensional Gaussian model potential as well as atomic hydrogen and helium, both in full dimensionality. The second process is the decoherence of two-electron quantum states when a time independent perturbation is applied to a planar two-electron quantum dot where both electrons are confined in an anharmonic potential. Even though the Hamiltonian of this system is time independent the corresponding differential equation shows a striking stiffness which makes the time integration extremely difficult. In the case of the one-dimensional Gaussian potential we discuss in detail the possibility of monitoring the time step for both explicit algorithms. In the other physical situations that are much more demanding in term of computations, we show that the accuracy of both algorithms depends strongly on the degree of stiffness of the problem.

I. INTRODUCTION
The numerical integration of the time-dependent Schrödinger equation (TDSE) has become the main theoretical approach for the quantitative study of a vast amount of phenomena, including strong field processes in atoms and molecules, quantum collisions and chemical reactions. In strong field physics, current light sources can create ultrashort pulses of very high intensity making the numerical solution of the TDSE unavoidable if accurate results are required. In the low frequency regime where the photon energy is much lower than the ionization potential, the advent of high-intensity lasers has allowed detailed investigations of phenomena such as above-threshold ionization [1], high order harmonic generation [2], multiphoton multiple ionization [3], attosecond pulse generation [4], molecular self-spectroscopy [5], etc. In the high frequency regime where the photon energy is of the order or larger than the ionization potential, very intense coherent X-ray sources are under development. They are based on the collective electronic response of a plasma to ultra intense laser fields [6] as well as the next generation free electron lasers (FEL) such as the European XFEL project. The latter is expected to boost the average photon flux by about two orders of magnitude in comparison with already existing FELs. The interaction of atoms or molecules with intense X-ray pulses with a duration in the femtosecond or subfemtosecond regime is expected to lead to highly non-linear processes which can no longer be described within perturbation theory as is currently the case.
In quantum collision theory, the interaction Hamiltonians do not usually depend explicitly on time. Time independent approaches such as the R-matrix [7] or the S-matrix methods [8] suffice. However, in many cases, the lack of knowledge of the asymptotic boundary conditions or the explicit introduction of a time in the interaction Hamiltonian through the classical description of a heavy projectile makes the numerical solution of the corresponding TDSE more convenient [9]. Methods such as the time-dependent close coupling [10] are particularly efficient. Nevertheless, when the quantum systems involved become more complex as is often the case for chemical reactions or in condensed matter physics, the numerical solution of the TDSE is no longer possible. Different approaches are necessary as for instance the time-dependent density functional theory (TDDFT) [11,12]. In that case, the Hamiltonian is replaced by the self-consistent Kohn-Sham Hamiltonian. In fact, TDDFT can be viewed as a reformulation of time-dependent quantum mechanics where the basic unknown is no longer the many-body wavefunction, but the time-dependent electron density [13]. This density can be obtained from the solution of a set of one-body equations, the so-called Kohn-Sham equations that have the same form as the usual TDSE.
The necessity of integrating numerically the TDSE motivates the development of efficient and accurate time propagators. Once the TDSE has been discretized in the spatial or/and energy domain by means of either a finite difference grid method or an approach based on spectral or finite element methods, the time integration of the TDSE reduces to the solution of a system of first order differential equations which may be written as follows: where H(t) is a matrix that depends explicitly on time. The main difficulty we have to face in solving such a system of ordinary differential equations is the fact that the spectrum of matrix H is not bound. In general, the matrix −iH has very high, purely imaginary, eigenvalues. These very high eigenvalues give rise to extremely fast oscillations of the true solution and usually determine the time step of the numerical time propagator. Another way to describe this problem is to say that the system behaves as a stiff system. Although there is no rigorous mathematical definition of the stiffness, a system is said to be stiff in a given interval of integration if the numerical method is forced to use a step length which is excessively small in relation to the smoothness of the exact solution [14]. In addition, increasing the size of the system generates eigenvalues each time higher thereby increasing its stiffness.
The problem associated with stiff systems is twofold: stability and accuracy. To each numerical method is associated a function named the stability function that determines the stability properties of the method and the range of time steps for which the numerical solution is stable and remains bounded. In the case of stiff systems, it can be shown, for instance, that none of the explicit methods of Runge Kutta (R-K) type is stable. In that case, it is necessary to use an implicit R-K scheme. Note that, by contrast to an explicit method which only requires matrix-vector products, all implicit schemes require solving systems of algebraic equations at each time step. For systems of considerable size, the computer time becomes, in these conditions, rapidly excessive. Fortunately, explicit schemes exist that are not of R-K type but having the stability properties required for dealing with such stiffness problems. The accuracy problem is more delicate. If an appropriate integration method is used, the stability problem may be avoided but, for a reasonable step length, the solution components corresponding to the largest eigenvalues are approximated very inaccurately [15]. However, there is no mathematical tool which allows one to predict whether the numerical solution of a stiff system will be accurate or not. Very often, the highest eigenvalues that correspond to very high energies do not play any physical role but, this does not imply that the error made in calculating the corresponding high energy components of the full numerical solution will not affect the final result. In fact, it is important to proceed on a case by case basis.
In this contribution, we analyze in detail two explicit one-step integration schemes that have the required stability properties for dealing with stiff systems. The first method is due to Fatunla [16,17] and the second one is a Krylov subspace method usually called the Arnoldi algorithm. The Arnoldi algorithm has already been used in many different contexts: strong field physics [18], condensed matter physics [13], etc... However, as far as we know, no systematic study of its stability and accuracy properties exists so far. In order to test the accuracy of both methods, we use a predictor-corrector scheme in which the predictor is either Fatunla's method or Arnoldi's algorithm while the corrector is a four-stage diagonally implicit Radau IIA method of order seven. Here, we consider the interaction of a quantum system with a strong and ultrashort electromagnetic pulse and test the three methods in the case of three different quantum systems: a model potential, atomic hydrogen and helium. We also examine the performances of these explicit schemes in a completely different context namely the calculation of a fidelity function that measures the decoherence of two-electron quantum states when a time independent perturbation is applied to a planar two-electron quantum dot where both electrons are confined in an anharmonic potential. In fact this is a difficult problem, which exhibits a strong degree of stiffness although its Hamiltonian is time independent. Finally, let us mention that a comparison of different time propagation algorithms for the time dependent Schrödinger equation may be found in [19].
This article is organized as follows. Section II is devoted to the general formulation of the TDSE. After some preliminary remarks, we give and discuss the general spectral representation of the TDSE and finally define the stability function of a given algorithm. In section III, we introduce the various algorithms (Fatunla's method, Arnoldi's algorithm and the predictor-corrector scheme) in the context of our model potential. For both explicit schemes, we give their stability function and analyze in detail their accuracy in the case of the model potential. Section IV is devoted to the results obtained with Fatunla's and Arnoldi's methods for the model potential, the interaction of atomic hydrogen with both a high and a very low frequency strong laser field and single ionization of helium. Finally, we consider the problem of the planar two-electron quantum dot. Unless otherwise stated, atomic units are used throughout this paper.

A. Preliminary remarks
Our aim is to study the interaction of a quantum system with an external time-dependent field. Solving numerically the corresponding TDSE proceeds in two steps: the discretization in the spatial or/and energy domain of the equation and the time propagation of the solution. There are typically three ways of discretizing the TDSE: the finite difference grid (FDG) methods and the approaches based on spectral or finite element methods. The simplest approaches are the FDG methods. These methods based on a spatial discretization are essentially local. They are very often used because the subsequent time propagation involves solving very sparse systems of algebraic equations. However, it is often tricky to extract information on how these methods account for the electronic structure of the quantum system under consideration and some observables are sometimes difficult to calculate. Furthermore, these methods yield finite-order rates of convergence in terms of the number of spatial grid points. In other words, the errors go to zero as the inverse of this number at a power given by the order of the method.
The spectral methods based on an energy discretization are non-local. They consist in writing the solution as a truncated expansion in terms of L 2 integrable functions. These functions form a complete basis set. Different choices of basis sets are possible, which usually depend on the physics of the problem. The most commonly used functions are the Hermite functions, the Coulomb Sturmian functions [20] and orthogonal polynomials. There are essentially two types of spectral methods: Galerkin and collocation [21]. In addition to the energy discretization as is the case in the Galerkin method, the collocation method involves also a spatial discretization. However, by contrast to the FDG methods, the grid mesh points are not arbitrarily chosen. They are the abscissae of the Gaussian quadrature associated with the basis functions. The spectral approaches are very appropriate for a very accurate description of the bound and resonance states of the quantum system under consideration. This is particularly true for resonance states very close to the ionization thresholds [22]. The convergence of the spectral methods in terms of the number of basis functions depends on the analytical properties of the solution. If the successive spatial derivatives of this solution do not exhibit singularities, the convergence is exponential. This means that the errors go to zero faster than any finite power of the inverse of the number of basis functions. On the other hand, if the solution exhibits singularities in its successive derivatives and if the basis wavefunctions do not account for these singularities, the convergence is much slower. Typical examples of such singularities are the Kato cusps present in many-particle system wavefunctions [23]. Another drawback of the spectral methods is the fact that the matrix associated with the Hamiltonian is, in most cases, not sparse.
The finite element methods which are based on a subdivision of the whole spatial domain of integration into simple subdomains are in fact closely related to the spectral methods. They differ, however, by the fact that the basis functions have bounded support, being therefore piecewise regular. In addition, these methods yield also finite-order rates of convergence like in the case of the FDG methods. Piecewise Lagrange polynomials or B-splines are very often used as basis functions. In general, these methods are particularly efficient in describing the electronic continuum states of the system under consideration. In addition, singularities or large gradients in the solution can be treated by considering non-regular subdomains. These methods are very often used, especially those based on B-splines [24] because the subsequent time propagation involves relatively sparse systems of equations to solve like in the case of the FDG methods. In the present contribution, we use spectral or/and B-spline based methods in all the cases treated.

B. The spectral representation of the TDSE
The TDSE for a quantum system interacting with an external field can be written as where Ψ (r, t) is the wavefunction of the system, r represents any set of n spatial coordinates and t is the time. The total Hamiltonian H (r, t), which depends explicitly on time, is given by with H 0 (r) the unperturbed Hamiltonian and V (r, t) the time-dependent interaction potential (velocity form in all the cases treated here). Using a complete basis set {f i (r)} of square integrable functions, we write the wavefunction Ψ(r, t), the solution of equation (2), as the following truncated expansion, where the expansion coefficients ψ i (t) are time-dependent. N represents the number of terms in the expansion and is taken sufficiently large to represent the wavefunction to the desired accuracy. As a result, the TDSE is transformed into a matrix equation for the vector Ψ(t) = {ψ i (t)} N , given by For a non-orthonormal basis, the overlap matrix B and the Hamiltonian H(t) have elements defined by The time evolution of the wavepacket is then given by the solution of the following N -dimensional system of first order differential equations: There is actually no need to evaluate explicitly the inverse of the overlap matrix B. This matrix is always symmetric and positive definite, which allows a numerically stable and fast Cholesky decomposition. In that case the action of B −1 on a vector can be calculated straightforwardly by solving a very sparse system of algebraic equations. The vector Ψ(t) is said to be B-orthogonal and its norm is given by Note that in the case of the FDG methods, a system of equations similar to the system (5) has to be solved but the matrix H is no longer associated with the Hamiltonian.

C. The boundary and asymptotic conditions
In solving numerically the TDSE, the discretization method has to account correctly for the non-trivial problems of the boundary and asymptotic conditions. By way of illustration, let us consider the case of the ionization of atomic hydrogen by an intense low frequency laser field. The amplitude of the electron quiver motion determines the minimum spatial grid size or the minimum number of basis functions to be included. For high intensities and very low frequencies, this amplitude may become of the order of thousands of atomic units thereby requiring excessively long computational times. In addition, during the interaction process, ionization takes place and fast emitted electrons will rapidly reach the boundaries of the computational domain. It is therefore important to choose appropriate boundary conditions to avoid spurious reflections of the wavefunction at these boundaries. Such reflections can be avoided by further increasing the size of the computational domain, but this becomes rapidly untractable. Instead, reflection problems can be overcome by introducing complex absorbing potentials [25,26]. Those potentials however are usually not completely reflection free. A better approach is exterior complex scaling (ECS) in which the outgoing electron coordinate becomes complex beyond a certain distance from the nucleus which is larger than the amplitude of the quiver motion [27,28].
For single electron systems, the extraction of the information on the differential probability densities does not cause any problem since the asymptotic behavior of the field-free continuum states is known. This contrasts with the multi-electron systems where the asymptotic behavior of the multiple continuum wavefunctions is unknown. In that case, one can either develop approximate expressions for these continuum wavefunctions or use more sophisticated time-dependent methods that circumvent the problem. When the outgoing electrons are sufficiently far from each other so that their interaction becomes negligible, multiple continuum wavefunctions are usually approximated by a product of Coulomb functions [29][30][31][32]. The validity of this approximation which gives reliable results is discussed in [33]. More sophisticated methods that avoid any projection of the final wavepacket on approximated multiple continuum wavefunctions have been developed. Palacios et al. [34] have derived a time-dependent method where the extraction of the information from the wavepacket is based on ECS. Malegat et al. extract the information from the total wavepacket after propagating semiclassically its Fourier components in space over very large distances [33,35]. Scrinzi has extended the time-dependent surface flux method to single and double ionization of two-electron systems [36]. Hutchinson et al. [37] are developing a time-dependent R-matrix approach that can describe the interaction of any (light) atomic systems with short electromagnetic pulses. More recently, Hamido et al. have developed the so-called time scaled coordinate (TSC) method [38]. This latter method which is used in some of the cases treated in this contribution, consists in performing a time-dependent scaling of the radial coordinates of the electrons together with a phase transformation of the wavefunction. As a result, an harmonic potential appears in the scaled Hamiltonian, which confines the wavefunction in configuration space. It can be shown that a relatively long time after the interaction, the wavefunction becomes stationary and its modulus gives directly the momentum distribution of the particles resulting from the fragmentation of the system. Consequently this method clearly circumvents the above mentioned difficulties. It however introduces different length scales that need to be treated with multiresolution techniques and that influence the stability of the numerical time propagation scheme.

D. The stability function
In order to analyze the stability of a one-step numerical time propagation scheme, it is convenient to consider the following standard test problem (Dahlquist's equation): where λ is a constant. If we assume that y(0) = η, the solution of this equation is y(t) = η exp(λt). Usually, a system of equations is said to be stiff when its Jacobian matrix has some eigenvalues with a very large negative real part.
In the case of Eq.(10), assuming that the real part of λ is very large and negative leads to a solution that tends extremely rapidly to zero. We have therefore to look for the conditions that have to be imposed on the numerical time propagation scheme in order that the numerical solution y n = y(nδt) → 0 as n → ∞ where δt is the time step. By applying the one-step numerical time propagation scheme to Eq.(10), we obtain where R(z) is the so-called stability function. In order that y n tends to zero as n → ∞, we must impose R(λδt) < 1 thereby implying some constraints on the time step δt. The set S = {z = λδt ∈ C; |R(z)| ≤ 1} is called the stability domain of the numerical scheme. This latter one is said to be A-stable if its stability domain is included in C − = {z; Re z ≤ 0}. It is L-stable if, apart from being A-stable, the stability function has the property lim Re(λδt)→−∞ |R(λδt)| = 0. L-stable methods are the most stable ones [14].
In the present case, the Jacobian of the system of equations we are interested in has large purely imaginary eigenvalues. Although such systems behave like a stiff system, the analysis of the stability of the numerical scheme is more delicate. Suppose for instance that the numerical scheme we use is L-stable and that its stability domain covers the half-plane C − as well as large parts of the right half-plane C + . In these conditions, uninteresting high oscillations of the true solution may be damped by the numerical scheme. However, the norm of the solution will not be necessarily preserved since |R(λδt)| ≤ 1. We must impose, as an additional constraint, that |R(λδt)| = 1. This means that if λ is purely imaginary, R(λδt) = exp(λδt). Following Fatunla [17], a numerical time propagation scheme is said to be exponentially fitted at a complex value λ = λ 0 if the stability function R(λδt) satisfies the relation

III. TIME PROPAGATION ALGORITHMS
In this section we describe and compare the performance of two explicit one-step time propagation schemes, namely Fatunla's method and Arnoldi's algorithm in terms of stability and accuracy. To test the accuracy of both schemes we use an implicit predictor-corrector (P-C) method. The predictor is either Fatunla's method or Arnoldi's one and the corrector is a four-stage Radau II-A implicit method which is of Runge-Kutta type. Monitoring the time step during the time propagation using both explicit schemes is a key point which will be addressed first within a simple one-dimensional model of one electron in a Gaussian potential of the form V (x) = −V 0 e −βx 2 , where V 0 and β are constants and exposed to a cosine square pulse.
In the following sections, these two algorithms will be also tested in two more demanding physical situations. The first situation is the interaction of a quantum system with a strong electromagnetic pulse. The quantum systems we shall be studying in that case are atomic hydrogen and helium, both treated in full dimensions. The second physical situation is the time evolution of a two-electron wavepacket in a two-dimensional quantum dot.

A. Fatunla's method
The idea behind Fatunla's method is to take into account the intrinsic frequencies of the atom-field system by introducing interpolating oscillatory functions that approximate the solution of the TDSE. This allows one to deal with problems displaying eigenvalues that differ by many orders of magnitude. That explains why Fatunla's method has the capability to solve stiff equations, while requiring only matrix vector products. More precisely, we write the first order differential equation (8) as The solution Ψ(t) over a subinterval [t n , t n + δt = t n+1 ] is approximated by the interpolating oscillatory function with I being the identity matrix. Ω 1 and Ω 2 are diagonal matrices, usually called the stiffness matrices, and a, b, c are constant vectors. By demanding that the interpolating function (14) coincides with the theoretical solution at the endpoints of the interval [t n , t n+1 ], and that it satisfies the differential equation at t = t n , we arrive at the recursion formula, where we use the notation f n = f (t n , ψ n ), f t=tn . R and S represent diagonal matrices defined as Φ and Ξ are diagonal matrices with non-zero entries given by [16,17], and The recursion formula (15) depends on the so far unknown stiffness matrices Ω 1 and Ω 2 . These matrices can be written in terms of the function f n and its derivatives up to third orther in t n . The use of the Taylor expansion of Ψ n+1 = Ψ(t n + δt) and of the Maclaurin series of in the recursion relation (15), leads to a simple system of algebraic equations for Ω 1 and Ω 2 . The components of the stiffness matrices obtained after solving these equations read as [17], where D j and E j (j = 1, ...., N ) are given in terms of the components of the derivatives f , provided that the denominator in Eq. (20) is not zero.
Fatunla [17] has established that his method is L-stable and exponentially fitted to any complex value λ. This means that the corresponding stability function R(λδt) = exp(λδt), gives the optimum stability properties. Furthermore, it can be shown that the jth component of the local truncation error at t = t n+1 is given by The implementation of Eq.(15) to calculate Ψ n+1 requires the calculation of the function f n and its first derivatives f (1) n at each value of t n , and also the stiffness matrices Ω 1 and Ω 2 to obtain the matrices R and S. We also calculate the truncation error T n+1 to control the size of the integration step imposing a boundary criterion for |T n+1 |. Note that to calculate the truncation error, we also need to evaluate f (4) n . The stiffness parameters carry the intrinsic information on the natural oscillations of the system. Due to this fact, Fatunla's scheme can afford larger values of the time step compared with other explicit methods of Runge-Kutta type [39].
In Fig.1, we show the evolution of the time step in Fatunla's propagation for our Gaussian model problem. The pulse envelope is also plotted in arbitrary units to illustrate the duration of the pulse. We see that the time step becomes increasingly large after the end of the pulse, reaching values of around 2 at the end of the total propagation (500 a.u. of time). It is clear that the most demanding part of the propagation, and therefore the most time consuming one, is during the interaction of the pulse with the system. This observation is important when it is necessary to propagate the wavefunction up to large distances after the end of the pulse, as is the case when the TSC method is used. In the results for the time step shown in Fig.1, the latter one is adapted according to the condition 10 −14 < |T n+1 | < 10 −9 , that is, if the truncation error is lower than the lower bound 10 −14 then we increase the time step, and if it is higher than the upper bound 10 −9 it is decreased. With this choice, the overall conservation of the norm is about 10 −5 , which is enough for the model problem case we are studying. For many physical problems, this level of accuracy in the norm is sufficient but, if a higher accuracy is needed, then we might expect that it is sufficient to shift the bounds of the truncation error. However, as shown in Fig.2, such a conclusion is not correct. In Fig.2, we consider three different constraints on the truncated error and calculate on a logarithmic scale, the absolute error on the norm denoted by ∆ as a function of time. This error is defined as the absolute value of the difference between 1 and the norm at time t . In these three cases, the time propagation is started with the same time step namely 10 −3 a.u. This time step always increases while the truncated error is smaller than the prescribed lower bound and decreases if the truncated error is above the upper bound. In all three cases, we observe a significant loss of accuracy in ∆ at the very beginning of the time propagation. As described by Madroñero and Piraux [39], this is due to initially very small values of the denominators in Eq. (20) which leads to inaccurate values of the stiffness matrix elements and of the truncated error. This problem is therefore intrinsically related to Fatunla's method and leads to difficulties in correctly controlling the time step. In fact , if we keep the time step constant from the beginning, we have a much better control of ∆.
We have also checked that this is true even in the field free case. On the other hand, in general, we see from the inset in Fig.2 that we maintain a higher accuracy when the constraint on the truncated error is more severe. In addition, we also observe several small jumps in ∆ the magnitude of which are much smaller than the jump in ∆ at the beginning of the time propagation. We attribute these jumps to an accumulation of roundoff errors. Indeed, we expect more roundoff errors in the case the constraint on the truncated error is the strongest since a smaller time step leads to a larger amount of calculations. Note that the jump observed in the red continuous line corresponds to a change of only one digit in the accuracy of the norm. The overall relative accuracy we obtain even for the most severe constraint we use on the truncation error is of the order of 10 −5 . To achieve a greater accuracy, it is necessary to use a fixed and very small time step. These results show clearly that the achievable accuracy for the adaptive time step approach in Fatunla's method has a lower bound for a given initial time step. As a result, the use of Fatunla's method rests on a compromise between the computer time required and the accuracy needed. In the following, we consider the interaction of helium with a strong laser pulse. In that case, the accuracy on the norm reduces to about 4 significant digits when Fatunla's method is used. This prevents us to calculate the probability of single ionization in various channels where the latter one is less than 10 −4 a.u. for field intensities currently used in the experiments.
In conclusion Fatunla's method allows one to treat stiff problems while fully exploiting the advantages of explicit schemes, namely that it only involves matrix vector multiplications. However it has its own limitations.

B. Krylov subspace method
In this section we consider a powerful method to propagate the TDSE solution, which provides accuracy of solutions and stability of propagation. It uses projection techniques on Krylov subspaces [40]. This approach was proposed by Arnoldi [41] in the calculation of the eigenstates of a matrix. Here we briefly recall the method used by Arnoldi as a time propagator [42], to solve the differential equation (5). Since the overlap matrix is positive-definite, we can use the Cholesky decomposition B = U † U to form an orthonormal basis defining the new coefficients Φ = UΨ. The TDSE for these coefficients is written in the form, where H = (U † ) −1 HU −1 . If we assume that the time interval is sufficiently small that the Hamiltonian may be treated as constant in time over a time step δt, it is trivial to demonstrate that Eq.(22) has a solution given by If H is diagonalizable and can be written as H = EΛE −1 , where Λ is a diagonal matrix with the eigenvalues λ i of H on the main diagonal and E is the matrix with the corresponding eigenvectors of H as its columns, then Eq.(23) can be reexpressed as follows However, for very large N this may be unnecessary and computationally very demanding. Instead, we can define the exponential in Eq.(23) using a Taylor expansion of the form, We use the successive matrix products as a basis set forming a Krylov subspace spanned by (m+1) linearly independent vectors, denoted by To build the Krylov subspace, we first use Gram-Schmidt orthogonalization of the initial vectors {Φ, HΦ, ..., H m Φ}, to obtain an orthonormal basis {q 0 , q 1 , ..., q m }. The procedure starts with q 0 = Φ/|Φ|, where the norm is defined as |Φ| = √ Φ † · Φ. The q k are obtained by calculating Hq k−1 and then orthonormalizing each vector with respect to q 0 , ..., q k−1 . If we define Q to be a matrix formed by the m + 1 column vectors (q 0 , ..., q m ), we finally get giving h = Q † HQ.
We see here that h is the Krylov subspace representation of the full Hamiltonian H, and that in this procedure, we obtain simultaneoulsy the Krylov vectors q 0 , ..., q m . Arnoldi's algorithm is general and applies to non-hermitian matrices. It reduces the dense matrix h to an upper Hessenberg form, and in the particular case of hermitian matrices, to a symmetric tridiagonal form. In this latter case, Lanczos has shown that this matrix can be obtained by means of a simple recursion formula. However, this formula is known to be problematic when the size of the Krylov subspace is large because the orthogonality of the Krylov vectors is rapidly lost [40]. It is the reason why we do not use this algorithm in the present case. Once we obtain the orthonormal Krylov subspace Q and the representation h of the Hamiltonian, it can be easily shown that Eq.(23) can be written as The matrix h for all our case studies is tridiagonal, and its size is never bigger than 100 × 100, so the calculation of the exponential through direct diagonalization, as in Eq. (24), is straightforward.
In actual numerical calculations, the Arnoldi algorithm [40] requires some modifications. After a first calculation of a new Krylov vector q j+1 , we ensure that the norm is equal to one, by re-checking the orthogonality against the previously calculated vectors, and perform again the Gram-Schmidt procedure if necessary. In principle the orthogonality condition determines the maximum size of the Krylov subspace and the algorithm can be used with a number m − 1 of vectors. Also, if we start generating the Krylov vectors from the ground state of the system, then Φ(t = t 0 ) is an eigenstate of the Hamiltonian, making it impossible to build a linearly independent set of Krylov vectors. To solve this problem, instead of using the vector Φ(t = t 0 ) as a starting point, we use a modified vector Φ(t = t 0 ) + , with a vector of random entries no larger than 10 −10 .
By construction (see Eq.23), the stability function associated to the numerical time propagator based on Arnoldi's algorithm is given by R(λδt) = exp(λδt). As a result, it has exactly the same stability properties as Fatunla's algorithm. However, it is worth remembering that Eq.(23) is only valid if the Hamiltonian is time independent. It is therefore a good approximation only for small values of δt. In the present case, there are two types of errors. The first one is directly related to Arnoldi's algorithm for the calculation of the exponential of a matrix. This type of error has been discussed in detail by Saad [43] and later on by Hochbruck and Lubich [44]. We have checked that this type of error is always negligible and does not depend on the time step. The second type of error is due to assuming that the Hamiltonian does not depend explicitly on time over the time step δt. We estimated this type of error by calculating dΨ/dt + iHΨ and checked that, as expected, it is of the order δt 2 . Another way to estimate this type of error is to compare our results with those obtained with an Arnoldi based method that takes explicitly into account the time dependence of the Hamiltonian. This can be done by using a Magnus expansion of the time evolution operator [45,46]. However, this method requires very time consuming calculations beyond the scope of this contribution. On the other hand, as it was already noted by other authors [42], enlarging the size of the Krylov space allows for larger time steps to be considered. In Fig.3, we give the number of Krylov vectors necessary to obtain convergence of our results as a function of the time step used in the calculations for the case of the one-dimensional Gaussian model potential with the same parameters as in Fig.1. In our calculations, the time step δt is kept constant during the propagation. The choice of the optimal value of the time step and of the corresponding dimension of the Krylov space is therefore the result of a compromise while trying to reduce the computer time.
The innovative use of the Arnoldi method as an explicit approach offers then the convenience that we only require matrix-vector and scalar products, which then transforms the method in a time-efficient approach as is the case for Fatunla's method. Furthermore, this particular scheme is norm-conserving with the advantage of providing a check for the method, even though it also means that it is not easy to quantify the error in the calculation of the norm. In the following sections the accuracy of both methods is tested in various situations by using a high order predictor-corrector method.

C. A predictor-corrector method.
In this subsection, we briefly describe the predictor-corrector (P-C) scheme we use to test the accuracy of both explicit methods described above. The predictor is either Fatunla's or Arnoldi's algorithm. The corrector is a fully implicit method of Runge-Kutta type which, here, is a four-stage Radau IIA method of order 7.
In a general Runge-Kutta method, the numerical solution Ψ n+1 of Eq.(8) at a given time t = t n+1 is obtained from the solution Ψ n at time t = t n as where δt is the time step and f (t i , Y i ) = −iB −1 H(t i )Y i with t i = t n + c i δt. b i and c i are coefficients defining the Runge-Kutta method for a number s of stages. We assume that the solution vector Ψ is of dimension N . The quantities Y i estimate the solution Ψ at the intermediate time t i . They are obtained by solving the following linear (sN × sN ) system of equations where the a ik are again given by the method. Solving such system represents the main difficulty of an implicit Runge-Kutta scheme. If this scheme is used for the corrector, we could in principle avoid solving such system of equations by using an iterative procedure in which we replace the vector Y k in the right hand side of Eq.(31) by Y where j gives the order of the iteration process. At the order 0, Y (0) k is provided by the predictor. However, such an iterative procedure is not stable. Instead, we follow a different iterative procedure that has been developed by van der Houwen and Sommeijer [47,48]. By introducing a diagonal matrix whose entries are calculated to guarantee optimum stability properties, they transform the (sN × sN ) system (31) into a set of uncoupled (N × N ) systems of equations that can be solved in parallel at each iteration. More precisely, they rewrite Eq.(31) as follows where the d ik are the entries of the diagonal matrix. The iterations in Eq.(32) start with Y (0) i , provided by the predictor. The iteration scheme is performed until a value j =max cor for which we have convergence. Then we can in Eq. (30) to obtain the solution at t = t n+1 . Once we have calculated Ψ n+1 , we can evaluate its norm and use its conservation as a criterion to monitor the size of the time step. Using the P-C method requires solving a large number of (N × N ) systems of equations. To solve these systems, we use an iterative method known as the biconjugate gradient stabilized method (Bi-CGSTAB) [49]. In order to reduce drastically the number of iterations, we use a pre-conditioner based on an incomplete LU factorization. In Fig.4, we consider the case of our one-dimensional Gaussian model potential with V 0 = 4 a.u. and β = 0.1 a.u. and a pulse of frequency ω = 0.5 a.u., duration of 8 optical cycles and peak intensity I = 10 16 Watt/cm 2 . We show, in this Fig.4, the multiplicity as a function of the number of iterations in the Bi-CGSTAB during the interaction. By multiplicity, we mean the number of times a given number of iterations is repeated during the whole propagation. It can be seen that without pre-conditioner, the number of iterations can grow significantly before reaching convergence, while using the preconditioned Bi-CGSTAB, the number of iterations is maintained below five. This reduces the computational time needed by 25%. However, care must be taken when including a pre-conditioner, since, by increasing the number of operations, it may increase the computational times even though it accelerates convergence. As mentioned above, the corrector scheme is iterated up to j =max cor where convergence is achieved. In Fig.5 we plot the time evolution of the time step during the propagation for different values of the maximum of iterations max cor in the corrector . We see here that, as we increase this maximum number of iterations, the value of the time step becomes larger. The relative error in the norm which is the same for all the calculations is of the order of 10 −11 . Moreover, the computational time with max cor =100 is half the time consumed for max cor =10 because it allows to use a much larger time step. It is therefore advisable to use large values of max cor to speed up the calculations.

A. Model potential
In this section we first present results for our case study of the one-dimensional Gaussian potential taking V 0 = 1 and β = 1. The electron wavepacket is developed in a basis of 200 B-splines and we use the time scaled coordinate method during the propagation [38]. We run our codes on an INTEL XEON 2.33 GHz Processor 51.40 (32 GB Ram). We choose a pulse of frequency ω = 0.7 a.u. and a full duration of 90 a.u. of time which corresponds to 10 optical cycles. The peak intensity I = 10 14 Watt/cm 2 . In Fig.6, the energy distribution is calculated by propagating the scaled wavepacket to a stationary state until a time of 1500 a.u. when convergence is achieved. The results shown are obtained using the two explicit propagators. Both methods converge to the same result but Fatunla's propagator uses 2.3 s of computer time with an adaptive time step while Arnoldi's propagator using five Krylov vectors and a fixed time step δt = 0.3 a.u. takes 6.2 s. For these values of intensity and frequency both methods give easily the correct result. However Arnoldi's method performs poorly from a computer time point of view. This can be understood by referring to Fig.1 where we show that Fatunla's propagator allows the use of ever larger time steps, particularly during the propagation after the end of the pulse, while Arnoldi's propagator keeps the same time-step throughout  To check how these methods behave in a more challenging case, we consider the same model potential with a pulse of frequency 0.1 a.u. with the same number of optical cycles and peak intensity. In this case, the total pulse duration is equal to 630 a.u. We see in Fig.7 that both methods give identical results. These results are obtained after propagating the wavepacket up to a time of 2500 a.u. The running time with Fatunla's propagator is equal to 958.27 s while in this case, Arnoldi's propagator performs better using 553.69 s for a subspace of 20 Krylov vectors and a fixed time step δt = 0.1 a.u. It can be seen that in general, Arnoldi's propagator performs better than Fatunla's propagator for long pulses.
To further probe these methods we increase the number of bound states supported by our potential by choosing V 0 = 4 and β = 0.1. The pulse has a frequency ω = 0.5 a.u. with a duration of 100.53 a.u. that corresponds to 8 optical cycles. The peak intensity I = 10 16 Watt/cm −2 . In this case, we use 1700 B-splines to propagate the wavepacket up to a time of 5000 a.u. Fig.8 shows the energy distribution obtained using Fatunla's propagator (straight line) and the predictor-corrector scheme (squares), which is used to test the accuracy of Fatunla's method.
Comparison of these two methods shows that Fatunla keeps the accuracy in the results down to a value of 10 −5 a.u. for the energy distribution. The TSC approach is used with an asymptotic scaling factor of 0.1. Fatunla's propagator takes 379.36 s while the P-C method with adaptive time step takes 1801.49 s. It is clear that Fatunla uses remarkably less computer time and works as long as the accuracy required is up to six digits. Fig.9 shows the energy distribution obtained with Arnoldi's propagator for the same parameters as in Fig.8. We show the results obtained with Arnoldi's approach and two different fixed time steps and compare these results with those obtained with the P-C scheme. The Krylov subspace contains 20 vectors and the wavepacket is again propagated up to 5000 a.u. The circles show the results for a time step δt = 0.03 a.u. and the stars for δt = 0.3 a.u. We note that increasing the time step leads to less accurate results by comparison with the P-C method. Arnoldi scheme takes 5957.19 s for a time step of 0.03 a.u. and 598.11 s for a time step of 0.3 a.u.

B. Hydrogen Atom
We now apply these methods to the more complex case of the interaction of hydrogen atom with a cosine square laser pulse. We use a spectral method based on the expansion of the wavefunction in a basis of Coulomb Sturmian functions [39], without implementing the TSC method. Unless otherwise stated, we performed all calculations on a laptop (with an INTEL core 2 duo processor of 2.4 GHz). The first pulse we use has a frequency of 0.7 a.u., a duration of 10 optical cycles and an intensity I=10 14 Watt/cm 2 , as in the case of Fig.6. In these rather simple conditions, we use 10 angular momenta. The non-linear parameter κ of the Coulomb Sturmian functions is taken equal to 0.3 a.u. Fatunla's and Arnoldi's algorithms produce the converged energy distribution as shown in Fig.10. The calculations carried on with Fatunla's propagator and an adaptive time step take 10.50 s of computer time. The integration performed with Arnoldi's method takes 13.72 s. For a time step of δt = 0.05 a.u. and 100 Coulomb Sturmian functions per angular momentum, it needs only 5 Krylov vectors. In this case Arnoldi's method is slower than Fatunla's method. This is related to the number of basis-set functions used. As this number increases, higher eigenvalues are generated in the Hamiltonian spectrum thereby increasing the stiff character of the system of equations to solve. In that case, more Krylov vectors have to be included to maintain the accuracy of the results. In Fig.11 we illustrate the effect of reducing the number of Krylov vectors n k from 5 to 4. It is surprising to see that the propagator gives a completely flat spectrum when the dimension of the Krylov space is insufficient. Fig.11 shows that for a basis set of 100 Coulomb Sturmian functions per angular momentum, accurate results for the energy distribution require a minimum of 5 Krylov vectors. It is interesting to note that to get an accurate spectrum, one of the eigenvalues of h must converge to 0.2 which corresponds to the position of the maximum of the spectrum which is what we expect from energy conservation (0.2 = −0.5+ω). If none of the eigenvalues converges to 0.2, the spectrum is completely flat because all the eigenvalues of h which are usually very high except the first one, do not contribute significantly to the spectrum. In addition, it is important to stress that decreasing the time step does not modify the minimum number of Krylov vectors to be used.
In Fig.12 we compare the performance of Arnoldi's and Fatunla's methods for a more difficult case. We consider a pulse of frequency ω = 0.114 a.u. and a duration of 20 optical cycles. The pulse intensity is the same as before, In Fig.13 we show results obtained for the challenging case of a pulse of very low frequency ω = 0.0228 a.u. and a duration of 4 optical cycles for the same intensity as before. To reproduce the energy distribution we need to use 1200 Coulomb Sturmian functions per angular momentum. 80 angular momenta are included in the calculations and the non-linear parameter κ of the Coulomb Sturmian functions is equal to 0.3. For this rather stiff problem Arnoldi's algorithm has to include a minimum of 70 Krylov vectors for a time step δt = 0.05 a.u. The calculation takes 24 hours on an 8 processor cluster using OpenMP. Fatunla's algorithm also reproduces the same energy distribution but the computer time used is more than four times larger. In fact, we observe that for larger scale problems where the degree of stiffness is important, Fatunla's method requires time steps that become prohibitively small thereby increasing the computational time.

C. Helium Atom
In this subsection, we show briefly results for the single ionization of helium by an intense electromagnetic pulse as an example of a more challenging problem. Following the remarks above, we perform here the calculations using only Arnoldi's algorithm. As mentioned before, Fatunla's algorithm is not accurate enough to calculate cross sections in various single ionization channels. The pulse has a peak intensity I=10 14 Watt/cm 2 , a frequency ω = 2.1 a.u. and a duration of 6 optical cycles. The wavefunction is expanded in a basis set that uses 140 B-spline functions of order 7 per electron angular momentum [24]. Three values of the total angular momentum (L=0,1,2) are taken into account and the maximum value of the individual electron angular momentum is three. The box size is 150 a.u. The step size during the interaction with the pulse is fixed at 0.01 a.u., while after the interaction the propagation used a step size of 1 a.u. The calculations are performed with 40 Krylov vectors. It takes 31 hours to run on a cluster with 10 Intel Xeon L5520 2.26 GHz processors using MPI (Message Passing Interface) and 3 GB of RAM per processor. Fig.14 shows the results obtained for the energy distribution of the single ionization of helium. As expected we observe a dominant peak at 1.2 a.u. which corresponds to the energy conservation. The spectrum is obtained by projecting the wave packet after the end of the pulse on a product of a Coulomb wave of the screened nucleus times a bound state FIG. 14. Single ionisation spectrum resulting from the interaction of an helium atom with a short cosine square pulse. Arnoldi's propagator is used. The pulse has a peak intensity of I=10 14 Watt/cm 2 , a frequency ω = 2.1 a.u. and a duration of 6 optical cycles. The basis-set of functions used is a set of 140 B-spline functions of order 7 per electron angular momentum. The total angular momentum L=0,1,2 and the maximum value of the individual electron angular momentum is three. The box size is is solved to obtain Ψ (r 1 , r 2 , t), with H given in Eq. (33). The question of decoherence of these quantum states can be studied through the quantum fidelity, which gives the overlap of the solutions of the TDSE, with and without the potential V anharmonic = ε r 2 1 2 + r 2 2 2 . The perturbation potential V p = r 2 1 2 + r 2 2 2 breaks the separability of Hooke's atom. We note that the Hamiltonian H ε in this case is not explicitly dependent on time and so it is different in nature to the Hamiltonians that we treated in previous examples. For a general Hamiltonian H 0 and a small real parameter ε that represents the strength of the perturbation, we write The quantum fidelity F ε at time t is defined as, where Ψ ε and Ψ 0 are the quantum states propagated with Eq.(35) for a perturbed and non-perturbed Hamiltonian, respectively. We can expand the quantum fidelity in terms of the perturbation parameter ε [51], as, with χ (t) being the quantum susceptibility. Taking the two first terms, we evaluate F ε (t) up to order ε 2 , valid near unity. For our particular case H 0 = H ε=0 and consequently Vp = r 2 1 2 + r 2 2 2 . The observable calculated in this problem is the susceptibility χ (t) and we take the harmonic frequency to be ω =1.0 a.u. and the perturbation parameter to be ε = 10 −5 . We study the evolution of the initial bound state of energy E=7 a.u. and vanishing angular momentum, singlet state with even parity [52]. The total number of functions in the basis set is 2370. The integration of the TDSE was first attempted using Fatunla's method. The stiffness of this problem forces the adaptive time step to become excessively small (of the order of 10 −5 ) so that the computer time needed by the method becomes of the order of several days instead of seconds. Furthermore the accuracy necessary to represent the effect of very small perturbations on the system could not be achieved. As a consequence we used Arnoldi's integrator, testing different combinations of the values of the time step and of the dimension of the Krylov subspace. It is worth stressing that the time evolution operator calculated within Arnoldi's method is essentially exact since the total Hamiltonian is time independent. However, the stiffness of the problem which is very strong because of the anharmonic character of the potential is expected to impose important constraints on the time step. In Fig.15 we show results for the quantum susceptibility using 5 Krylov vectors and two different time steps. In order to get converged results, this shows that we need a time step of at least δt = 10 −4 a.u., leading to a computational time of 4 hours. The same calculation performed with the P-C method took 17 days, 8 hours and 29 minutes. Fig.16 shows results for the observable χ (t) under the same conditions as in Fig.15 but using Krylov subspaces of higher dimension (n k = 7 and n k = 9). For n k = 7 converged results were obtained with a step size of δt = 5 × 10 −4 a.u. leading to a computational time of 45 minutes. However this figure illustrates the compromise to be achieved between the size of the time step used and the dimension of the Krylov subspace. For n k = 9, a time step of δt = 10 −3 a.u. leads to a calculation taking 30 minutes of computer time only. The choice of the optimal value of time step and of the Krylov subspace dimension needs to be balanced. This means to search for the optimal larger value of the time step for which the propagation will take less iterations. These calculations performed with the P-C method take 9 days, 17 hours and 14 minutes for n k = 7 and δt = 5 × 10 −4 and 8 days, 5 hours and 43 minutes for a.u. n k = 9 and δt = 10 −3 a.u. The computer used in these calculations was a single core of a Intel(R) Core (TM) 2 Quad CPU Q 9400(2.66 GHz) with 8 GB main memory.

V. CONCLUSIONS
In this contribution, we addressed the problem of the numerical integration of the time-dependent Schrödinger equation describing physical processes whose complexity requires the use of state of the art methods. The problem can be reduced to the solution of a system of first order differential equations. The main difficulties we have to face are the size of the system and its stiff character which results from the presence of very high energy eigenvalues in the Hamiltonian spectrum. These difficulties impose important constraints on the choice of the time propagator. Given the size of the system, this time propagator must be explicit. This means that it involves only matrix-vector products instead of solving large system of algebraic equations at each time step as is the case for implicit methods. In addition, this propagator must have optimum stability and accuracy properties to cope with the stiffness of the system. We have analyzed and compared the performance of two one-step explicit time propagators namely Fatunla's and Arnoldi's algorithms. It turns out that both of these methods share the same optimum stability properties. Nevertheless, we show that their accuracy properties differ significantly in most of the problems that we treat here. As a matter of fact, the accuracy of the method depends essentially on the stiffness of the system to solve which determines the appropriate choice of the propagator.
In all the problems considered here, the relative accuracy of Fatunla's method is always limited to about 10 −6 . In some cases, this might be sufficient but we should not forget that when the degree of stiffness increases, the adaptive time step becomes excessively small making the method inapplicable. By contrast, highly accurate results are obtained with Arnoldi's algorithm in all cases treated here. However, for a given time step, there is a minimal number of Krylov vectors to take into account. If the actual number used is smaller than this minimal number, generally there is an abrupt transition and the results are wrong giving a flat spectrum (in some cases this transition is not so abrupt but is rapid nevertheless.) On the other hand, when the degree of stiffness is high, this minimal number may become very large thereby imposing strong limitations on the applicability of the method. This is the case when the spacing between grid points becomes very small or, for spectral methods, when the size of the basis set is very large. In applying Arnoldi's scheme, it is therefore important to try to reduce the stiffness as much as possible. An obvious way to achieve this is to move to the atomic basis in which the Hamiltonian is diagonal and to eliminate the highest energy eigenvalues which, in principle do not play any physical role. In that case however, the ac-Stark shift of the levels will not be evaluated accurately. In addition, our calculations in the case of the Gaussian potential model clearly show that the energy electron spectrum calculated with Arnoldi's algorithm deteriorates.