High-order commutator-free quasi-Magnus exponential integrators for non-autonomous linear evolution equations

The class of commutator-free quasi-Magnus (CFQM) exponential integrators provides a favourable alternative to standard Magnus integrators, in particular for large-scale applications arising in the time integration of non-autonomous linear evolution equations. The schemes are given by compositions of several exponentials that comprise certain linear combinations of the values of the deﬁning operator at speciﬁed nodes. Due to the fact that previously proposed CFQM exponential integrators of order ﬁve or higher involve negative coeﬃcients in the linear combinations, severe instabilities are observed for spatially semi-discretised parabolic equations or for master equations describing dissipative quantum systems. In order to remedy this issue, two diﬀerent approaches for the design of eﬃcient time integrators of orders four, ﬁve, and six are pursued: (i) the study of CFQM exponential integrators involving complex coeﬃcients that satisfy a positivity condition, and (ii) the study of unconventional methods in the sense that an additional exponential involving a commutator of higher order with respect to the time stepsize occurs. Numerical experiments conﬁrm that the identiﬁed novel time integrators are superior to other integrators of the same family previously proposed in the literature.


Introduction
Simulation of quantum systems. In order to describe the evolution of a quantum system subjected to a time-dependent potential, one must deal with the time-dependent Schrödinger equation i ∂ t ψ(t) = H(t) ψ(t) ; (1) here, H denotes the time-dependent Hamiltonian operator, and the wave function ψ(t) represents the state of the system at time t. For convenience, we employ a compact formulation as evolution equation and the normalisation = 1.
Dynamical processes involving dissipation are better described by a master equation defined by a time-dependent Liouville operator L , see [19]; this in particular applies to quantum systems with dissipation caused by cavity losses and spontaneous emission of energy input from external pumping [3,30] and also when a complex absorbing potential is present [36].
Due to the fact that analytical solutions to time-dependent partial differential equations such as Schrödinger equations and master equations are available only in academic examples, numerical methods are inevitable tools in their simulation. A common approach is to first discretise in space, for which procedures such as finite differences, pseudo-spectral methods based on collocation with trigonometric polynomials, or Galerkin methods with a Hermite basis are used, see [31] and references therein. The (large) system matrix of the resulting linear ordinary differential equation inherites fundamental properties of the Hamiltonian and the Liouville operator, respectively. In the first case, it is skew-Hermitian reflecting (high) oscillations in the solution; in the second case, it comprises eigenvalues with large negative real part reflecting dissipation.
Evolution equation. In the simulation of (dissipative) quantum systems, one is thus led to numerically integrate a non-autonomous evolution equation of the form defined by a family of time-dependent linear operators (A(t)) t∈[t 0 ,T ] . In view of a spatial semi-discretisation, we assume in addition that A(t) is a family of real or complex matrices which have large dimension and large operator norm, in general.
In addition to (1) and (2) leading in a natural way to a non-autonomous linear evolution equation (3), we mention the variational equations related to periodic solutions in nonlinear parabolic equations, see [39].
Exponential time integrators. Compared to standard time integration methods such as explicit or implicit Runge-Kutta methods, exponential integrators have shown to be particularly useful in approximating the exact solution of (oscillatory) stiff differential equations in an efficient manner whilst preserving qualitative properties of the system, see [13] and references therein. As a matter of fact, there is a considerable research interest in the study of different methods involving exponentials and the design of novel schemes; we mention time-splitting methods, Magnus and Fer integrators, Crouch-Grossman methods, and commutator-free exponential integrators, see [6,10,13,23,26,27,33]. Two relevant classes, time-splitting methods and interpolatory Magnus integrators, will be reviewed in Section 2.
Exponential midpoint rule. An elementary time integrator of exponential type requiring the computation of one single matrix-exponential per time step is the so-called exponential midpoint rule. For certain equidistant grid points the time-discrete solution of (3) is defined by the recurrence u n+1 = CF [2] 1 (τ, t n ) u n = e τ A(tn+ 1 2 τ ) u n ≈ u(t n+1 ) , n ∈ {0, 1, . . . , N − 1} .
The exponential midpoint rule provides an approximation that has non-stiff order two in the time stepsize. When the arising matrices are of large dimension, the scheme is frequently implemented by means of an efficient algorithm that approximates the action of a matrix-exponential on a vector, instead of computing the matrix-exponential itself, see [1,26,35,38,41] and references therein; in this respect, procedures that involve only matrix-vector products such as Chebyshev or Krylov methods are usually preferred, see [25,31].
Higher order methods can be achieved by applying different strategies, in particular by constructing the scheme as a single exponential (Magnus integrators) or as a product of exponentials (Fer integrators, splitting methods, commutator-free methods, etc.). In any case, the presence of commutators and/or negative coefficients in the methods lead to serious issues concerning well-definedness and stability. The reason is easy to grasp making use of backward error analysis [23]. Suppose we have a numerical map provided by a given integrator of order, say, r, u 1 = Φ h u 0 . Then u 1 =ũ(h) = u(h) + O(h r+1 ) andũ(h) is the exact solution of the modified equation where A = A 1 (t) + hA 2 (t) + h 2 A 3 (t) + . . . is an asymptotic series in general. Whereas low order exponential methods usually show excellent stability and convergence properties, for higher order methods the corresponding modified equation is not well conditioned, the norm of the perturbation term A takes exceedingly large values and thus only very small time steps provide meaningful results.
One of the goals of the present paper is indeed to develop high-order integrators avoiding these drawbacks. This can be achieved, in particular, by considering commutator-free methods of exponential type involving complex coefficients. Splitting methods with complex coefficients have been proposed for Hamiltonian systems, see [21], in the context of the time-dependent Schrödinger equation in quantum mechanics, see [5,7,8], and also in the more abstract setting of evolution equations with unbounded operators generating analytic semigroups, see [11,20,24].
Commutator-free quasi-Magnus exponential integrators. To be more specific, in this work we primarily focus on a class of high-order commutator-free integrators which can be considered as a generalisation of (4): they are given by a composition of several exponentials that comprise certain linear combinations of the values of the defining operator at specified quadrature nodes. We call them commutator-free quasi-Magnus (CFQM) exponential integrators to distinguish themselves from other classes of commutator-free integrators and to emphasise the close relationship with interpolatory Magnus methods, in the sense that both make use of the Magnus expansion in their derivation.
For our purposes, it suffices to study the initial step of a commutator-free quasi-Magnus exponential integrator of order r for notational simplicity, here and in the sequel, we only indicate the dependence on the decisive quantity τ > 0.
Integrators of this class with real coefficients have been designed and analysed in [18,43]; we in particular mention [2], where optimised methods of orders four, six, and eight have been proposed and tested.
Evidently, for the special case of a time-independent operator A, the composition (5a) reads CF J (τ ) = e τ a J A · · · e τ a 1 A , where As shown in [16], the positivity condition ensures that a CFQM exponential integrator is well-defined within the framework of sectorial operators and analytic semi-groups.
Depending on the prescribed order, condition (5f) can be satisfied by a commutator-free quasi-Magnus exponential integrator with real coefficients or demands the consideration of complex coefficients. For instance, for order four, there exist schemes comprising two exponentials that satisfy property (5f). On the contrary, it has been noticed in [3] that for all schemes with real coefficients of order six and higher the positivity condition (5f) is violated for at least one index; numerical experiments confirm that these schemes exhibit poor stability when applied to spatially semi-discretised partial differential equations of parabolic type or to equations modelling driven open quantum systems, respectively. Since we are interested in designing high-order schemes that are also suitable for the time integration of parabolic equations and for master equations describing dissipative quantum systems, we explore the more general approach of complex coefficients.
In the present work, we trace back this unfavourable feature to the non-stiff order conditions and show that, at least if the CFQM exponential integrator does not involve an exceedingly large number of stages, this property is characteristic for schemes with real coefficients of order five or higher. To remedy the issue, we explore the following approaches: The design of novel commutator-free quasi-Magnus exponential integrators that involve more stages than strictly required to satisfy as many conditions at order five as possible whereas keeping positive coefficients, and the design of novel fifth-and sixth-order CFQM exponential integrators with complex coefficients satisfying the positivity condition. As confirmed by numerical experiments, the proposed schemes can be safely applied to problems related to parabolic equations and to dissipative quantum systems. Moreover, we study unconventional methods that comprise in addition a single exponential involving a commutator which is of higher order with respect to the time stepsize.

