{'fraction_non_alphanumeric': 0.037665477149764016, 'fraction_numerical': 0.010268216875791412, 'mean_word_length': 4.5381869182710695, 'pattern_counts': {'":': 0, '<': 4, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Propagators for the time-dependent Kohn-Sham equations: multistep, Runge-Kutta, exponential Runge-Kutta, and commutator free Magnus methods

We examine various integration schemes for the time-dependent Kohn-Sham equations. Contrary to the time-dependent Schr\"odinger's equation, this set of equations is non-linear, due to the dependence of the Hamiltonian on the electronic density. We discuss some of their exact properties, and in particular their symplectic structure. Four different families of propagators are considered, specifically the linear multistep, Runge-Kutta, exponential Runge-Kutta, and the commutator-free Magnus schemes. These have been chosen because they have been largely ignored in the past for time-dependent electronic structure calculations. The performance is analyzed in terms of cost-versus-accuracy. The clear winner, in terms of robustness, simplicity, and efficiency is a simplified version of a fourth-order commutator-free Magnus integrator. However, in some specific cases, other propagators, such as some implicit versions of the multistep methods, may be useful.


Introduction
In 1984, Runge and Gross 1 extended the fundamental theorems of density-functional theory to the time-dependent case, thereby founding time-dependent density-functional theory (TDDFT). 2 Over the years, TDDFT has become a very popular tool for the calculation of properties of atoms, molecules, nanostructures, or bulk materials thanks to its favorable accuracy/computational cost relation. It can also be used for a wide range of applications, e.g. to calculate optical properties, 3 to study nuclear dynamics, 4 charge transfer processes, 5 electronic excitations 6 and ultrafast interaction of electrons with strong laser fields, 7 to name a few. At the core of these simulations are the time-dependent Kohn-Sham equations (2) The Kohn-Sham Hamiltonian is a linear Hermitian operator, that can have an explicit timedependence (e.g., if the atoms are moving, or in the presence of a laser field), and an implicit time-dependence though the density. As the density (2) is written in terms of the Kohn-Sham orbitals, Eqs. (1) are indeed a set of nonlinear equations. Moreover, the KS Hamiltonian at time t depends on the full history of the density at all times t ′ ≤ t, and not only on its value at time t. These "memory" effects are rather important in several circumstances (for example, for multiple excitations), and have been extensively studied. 8 Unfortunately, there is a lack of memory-including exchangecorrelation functionals, and the dependence on the full history make the solution of the TDKS equations rather complex. Therefore, almost all applications of real-time TDDFT invoke the adiabatic approximation, that states that the KS Hamiltonian at time t only depends on the instantaneous density at the same time, as we already assumed in Eqs. (1). Upon discretization of the electronic Hilbert space, the TDKS equations fall into the category of systems of initial-value first-order ordinary differential equations (ODEs), i.e. they have the general form: Note that if history were to be considered, the TDKS equations would no longer be an ODE system: they would belong instead to the more general family of delay differential equations, or time-delay systems. 9 Centuries of research since the early days of Newton, Euler, etc. have produced a wide variety of numerical methods to solve ODEs. [10][11][12] Any of those can in principle be applied to the TDKS equations, but finding the most efficient one is a difficult task. [13][14][15][16][17][18][19][20][21][22] In the following paragraphs, we make a necessarily non-exhaustive recap of the ODE schemes that have, or have not, been tried for TDDFT problems.
A first division can be established between one-step and multi-step methods. The former provide a recipe to compute an approximation to the solution at some time t from the knowledge of the solution at a single previous time t − ∆t. The latter use information from a number of previous steps t − ∆t, t − 2∆t, etc. Multistep formulas have been scarcely used in the quantum chemistry or electronic structure community, and to our knowledge never for TDDFT calculations. Perhaps the reason is the need to store the information about a number of previous steps, a large amount of data for this type of problems. The most common alternatives are the implicit and explicit formulas of Adams, and the backward-differentiation formulas (BDFs).
For what concerns single-step methods, arguably the most used and studied one is the family of Runge-Kutta (RK) integrators. 23 This includes the implicit and explicit Euler formulas, the trapezoidal (also known as Crank-Nicolson 24 ) and implicit midpoint rules, the explicit RK4 formula (considered "the" RK formula since it is perhaps the most common), the Gauss-Legendre collocation methods, the Radau and Lobatto families, etc. Moreover, numerous possible extensions and variations are possible: partitioned RK, embedded formulas, use of variable time-step, extrapolation methods on top of the RK schemes (e.g. the Gragg-Bulirsch-Stoer algorithm 25 ), composition techniques, the linearly implicit Rosenbrock methods, etc. (see Refs. 10-12 for a description of these and other ideas). Once again, many of these options have never been tested for TDDFT problems.
Linear autonomous ODE systems can also be solved directly by acting on the initial state with the exponential of the operator that de-fines the system. Quantum problems without an explicitly time-dependent Hamiltonian belong to this class. The problem of the quantum propagator can therefore be reduced to finding a good algorithm to compute the action of the exponential of an operator. Various alternatives exist: a truncation of the Taylor expansion, 26 the Chebychev 27 and Krylov polynomial expansions, 28,29 Leja and Padé interpolations, 30 etc.
For non-autonomous linear systems (e.g. quantum problems with time-dependent Hamiltonians), a time-ordered exponential must substitute the simple one. By using short-enough time-steps, however, a constant Hamiltonian can be applied within each interval, and the simple exponential methods mentioned above may suffice. Otherwise, one can resort to Magnus expansions. 31 Perhaps the most used one is also the simplest: the second-order Magnus expansion, also known as the exponential midpoint rule. More sophisticated (higher order) expansions require the computation of commutators of the Hamiltonian at different times, a costly operation. Recently, commutator-free Magnus expansions have also been proposed. 32 Other recent options essentially based on the exponential (and tested for TDDFT) are the non-recursive Chebychev expansion of William Young et al, 33 or the three-term recurrence of Akama et al. 17 An old-time favourite in condensed matter physics is the split-operator formula. 34 It belongs to the wide class of splitting techniques, whose simplest members are the Lie-Trotter 35 and Strang 36 splittings. In chemistry and physics, these use the usual division of the Hamiltonian into a kinetic and a potential part, as both can be treated exactly in the proper representation -the main computational problem is then reduced to the transformation to and from real and Fourier space. More sophisticated splitting formulas have also been developed (see e.g. Refs. [37][38][39][40]. The TDDFT Hamiltonian may also be divided into a linear and a non-linear part. The non-linear part must of course include the Hartree, exchange, and correlation potentials. The kinetic term is almost always included in the linear term. This is considered to be the term responsible for the possible stiffness of the equations. It is difficult to give a precise definition of stiffness, and a pragmatic one is generally accepted: "stiff equations are equations where certain implicit methods perform better, usually tremendously better, than explicit ones". 41 Implicit methods require the solution of nonlinear algebraic equations. Besides outperforming explicit methods for stiff cases, some of them may also have the advantage of preserving structural properties of the system, such as symplecticity -a topic that we will discuss later on. For cases in which one part of the equation requires an implicit method, but another part does not, the implicit-explicit (IMEX) methods were invented. 42,43 Another recent approach that relies on the separation of a linear and a non-linear part are the exponential integrators. [44][45][46][47] There are various subfamilies: "integrating factor" (IF), "exponential time-differencing" (ETD) formulas, exponential RK, etc. These techniques have not been tested for non-equilibrium electron dynamics in general, or TDDFT in particular, until very recently. 15 An alternative that has been followed by several groups is the transformation of the system to the adiabatic eigenbasis, or to a closely related one (a "spectral basis", generally speaking). In that appropriately chosen basis, the dimension of the system is small, and any method can do the job. The burden of the task is then transferred to the construction and update of the basis along the time evolution, an operation that involves diagonalization. Refs. 14,16,48-50 are some recent examples, some of them based on Houston states, 51 that report notable speed-ups over conventional methods. This result seems, however, to depend on the type of problem, implementation details, etc.
The former list of algorithms, though long, was not exhaustive: for example, we can also mention Fatunla's algorithm, 29,52,53 or the very recent semi-global approach of Schaefer et al. 54 based on the Chebychev propagator. It becomes evident that the list of options is extensive, making the identification of the most efficient, accurate, or reliable algorithm a hard task. Some of us presented in 2003 a performance analysis of various propagation methods for the TDKS equations; 13 it is the purpose of this Article to continue along those lines, by investigating other promising propagation schemes and by providing several benchmarks in order to assert their efficiency in realworld applications. In particular, we look here at multi-step based propagators, exponential RK integrators (along with the standard RK), and a commutator-free version of the Magnus propagator. We implemented these propagation schemes in our code octopus, 55,56 a general purpose pseudopotential, real-space and realtime code.
The remaining of this article is organized in the following way: first we study in Section 2 the theory regarding the propagation schemes and its relation with the KS equations, paying special attention at the issue of symplecticity. Then in Section 3 we show the benchmarks obtained for the different propagation schemes. Finally, in Section 4 we state our conclusions.
for some discrete time step ∆t (that we will consider to be constant along the evolution in this work). The evolution operator is given bŷ i.e. it is the time-ordered evolution operator. The non-linearity, however, implies that a linear evolution operator linking ϕ m (t − ∆t) to ϕ m (t) does not exist. We may however still assume the existence of a nonlinear evolution operator, that is usually called a flow in the general case, [Eqs. (3) and (4)]; it is defined as This is the object that must be approximated through some algorithm -an algorithm that of course takes the form of a linear operator whenever employed on linear systems. To choose a numerical method to propagate the TDKS equations one is usually concerned by its performance and stability. Performance is loosely speaking related to the computer time required to propagate the equations for a certain amount of time. Stability, on the other hand, is a measure on the quality of the solution after a certain time. For linear systems (or for propagators applied to linear systems), it is possible to give a simple mathematical definition of stability. A propagator is stable below ∆t max if, for any ∆t < ∆t max and n > 0,Û n (t + ∆t, t) is uniformly bounded. One way to assure that the algorithm is stable is by making it "contractive", which means that ||Û(t + ∆t)|| ≤ 1. Of course, if the algorithm is unitary, it is also contractive and hence stable; but if the algorithm is only approximately unitary, it is better if it is contractive. We can also talk about unconditionally stable algorithms if their stability is independent of ∆t and of the spectral characteristics ofĤ.
Clearly, in many cases stability can be enhanced by decreasing the time-step of the algorithm, i.e., by decreasing its numerical performance. In other cases, however, long-time stability is almost impossible to achieve for some methods.
A common strategy to develop stable numerical methods is to request that these obey a number of exact features (although this does not ensure the stability or the performance). There are a series of exact conditions that can be easily derived. For example, it is well known that, for linear systems with Hermitian Hamiltonians, the propagator is unitarŷ This property ensures that the KS wavefunctions remain orthonormal during the timeevolution. Algorithms that severely violate (8) will have to orthogonalize regularly the wavefunctions, a rather expensive (N 3 ) operation, especially for large systems. Note that for the TDKS equationss it is not, strictly speaking, correct to speak of unitarity due to the nonlinear character of the propagators even if the orthogonality condition still holds among the KS orbitals (see the discussion in Ref. 57). For systems that do not contain a magneticfield or a spin-orbit coupling term (or any other term that breaks time-reversal symmetry), the evolution operator fullfilŝ This relation is rather important in order to ensure stability of the numerical propagator, and it is often violated by many explicit methods.

Symplecticity
The geometrical structure of an ODE system, as well as that of its numerical representation (i.e. the propagator), is another important issue to consider. 12 In this context, an important property is symplecticity.
Given any system of ODEs, the flow is a differentiable map. The first requirement for a flow to be symplectic is that the system is formed by an even number of real equations. Any complex system, however, may be split into its real and imaginary parts, and is equivalent to a system with an even number of real equations.
The system of equations is also required to be Hamiltonian: a system is Hamiltonian if it follows the equation of motioṅ where y ∈ R 2n , and H is some scalar function of y. It is usual to decompose y = (q, p) T , leading to the well-known Hamilton equations of motionq with q i and p i elements of the vectors q and p. The flow of a Hamiltonian system is symplectic. Roughly speaking, the inverse is also true. 12,58 One can easily prove that the (usual) Schrödinger equation forms a Hamiltonian system, 59 and is therefore symplectic. It is possible to perform the derivation in coordinate space, but the proof is somewhat simpler if we expand the wave-function in a given basis set where {|Ψ i } forms an orthonormal basis and c i (t) = Ψ i |Ψ(t) are the time-dependent expansion coefficients. The Schrödinger equation is thus transformed into: where the Hamiltonian matrix H defined by H ij = Ψ i |H|Ψ j , and c is the vector of the coefficients. We now split the coefficients c i of the wave-function into their real and imaginary parts i.e., q i = √ 2 ℜc i , p i = √ 2 ℑc i . We can now define a Hamiltonian function of the vectors q and p H(q, p) = Ψ(q, p)|Ĥ|Ψ(q, p) where we have used the fact that, if H is Hermitian, then H(q, p) must be real, and the real part of H is symmetric (ℜH T = ℜH) while its imaginary part is anti-symmetric (ℑH T = −ℑH). We can now calculate the partial derivatives ∂H(q, p) In order to find the equations of motion for the q and p coordinates, we rewrite Eq. (15a) as: The proof follows by separating the real and imaginary parts of (18) and comparing them to Eqs. (17). The Schrödinger's equation forms a Hamiltonian, symplectic system.
Whether or not an ODE system is symplectic has important theoretical consequences (the flow preserves the volume in phase space, for example). Numerically, the algorithm that we choose to approximate the real flow defines a numerical flow that may or may not be symplectic. It is of course convenient for it to be: for example, one can demonstrate 60 that symplectic numerical flows lead to long term stability of the energy, that typically oscillates around its true value. Usually, the error in the conservation of other constants of motion also behaves better when symplectic algorithms are used. In the following, we shall prove, following a similar procedure to the one above for the Schrödinger equation, that the TDKS equations, in the adiabatic approximation, form a symplectic, Hamiltonian system. Therefore, it is convenient (although not strictly necessary) to choose symplectic algorithms to approximate the TDKS propagator.

Symplecticity and the TDKS equations
For the TDKS equations, Eqs. (1), the Hamiltonian operator can be written aŝ where the terms represent the kinetic energy operator, the external potential, and the Hartreeexchange-correlation (Hxc) potential. In the coordinate representation, we have We now expand the KS orbitals in a one- The TDKS equations are thus transformed into the initial value probleṁ where the matrix H[c] is given by Note that the dependence on the (instantaneous) density is rewritten as a dependence on the full set of coefficients c. We again split the coefficients into their real and imaginary parts The TDKS equation may then be rewritten aṡ and separating into real and imaginary partṡ In order to rewrite the TDKS system as a classical Hamiltonian system, we need to find a Hamiltonian function H(q, p). It can be easily seen that the non-interacting energy of the KS system does not work. However, we can use the ground-state energy functional, which is given by evaluated adiabatically with the timedependent density. Remembering that the density is evaluated from the KS orbitals, we can write the energy as a functional of these Representing the orbitals by the new variables (q, p), we define a Hamiltonian function The first two terms can be treated exactly in the same way as for the standard Schrödinger equation. The non-interacting kinetic energy function reads while the external potential is Calculating the partial derivatives of the previous expressions, we arrive at and with a similar expression for ∂V (q, p)/∂q mi and ∂V (q, p)/∂p mi . Using (26), we see that these two terms verify the necessary conditions for a Hamiltonian system. We are left with the term E Hxc (q, p). Its partial derivatives can be computed with the help of the chain rule ∂E Hxc (p, q) ∂q mi = d 3 r δE Hxc δn(q, p; r) The density n(q, p, r) is the one that corresponds to the set of Kohn-Sham orbitals defined by the (q, p) coordinates. The functional derivative of E Hxc is the Hartree, exchange, and correlation potential In order to compute the partial derivatives of the density with respect to q and p, one needs to write it in terms of those variables Plugging these expressions into Eqs. (33) results in where the matrix V Hxc [q, p] is given by: Therefore, the partial derivatives of E Hxc also have the right structure, which concludes the proof that the TDKS equations form a Hamiltonian system.

Results
In order to analyze the performance of the integration schemes we used a "real world" benchmark based on the propagation of a benzene molecule. We placed the molecule in a spherical simulation box of radius r = 12 a.u., with a grid spacing of a = 0.4 a.u. At time zero, the system is subject to an instantaneous perturbation: i.e. each KS orbital, initially at its ground-state equilibrium value ϕ GS j is transformed at time zero into a slightly perturbed orbital ϕ j (t = 0 + ), corresponding to a sudden application of an electric field with strength k = 0.1 a.u. in the z-direction. Then it evolves freely for a total propagation time T = 2π a.u. We compared both the wave-function and the energy obtained at the end of the run with a reference "exact" calculation, performed with a very small time step and the explicit RK4 propagator. The error in the wave-function is then defined as and the error in the energy is defined as where ϕ exact m and E exact are the KS orbitals and the energy obtained from the "exact" calculation.

Exponential midpoint rule
We used the exponential midpoint rule (EMR), one of the propagators studied in Ref. 13, as a base for comparison with the new schemes. The EMR prescribes: where ϕ is the average wavefunction: The EMR is second order in ∆t, symplectic, and preserves time-reversal symmetry. It is also an implicit scheme as it requires the Hamiltonian calculated with the average wave-function. The non-linear equations (42) and (43) can be solved, e.g., by iteration until self-consistence is achieved. The first iteration can be started by making use of an extrapolated Hamiltonian We will use the shorthand notationĤ extr for a Hamiltonian that is obtained via extrapolation or interpolation to time τ from a number p of known Hamiltonians: H t−∆t ,Ĥ t−2∆t , . . . ,Ĥ t−p∆t . We will also use the notationĤ In practice, most of the times one does not iterate the self-consistent procedure, but uses Eq. (44) directly. This leads to an explicit EMR, that is the method used in the remainder of this work. Of course, this approximated method no longer fulfills the exact properties stated above.
The definition of the algorithm must be complemented with a recipe to compute the action of the exponential of an operator on a vector. There are a variety of possibilities, passing by a truncated Taylor expansion, a Lanczos expansion, the split-operator scheme (as well as any of the higher-order variants of this), etc. For our purposes we decided to use the first, namely a Taylor expansion truncated to fourth order.
One may also design other exponential-based methods that can be considered variations of the EMR. For example, in Ref. 13 we defined the "enforced time-reversal symmetry" (ETRS) scheme This algorithm was designed to improve on the preservation of time-reversal symmetry. It is also an implicit method, and the non-linear equation (46) can be solved iteratively. Alternatively, one can use an extrapolated Hamilto-nianĤ extr (t) in Eq. (46), leading to the approximate ETRS (AETRS) algorithm.

Commutator-free
Magnus expansions Restricting the discussion momentarily to linear systems, for time-dependent Hamiltonians, the evolution operator (6) has a rather complicated form involving an integral over time and a time-ordering operator. It is natural to wonder, however, if there exists an operatorΩ(t, t − ∆t) that makes the following expression exact In 1954 Magnus 31 showed that, for some neighborhood of t, there exists an infinite series such thatΩ and provided a recursive relation to find the operatorsΩ m . This recursive relation involves nested commutators of the Hamiltonian at different times. To obtain a Magnus propagator of order 2n,Û M (2n) , one truncates the series at the n-th term, and approximates each term with some n-th order quadrature formula. As an example of this procedure, the aforementioned EMR can be regarded as a Magnus expansion of order two (although, strictly speaking, only for linear systems): This second order formula is unusual as it does not involve commutators. For higher orders, the main difficulty arises from the evaluation of the commutators. To circumvent this problem, Blanes et al. 32 developed a series of alternative Magnus expansions that are free of commutators, and also address the nonlinear case.
We implemented the fourth order commutatorfree version of the Magnus expansion, presented in Eq. (43) of Ref. 32, and labelled as the "Method 4" in page 6 of Ref. 61. This method, that we will refer to in the following by "CFM4", is given by: for some carefully chosen constants α 1 , α 2 and intermediate times t 1 , t 2 . The application of this method to the nonlinear TDKS equations leads again to an implicit rule, as we needĤ (t 1 ) andĤ (t 2 ) . Therefore, we have implemented an approximate version, again relying on extrapolated Hamiltonians. If this extrapolation is peformed at fourth order (i.e. using at least four previous steps), then the order of the method is preserved. Figures. 1 and 2 depict the results obtained with the CFM4, EMR, ETRS and AETRS methods. The top (botton) panel of Fig. 1 shows the error in the energy (wave-function) as a function of the time-step. We used logarithmic scales in both axes, so that the curves become straight lines in the small ∆t limit (until numerical precision is reached). The slope of those lines is given by the order of each method -at least for the error in the wave-function. For larger values of the time-step, the curves are no longer straight lines, and may actually exhibit a faster behavior: for example, the EMR, ETRS and AETRS methods behave as fourth order propagators for larger ∆t, whereas their order is actually two. As we can see in Fig. 1, for the largest time-steps (up to 10 −2 a.u.) all the methods have similar precision, except the EMR, which becomes unstable (this is the reason why this data point is missing). When the time-step decreases, EMR, ETRS and AETRS behave as order-two methods while CFM4 maintains its fourth order throughout the whole range of ∆t. This makes CFM4 significantly more precise than the other propaga-tors for ∆t < 10 −2 a.u.
In Fig. 2 we show the cost (measured in seconds) of the propagation as a function of the error in the wave-function, again in logarithmic scale. From these kind of plots one can identify the best performing method for a given required precision. This required precision must be decided a priori by the user, and it is problem dependent. For the largest values of the error the performance of all the integrators is very similar. For smaller values, EMR, ETRS and AETRS have similar cost, but CFM4 is significantly faster. This makes CFM4 the best method overall.

Multistep methods
In 1883 J. C. Adams and F. Bashforth proposed multistep methods in the context of fluid mechanics. 62 These methods use s > 1 previous steps in order to calculate the following one. They require a starting procedure to provide those first s steps. The simplest procedure consists in using a single-step method. In our case we used the standard explicit fourth-order RK (described below).
We examined linear multistep formulas given by (51) where {a k } s−1 k=0 and {b k } s k=0 are the coefficients that determine the method. If b s = 0, then the method is explicit, since the equation is an explicit formula for ϕ(t). If b s = 0 then the method is implicit, as it provides a relation between ϕ(t) and f (ϕ(t), t). If we consider the dynamical function relevant for TDDFT and we define the shorthand notation we finally arrive at The first multistep integrators that we studied belong to the family of explicit Adams methods, also known as Adams-Bashforth (AB) methods. They are explicit (b s = 0) and the coefficients a are: a s−1 = −1 , and a s−2 = · · · = a 0 = 0. The remaining b k are chosen such that the methods have order s, which determines them uniquely. The method then reads: The implicit Adams, or Adams-Moulton (AM) family is similar to the Adams-Bashforth methods in that they also have a s−1 = −1 and a s−2 = ... = a 0 = 0: Again, the b coefficients are chosen to obtain the highest possible order. The Adams-Moulton methods are implicit methods, since the restriction b s = 0 is removed. This fact permits the increase of the order of the error: an s-step Adams-Moulton method is of order s + 1, while an s-step Adams-Bashforth method is only of order s. Equation (56) was solved iteratively. We also implemented a "linearized" version of the Adams-Moulton formula (lAM), where we used an extrapolation of the Hamiltonian at time t, thereby transforming Eq. (56) into a linear equation. Another possible simplification of the Adams-Moulton formula regards the use of the so-called "predictor-corrector" schemes, which avoid the linear system solution altogether by turning the implicit method into an explicit one. In our implementation, it consists of using Adams-Bashforth to get an approximated ("predictor")φ(t), and use this to ob-tain the Hamiltonian in the left-hand side of Eq. (56). We named this procedure the Adams-Bashforth-Moulton (ABM) method.
The backward differentiation formulas (BDF) are implicit methods with b s−1 = · · · = b 0 = 0 and the other coefficients chosen such that the method has order s (the maximum possible). These methods are especially suited for the solution of stiff differential equations. The general formula for a BDF can be written as:  We implemented these families of integrators in octopus, ran these five methods with steps s = 1, . . . , 5, and compared them among each other and with the EMR. Figures. 3 and 4 show the best candidate from each family. The number accompanying the name of the propagator indicates the number of previous steps s used in the calculation. As we can see, the EMR is more stable than any of the multistep methods for large time-steps (especially AB4, which is the most unstable), but is outclassed in precision by every other propagator. This is not surprising, as they are methods of order 4 (AB4 and ABM4), 5 (AM4) or 6 (linearized AM5). The most precise method for a given time-step is the lAM5, reaching the numerical precision of our machines for the smallest time-steps. In Fig. 4 we plot the cost of the methods as a function of the error. AB4 cannot compete in precision, stability or performance with the EMR. For error values larger than 10 −7 the EMR is the fastest propagator, while for smaller values it is overcome by AMB4 and linearized AM5. AM4 is, as expected, the most computationally expensive method, with the linearization procedure dramatically improving its speed.
In Fig. 5 we represent the BDF results for s = 1, . . . , 5. Our aim here is to illustrate one important characteristic of the multistep methods, namely that the cost does not increase significantly with the number of previous steps s. This can be clearly seen on the left panel of Fig. 5. Furthermore, in the right panel we can see that each extra step included in the method increases its order by one. Then, why not increasing the number of steps to a very larger number? First, there is a memory issue, as the previous s steps have to be stored in memory. But more importantly, as the number of previous steps increases, the stability of the method decreases. This can be seen in both panels of Fig 5. Both BDF1 and BDF2 have better stability properties than EMR, but as soon as we make s ≥ 3 we need to reduce the time-step by a factor of 16 to avoid the breaking down of the method. This reduction of the stability region with the number of steps seems to hold for all linear multistep methods. Moreover, for BDF there is a mathematical proof that states that for s ≥ 7 these methods are unstable (check  section III.3 from Ref. 11 for a more detailed  explanation).
Finally, one important caveat of multistep methods is that they cannot be symplectic. In fact, the definition does not even apply, as a multistep algorithm is a map from several previous steps into the next one, and one cannot speak of a flow in the usual way. There are however some ways to understand symplecticity also for these methods, 12 but the conclusion is in any case negative, and the long-term stability properties of these methods is disappointing.
where the functions Y i are defined as at the time-steps is called an s-stage RK scheme. To specify a particular method, one needs to provide the integer s (the number of stages), and the coefficients a ij , b i , and c i (for i = 1, 2, . . . , s). These are usually arranged in a mnemonic device, known as a Butcher tableau: When a ij = 0 for i ≤ j the method is explicit, whereas in all other cases the method is implicit. Explicit RK methods are generally unsuitable for the solution of stiff equations because their region of absolute stability is small. These shortcomings motivated the development of implicit methods. They are visually easy to identify looking at their tableaux, as they include non-zero entries in the upper triangle.
For the implicit methods we need to solve a system of algebraic equations, the dimension of which grows with the number of stages: For a method with s stages, the equation has m × s unknowns, where m is the dimension of the original system. In contrast, linear multistep methods only require the solution of mdimensional algebraic equations.
For a RK scheme to be symplectic, one can prove 64 that the s × s-matrix M with coefficients has to satisfy M = 0. This implies that no explicit RK scheme can be symplectic. We studied the RK propagators up to order four. The reason behind this choice is that, up to this order, the required number of stages s for explicit methods is equal to the desired order of the method. From order five onward, however, s is strictly greater than the desired order. 12,23,63 Therefore, the precision gained by increasing the order does not compensate for the increase in the computational cost. Regarding explicit methods, the most widely known RK scheme is the fourth order explicit RK method, also known as "RK4" or simply "the" RK method. Its Butcher tableau is (63) Unfortunately, this method, as any other explicit one, is not symplectic.
A particularly relevant branch of the RK family is the Gauss-collocation scheme. Gausscollocation methods of s stages have order 2s, and they are both symplectic and symmetric. We chose two of these methods for our benchmarks, specifically the second-order "implicit midpoint rule" (that we will call imRK2), and the fourth-order method (imRK4).
The tableaux for imRK2 is and leads to the non-linear equation for ϕ(t) where ϕ = 1 2 [ϕ(t) + ϕ(t − ∆t)]. Note that this equation is similar, but not identical, to the trapezoidal or Crank-Nicolson rule. These two methods are in fact conjugate, 10 and the name "Crank-Nicolson" is sometimes used indistinctly for both.
The Butcher tableux for imRK4 is Once again, we face non-linear equations, that we implemented through self-consistent iterative procedures similar to the ones described for the AM formulas. Each iteration requires the solution of a linear system. We also define "linearized" variants of RK as the simplified versions in which we just perform the first step of the self-consistent cycle with an extrapolated Hamiltonian -a strategy that always seems to produce the best performing algorithm.
We plotted the errors in the energy and wavefunction as a function of the time-step in Fig. 6. The points for imRK4 and lRK4 in the energy panel that do not appear for time-steps smaller than ∼ 10 −1.8 are those that reached the precision of our machines. From the bottom panel of this figure we can see that the EMR is significantly more precise than the second order RK methods for the wave-function, and it can even compete with the fourth-order methods for time-steps smaller than 10 −2.5 . On the other hand, as far as the energy is concerned, we can see that EMR is outclassed by every RK method, and especially by the implicit methods and their linearized versions. The EMR also breaks down for the larger time-steps values, like the RK4 method.
We also found that the linearized versions of the implicit methods behave similarly to their full counterpart as far as the wave-function is concerned, but that there is a significant difference in the error of the energy (see the curves for imRK2 and lRK2). The explicit RK4 method performs worse both for the wave-function and the energy when compared with the implicit methods.
In Fig. 7 we show the cost as a function of the error in the wave-function. Here the explicit methods have the advantage, with both the EMR and RK4 performing around an order of magnitude faster than the implicit propagators. The EMR has the best performance up to an error of 10 −8 . The lRK2 performs better than the imRK2, while imRK4 and the lRK4 have the same cost (implying that the self-consistent cycle converged in one iteration). Among the implicit methods, lRK2 is the best performing method up to an error of 10 −4 , but for smaller values either the lRK4 or the imRK4 propagators are the best choice.
Finally, a word of caution regarding these comparisons between explicit and implicit methods: the latter require the solution of linear systems, and their performance will depend on the performance of the linear solvers. The existence or not of preconditioners, for example, make these comparisons very system and implementation dependent.   Figure 7: Cost of the method, as a function of the error obtained (in the wave-function), for the various RK methods (implicit and linearized RK2 and RK4 and explicit RK4) and for the EMR propagator.

Exponential RK Schemes
Recently, we saw the appearance of the socalled "exponential Runge-Kutta" (ERK) schemes. 45,65,66 The main appeal in this family of propagators lies in its ability to tackle stiff problems. The key idea is solving the stiff part of the equation precisely, and approximating the remaining part by a quadrature formula. Let us rewrite our non-linear TDKS equation as where T is the kinetic operator (the stiff part), and the last term is the Kohn-Sham potential acting on the orbitals. An ERK scheme for this equation has the form with the definition Equations (69) are in general a set of s nonlinear equations. The constants c i and the operator functions a ij and b i fully determine the algorithm. These constants reduce to an underlying RK scheme at T = 0, so that a ij (0) = a ij , and b ij (0) = b ij . Just as with normal RK schemes, the methods can be explicit or implicit.
For the explicit ERK schemes, we have to compute some auxiliary functions Y i (Y j , t i ), with j < i and i = 1, ..., s, where s is the number of stages of the method. Here, the coefficients a ij and b i are linear combinations of the so-called φ k functions, defined by the recurrence relation For the evaluations of the φ k functions we used this Taylor expansion. This allows us to compute both the regular exponential function and these φ k , and any linear combination of them in a simple subroutine, simplifying the implementation of the generalization of the explicit ERK methods.
The simplest example of this family is the exponential version of the Euler method, given by which is an order 1 method (we call it ERK1). We implemented a general algorithm for a broad family of ERK schemes of any order described by Hochbruck. 46 We show results for the best performing methods for orders 2, 3 and 4: method (5.4) for order 2 (we name it ERK2), method (5.8) for order 3 (we name it ERK3) and method (5.17) for order 4 (we name it ERK4).
Explicit exponential RK schemes cannot be symplectic (just as explicit "normal" RK ones), but implicit ones can. 66 This is achieved if the underlying RK method is symplectic and if the functions a ij and b i obey Following this recipe, we implemented the exponential version of RK2 (labeled imERK2 in the figures), characterized by s = 1, c 1 = a 11 = 1 2 , and b 1 = 1, resulting in the equations Figs. 8 and 9 display the numerical results obtained for this implicit method, and for the four explicit methods mentioned above, compared to the EMR scheme. As we can see, all the exponential methods keep their order for the values of the time-step where they do not break down. Furthermore, the implicit version of the exponential RK2 propagator has a wider range of stability than the explicit versions. From the top panel we can see that imERK2 has slightly smaller errors for the energy than EMR2 up to ∆t ∼ 10 −2.7 , and it is always better than the EMR as far as the error in the energy is concerned. The ERK4 beats every other propagator in the top panel. On the other hand, for the wave-function the EMR is on par with ERK4, being slightly more precise for the largest values of the time-step, and only being overtaken by ERK4 when ∆t < 10 2.5 . If we compare the explicit and implicit ERK2, we can see that imERK2 has a smaller error in the wave-function than ERK2, but in the energy comparison ERK2 is better for time-step values below 10 −2.5 . As we have mentioned before, this figure clearly shows that the methods behave as expected from the theoretical formulas, maintaining their order during the whole range of time-steps studied.
From Fig. 9 we can see that these methods are computationally expensive, with none of the exponential RK methods coming close to the EMR cost. Among the ERK family, ERK4 has the best performance for values of the error in the wave-function below 10 −2.2 , making it the best overall ERK method from the ones we tested.

Conclusions
We implemented and analyzed four families of numerical integrators for the Kohn-Sham equations in our code octopus, specifically commutator-free Magnus expansions, multistep methods, Runge-Kutta propagators and exponential Runge-Kutta integrators. These were compared to the previously studied exponential mid-point rule, enforced time-reversal symmetry and approximately enforced time-reversal symmetry propagators. For each method we evaluated the error in the wave-function and the energy as a function of the time-step, together with the cost in computational time as a function of the error. ERK2 ERK3 ERK4 Figure 9: Cost of the method, as a function of the error obtained (in the wave-function), for the various exponential RK methods (exponential Euler method, implicit RK2 and explicit RK2, RK3 and RK4) and for the EMR propagator.
Among the new families of propagators studied in this paper, the fourth-order commutatorfree Magnus expansion beats every other propagator in terms of cost/accuracy, making it the recommended method for TDDFT. The multistep integrators main advantage is that the computational cost remains constant independently of the number of previous steps considered. The exponential Runge-Kutta propagators do not show any clear advantage over the regular Runge-Kutta methods, with the explicit Runge-Kutta method of fourth-order being usually the best choice. The exception are stiff problems or in situations where a high degree of conservation of some quantity is required. In such cases, the symplecticity of the implicit versions of Runge-Kutta comes into play.
We have shown how the TDKS equations, in the adiabatic approximation, form a Hamiltonian and therefore a symplectic ODE system. Therefore, for long time propagations one should benefit from the use of structure preserving algorithms. This fact discourages the use of multistep schemes, for example, and favors implicit schemes that are unfortunately less costeffective.
The numerical integration of first-order ordinary differential equations is a very active field of research, with new schemes being proposed and old ones refined regularly, and many other methods still untested. We can therefore still expect new developments in the numerical propagation of the time-dependent Kohn-Sham equations, opening the way for the study of larger systems for longer periods of time.