Solving the Schr\"odinger eigenvalue problem by the imaginary time propagation technique using splitting methods with complex coefficients

The Schr\"odinger eigenvalue problem is solved with the imaginary time propagation technique. The separability of the Hamiltonian makes the problem suitable for the application of splitting methods. High order fractional time steps of order greater than two necessarily have negative steps and can not be used for this class of diffusive problems. However, there exist methods which use fractional complex time steps with positive real parts which can be used with only a moderate increase in the computational cost. We analyze the performance of this class of schemes and propose new methods which outperform the existing ones in most cases. On the other hand, if the gradient of the potential is available, methods up to fourth order with real and positive coefficients exist. We also explore this case and propose new methods as well as sixth-order methods with complex coefficients. In particular, highly optimized sixth-order schemes for near integrable systems using positive real part complex coefficients with and without modified potentials are presented. A time-stepping variable order algorithm is proposed and numerical results show the enhanced efficiency of the new methods.


I. INTRODUCTION
We consider the eigenvalue problem for the stationary Schrödinger equation (SE) ( = m = 1), where V (x) denotes the interaction potential and ∆ is the Laplacian operator. Since the Hamiltonian H is a Hermitian operator, its eigenvalues E i are real valued, and its corresponding real eigenfunctions φ i (x) form a basis of the underlying Hilbert space. This particular problem has attracted great interest among theorists and practitioners 1-3 due to its relevance for the understanding of the atomic and molecular structure of matter. A widely used approach to solve this problem is based on using the corresponding time-dependent Schrödinger equation in imaginary time (t = −iτ ), whose formal solution is given by the evolution operator exp(−τ H). In this way, in general, any initial condition, under the action of exp(−τ H), converges asymptotically to the ground state solution when τ → ∞. Notice that the evolution operator exp(−τ H) has the same eigenfunctions as the problem (1)- (2). This technique is usually referred to as the imaginary time propagation method (ITP for short). In this setting, only the action of exp(−τ H) on a wave function has to be computed 4,5 .
The ITP method can be regarded as an analog of the well-known power method in numerical linear algebra 6 .
In this sense, one may also consider the inverse power method: instead of the iterative application of the exponential operator exp(−τ H), the scheme v n+1 = (H − E i ) −1 v n , n = 0, 1, 2, . . . is used for some givenẼ i . This iteration is known to converge after normalization to the eigenvector with eigenvalue closest toẼ i . Although faster convergence than for the ITP method can be observed for an accurate initial guessẼ i ≈ E i , in general, the algorithm needs more iterations until convergence 7 .
Since the operators e −τ V and e −τ T can be exactly computed in the coordinate and momentum space, respectively, the operator splitting technique involving a composition of these exponential operators with appropriate coefficients can be used to approximate e −τ H . The computational cost depends on the number of changes between these coordinates which are cheaply performed by Fast Fourier transforms (FFT).
However, the operator splitting technique has some limitations. In particular, splitting methods of order p > 2 require negative time-steps 8,9 and the instabilities caused thereof are analogous to the ones for the integration of a diffusion equation backwards in time. If it is feasible to compute the gradient of the potential V , generalized splitting methods allow to build methods with positive coefficients up to fourth order [10][11][12] , but higher order methods also use negative time-steps. In this paper we propose methods to overcome the order barriers for both cases by using complex time-steps. Splitting methods can be tailored to particular equations to achieve better performances and we present criteria based on near-integrability that apply to a wide range of problems and thus yield highly efficient high order schemes. The obtained methods outperform the existing splitting schemes when high accuracy is desired and could be appropriate for elaborating a variable order algorithm. We also report some numerical experiments illustrating the efficiency of the new methods.

