Symmetric-conjugate splitting methods for linear unitary problems

We analyze the preservation properties of a family of reversible splitting methods when they are applied to the numerical time integration of linear differential equations defined in the unitary group. The schemes involve complex coefficients and are conjugated to unitary transformations for sufficiently small values of the time step-size. New and efficient methods up to order six are constructed and tested on the linear Schrödinger equation.


Introduction
We are concerned in this work with the numerical integration of the linear ordinary differential equation where u ∈ C N and H ∈ R N ×N is a real matrix.A particular example of paramount importance leading to eq. (1.1) is the time-dependent Schrödinger equation once it is discretized in space.In that case H (related to the Hamiltonian of the system) can be typically split into two parts, H = A + B .The equation with y ∈ R d , K ∈ R d×d can also be recast in the form (1.1) if the matrix K satisfy certain conditions [6].
Although the solution of (1.1) is given by u(t) = e itH u 0 , very often the dimension of H is so large that evaluating directly the action of the matrix exponential on u 0 is computationally very expensive, and so other approximation techniques are desirable.When H = A + B and e itA u 0 , e itB u 0 can be efficiently evaluated, then splitting methods constitute a natural option [17].They are of the form for a time step h .Here a j , b j are coefficients chosen in such a way that S h = e ihH + O(h p+1 ) when h → 0 for a given p ≥ 1 .After applying the Baker-Campbell-Hausdorff (BCH) formula, S h can be formally expressed as S h = exp (ihH h ) , with iH h = iH o h + H e h and Here [A, B] := AB − BA , g k,j are polynomials of degree k in the coefficients a i , b i verifying g 1,1 = g 1,2 = 1 (for consistency), and g k,j = 0, k = 1, 2, . . ., p, ∀j for achieving order p .
If A and B are real symmetric matrices, then [A, B] is skew-symmetric and [A, [A, B]] is symmetric.In general, all nested commutators with an even number of matrices A, B are skew-symmetric and those containing an odd number are symmetric, so that (H o h ) T = H o h and (H e h ) T = −H e h .When the coefficients a j , b j are real, then g k,j are also real and therefore S h = e ihH h is a unitary matrix.In addition, if the composition (1.2) is palindromic, i.e., a 2n−j = a j , b 2n−1−j = b j , j = 1, 2, . . ., then g 2k,j = 0 and H −h = H h , thus leading to a time-reversible method, S −h = S −1 h .In other words, if u n denotes the approximation at time t = nh , then S −h (u n+1 ) = u n .As a result, one gets a very favorable long-time behavior of the error for this type of integrators [16].Thus, in particular, M(u) := |u| 2  (norm) and H(u) := ūT Hu (expected value of the energy) are almost globally preserved.Recently, some preliminary results obtained with a different class of splitting methods (1.2) have been reported when they are applied to the semi-discretized Schrödinger equation [4].These schemes are characterized by the fact that the coefficients in (1.2) are complex numbers.Notice, however, that in this case the polynomials g k,j ∈ C , so that S h = e ihH h is not unitary in general.This is so even for palindromic compositions, since g 2ℓ+1,j are complex anyway.
There is nevertheless a special symmetry in the coefficients, namely a 2n−j = a j and b 2n−1−j = b j , j = 1, 2, . . ., (1.3) worth to be considered.Methods of this class can be properly called symmetric-conjugate compositions.In that case, a straightforward computation shows that the resulting composition satisfies for real matrices A and B , and in addition (S h ) T = S −h (1.5) if A and B are real symmetric.In consequence, for certain real matrices Ĥo h (symmetric), and Ĥe h (skew-symmetric).Since i Ĥe h is not real, then unitarity is lost.In spite of that, the examples collected in [4] seem to indicate that this class of schemes behave as compositions with real coefficients regarding preservation properties, at least for sufficiently small values of h .Intuitively, this can be traced back to the fact that i Ĥe h = O(h p ) and is purely imaginary.One of the purposes of this paper is to provide a rigorous justification of this behavior by generalizing the treatment done in [4] for the problem (1.1) defined in the group SU(2), i.e., when H is a linear combination of Pauli matrices.In particular, we prove here that, typically, any consistent symmetric-conjugate splitting method applied to (1.1) when H is real symmetric, is conjugated to a unitary method for sufficiently small values of h .In fact, this property can be related to the reversibility of the map S h with respect to complex conjugation, as specified next.
Let C be the linear transformation defined by C(u) = u for all u ∈ C N .Then, the differential equation In other words, the map S h (u) is C -reversible [12] (or reversible for short).Notice that this also holds for palindromic compositions (1.2) with real coefficients.
In the sequel we will refer to compositions verifying (1.3) as symmetric-conjugate or reversible methods.Splitting and composition methods with complex coefficients have also interesting properties concerning the magnitude of the successive terms in the asymptotic expansion of the local truncation error.Contrarily to methods with real coefficients, higher order error terms in the expansion of a given method have essentially a similar size as lower order terms [3].In addition, an integrator of a given order with the minimum number of flows typically achieves a good efficiency, whereas with real coefficients one has to introduce additional parameters (and therefore more flows in the composition) for optimization purposes.It makes sense, then, to apply this class of schemes to equation (1.1) and eventually compare their performance with splitting methods involving real coefficients, since in any case the presence of complex coefficients does not lead to an increment in the overall computational cost.
The structure of the paper goes as follows.In section 2 we provide further experimental evidence of the preservation properties exhibited by C -reversible splitting methods applied to different classes of matrices H by considering several illustrative numerical examples.In section 3 we analyze in detail this type of methods and validate theoretically the observed results by stating two theorems concerning consistent reversible maps.Then, in section 4 we present new symmetric-conjugate schemes up to order 6 specifically designed for the semi-discretized Schrödinger equation and other problems with the same algebraic structure.Finally, these new methods are tested in section 5 for a specific potential.