Convergence analysis.
A rigorous stability and error analysis of high-order commutator-free quasi-Magnus exponential integrators with complex coefficients applied to non-autonomous linear evolution equations of parabolic type is carried out in [16]; the stated convergence estimate in particular implies that a CFQM exponential integrator retains its non-stiff order under suitable regularity and compatibility requirements on the exact solution.
Extension. The time integration of nonlinear evolution equations of the form is often advantageously based on the previous accurate numerical solution of the nonautonomous linear equation (3); the equation associated to the perturbation, u (t) = ε g t, u(t) , is solved separately with the time frozen. The techniques provided in this work and the obtained novel schemes are thus of use in this setting.
Outline and overview. The present work is organised as follows. In Section 2, we introduce two classes of numerical methods based on exponentials, to put commutatorfree quasi-Magnus exponential integrators into perspective with them. In Section 3, we in particular state a set of independent order conditions to be satisfied by time-symmetric CFQM exponential integrators of order six. In Section 4, we solve the order conditions and collect several families of commutator-free quasi-Magnus exponential integrators and related methods. Finally, in Section 5, efficient novel schemes are tested on model equations.
An overview of the different time-splitting methods, interpolatory Magnus integrators, and commutator-free quasi-Magnus exponential integrators that are discussed in this work is given in Table 1. In particular, for sixth-order methods, it is seen that CFQM exponential integrators involve a significantly reduced number of exponentials, which is a first indicator of the computational cost. In the design of novel schemes, an additional objective is to keep the quantity ρ defined in (18) as small as possible.