II. THE IMAGINARY TIME INTEGRATION METHOD FOR THE SCHRÖDINGER EQUATION
An important property of the Hermitian operator H is that (choosing properly the origin of the potential) its eigenvalues 0 ≤ E 0 ≤ E 1 ≤ . . . are real and nonnegative, and the corresponding eigenfunctions φ i can be chosen to form a real orthonormal basis on its domain. The problem (1) originates from the time-dependent SE i ∂ ∂t ψ(x, t) = Hψ(x, t), ψ(x, 0) = ψ 0 (x).
A Wick rotation of the time coordinate, t = −iτ , transforms (3) into a diffusion type equation with formal solution ψ(x, τ ) = e −τ H ψ(x, 0). After expanding the initial condition ψ 0 in the basis of eigen- where · | · is the usual L 2 scalar product, the time evolution of (4) is given by Asymptotically, for a sufficiently long time integration, we get ψ(x, τ ) → e −τ E0 c 0 φ 0 since the other exponentials decay more rapidly. The convergence rate depends of course on the separation of the eigenvalues. For simplicity, we restrict ourselves to the non-degenerate case E 0 < E 1 . If there is degeneracy, it converges to a linear combination of eigenfunctions, and repeating this process with different initial conditions one can obtain a complete set of independent vectors of the subspace which can be orthonormalized. Normalization of the asymptotic value yields the eigenfunction φ 0 and the corresponding eigenvalue is computed via E 0 = φ 0 |Hφ 0 . Excited states can be obtained by propagating different wave functions simultaneously (or successively) in time and using, for example, the Gram-Schmidt orthonormalization or diagonalizing the overlap matrix 7 .
For simplicity in the presentation, the spatial dimension is set to one unless it is explicitly stated, but our results also apply to higher dimensions.
The problem is further simplified by assuming x ∈ [a, b] with the interval [a, b] sufficiently large such that the wave function and all its derivatives of interest vanish at the boundaries. For numerical computations, the infinite dimensional domain of H has to be truncated, which is done by discretizing the spatial coordinate x: we fix N equally spaced grid points x i = x 0 + k∆x, k = 0, 1, 2, . . . , N − 1, with a = x 0 and b = x N . In this way, the interval is divided into N subintervals of size ∆x = (b − a)/N .
The potential V is represented in this grid by a diagonal matrix and the periodicity of the system (ψ (n) (a) = ψ (n) (b) = 0, n = 0, 1, 2, . . .) allows for the use of spectral methods (in space) for the calculation of T , namely the Fast Fourier Transform after which the matrix representation of T also becomes diagonal. The computational costs for the application of V and T to a vector are thus proportional to N and N log N operations, respectively. In a d-dimensional space with N mesh points on each dimension, their costs are proportional to N d and N d log N , respectively.