Symmetric-conjugate splitting methods in practice: some illustrative examples
To illustrate the preservation properties exhibited by symmetric-conjugate (or reversible) methods when applied to (1.1) with H = A + B , we consider some low order compositions of this type.Specifically, the tests will be carried out with the following schemes: Order 3. The simplest symmetric-conjugate method corresponds to with and was first obtained in [1].In addition, and as a representative of the schemes considered in section 4, we also use the following method, with a j > 0 and b j ∈ C , ℜ(b j ) > 0 : where h = e ihb 0 B e iha 1 A e ihb 1 B e iha 2 A e ihb 1 B e iha 1 A e ihb 0 B , ( but now When the matrix H results from a space discretization of the time-dependent Schrödinger equation (for instance, by means of a pseudo-spectral method), then it is real symmetric and A , B are also symmetric (in fact, B is diagonal).It makes sense, then, start analyzing this situation, where, in addition, all the eigenvalues of H are simple.To proceed, we generate a N × N real matrix with N = 10 and uniformly distributed elements in the interval (0, 1) , and take H as its symmetric part.The symmetric matrix A is generated analogously, and finally we fix B = H − A .Next we compute the approximations obtained by S and S [4] h for different values of h , determine their eigenvalues ω j and compute the quantity for each h .Finally, we depict D h as a function of h .
Figure 1 (left) is representative of the results obtained in all cases we have tested: all |ω j | are 1 (except round-off) for some interval 0 < h < h * , and then there is always some ω ℓ such that |ω ℓ | > 1 .In other words, S and S [4] h behave as unitary maps in this interval.This is precisely what happens in the group SU(2), as shown in [4].
The right panel of Figure 1 is obtained in the same situation (i.e., H real symmetric with simple eigenvalues), but now both A and B are no longer symmetric: essentially the same behavior as before is observed.Of course, when h < h * , both the norm of u , M(u) , and the expected value of the energy, H(u) are preserved for long times, as shown in [4].
Our next simulation concerns a real (but not symmetric) matrix H with all its eigenvalues real and simple.Again, there exists a threshold h * > 0 such that for h < h * the schemes render unitary approximations.This is clearly visible in Figure 2 (left panel).If we consider instead a completely arbitrary real matrix H , then the outcome is rather different: D h > 0 for any h > 0 (right panel; for this example D h = 9.79 • 10 −4 already for h = 0.001 ).
Next we illustrate the situation when the real matrix H has multiple eigenvalues but is still diagonalizable.As before, we consider first the analogue of Figure 1  symmetric matrices (Figure 3, left panel) and A and B are real, but not symmetric (right panel).In the first case we notice that, whereas all the eigenvalues of the approximations rendered by S and S [4] h still have absolute value 1 for some interval 0 < h < h * , this is clearly not the case of S [3,2] h .If, on the other hand, the splitting is done is such a way that A and B are not symmetric (but still real), then D h > 0 even for very small values of h .The same behavior is observed when H is taken as a real (but not symmetric), diagonalizable matrix with multiple real eigenvalues.
The different phenomena exhibited by these examples require then a detailed numerical analysis of the class of schemes involved, trying to explain in particular the role played by the eigenvalues of the matrix H in the final outcome, as well as the different behavior of S [3,1] h and S [3,2] h .This will be the subject of the  next section.
3 Numerical analysis of reversible integration schemes

Main results
We next state two theorems and two additional corollaries that, generally speaking, justify the previous experiments and explain the good behavior exhibited by reversible methods.
Theorem 3.1 Let H ∈ R N ×N be a real matrix and let S h ∈ C N ×N be a family of complex matrices depending smoothly on h ∈ R such that • S h is a reversible map in the previous sense, so that • S h is consistent with exp(ihH) , i.e. there exists p ≥ 1 such that • the eigenvalues of H are real and simple.
Then there exist • D h , a family of real diagonal matrices depending smoothly on h , • P h , a family of real invertible matrices depending smoothly on h , such that ) and, provided that |h| is small enough, Corollary 3.2 In the setting of Theorem 3.1, there exists a constant C > 0 such that, provided that |h| is small enough, for all u ∈ C N and all eigenvalues ω ∈ σ(H) , one has where Π ω denotes the spectral projector onto Ker(H − ωI N ) .Moreover, if H is symmetric, the norm and the energy are almost conserved, in the sense that, for all u ∈ C N , it holds that where M(u) = |u| 2 and H(u) = u T Hu .
Proof:[Proof of Corollary 3.2] First, we focus on (3.3).We note that by consistency, we have Since the eigenvalues of H are simple, it follows that the spectral projectors are all of the form where e 1 , . . ., e N denotes the canonical basis of R N .Then, we note that for all n ∈ Z , we have Therefore, since e inhD h is uniformly bounded with respect to h and n (because D h is a real diagonal matrix) and P h = P 0 + O(h p ) , it follows that where the implicit constant in O term does not depend on n (here and later).Therefore, it is enough to use the explicit formula (3.5) to prove that As a consequence, the estimate (3.3) follows directly by the triangular inequality : Now, we focus on (3.4).Here, since H is assumed to be symmetric, its eigenspaces are orthogonal.Therefore by the Pythagorean theorem, we have As a consequence, (3.4) follows directly of (3.3). 2 The main limitation of Theorem 3.1 is the assumption on the simplicity of the eigenvalues of H . Indeed, even if this assumption is typically satisfied, it depends only on the equation we aim at solving and not of the numerical method one uses.The following theorem, which is a refinement of Theorem 3.1, remedies this point by making an assumption on the leading term of the consistency error (which is typically satisfied for generic choices of numerical integrators).Theorem 3.3 Let H ∈ R N ×N be a real matrix and let S h ∈ C N ×N be a family of complex matrices depending smoothly on h such that • S h is consistent with exp(ihH) , i.e.
where p ≥ 1 is the order of consistency and R is a real matrix1 ; • H is diagonalizable and its eigenvalues are real; • for all ω ∈ σ(H) , the eigenvalues of Π ω R |Eω(H) are real and simple, where Π ω denotes the spectral projector on E ω (H) := Ker(H − ωI N ) .
Then there exist • D h , a family of real diagonal matrices depending smoothly on h , • P h , a family of real invertible matrices depending smoothly on h , such that, both P −1 0 BP 0 and P −1 0 HP 0 are diagonal, where B := ω∈σ(H) Π ω R Π ω , and provided that |h| is small enough, it holds that Corollary 3.4 In the setting of Theorem 3.3, there exists a constant C > 0 such that, provided that |h| is small enough, for all u ∈ C N , all ω ∈ σ(H) and all λ ∈ σ(Π ω R |Eω(H) ) , we have where P λ,ω denotes the projector along ) .Moreover, if H and R are symmetric, for all ω ∈ σ(H) , one gets and the mass and the energy are almost conserved, i.e. for all u ∈ C N , it holds that where, as before, M(u) = |u| 2 and H(u) = u T Hu .
Proof:[Proof of Corollary 3.4] The proof is almost identical to the one of Corollary 3.2.The key point is that, since both P −1 0 BP 0 and P −1 0 HP 0 are diagonal, then the projectors P λ,ω are exactly the projectors Π (j) , 1 ≤ j ≤ N (given by (3.5)).Note that, contrary to Theorem 3.1, in Theorem 3.3 one does not claim that P h = P 0 +O(h p ) .A priori, here, in general, the best estimate we expect is P h = P 0 +O(h) (which follows directly from the smoothness of P h with respect to h ).It is this loss which explains why, in Corollary 3.4, the error terms are of order O(h) whereas they are of order O(h p ) in Corollary 3.2. 2 Remark.Before starting the proof of these theorems, let us provide some comments about the context and the ideas involved.
• In Theorem 3.1 and its proof, we are just putting S h in Birkhoff normal form.The fact that S h can be diagonalized is due to the simplicity of the eigenvalues of H while the fact that its eigenvalues are complex numbers of modulus 1 is due the reversibility of S h .This approach is robust and well known, in particular it can be extended to the nonlinear setting (see e.g.[12, section V.1]).Note that here, we reach convergence of the Birkhoff normal form because the system is linear.
• Theorem 3.3 is a refinement of Theorem 3.1.To prove the absence of resonances due to the multiplicity of the eigenvalues of H , we use the first correction to the frequencies generated by the perturbation of H (i.e., the projections of R in Theorem 3.3).This approach is typical of what one does in the proof of Nekhoroshev theorems or KAM theorems (see also [12]).
• In order to give some intuition about the proof and the assumptions of Theorem 3.1, let us prove simply that, provided h is small enough, S h is conjugated to a unitary matrix.Indeed, since S h is reversible it writes as where ) is a real matrix (provided that h is small enough).Now, since the set of the real matrices whose eigenvalues are simple and real is open in the space of the real matrices (by continuity of the eigenvalues) and H h is a real perturbation of such a matrix ( H by assumption), we deduce that, provided h is small enough, its eigenvalues are simple and real.This implies that H h is conjugated to a real diagonal matrix and so that S h is conjugated to a unitary matrix.