Splitting methods and Magnus integrators
In this section, we briefly review splitting methods and interpolatory Magnus integrators for the time integration of non-autonomous linear evolution equations, paying special attention to stability issues and computational effort. In conclusion, we relate commutator-free quasi-Magnus exponential integrators to them.
Estimating the computational cost. In view of the design of novel time integration methods, it is essential to compare efficiency. To this end, we introduce the computational cost of evaluating the action of a matrix-exponential of the form e a 0 A(t) with real or complex coefficient a 0 on a vector u by polynomial approximations such as Krylov-type methods, as a basic measure.
Working with complex coefficients instead of real coefficients potentially makes an algorithm about four times more costly; however, this is not necessarily what happens in our

Coefficients
Order Indicator for computational cost S (see [17]) real r = 4 J = 6 S (see [20]) complex r = 4 J = 4 S (see [17]) real r = 6 J = 10 S (see [11]) complex r = 6 J = 16 M (see (9)) real r = 4 C e Ω [4] (τ ) u 0 5 C e A(t) u 0 M real r = 6 C e Ω [6] (τ ) u 0 13 C e A(t) u 0 CFQM (see (12)) real r = 4 = 2 s ρ(CF [4] 2 ) = 1 CFQM (see [2,Eq. (43) 5C , additional commutator Table 1 Overview of various higher-order time-splitting methods (S), Magnus integrators (M), and commutator-free quasi-Magnus exponential integrators (CFQM). The number of Gaussian quadrature nodes s corresponds to the minimum number of nodes. Indicators for the computational cost are the number of exponentials J and the quantity (18). For Magnus integrators, the cost comparable to the computation of e A(t) u 0 is given. We recall that the abbreviation CF situation. For instance, in applications from quantum mechanics, where A ∈ C d×d and the solution of (3) evolves in the complex space, the evaluation of e a 0 A(t) u by a Krylov method requires the computation of certain products (A(t)) m u and the computational cost is comparable for real and complex coefficients; if A(t) ∈ R d×d and u ∈ R d , the extra cost to compute U = e a 0 A(t) u with a 0 , ∈ C corresponds to the computation of the exponential of a complex matrix of relatively small dimension. However, if one has next to compute, say, U = e a 0 A(t) U , the products (A(t)) m U duplicate their computational cost, since U is a complex vector.
Due to the fact that a matrix-vector product of the form (A 1 + A 2 ) u commonly has the same complexity as either A 1 u or A 2 u, first considerations given below for time-splitting methods essentially extend to commutator-free quasi-Magnus exponential integrators.

Time-splitting methods
Splitting methods. Given two non-commuting square matrices X and Y , a timesplitting method of order r for approximating the solution of the linear differential equation or equivalently the corresponding matrix-exponential is given by a composition of the form Simply by adding the trivial relation d dt t = 1, this splitting method applied to (3) yields Efficient splitting methods are found in [17,33]; for instance, the coefficients of a fourthorder scheme comprising six exponentials and of a sixth-order scheme comprising ten exponentials are given there.

Stability.
A well-established result states that time-splitting methods of order higher than two necessarily involve negative coefficients and hence are not suitable for non-reversible systems [9,22,40,42]; in particular, such schemes have poor stability properties for nonautonomous linear evolution equations of parabolic type and for dissipative quantum systems.
An alternative approach is to employ time-splitting methods that are defined by complex coefficients a ∈ C + and real coefficients b ∈ R for ∈ {1, . . . , L}; we note that this ensures t 0 +c τ ∈ R for ∈ {1, . . . , L}. A fourth-order scheme comprising four exponentials is given in [20], and a sixth-order scheme comprising 16 exponentials is found in [11].
Computational effort. In situations where the defining operator family and the timediscrete solution are complex-valued, the presence of complex coefficients in the composition (7) does not lead to an increment in the computational cost, compared to splitting methods with real coefficients; however, if A(t) ∈ R d×d and u 0 ∈ R d , the computational effort of schemes with complex coefficients will increase, in general. In this situation, the computation from time t 0 to t 1 proceeds as and the projection on the real part defines the new time-discrete solution value Re(U L ) ≈ u(t 1 ) ; the computational effort for evaluating U 1 by means of polynomial approximations such as Krylov methods is basically the same for a 1 ∈ C + and a 1 ∈ R, whereas the effort for computing the remaining vectors doubles for complex coefficients.

Interpolatory Magnus integrators
Magnus expansion. The Magnus expansion provides a formal solution representation of (3) as the exponential of an infinite series which involves multiple integrals of nested matrix-commutators such as here, as standard, the commutator of two square matrices is defined by Recursive procedures to generate the Magnus expansion are found in [13,27,32].
Magnus integrators. The natural approach of truncating the infinite series and applying quadrature approximations of the integrals leads to the class of interpolatory Magnus integrators. For instance, from considering the two principal terms in (8), Ω 1 + Ω 2 , and applying a quadrature formula of order four or higher, defined by K nodes and corresponding weights, a fourth-order Magnus integrator results [14,28] in comparison with (5); the particular choice of the fourth-order Gauss-Legendre quadrature formula yields In the context of partial differential equations, however, higher-order interpolatory Magnus integrators exhibit serious difficulties, due to the presence of nested commutators.
Well-definedness and stability. A first issue is that commutators in general do not inherit the properties of the operator family (A(t)) t∈[t 0 ,T ] ; as a consequence, characterising well-definedness and stability of higher-order interpolatory Magnus integrators for partial differential equations is a delicate question.
Computational effort. A second issue is the potential exceedingly large computational cost of higher-order interpolatory Magnus integrators. Suppose for instance that the fourth-order Magnus integrator (9) is implemented by Krylov techniques. The computation of e Ω [4] (τ ) u 0 requires to evaluate the action of certain powers of Ω [4] (τ ) at the cost of five matrix-vector multiplications per product; specifically, the relation in general implies the computational cost C e Ω [4] (τ ) u 0 5 C e A(t) u 0 .
In some cases, however, the evaluation of the commutators acting on the vector can be efficiently carried using an appropriate approach, e.g. the linear Schrödinger equation in the semiclassical regime with time-dependent potential [6].
For a sixth-order Magnus integrator at least 13 matrix-vector multiplications are needed [15]; consequently, computing the action of the matrix-exponential is about 13 times more expensive than the cost of evaluating a single matrix-exponential C e Ω [6] (τ ) u 0 13 C e A(t) u 0 .

Commutator-free quasi-Magnus exponential integrators
Favourable conjunction of method classes. The class of commutator-free quasi-Magnus exponential integrators combines the favourable properties of time-splitting meth-ods and standard interpolatory Magnus integrators. On the one hand, CFQM exponential integrators can be considered as generalisations of time-splitting methods (7), but with a significantly reduced number of exponentials, see also Sections 3 and 4. On the other hand, they circumvent difficulties exhibited by higher-order interpolatory Magnus integrators, especially the presence of the often cumbersome (from the stability point of view) and in many cases computationally costly nested commutators; the single exponential involving commutators is replaced by a composition of several exponentials that comprise linear combinations of the values of the defining operator at certain quadrature nodes, which leads to the general format (5). As for interpolatory Magnus integrators, the coefficients of a commutator-free quasi-Magnus exponential integrator can be easily recalculated for each particular choice of the nodes.

Derivation and analysis of order conditions
In this section, we deduce the non-stiff order conditions for CFQM exponential integrators; we thus may restrict ourselves to the case where the operator norm of A(t) ∈ K d×d is of moderate size, see (3). We first indicate the general approach; as illustration, the expansions are rendered more precisely for the special case of a fourth-order CFQM exponential integrator involving two exponentials. Following up these considerations, and employing additional algebraic means, we attain a set of independent order conditions that is suitable for a detailed analysis of its solutions.

Expansions
The requirement that a commutator-free quasi-Magnus exponential integrator of the form (5a)-(5c) satisfies relation (5d) for a given r, implies certain conditions on the coefficients. Among a number of ways, suitable for ordinary differential equations, a possible approach for the derivation of these conditions consists in representing the exact solution by the Magnus expansion (8) and reproducing it by the time-discrete solution up to the prescribed order CF For simplicity, in the presentation, we take t 0 = 0. Given (w i , c i ), i = 1, . . . , K, the weights and nodes of a quadrature rule of order r, it is shown in [28] that the evaluations A(c i τ ), i = 1, . . . , K suffice to get an approximation of the solution of (3) up to order O τ r (see also [34,37,44]). This is so because the solution of the differential equation where A(t) is the polynomial of degree K − 1 in t that interpolates A(t) on the interval [0, τ ] at the nodes c i τ, i = 1, . . . , K, i.e.
On the other hand, the Magnus expansion can be used to solve the initial value problem (10), so that Now, since A(t) is a polynomial of degree K − 1 in t, it is clear that all the terms in the series (11) can be evaluated analytically. An approximation to u(τ ) up to order p can then be achieved by taking into account the first terms of the series, since it holds that Let us take in particular the K nodes of a Gauss-Legendre quadrature rule, so that p = 2K. Since these nodes are symmetrically distributed with respect to the midpoint, it is advantageous to consider a Taylor expansion of A(t) around τ 2 , , etc., and Ω(τ ) can be written as an element of the free Lie algebra generated by α 1 , . . . , α K . Since the exact solution is time-symmetric and α i are derivatives with respect to the midpoint, no even terms will appear in the expansion [15].
Specifically, with r = 2K = 4 (order four) one has The approximation e Ω 1 (τ )+ Ω 2 (τ ) can also be reproduced up to order O(τ 5 ) by a product of two exponentials. Specifically, Finally, expressing α 1 , α 2 in terms of the evaluations A 1 , A 2 , results in the fourth-order CFQM exponential integrator for the benefit of a compact formulation, the coefficients are given in matrix-form. We point out that although negative coefficients necessarily arise, the positivity condition (5f) is satisfied, since Additional considerations (r = 2K = 6). For order six we consider the nodes and weights of the sixth-order Gauss-Legendre quadrature formula, and The exact solution can be approximated by a composition of the form with r = 6 and appropriately chosen coefficients in the linear combinations. Once such coefficients are determined by solving certain order conditions, see Section 3.2, this product defines a CFQM exponential integrator (5) with K = 3 and besides, the positivity condition (5f), where α 1 = τ A, α 2 = α 3 = 0, corresponds to the requirement ∀ j ∈ {1, . . . , J} :

Order conditions
Approach based on BCH formula (r = 4). Different approaches can be used to determine the coefficients that arise in the factorization (14). For instance, we can apply recursively the Baker-Campbell-Hausdorff (BCH) formula to get a series expansion with respect to the time stepsize τ ; for (12), this yields log Ψ [r] where w 1 , w 2 , w 3 , w 12 , w 13 are polynomials in terms of x j1 , x j2 , x j3 . The order conditions are then obtained by requiring that this expansion agrees with the corresponding expansion of Ω(τ ) up to the desired order, that is Order conditions and time-symmetry (r = 6). A more elegant approach to deduce a set of necessary and sufficient independent order conditions for splitting methods has been proposed in [12] (see also [4]) and makes use of Lyndon words. Roughly speaking, one computes the series expansion of the product of exponentials defining Ψ [r] J (τ ), but only takes into account those terms in the expansion corresponding to a Lyndon word. Then these terms are equating to those obtained by applying the same procedure to e Ω(τ ) . The order conditions thus derived are shown to be equivalent to those given by applying the BCH formula.
More specifically, let us consider the set of all multi-indices comprising positive integers with lexicographical order <; a tuple (i 1 , . . . , i m ) ∈ N m is called a Lyndon multi-index if for each integer ∈ {1, . . . , m} the condition (i 1 , . . . , i ) < (i +1 , . . . , i m ) is satisfied, and for each Lyndon multi-index there is a Lyndon word on the alphabet {α 1 , α 2 , . . . , α m }. In view of (14), anticipating that the additional requirement of time-symmetry will raise the order to six, we list all Lyndon multi-indices with cross sum at most five the word α 1 α 3 corresponds to (1, 3), the word α 1 α 2 α 2 corresponds to (1, 2, 2), and so on.
Using that the composition in (14) agrees up to order five or higher with the power series corresponding to the Magnus expansion (8) if and only if the coefficients multiplying the Lyndon words on the alphabet {α 1 , α 2 , α 3 } are equal in both expansions, gives us ten independent polynomial equations to be satisfied by the coefficients. Setting explicit expressions of these order conditions are (1) : : x j3 x j1 + 2 (1 − y j ) = 1 12 , x j3 x j2 − 2 z j = 1 120 , (1, 2, 2) : When dealing with a non-autonomous linear evolution equation (3), time-symmetric meth-ods (14) satisfying the condition are usually preferred. As shown in [18], this corresponds to the requirement due to the fact that the order conditions at odd orders are then automatically satisfied, a time-symmetric scheme is of even order. In particular, to achieve order six, the seven equations in (17b) associated with with x j3 = 0 for j ∈ {1, . . . , J} need to be solved; setting J = 2 and employing the requirement of time-symmetry, we retain (12)

Families of novel commutator-free quasi-Magnus exponential integrators
In this section, we design efficient commutator-free quasi-Magnus exponential integrators of orders four, five, and six and related methods of order six.
Efficient novel schemes. From the analysis provided in Section 3, we conclude that for CFQM exponential integrators, cast into the form (14) with real coefficients and satisfying the positivity condition x j1 > 0 for j ∈ {1, . . . , J}, there is an order barrier at order four, at least from a practical point of view. This is also confirmed by [2], where an exhaustive study of fourth-, sixth-and eighth-order time-symmetric CFQM exponential integrators with up to three, six and eleven exponentials, respectively, was carried out; when more exponentials than strictly necessary were included in the composition, an optimisation criterion was adopted. Any of the sixth-and eighth-order schemes considered in [2] contain at least one index j ∈ {1, . . . , J} such that x j1 < 0, and hence their applicability to nonreversible systems is problematic. Our objective is to design CFQM exponential integrators which do not present this limitation; in doing so, we eventually consider schemes including compositions with complex coefficients. Moreover, we study unconventional methods comprising an additional exponential that involves a commutator which is of higher order with respect to the time stepsize.
Short notation. For the sake of a compact formulation, instead of see (5a)-(5c) and (14), we henceforth write for short. We recall that the coefficients of a sixth-order time-symmetric CFQM exponential integrator are determined by seven conditions associated with the multi-indices  (17); we collect these coefficients in matrix-form.
In order to measure the computational effort needed to achieve a given accuracy 1 , we introduce, as a reference, the quantity ρ CF Thus, in particular, for the method (12) one has ρ(CF [4] 2 ) = 1.
Notice that ρ(CF

[r]
J ) is independent of the quadrature nodes used.

Fourth-order schemes with real coefficients
First, we design favourable fourth-order time-symmetric commutator-free quasi-Magnus exponential integrators with real coefficients satisfying the positivity condition x j1 > 0 for any j ∈ {1, . . . , J}, see (14) and (17c). We use the extra degrees of freedom due to the inclusion of sixth-order Gaussian quadrature nodes and additional exponentials to verify certain conditions at order five and to minimise the deviation of the remaining fifth-order conditions without increasing the overall computational cost.
Scheme (r = 4, s = 3, J = 3). For three exponentials, we have not found fourth-order time-symmetric CFQM exponential integrators that are substantially more efficient than the optimised scheme proposed in [2, Eq. (43)]. For this reason, we take this scheme as representative. it turns out that the condition related to (1, 2, 2) is automatically satisfied, and so one ends up with a sixth-order method; this feature was already noticed in [2]. However, as the resulting scheme does not fulfill the positivity condition, we discard it for non-reversible systems. (ii) Alternatively, we solve the order conditions that correspond to in terms of x 11 ; moreover, we minimise the amount by which the remaining order conditions (1, 2, 2), (1, 1, 1, 2) are not satisfied. This yields the solution x 23 by (15) we first solve the order conditions associated with in terms of x 11 and x 31 ; then, we minimise the amount by which the conditions related to (1, 1, 1, 2) at order five and to (1, 1, 1, 1, 1, 2) at order seven are not satisfied. This results in the coefficients

Schemes with complex coefficients
Complex coefficients. In the following, we make use of the fact that the order conditions for higher-order commutator-free quasi-Magnus exponential integrators admit complex solutions such that Re(x j1 ) > 0 for any j ∈ {1, . . . , J}. Irrespective of their favourable stability behaviour for non-reversible systems, schemes with complex coefficients generically involve smaller values and exhibit better accuracy, compared to related schemes with real coefficients; in certain situations, these benefits compensate the additional computational cost. Besides, contrary to the real case, the optimisation strategy of considering additional parameters to design more efficient schemes is usually not needed, since solving the order conditions with the minimum number of degrees of freedom often yields more efficient schemes.
Novel scheme (r = 5, s = 3, J = 3). Remarkably, the ten conditions (17b) for order five can be verified by a composition of three exponentials which only contains nine coefficients. The resulting scheme is not time-symmetric but exhibits a special symmetry the coefficients are given by Novel scheme (r = 6, s = 3, J = 4). As alternative to a sixth-order time-symmetric scheme with real coefficients that violates the positivity condition, since x 21 < 0, we consider the composition defined by eight coefficients, we can choose one free parameter to fulfill the seven order conditions (17d); requiring in addition the condition (1, 1, 1, 1, 1, 2) at order seven to be satisfied, we obtain  Re(a 13 ) Re(a 12 ) Re(a 11 ) Im(a 13 ) Im(a 12 ) Im(a 11 ) which yields ρ CF [6] 5 = 1.29727 . It turns out that for this scheme also the condition at order seven related to (1, 1, 1, 1, 3) is verified.

Sixth-order schemes with one commutator
Unconventional schemes. The reason for the high computational cost of interpolatory Magnus integrators can be traced back to the large number of matrix-vector products that are required to approximate e Ω [r] u, since Ω [r] (τ ) = O(τ ). However, if the considered matrix has a special form say, D = [B 1 , B 2 ], such that the evaluation of D u requires only four matrix-vector products and hence the action of the matrix-exponential can be approximated to a certain accuracy with 4 m products thus, in practice, the computational effort of e D u is comparable to that of a single matrixexponential C e D u ∼ C e τ A(t) u .
This observation opens the door to consider yet another unconventional possibility to design sixth-order schemes satisfying the positivity condition by including a single commutator of the form in the composition (14).
Novel scheme involving commutator (r = 6, s = 3, J = 5). Based on our considerations, we study a time-symmetric composition of the form CF [6] involving eight coefficients, required to fulfill the seven order conditions (17d); the free parameter is chosen such that the condition related with (1, 1, 1, 1, 1, 2) at order seven is verified. The uniquely determined solution is given by

Numerical experiments
In this section, we present numerical results for two representative examples, a diffusionadvection-reaction equation and a high-dimensional (dissipative) Rosen-Zener model. In particular, we illustrate the behaviour of the novel schemes CF [4] 4 , CF [4] 5 , CF [5] 3 , CF [6] 4 , CF [6] 5 , CF [6] 5C , compared to the optimised fourth-and sixth-order commutator-free quasi-Magnus exponential integrators proposed in [2,Eq. (43), We note that the sixth-order scheme does not satisfy the positivity condition (5f) resulting in poor stability properties for evolution equations of parabolic type and dissipative quantum systems; besides, as the value of ρ(CF [6] 6 ) is considerably larger than the corresponding values for the novel schemes, the computational effort for the evaluation of the matrix-exponential will in general be higher.

Test equation of parabolic type
Variational equation. Following [16], we consider a nonlinear diffusion-advectionreaction equation, prescribing a smooth solution U : The associated variational equation is of the form (3) We in particular consider the initial state and impose periodic boundary conditions; the solution profile is displayed in Figure 1.
Numerical results. For the space discretisation, we choose M = 100 equidistant grid points and apply standard symmetric finite differences. The time integration is performed by the schemes introduced in the previous section with time stepsizes 2 − for ∈ {1, . . . , 10}; the numerical approximation obtained for the finest time stepsize serves as reference solution. The global errors at time T = 1 versus the time stepsizes are displayed in Figure 2; in addition, as a rough estimate of the computational cost, we include the errors versus the numbers of exponentials. In accordance with the convergence result given in [16], the nonstiff orders are retained for smaller time stepsizes. The numerical comparisons also confirm the poor stability behaviour of the sixth-order scheme CF    Effect of space discretisation. For completeness, we include the results obtained for M = 150, see Figure 3. Here, the instabilities are more severe; otherwise, we observe that the spatial error is negligible and that the numerical results are independent of the considered space discretisation method.

High-dimensional Rosen-Zener model with dissipation
Non-autonomous linear system. As a further application, we consider a quantum system that is closely related to the model analysed in [29]. For a Rosen-Zener model with dissipation, the associated Schrödinger equation in the interaction picture has the following form see (3). After normalisation, the time-dependent Hamiltonian reads with identity matrix I ∈ R k×k , Pauli matrices and dissipation parameter δ > 0 with Evidently, the components of the matrix D increase in modulus for larger dimensions. We in particular set   [4] CF 4 [4] CF 5 [4] V 0 =2, δ=0 Time integration. In principle, commutator-free quasi-Magnus exponential integrators (5a)-(5c) that do not satisfy the positivity condition (5f) are suitable for the time integration of a Rosen-Zener model with dissipation (21). However, if the dimension of the system is large, the dissipative term forces the time stepsize to be sufficiently small, in order to avoid instabilities; in this situation, the proposed novel schemes involving complex coefficients provide a favourable alternative.
Global error. In the sequel, we choose t 0 = − 4 T 0 for some T 0 > 0 as initial time and determine numerical approximations u app (T ) at final time T = 4 T 0 for different time stepsizes τ = T −t 0 N ; in addition, a reference solution u ref (T ) is computed numerically to high accuracy. In order to obtain numerical results that are independent of the prescribed initial value, we determine the global errors in the corresponding fundamental matrix solutions  chosen parameter regimes, it suffices to use a relatively low order truncated Taylor series approximation for the computation of the matrix-exponentials to illustrate the performance of the novel schemes, whereas other approximations to the matrix-exponentials might be more appropriate for longer times or more challenging models, respectively. In fact, for the novel schemes, it turned out that the Taylor method is nearly optimal, since polynomials of degrees four to eight suffice in most cases; contrary, the optimal truncation order for the schemes proposed in [2], resulting in larger values (19), is usually higher. In order to measure the computational cost of a commutator-free quasi-Magnus exponential integrator by an M th-order Taylor approximation, we count the number of matrix-vector products per time step cost Numerical results. In a first test comparing the performance of the novel schemes, we fix the frequency, the time interval, and the dimension of the system The order of the Taylor series approximation is set to M = 6; for the computation of the action of the matrix-exponential involving the commutator, a second-order Taylor approximation is applied, and a first-order approximation essentially yields the same results. The  amplitude of the initial value is set to V 0 = 2 or V 0 = 1 2 , respectively. Besides, we contrast the special case of vanishing dissipation to the case δ = 10 −1 , respectively.
(i) Fourth-order schemes. The numerical results obtained for the optimised CFQM exponential integrator with three exponentials proposed in [2] and for the novel schemes with J = 4 and J = 5, respectively, are displayed in Figure 4. We observe that for medium to high accuracies the novel schemes perform better. (ii) Sixth-order schemes with real coefficients. The numerical results obtained for the optimised CFQM exponential integrator with six exponentials proposed in [2] and for the novel method involving a single commutator are displayed in Figure 5. We observe that the novel scheme performs better; in addition, in the presence of a dissipative term, its stability properties are clearly superior. (iii) Fifth-and sixth-order schemes with complex coefficients. The numerical results obtained for the novel fifth-and sixth-order schemes with complex coefficients are displayed in Figure 6. The schemes are appropriate for dissipative quantum systems and show a good performance, whenever high accuracy is desired, see also Figure 5. In particular, the non-time-symmetric fifth-order scheme involving three exponentials is favourable in efficiency.  [5] CF 5C [6] CF 4 [4] RK 7 [6] M 13 [6] log 10  Comparison with standard integrators. Next, we take into account that the order of accuracy when computing the action of the arising matrix-exponentials on vectors influences the performance of the novel schemes; in our case, this corresponds to the order of the Taylor series approximation. In order to illustrate this influence, we consider the novel schemes CF [4] 4 , CF [5] 3 , CF [6] 5C , which performed well in the previous test, and compare them with the standard sixth-order interpolatory Magnus integrator involving a single exponential and a standard explicit sixth-order Runge-Kutta method with seven stages RK [6] 7 , M [6] 13 ; for the Magnus integrator, each product of the exponent with a vector requires 13 products, whereas the explicit Runge-Kutta method only requires seven products per time step and thus is considerably cheaper than all other exponential integrators. Consequently, for our simple test problem, this would be in favour of the Runge-Kutta method.
For most accuracies, a Taylor series approximation of order M = 6 for the fourth-order scheme and M = 8 for the fifth-and 6th-order schemes is nearly optimal. However, due to the fact that the matrix-exponentials are replaced by polynomial approximations, the U app (T ) is not exactly preserved.
In Figure 7, we display the global errors as well as the errors in the norm preservation | U app (T ) − 1|. We observe that the standard Magnus integrator is not appropriate for the considered situation; its performance would improve for low dimensional problems, where matrix-matrix products are feasible. In spite of its very low computational cost, the Comparison of optimised sixth-order scheme CF [6] 6 (dashed lines), see [2], and novel scheme CF standard explicit Runge-Kutta method shows a poor performance; this is more manifest when qualitative properties are of interest or for longer time integrations.
A benefit of (commutator-free quasi-) Magnus integrators is that different underlying quadrature rules as well as different approximations to the exponentials can be used; their optimal choice will depend on the considered problem and the sizes of the parameters. In order to illustrate this, we set δ = 10 −2 , (V 0 , ω) ∈ {(5, 2), (2, 5)} , T 0 = 5 , d = 10 .
In Figure 8, we display the results obtained for different orders of Taylor series approximations and the optimised sixth-order CFQM exponential integrator proposed in [2] and for the novel fifth-order scheme with complex coefficients. We conclude that the optimal choice of Taylor series approximations indeed depends on the accuracy desired and on the particular problem considered; in all cases, the novel scheme is superior.

A.1 Magnus integrators in terms of arbitrary quadrature rules
To build methods of order 2K that can easily be used with any quadrature rule of order 2K or higher we introduce the one-dimensional momentum matrices with respect to the interpolating matrix A in (10) 1 − (−1) i+j (i + j)2 i+j α j , i = 0, . . . , K − 1 .
For example, at order four and six we have , and the Magnus expansion can be written in terms of A (i) . If we take into account that we can safely replace A (i) by A (i) in the Magnus integrators and then to use any quadrature rule of order p ≥ 2K to approximate A (i) , i = 0, . . . , K − 1. If { w i , c i } K i=1 correspond to the weights and nodes of a quadrature rule of order p ≥ 2K, then we have that where A j = A( c j τ ). This relation allows to change from the Gauss-Legendre quadrature rule to any other quadrature rule, e.g. we can replace that keeps accuracy up to order 2K so, in (5) we can take with coefficients a = a G −1 Q to obtain another approximation in terms of the evaluations at the new quadrature rule that is correct up to terms of order τ 2K .