III. SPLITTING METHODS FOR THE SCHRÖDINGER EQUATION
To approximate the time evolution (5), i.e., the computation of e −τ H acting on a vector, we propose to use compositions of the operators e −τ V and e −τ T evaluated at different times. A first example is provided by the well-known Strang splitting verifying Ψ [2] h = e −hH + O(h 3 ) with h ≡ ∆τ . Higher order approximations can be obtained by a more general composition where Ψ [p] h = e −hH + O(h p+1 ) if the coefficients a i , b i are chosen such that they satisfy a number of order conditions (with m sufficiently large).
It is well-known, however, that methods of order greater than two (p > 2) necessarily have negative coefficients 8,9,13 (a simple proof can be found in Ref. 14 ). While this is usually not a problem for the coefficients b i , having negative a i coefficients makes the algorithm badly conditioned (in the limit N → ∞).
Composition methods with coefficients b i positive are also convenient for unbounded potentials, e.g., V (x) = x 2 , since negative values of b i can generate large roundoff errors in the exponential e −biV at the boundaries if the interval-size of the spatial discretization is not appropriately chosen and the potential takes exceedingly large values.
Splitting methods are particularly appropriate for the numerical integration of this problem since the choice of the time step, h, is not affected by the mesh size. Taking a finer mesh (i.e., a larger value of N ) does not necessarily lead to a smaller time step, and the extra computational effort originates only from the FFTs, whose cost is N log(N ) (or N d log(N ) in a d-dimensional problem with N points on each coordinate).
One possible approach to derive the order conditions to be satisfied by the coefficients a i , b i consists in applying the Baker-Campbell-Hausdorff formula to the composition (7), which we assume consistent ( i a i = i b i = 1) 15 . Thus we get Ψ where f i,j are polynomials of degree i in the coefficients a k , b k and the symbol [T, V ] stands for the commutator of the operators T and V . Condition f 2,1 = 0 leads to second order methods, and this can always be achieved by taking a left-right symmetric composition in (7) because all even terms automatically vanish. Methods of higher orders require in addition f 3,1 = f 3,2 = 0. Taking into account consistency, these equations can be written as 16 These two conditions imply that at least one of the a i as well as one of the b i become negative (see 14 and references therein), so that only methods of order two can be used for this problem. There are several possibilities to circumvent this limitation, and in the following, we enumerate some of them.
a. Modified potentials. If the kinetic energy operator in (4) is quadratic in momenta, then the nested commutator is diagonal in coordinate space. For this reason, (11) is usually called modifying potential. In consequence, involving two parameters. As a result, condition (10) becomes This equation can always be satisfied with a proper choice of the coefficients c i , so that the constraints on the coefficients a i , b i reduce to the single condition f 3,1 = 0, allowing for positive coefficients. In addition, solutions with positive c i coefficients also exist. A first example is the 4th-order composition 10,17 It turns out, however, that 6th-order methods using the operator (11) necessarily have some negative coefficients a i 18 . b. Near-integrable systems. When the Hamiltonian can be considered as a perturbed system, i.e., H = H 0 + εV ε (x) with an exactly solvable part H 0 = T + V 0 (x) and a small perturbation εV ε (x), it is advantageous to split the Hamiltonian into the dominant part H 0 and its perturbation εV ε . For example, if one is interested in the lower excited states, which evolve near the minimum of the potential, it can be useful to separate the quadratic part and to treat the remainder as a perturbation since the harmonic oscillator has a simple and fast solution using FFTs 19,20 .
Notice that in this case, the commutator depends only on the coordinates and modified potentials can also be applied as before. Then, all compositions remain the same except for replacing T by H 0 and V by εV ε .
With the additional information that one part of the operator is significantly smaller than the other, it is clear that the error expansion for a consistent method Ψ h can be asymptotically expressed as where the s i start from the first non-vanishing error coefficient e si,k . We say that Ψ h is of generalized order (s 1 , s 2 , . . . , s m ) (where s 1 ≥ s 2 ≥ · · · ≥ s m ) if the local error satisfies that Thus, for a method of generalized order (8, 2), denoted by Ψ This class of schemes can also be applied in several other situations. For instance, suppose one takes a sufficiently fine mesh. Then T ≫ V and the previous consider- c. Complex coefficients. A third possibility consists of considering complex coefficients in the composition (7) (with or without modified potentials). In other problems where the presence of negative real coefficients is unacceptable, the use of high-order splitting methods with complex coefficients having positive real part has shown to possess some advantages. In recent years a systematic search for new methods with complex coefficients has been carried out and the resulting schemes have been tested in different settings: Hamiltonian systems in celestial mechanics 21 , the time-dependent Schrödinger equation in quantum mechanics 22 and also in the more abstract setting of evolution equations with unbounded operators generating analytic semigroups 23,24 . It is worth noticing that the propagator exp(z∆) (z ∈ C) associated with the Laplacian is well-defined (in a reasonable distributional sense) if and only if Re(z) ≥ 0 23 , which is the case for the presented methods.
Many of the existing splitting methods with complex coefficients have been constructed by applying the composition technique to the symmetric second-order leapfrog scheme (6). For example, a fourth-order integrator can be obtained with the symmetric composition Ψ [4] h = Ψ [2] αh Ψ [2] βh Ψ [2] αh , where and k = 1, 2. In both cases, one has Re(α), Re(β) > 0. Higher order composition methods with complex coefficients and positive real part can be found in Refs. 23, 24, and 26, where several numerical examples are also reported.

IV. NEW SPLITTING METHODS FOR THE ITP PROBLEM
In this section, we carry out a systematic search of methods within the classes (a)-(c) above enumerated. The best methods for each subclass can be found online 25 with 25 digits of accuracy whereas the methods used in the numerical examples (Sec. V) are given in the subsequent tables with 18 digits for simplicity.
We only consider symmetric methods and, since T and V have qualitatively different properties, we analyze both TVT-and VTV-type compositions, defined as and respectively. In principle, both compositions have the same computational cost for the same number of exponentials. Nevertheless, due to a projection step to the real part after each full time-step, only in the VTV composition we can concatenate the last map in the current step with the first stage in the next one. The TVT compositions thus require two additional FFTs in comparison with the VTV composition, and this is accounted for in the numerical experiments. The methods we obtain are classified into two families: (I) methods without modified potentials and (II) methods with modified potentials. For each class we distinguish between methods for general problems (with the unique constraint that [V, [V, [T, V ]]] = 0) and methods for near-integrable problems (when the main dominant part contains the kinetic energy).
We have explored both TVT and VTV compositions with different number of stages. In some cases we consider extra stages to have free parameters for optimization. When the number and complexity of the order conditions is relatively low, we get all solutions. We then select the solutions having all of their coefficients with positive real part. Finally, we choose the solution which minimizes and/or minimizes the absolute value of the real part of the coefficients appearing at the leading error terms. These methods are subsequently tested on several numerical examples. After this process, we collect a number of schemes offering the best performance for most of the problems considered. In practice, however, one has to bear in mind that the relative performance between different methods depends eventually on the particular problem considered, the desired accuracy, the initial conditions, etc.