Technical lemmas
In the proof of the previous theorems we will make use of the following three lemmas.
Lemma 3.5 Let M be a complex matrix and let P be a complex invertible matrix.Then ad P −1 M P and ad M are similar.More precisely, where int P M := P −1 M P .Here ad M stands for the adjoint operator: ad M X := [M, X] = M X−XM , for any matrix X .
Proof: A straightforward calculation shows that, for any X , If M 0 is diagonalizable on C , then there exists a family of real matrices χ h , depending smoothly on h , such that if |h| is small enough, e −h p χ h M h e h p χ h commutes with M 0 , i.e.
[e h p χ h M h e −h p χ h , M 0 ] = 0.
Proof: We aim at designing the family χ h as solution of the equation Thanks to the well known identity e A B e −A = e ad A B , this equation rewrites as Next we write the Taylor expansion of M h at order p as where R h is a family of real matrices depending smoothly on h .Then, isolating the terms of order 0 (and dividing by h p ), the equation (3.9) leads to where φ 1 (z) := e z −1 z .We restrict ourselves to χ h in Im R ad M 0 and consider f as a smooth map from R × Im R ad M 0 to Im R ad M 0 .To solve the equation f (h, χ h ) = 0 using the implicit function theorem, we just have to design χ 0 so that Actually, these properties are clear because the first one is a consequence of the second one, whereas the second follows directly from Lemma 3.6.2