A. Methods without modified potentials
TVT and VTV compositions with 3 up to 9 stages have been analyzed. To simplify the notation, we denote compositions (16) and (17) as respectively. Here n indicates the order (or generalized order) of the method and m corresponds to the number of stages, i.e., the number of b i coefficients in the TVT composition or the number of a i coefficients in the VTV composition. The coefficients of the selected TVT methods are collected in Table I, whereas those corresponding to the TVT methods are displayed in Table II.

Methods for general problems
Analogously to (8), the symmetric compositions (16) and (17) can be formally expressed as a single exponential Ψ where the E i,j are chosen to form a basis of the algebra of commutators of length i. The chosen basis elements relevant for our exposition are Here we summarize some of the methods which have been analyzed: a. 3-stage compositions. A 3-stage composition has sufficient parameters to build 4th-order methods. There is one real solution and two complex solutions (conjugate to each other). For example, the VTV method corresponds to the composition (14) when Ψ [2] h is given by (6). The TVT version is obtained by interchanging T and V .
b. 5-stage compositions. Fourth-order methods with two free parameters can be obtained using 5-stage symmetric compositions. These two parameters can be used to build methods of effective order 6 (i.e., 4th-order methods that are conjugate to 6th-order methods by a near-identity change of variables). This requires to impose some additional constraints on the leading error terms, f 5,j , j = 1, 2, 3, 4. Specifically, these are f 5,1 − f 5,2 = 0 and f 5,3 + f 5,4 = 0 27 . We have found six solutions for the TVT composition and three solutions for the VTV composition with coefficients having positive real part. The solutions with smallest error terms at order 5 are denoted by T4 5 and V4 5 25 . c. 7-stage compositions. In principle, there are sufficient parameters to build 6th-order methods with 7 stages. For the TVT composition there are 11 solutions with all coefficients having positive real parts. The solution leading to a minimum value of the norm of the error at order 7 can be found online 25 .
With respect to the VTV composition, the best method we have found is identical with the most efficient sixth-order method obtained by Chambers 21 , where it has been presented as a symmetric composition similar to (14) but with 7 stages instead of 3, and with Ψ [2] h given by (6).

Methods for near-integrable problems
Proceeding analogously as before, we arrive at the following methods. We recall that in all compositions one should replace T by H 0 and V by εV ε . a. n-stage compositions of generalized order (2n, 2). This class of compositions has real and positive coefficients 28,29 . A 4-stage VTV composition of generalized order (8, 2) is given by scheme V84M LR 4 in Table IV with c 1 = 0.
b. 5-stage compositions. To build a method of generalized order (8,4) the following conditions must be satisfied by a consistent and symmetric method: f 3,1 = f 3,2 = f 5,1 = f 7,1 = 0. It requires at least 5 stages, and in this case only one solution with all coefficients having positive real part is found both for the TVT and VTV compositions. The coefficients of these methods, denoted by T84 5 and V84 5 , are collected in Table I and  Table II, respectively.
d. (8,6) methods. Increasing the number of stages to 9 we have two free parameters, which are used to satisfy in addition the following conditions: f 5,4 = f 7,2 = 0. In this way, methods of generalized order (8,6) and effective order (10,8,6) are obtained. Two efficient schemes correspond to T86 9 and V86 9 in Table I and Table II, respectively 30 .

B. Methods with modified potentials
Fourth-order methods incorporating modified potentials do exist with real and positive coefficients. In fact, 2-and 3-stage schemes have been extensively studied 11,12,18 . Methods of generalized order (n, 4) also exist with positive real coefficients 29 . Here we construct new methods of generalized order (6,4) and (8,4) with this property and generalize the treatment to 6th-order schemes with complex coefficients. In all cases, we take compositions TVT and VTV with up to 5 stages and denote them as Here, the parenthesis is used to help counting of the number of exponentials, and the letter M indicates that the methods use modified potentials. Notice that the number of free parameters can differ for the TVT and VTV sequences with the same number of exponentials because the exponent of a modified potential contains two parameters. The coefficients of the selected methods are collected in Table III and Table IV for the TVT and VTV compositions, respectively. The VTV composition allows one to build 6th-order methods, whereas the TVT needs an extra stage. There is only one solution (and its complex conjugate) with all coefficients having positive real part. It is denoted by V6M 4 and can be found online 25 .

Methods for near-integrable problems.
We first consider (n, 4) methods with real and positive coefficients. For schemes of generalized order (8,6) we collect only complex solutions with positive real part.
a. (6,4) methods They require at least 3 stages to satisfy the following order conditions: f 3,1 = f 3,2 = f 5,1 = 0. The coefficients a i and b i correspond to the methods (6,2) obtained in Ref. 28 (without modified potentials). We have also considered methods with 4 stages in order to have additional free parameters. As previously mentioned, there is the same number of order conditions as parameters to get a method of order 6 for the VTV sequence, but there are no solutions with coefficients being real and positive. To get a sixth-order method the following conditions must also to be satisfied: f 5,2 = f 5,3 = f 5,4 = 0. The coefficients c i only appear in f 5,3 and f 5,4 and can only be used to cancel these terms. The VTV sequence has three free parameters which can be used to annihilate f 5,3 and f 5,4 and to minimize the absolute value of f 5,2 under the constraint that all coefficients must be real and positive. The TVT sequence has only two free parameters which can be used to annihilate f 5,3 and to minimize the absolute value of the dominant term, f 5,2 , under the same constraint on the coefficients. The best methods we have obtained are denoted by T64M 4 and V64M 4 and are published online 25 .
b. (8,4) methods They require at least 4 stages. The coefficients a i and b i correspond to the methods (8,2) without using modified potentials and obtained in 28 . There is one coefficients c i in the TVT composition which can be used to cancel f 5,3 , and two coefficients c i in the VTV composition which can be used to annihilate f 5,3 and f 5,4 . The solution with c 2 = c 3 = 0 was already obtained in 29 . We have collected the corresponding coefficients for this method, V84M LR 4 , in Table IV. We have also considered methods with 5 stages in order to have an additional free parameters. There is the same number of order conditions as parameters to get a method of order (8,6) (which would be of order 6 for a general problem) but, obviously, there are no solutions with coefficients real and positive. As in the previous case, the term f 5,2 can not be zeroed using real positive coefficients. Then in both TVT and VTV compositions we have chosen the method which, while having real and positive coefficients,  c. (8,6) methods They require at least 5 stages and do not admit real and positive solutions for the coefficients and we are forced to consider complex solutions. We have found only one solution with positive real part in the coefficients for both TVT and VTV compositions. The coefficients for the methods denoted by T86M 5 and V86M 5 are given in Table III and Table IV, respectively.

A. Efficiency of the methods
As test bench for the numerical methods, we consider in the following two qualitatively different cases, the Pöschl-Teller potential and a perturbed harmonic oscillator, the latter being a classic example of a near-integrable system and of practical interest 3 . These two problems can be numerically integrated using modified potentials. However, we compare the relative performance of the methods (with and without modified potentials) separately in order to study the performance of the methods when it is not feasible to compute the gradient of the potential.  The numerical integration proceeds as follows: starting from random initial data, we iterate with fixed time-step until the sufficiently large final time T = 100 and compare the result with the exact solution, ψ(T ), which has been obtained by integrating with a much smaller time step. The spatial interval is fixed for all experiments to [−10, 10] and is discretized with N = 128 equidistant mesh points. Similar results are obtained for larger N = 256, 512, 1024. At each step, we project the obtained vector to its real part and normalize it to one in ℓ 2 (R), i.e., given the method Ψ [p] h and initial conditions, u n ∈ R N , we compute u n+1 as h u n ; then, sinceũ n+1 is a complex vector (but O(h p ) away from a real vector) we project on the real space by removing the imaginary part u n+1 = Re(ũ n+1 ) and then normalize the solution u n+1 =ū n+1 ū n+1 , where the norm is given by  This procedure will allow us to determine the efficiency of the new splitting methods, which will depend on the desired accuracy, and thereby choose the methods which are most appropriate for implementation with a more efficient algorithm that is based on variable time step and order. We distinguish two types of problems: on the one hand, methods that include modified potentials, the reference methods being Chin-4M (13), OMF-4M 11 and V84M LR 4 29 given in Table IV as well as a differently optimized scheme SCF-4M 31 and on the other hand, methods without modifying potentials with the reference methods V82 28 , the fourth order complex triple-jump scheme (14), referenced as Yoshida 4 and a 6th-order complex coefficient method by Chambers 21 . We remark that all relevant methods in the cited papers have been tested and the most efficient ones for this problem are included in the plots.