Proofs of the theorems
We are now in a position to prove Theorems 3.1 and 3.3.Without loss of generality, and to simplify notations, we assume that H is diagonal where Thanks to the consistency assumption (3.6) (which is equivalent to (3.1)), provided that |h| is small enough, S h rewrites as Moreover, the reversibility assumption S −1 h = S h implies that H h is a real matrix (provided that |h| is small enough).Note that, hence, we deduce that R is also a real matrix.Then, applying Lemma 3.7 to H h , we get a family of real matrices χ h such that, provided that |h| is small enough, We conclude that W h is block-diagonal (with the same structure of blocks as H ), i.e. there exists some n j × n j real matrices W (j) h such that As a consequence, if the eigenvalues of H are simple (i.e.d = N and n j = 1 for all j ) then W h is diagonal.Therefore, in this case, it is enough to set P h = e −h p χ h and W h = D h to conclude the proof of Theorem 3.1.
So, from now on, we only focus on the proof of Theorem 3.3.First, we aim at identifying the matrices on the blocks in (3.10).The Taylor expansion of W h is clearly However, since [W h , H] = 0 , we deduce that [B, H] = 0 and so that B is block-diagonal.Moreover, since the matrix [χ 0 , H] is identically equal to zero on the diagonal blocks, the diagonal blocks of B are exactly those of R .As a consequence, with a slight abuse of notations, we may write and Y (j) h is a family of real matrices depending smoothly on h .
Next we aim at diagonalizing these blocks.By assumption, the eigenvalues of each matrix B (j) are real and simple.Therefore, all B (j) are diagonalizable.As a consequence, and again by applying Lemma 3.7, we get a family of real matrices Υ (j) h such that if |h| is small enough, for all j ∈ 1, d we have e hΥ (j) h (B (j) + hY This means that the eigenspaces of B (j) are stable by the action of e hΥ (j) h (B (j) + hY h )e −hΥ (j) h .However, by assumption, these spaces are lines.Therefore, if Q (j) is a real invertible matrix such that Q h )e −hΥ (j) h (Q (j) ) −1 is also diagonal.Finally, as a consequence, setting (1)  h Q (1)  . . .
we have proven that D h := P −1 h H h P h is real diagonal, which concludes the proof of Theorem 3.3.