Pöschl-Teller potential
We have chosen the well-known one-dimensional Pöschl-Teller potential for the availability of analytic solutions of the eigenstates with λ(λ + 1) = 10. The results of our computation are shown in Figure 1a. The higher order of the complex coefficient methods outweighs their extra cost starting from moderate accuracy. The optimizations of the error terms can be clearly appreciated in the comparison with the 4th order triple-jump (14). When we consider methods with modified potentials, we observe that the new methods show only slight improvements with respect to the method OMF-4M since both parts of the splitting T and V are of comparable size. As the desired precision is increased, the new sixth order methods dominate in efficiency.

Perturbed harmonic oscillator
To illustrate the benefits of methods designed for near integrable systems, we use the Hamiltonian and split it in a large part H HO = − 1 2 ∂ 2 ∂x 2 + 1 2 ω 2 x 2 and a small part εV ε (x). The trap frequency is set to ω = 1 and the perturbation εV ε is given by the Pöschl-Teller potential in (19), with λ(λ + 1) = 2/5. The harmonic part H HO can be solved exactly via an exact splitting using Fourier transforms, cf. 19 , where it is shown that for |δω| < π and p 2 ≡ − ∂ 2 ∂x 2 . From the computational point of view, it is suggested 20 to consider the VTV split instead of the TVT split because it can be concatenated with the perturbation which only depends on the coordinates and no additional FFTs are necessary, i.e.
In 19 this decomposition is generalized to the twodimensional problem H = 1 2 (p 2 x + p 2 y ) + 1 2 (w 2 1 x 2 + w 2 2 y 2 ) − Ω(xp y − yp x ) and in 20 to the non homogeneous and possibly time-dependent one-dimensional problem H = for |Im(h)ω| < π and Re(h) > 0 (for numerical stability) and the perturbation part is easily propagated after discretization by the exponential of a diagonal matrix.
In this setting, the higher order in the small parameter is amplified and the efficiency plots in Figure 1b indicate that the new methods outperform the existing ones when high precision is sought and overall when modified potentials are allowed. We observe in both examples that, when modified potentials can be computed without exceedingly large computational cost, they should be used.
Further numerical experiments show that the efficiency curves are independent of the mesh size, i.e., the norm of T , and the cost only increases as N log(N ) as expected. The reason for this can be understood by following the evolution of the state vector along the iterations of the algorithm. Whereas in the beginning one has a non-smooth configuration u 0 , after a few steps the vector u i is close to an eigenstate and thus smoothened.
It is important to remark that the methods proposed in this work can be implemented in an algorithm which uses variable step, variable order, variable mesh size and variable simple-double precision. The best implementation can depend on the class of problems to be solved. For illustration, we present an implementation with variable time steps.

B. Variable step method
The previous examples show that for low accuracies and large time steps, the (8,2) method (with real coefficients) performs best. However, if we allow for variable time steps, as proposed in 5,7 , the computational cost is drastically reduced. We propose an improved timestepping algorithm that is based on two different estimators for the eigenvalue. In the first row, efficiency curves (error vs. number of FFTs) for methods without force evaluations are presented, with the new methods (triangles) performing best for high accuracies. The middle rows depicts methods based on modified potentials. In the right column, T86M5 intersects with V84M5 at precision 10 −13 , whereas it already improves on T84M5 at 10 −9 for the left column. SCB-4M overlaps with Chin-4M and has thus been omitted in the plot. In the bottom row, the random initial conditions (green), the ground states (black) and the potentials (dashed blue), scaled by 1/5 to fit the axis, are shown.
Recall that fixing the time-step and iterating to convergence will yield an eigenvector with the error being of the order of the method O(h p ) since we are computing exactly the spectrum of a perturbed Hamiltonian. Assume now that we are close to convergence, i.e, one has obtained an eigenvector u n = v 0 +O (h p ) and we consider the decomposition in the basis of exact eigenvectors v i of H, It is clear that d i = O(h p ), i > 1 and due to the normalization d 0 = 1 + O(h 2p) ). Then, an energy estimation is given by Alternatively, the energy can be estimated by the loss of norm in each time step, and then Combining both expressions yields an error estimate for the energy, The convergence in energy is measured by comparison with the previous time step, The time stepper then works as follows: starting from a large step size, the time step is decreased by a factor 1/2 whenever the actual reduction in energy of the iteration δE falls below the the maximally reachable precision ∆E, i.e., |δE| < (∆E) 2 and the iteration is terminated once the error estimate ∆E has reached a given tolerance. For the numerical experiments, we use the same configurations as for constant time step but terminate the algorithm when convergence in energy is reached at ∆E < 10 −10 . The iterations are initialized with random normalized data and a time step of τ = 10. The results are displayed in Figure 2a for the Pöschl-Teller potential and in Figure 2b for the perturbed harmonic oscillator with the same parameters as in the fixed-step size experiments. The error is measured as the ℓ 2 norm of the difference between the current value of the algorithm ψ(t) and the exact ground state φ(T ) as in the previous experiments, error = ψ(t) − φ(T ) .
As expected, it is apparent that lower order methods show better smoothing behavior for the first steps, when the wave function is still rough (recall that the algorithm is initialized by a worst-case wave function). For higher precisions, the new methods clearly outperform the existing ones, with the sole exception of the unperturbed setting with modified potentials, where the globally optimized OMF-4M method can hardly be improved unless extremely high precision is sought and the 6th order methods of Table IV and III become favorable (not  shown). Finally, if one is interested in very high accuracies, high order extrapolation methods 1,32 can be used for the last part of the time integration.
The results indicate that for low precision, i.e., for the first iterations, a lower order method should be used and then, after a certain precision is reached, e.g., when the higher order methods exhibit their superiority the algorithm should change to the optimal method, either V864 7 or V86M 5 until convergence. Further preliminary experiments on this adaptive order strategy have shown that there is plenty of room for optimization, e.g., by changing the initial step-size, adjusting the step-size by a different factor or by modifying the control criterion. Each of which has certain advantages and disadvantages, depending on the initial conditions and the range of precision.
For excited states, one expects an even better performance of the new methods since several states have to be computed to high precision in order to avoid error accumulation and the gains of the new methods are thus amplified. We have confirmed this conjecture by numerical experiments. The results thereof are omitted in the manuscript since they do not contribute insight beyond the presented experiments: they are qualitatively identical.

VI. CONCLUSIONS AND OUTLOOK
We have studied the Schrödinger eigenvalue problem by the imaginary time propagation method and proposed splitting schemes with positive real coefficients using modified potentials as well as with complex coefficients that can overcome the order barrier for parabolic problems since the coefficients have only positive real parts. The obtained sixth order methods are clearly superior to any classical ones for high precisions. On the other hand, when the gradient of the potential can be cheaply evaluated, the high order methods with complex coefficients are efficient only at very high accuracies due to the double cost caused by complex arithmetic.
We have proposed different high order methods to reach highly accurate results. An efficient implementation should take into account, for example, a preliminary time integration on a coarse mesh using simple precision arithmetic in order to get, as fast as possible, a smooth and relatively accurate solution from a random initial guess, and next consider a refined mesh using arithmetic in double precision. For simple precision arithmetic and low accuracies, it suffices to consider only low order methods, and when higher accuracies are desired we turn to double precision, variable time step and variable order methods. The best algorithm could depend on the class of problems to solve.
It is also important to remark that the form of the exponent allows that the techniques presented in this work can also be transferred to other areas whenever splitting is appropriate and the integration has to be performed forward in time, e.g., statistical mechanics of quantum systems, where one has to compute the Boltzman operator exp(−βH), with β = (kT ) −1 or quantum Monte-Carlo simulations 22 .
Finally, we would like to mention that real time integration with complex coefficients is under investigation.
To compute e −iaitT requires complex FFTs, and this is irrespective of the coefficients a i being real or complex. However, the constraint Im(a i ) ≤ 0 and the consistency condition, i a i = 1, necessarily requires a i ∈ R, while b i can be complex. A large number of new methods have been explored, but the superiority is not yet clear since there exist highly efficient methods with real coefficients for perturbed problems 27,33 and using modified potentials 34 .