Applications to reversible splitting and composition methods
Theorems 3.1 and 3.3 shed light on the behavior observed in the examples collected in Section 2. Thus, suppose H = A + B is a real symmetric matrix, with A , B also real.Furthermore, consider a splitting scheme S h of the form (1.2) with coefficients satisfying the symmetry conditions (1.3) and consistency, Clearly, S h is a reversible map and moreover, it is consistent with e ihH at least at order 1 , so that (3.1) holds with p ≥ 1 .Since H is real symmetric, it is diagonalizable.Therefore, if the eigenvalues of H are simple, the dynamics of (S n h ) n∈Z is given by Theorem 3.1: for sufficiently small h , there exist real matrices D h (diagonal) and P h (invertible) so that S n h = P h e inD h P −1 h , all the eigenvalues of S h verify |ω j | = 1 and M(u) and H(u) are almost preserved for long times.This corresponds to the examples of Figure 1.The same conclusions apply as long as H is a real matrix with all its eigenvalues real and simple (Figure 2, left), whereas the general case of complex eigenvalues is not covered by the theorem, and no preservation is ensured (Figure 2, right).
Suppose now that the real matrix H has multiple real eigenvalues, but is still diagonalizable, and that A and B are real and symmetric.In that case, a symmetric-conjugate splitting method satisfy both conditions (1.4) and (1.5), so that it can be written as where H h is a family of real matrices whose even terms in h are symmetric and odd terms are skewsymmetric.Suppose in addition that S h is of even order (i.e., p is even in (3.6)).In that case the matrix R in Theorem 3.3 is symmetric, and so its eigenvalues are real.Moreover, since R strongly depends on the coefficients a j , b j and the decomposition H = A + B , it is very likely that typically the eigenvalues of the operators Π ω R |Eω(H) are simple and so that the dynamics of (S n h ) n∈Z is given by Theorem 3.3 and is therefore similar to the one of (e inhH ) n∈Z .Notice that this does not necessarily hold if the scheme is of odd order and/or A and B are not symmetric.This phenomenon is clearly illustrated in the examples of Figure 3 by methods S [3,2] h and S [4] h .Notice, however, that method S [3,1] h , although of odd order, works in fact better than expected from the previous considerations.The reason for this behavior resides in the following Proposition 3.8 The 3th-order symmetric-conjugate splitting method , is indeed conjugate to a reversible integrator V h of order 4, i.e., there exists a real near-identity transformation ᾱh S [2] αh , where S [2] h is a time-symmetric 2nd-order method and α = a 1 .Specifically, Therefore, it can be written as for certain real matrices F 2j+1 .In consequence, by applying the BCH formula, one gets ψ h = e W (h) , with Here w 5,j are polynomials in α .Now let us consider for a given parameter λ .Then, clearly, ) and the stated result is obtained, with F h = e λh 3 F 3 .2 This result can be generalized as follows: given a time-symmetric method S [2k] h of order 2k , if α is chosen so that the composition αh is of order 2k + 1 , then ψ h is conjugate to a reversible method of order 2k + 2 .Theorems 3.1 and 3.3 also allow one to explain the good behavior shown by symmetric-conjugate composition methods for this type of problems.In fact, suppose H is a real symmetric matrix and Φ z H is a family of linear maps which are consistent with e izH at least at order 1 and satisfy If we define S h as the symmetric-conjugate composition where α j are some complex coefficients satisfying the symmetry condition

and the consistency condition
then S h is a reversible map.Moreover, it is consistent with e ihH at least at order 1 .Therefore, one can apply Theorem 3.1 and Theorem 3.3 also in this case.Notice, in particular, that even if the maps e iha j A and/or e ihb j B in the symmetric-conjugate splitting method (1.2) are not computed exactly, but only conveniently approximated (for instance, by the midpoint rule), the previous theorems still apply, so that one can expect good long term behavior from the resulting approximation.
4 Symmetric-conjugate splitting methods for the Schrödinger equation An important application of the previous results corresponds to the numerical integration of the time dependent Schrödinger equation ( where ψ : R 3 × R −→ C .The Hamiltonian operator Ĥ is the sum Ĥ = T + V of the kinetic energy operator T and the potential V .Specifically, In addition, a simple computation shows that Assuming d = 1 and periodic boundary conditions, the application of a pseudo-spectral method in space (with N points) leads to the N -dimensional system (1.1),where u(0) = u 0 ∈ C N and H represents the (real symmetric) N × N matrix associated with the operator − Ĥ [16].Now where A is the (minus) differentiation matrix corresponding to the discretization of T (a real and symmetric matrix) and B is the diagonal matrix associated to − V at the grid points.Since exp(tA) can be efficiently computed with the fast Fourier transform (FFT) algorithm, it is a common practice to use splitting methods of the form (1.2) to integrate this problem.In this respect, notice that property (4.2) will be inherited by the matrices A and B only if the number of discretization points N is sufficiently large to achieve spectral accuracy, i.e., Assuming this is satisfied, then there is a reduction in the number of conditions necessary to construct a method (1.2) of a given order p [12,2].Integrators of this class are sometimes called Runge-Kutta-Nyström (RKN) splitting methods [5].Two further points are worth remarking.First, the computational cost of evaluating (1.2) is not significantly increased by incorporating complex coefficients into the scheme, since one has to use complex arithmetic anyway.Second, since j a j = 1 for a consistent method, if a j ∈ C , then both positive and negative imaginary parts are present, and this can lead to severe instabilities due to the unboundedness of the Laplace operator [8,14].On the other hand, the spurious effects introduced by complex b j can be eliminated (at least for sufficiently small values of h ) by introducing an artificial cut-off bound in the potential when necessary.
In view of these considerations, we next limit our exploration to symmetric-conjugate splitting methods of the form (1.2) with 0 < a j < 1 and b j ∈ C with ℜ(b j ) > 0 to try to reduce the size of the error terms appearing in the asymptotic expansion of the modified Hamiltonian H h associated with the integrator.
For simplicity, we denote the symmetric-conjugate splitting schemes S h by their sequence of coefficients as As a matter of fact, since A and B are sought to verify (4.3), sequences starting with B may lead to schemes with a different efficiency, so that we also analyze methods of the form include integrators where the central exponential corresponds to A (when b r = 0 ) and B (when a r = 0 ), respectively.The method has s stages if the number of exponentials of A is precisely s for the scheme (4.5) or s + 1 for the scheme (4.4).
The construction process of methods within this class is detailed elsewhere (e.g.[7,5] and references therein), so that it is only summarized here.First, we get the order conditions a symmetric-conjugate scheme has to satisfy to achieve a given order p = 4, 5 and 6.These are polynomial equations depending on the coefficients a j , b j , and can be obtained by identifying a basis in the Lie algebra generated by {A, B} and using repeatedly the BCH formula to express the splitting method as S h = exp(hH h ) , with H h in terms of A , B and their nested commutators.The order conditions up to order p are obtained by requiring that H h = H + O(h) p+1 , and the number is 7, 11 and 16 for orders 4, 5 and 6, respectively.
Second, we take compositions (4.4) and (4.5) involving the minimum number of stages required to solve the order conditions and get eventually all possible solutions with the appropriate symmetry.Sometimes, one has to add parameters, because there are no such solutions.In particular, there are no 4th-order schemes with 4 stages with both a j > 0 and ℜ(b j ) > 0 .
Even when there are appropriate solutions, it may be convenient to explore compositions with additional stages to have free parameters for optimization.This strategy usually pays off when purely real coefficients are involved, and so it is worth to be explored also in this context.Of course, some optimization criterion related with the error terms and the computational effort has to be adopted.In our study we look at the error terms in the expansion of H h at successive orders and the size of the b j coefficients.Specifically, we compute for each method of order, say, p , the quantities Here s is the number of stages and E r+1 is the Euclidean norm of the vector of error coefficients in H h at higher orders than the method itself.In particular, for a method of order 6 , E for this method we get an idea of how the higher order error terms behave.It will be of interest, of course, to reduce these quantities as much as possible to get efficient schemes.
Solving the polynomial equations required to construct splitting methods with additional stages is not a trivial task, especially for orders 5 and 6.In these cases we have used the Python function fsolve of the SciPy library, with a large number of initial points in the space of parameters to start the procedure.From the total number of valid solutions thus obtained, we have selected those leading to reasonably small values of all quantities (4.6) and checked them on numerical examples.
The corresponding values for the most efficient methods we have found by following this approach have been collected in Table 1, where N A 1.000 1.267 0.400 0.821 0.704 1.082 1.012 1.000 1.141 0.352 0.698 0.559 0.913 0.789 1.000 1.416 0.322 0.766 0.666 1.025 0.866 1.000 These conditions restrict the coefficients a j , b j in (4.8) to be positive, however, and thus the method to be of second order at most.Nevertheless, it has been shown in [14,8] that, if in addition L , Â and B generate analytic semigroups on X defined in the sector Σ ϕ = {z ∈ C : | arg z| < ϕ} , for a given angle ϕ ∈ (0, π/2] and the operators Â and B verify ∥e z Â∥ ≤ e ω|z| , ∥e z B ∥ ≤ e ω|z| for some ω ≥ 0 and all z ∈ Σ ϕ , then a splitting method of the form (4.8) of classical order p with all its where C is a constant independent of n and h .

Numerical illustration: Modified Pöschl-Teller potential
The so-called modified Pöschl-Teller potential takes the form with λ > 1 , and admits an analytic treatment to compute explicitly the eigenvalues for negative energies [9].
For the simulations we take α = 1 , λ(λ − 1) = 10 and the initial condition ψ 0 (x) = σ e −x We first check how the errors in the norm M(u) and in the energy H(u) evolve with time according with each type of integrator.To this end we integrate numerically until the final time t f = 10 4 with three 6th-order compositions involving complex coefficients: (i) the new symmetric-conjugate scheme N B * [6] 11 collected in Table 2 (h = 100/909 ≈ 0.11) , (ii) the palindromic scheme denoted by B [6] 16 with all a j taking the same value a j = 1/16 , j = 1, . . ., 8 and complex b j with positive real part3 (h = 0.16) , and (iii) the method obtained by composing B [6] 16 with its complex conjugate (B [6] 16 ) * , resulting in a symmetricconjugate integrator (h = 0.32) .The step size is chosen in such a way that all the methods require the same number of FFTs.The results are depicted in Figure 4. We see that, according with the previous analysis, the error in both unitarity and energy furnished by the new scheme N B * [6] 11 does not grow with time, in contrast with palindromic compositions involving complex coefficients.Notice also that the composition of the palindromic scheme B [6] 16 with its complex conjugate leads to a new (symmetric-conjugate) integrator with good preservation properties.On the other hand, composing a symmetric-conjugate method with its complex conjugate results in a palindromic scheme showing a drift in the error of both the norm and the energy [4].
In our second experiment, we test the efficiency of the different schemes.To this end we integrate until the final time t f = 100 , compute the expectation value of the energy, H(u app (t)) , and measure the error as the maximum of the difference with respect to the exact value along the integration: (5. 2) The corresponding results are displayed as a function of the computational cost measured by the number of FFTs necessary to carry out the calculations (in log-log plots) in Figure 5. Notice how the new symmetricconjugate schemes offer a better efficiency than standard splitting methods for this problem.The improvement is particularly significant in the 6th-order case.

Figure 2 :
Figure 2: Same as Figure 1 when H = A + B is a real (but not symmetric) matrix.Left: the eigenvalues of H are real and simple.Right: the eigenvalues of H are arbitrary.

Figure 3 :
Figure 3: Same as Figure 1 when H = A + B is a real symmetric matrix with multiple eigenvalues.Left: A and B are real symmetric matrices.Right: A and B are real, but not symmetric.
fgives an estimate of the efficiency of the scheme by considering only the error at order 7.By computing E

Figure 4 :
Figure 4: Error in norm M(u) (left) and in energy H(u) (right) as a function of time for complex-conjugate and palindromic methods involving complex coefficients.
, namely: H is symmetric, with A and B (blue dashed line) for different values of h when H = A + B is a real symmetric matrix with simple eigenvalues.Left: A and B are also real symmetric.Right: A and B are real, but not symmetric. h Lemma 3.6Let M be a complex matrix.Then M is diagonalizable if and only if the kernel and the image of ad We can assume, in virtue of Lemma 3.5 and without loss of generality, that M is in Jordan normal form2.On the one hand, if M is diagonal, we have ad M A = ((m i,i − m j,j )A) i,j and so the support of the matrices in Ker C ad M and Im C ad M are clearly disjoint (which implies (3.8)).Conversely, doing calculations by blocks it is enough to consider the case where M = λI N + N is a Jordan matrix (i.e.λ ∈ C and N nilpotent).Then we just have to note that ad λI N +N = ad N and that since ad N is nilpotent necessarily we have Ker C ad N ∩ Im C ad N ̸ = {0} . 2 Lemma 3.7 Let M h be a family of real matrices depending smoothly on h and of the form M are supplementary, i.e.Ker C ad M ∩ Im C ad M = {0}.(3.8)Proof: Smoothness property: For any pair of multi-indices (i 1 , ..., i m ) and (j 1 , ..., j m ) with i1 + • • • + i m + j 1 + • • • + j m = p + 1 ,and for all t ∈ [0, T ] , ∥ Âi 1 Bj 1 . . .Âim Bjm e t Lu 0 ∥ ≤ C for a positive constant C .

Table 2 :
Coefficients of the most efficient symmetric-conjugate RKN splitting methods of order 4, 5 and 6.
coefficients a j , b j in the sector Σ 2 /2 , with σ a normalizing constant.We discretize the interval x ∈ [−8, 8] with N = 256 equispaced points and apply Fourier spectral methods.With this value of N it turns out that ∥([B,[B, [A, B]]])u 0 ∥ is sufficiently close to zero to be negligible, so that we can safely apply the schemes of Table2.If N is not sufficiently large, then the corresponding matrices A and B do not satisfy (4.3), and as a consequence, the schemes are only of order three.This can be indeed observed in practice.

Table 4 :
Coefficients of the most efficient splitting methods collected in Table3.