Numerical integration of quantum time evolution in a curved manifold

The numerical integration of the Schr\"odinger equation by discretization of time is explored for the curved manifolds arising from finite representations based on evolving basis states. In particular, the unitarity of the evolution is assessed, in the sense of the conservation of mutual scalar products in a set of evolving states, and with them the conservation of orthonormality and particle number. Although the adequately represented equation is known to give rise to unitary evolution in spite of curvature, discretized integrators easily break that conservation, thereby deteriorating their stability. The Crank Nicolson algorithm, which offers unitary evolution in Euclidian spaces independent of time-step size $\mathrm{d}t$, can be generalised to curved manifolds in different ways. Here we compare a previously proposed algorithm that is unitary by construction, albeit integrating the wrong equation, with a faithful generalisation of the algorithm, which is, however, not strictly unitary for finite $\mathrm{d}t$.

The numerical integration of the Schrödinger equation by discretization of time is explored for the curved manifolds arising from finite representations based on evolving basis states. In particular, the unitarity of the evolution is assessed, in the sense of the conservation of mutual scalar products in a set of evolving states, and with them the conservation of orthonormality and particle number. Although the adequately represented equation is known to give rise to unitary evolution in spite of curvature, discretized integrators easily break that conservation, thereby deteriorating their stability. The Crank Nicolson algorithm, which offers unitary evolution in Euclidian spaces independent of time-step size dt, can be generalised to curved manifolds in different ways. Here we compare a previously proposed algorithm that is unitary by construction, albeit integrating the wrong equation, with a faithful generalisation of the algorithm, which is, however, not strictly unitary for finite dt.

I. INTRODUCTION
Algorithms for the efficient numerical integration of the time-dependent Schrödinger equation in discretised real time for finite representations have been discussed at some length in the literature (see e.g. Refs [1][2][3][4]), normally in the context of the single-particle states in time-dependent Hartree-Fock or time-dependent densityfunctional theory (TD-DFT) [5]. The consideration of evolving basis sets complicates matters, and there is less knowledge accumulated on good integrators for them and for arbitrarily quick basis evolution [6][7][8][9][10]. Evolving basis sets are routinely encountered in electronic structure calculations for which atom-centered basis functions are used, and where atoms move, that is, any first-principles dynamical calculation method in quantum chemistry or condensed matter and materials physics using atomic orbitals as basis sets. There are many such software packages that are widely used in either or both communities. For a brief review and links to codes used in quantum chemistry see, e.g., Ref. [11]; for methods and programs using atomic orbitals in condensed matter see, e.g., Refs. [12][13][14][15][16][17][18].
The equation for the evolution of quantum states for a moving basis is easily obtained. For a basis set {|e µ , t , µ = 1 . . . N }, H|ψ n = ih∂ t |ψ n straightforwardly becomes with H µν = e µ |H|e ν , S µν = e µ |e ν , D µνt = e µ |∂ t |e ν , and ψ ν n the coefficients in the expansion It is known that the evolution of a set of states following this equation is unitary in the sense that it preserves their mutual scalar products (see e.g. [19]). Therefore, if the evolving states are, for instance, the occupied Kohn-Sham states in TD-DFT evolution, they preserve their orthonormality, and the number of particles is conserved.
It is not obvious, however, how to guarantee such unitarity for approximate algorithms based on time discretization. Notice that the D µνt matrix does not need to be anti-hermitian if the evolution of basis vectors |e µ and |e ν is arbitrary (think e.g. of one of them not evolving while the other does), and therefore the usual thinking in terms of unitary matrices from the exponential of hermitian matrices does not apply, at least directly. In this work we focus on how the unitarity of the time evolution is affected by the discretization of time for numerical integration.
The semiclassical description of atomic collisions has made use of travelling orbitals [20], defined as , is a time-independent atomic-like orbital, and v is the velocity at which it is travelling, normally attached to a nucleus. They are very well adapted to the situation in which the electrons themselves travel with the atoms and, therefore, with the basis states, which is frequently the case in atomic collisions [20]. They are not so useful beyond that realm, as in e.g. atoms moving in metals, where the electrons are pushed around by a moving nucleus, but do not necessarily accompany it.
A more general procedure based on a Löwdin orthonormalization was proposed by Tomfohr and Sankey (TS hernceforth) in Ref. [6], which was built strictly to preserve unitarity of evolution for finite time steps for any evolving basis. It has been quite successfully used for projectiles traversing solids at non-adiabatic but relatively low velocity [21][22][23][24][25]. The method was further discussed in Ref. [9], where it was shown to affect the equation being integrated, with the potential to lead to inaccurate arXiv:2108.12614v2 [physics.comp-ph] 1 Nov 2021 propagation for high velocity, even in the converged low time-step limit.
Finally, a different way of approaching the integration of Eq. (1) is by rewriting it as and taking H µν − ihD µνt as a modified Hamiltonian matrix (see e.g. Ref. [7]). The behavior for finite time step of this pragmatic approach is not easy to discern from general considerations, since the D µνt matrix is not expected to be anti-hermitian, as mentioned above. It does work reasonably well, however [7]. Here we explore it further, both formally, and by explicitly comparing its time-step convergence with the very stable TS integrator [6], while the accuracy of the latter is further scrutinized. For a better understanding of the evolution we use in the following the recent geometric interpretation [9] of Eq. (1). The same paper proposed a strictly unitary integrator for an evolving basis, as long as the spanned Hilbert space Ω were invariant at all times, which is approximately the case for a well converged basis set. Here we explore the situation for an evolving Hilbert space Ω(t), defining a curved fibre bundle [9], for the general situation in which its curvature cannot be neglected.
Ref. [9] identifies the D µνt matrix as a connection in differential geometry. The integration procedure of Ref. [7] can then be interpreted in this context as using the connection as a gauge potential in the Hamiltonian [9]. We will hence refer to this procedure as the gauge-potential (GP) integrator.
The set {|e µ , t , µ = 1 . . . N } is the dual basis of {|e µ , t }, also a basis of Ω(t), satisfying e µ , t|e ν , t = δ µ ν , ∀µ, ν at any time t. ð t represents the covariant time derivative [9] defined as where D µ νt = e µ |∂ t e ν gives the connection in the manifold (note the convention in the order of indices).
Similarly, a bra evolves following where we have made use of the hermiticity of the Hamiltonian operator. It is represented by [9] with ψ mν = ψ m |e ν , and the latter equality being due to ∂ t e ν |e µ = 0. Eq. (5) is the key for the unitarity of the propagation, replacing the conventional expectation of hermiticity of iD µνt .

A. Conservation of scalar products
We start by showing the expected [19] unitarity of the evolution in the manifold, defined here as conservation of scalar products among the propagating states {|ψ m , m = 1 . . . N e } at any time. The evolution of the coefficients for the ket and bra then is determined by Eqs. (3) and (4), as It is easy to check that the scalar products ψ m |ψ n = ψ mµ ψ µ n are preserved in time: In addition to the Hamiltonian operator being hermitian, the unitarity of the propagation is therefore direct consequence of Eq. (5). The natural representation does not recover the usual self-adjoint matrix shape, but offers quite transparent relations and derivations. If the dealings above seem a bit of a sleight of hand, Appendix A reproduces the result in the matrix representation.

III. FINITE TIME STEP dt
After time discretisation for numerical integration, we are interested in propagating the set of coefficients for finite dt, trying to maximise both the quality and the stability of whatever the algorithm we use. Preservation of the orthonormality of the propagating states is key for that purpose. Here we will explore the behaviour of the Crank Nicolson algorithm in the Ξ curved manifold.
For the finite time-step dt we will neglect henceforth the time evolution of both the hamiltonian and the connection between t and t + dt. This is compatible with various integration algorithms such as extrapolating the Hamiltonian to t + dt/2. Other algorithms, such as the self-consistent Crank Nicolson [27] or self-consistent predictor-corrector schemes [28] may need further consideration.

A. Unitary propagation for static basis
Let us start with a reminder of unitary propagation when the basis set may be non-orthogonal but does not evolve, and, consequently, the state manifold is a regular Hilbert space Ω. It means that D µ νt = D µ tν = 0, S µν and S µν are constant, and the solutions for Eqs. (7) and (8) are Scalar products among |ψ m 's are preserved as the states evolve, as expected. Appendices B 1 a and B 1 b show the same in the matrix representation, and that the evolved bra in Eq. (9) remains the bra of the evolved ket, respectively.

Crank-Nicolson for a static basis
The approximate evolution given by is the direct generalisation (in the natural representation) of the usual Crank-Nicolson evolution step for orthonormal bases [9]. We have abstracted the notation, with the circles standing for indices, to be contracted with contiguous ones if full, and always up with down. Eq. (11) can also be recast in the matrix representation as follows where we use the fact that In the first equation we find the conventional expression in terms of the S −1 H matrix product, while the last one shows a variant that does not require overlap inversion. Although approximate in the evolution, it can be shown to be strictly unitary for finite dt (see Appendix B 2).

B. Constant connection
For a situation in which both the basis and the (tangent) Hilbert space Ω(t) do change with time, we consider now the case in which dt is finite but small enough so that the connection D µ νt can be taken as constant (the situation for varying basis set but within an invariant or converged Ω was contemplated in Ref. [9]).

Metric tensor evolution under constant connection
First, let us see how the metric tensors evolve between t and t + dt. In general, and still exact, If both D µ νt and D ν µt are taken as constant, given the premise of this Section, and given the fact that D ν µt = D ν * µt , the solution of Eq. (12) is which will be exact as long as those connections are strictly constant, but will represent an approximate solution for small dt but varying connections. Analogously, since both D µ tσ and D ν λ t are constant again.

Calculation of the connection
The results in Eq. (13) and Eq. (14) are important, not only for further algebraic manipulations, but because they represent consistency conditions for the evolution, and, to some extent, they define the connection. In explicit calculations, the overlap matrix is defined extrinsically at any time step, i.e., it does not arise from evolution, but given the positions of atoms at a given time step and given the basis set definition in the larger ambient Hilbert space, as , and where R µ (t) represents the position of the centre of the orbital at that time.
The connection itself can be computed extrinsically, normally as There are therefore two possibilities for an approximate evolution for finite dt. (i) The connection can be calculated as in Eq. (16) neglecting the small discrepancy in the evolved overlap S µν (t+dt) between the actual one, as in Eq. (15), and the one that would result from evolving under the calculated connection, as in Eq. (13). Alternatively, (ii) a connection can be proposed, instead of via Eq. (16), by construction to satisfy Eq. (13) from the overlaps calculated using Eq. (15) both at t and t + dt.
In this work we test option (i), since it is the most attractive numerically. The calculation of the connection via Eq. (16) can be done in linear-scaling operations. We have implemented it in the Siesta program [12,29,30], which uses finite-support atomic orbitals as basis sets, and contains a real-time TD-DFT implementation [31]. The integrals in Eq. (16) represent two-centre integrals which are evaluated by a trivial extension of the method explained in Section 5 of Ref. [29]. The implemented connection is tested below. Option (ii) is further explored in Appendix D, where a possible alternative direction towards better orthonormality preservation is outlined, albeit with worse scaling in computational expense as far as we can see.
C. Unitarity and convergence with time step

Exponential evolution
For constant D µ νt and H µ ν the coefficients for the ket and the bra evolve as given by the solution of Eqs. (7) and (8), namely, Checking again for unitarity, Since e A e B = e A+B when [A, B] = 0, then and unitarity of the propagation in Eqs. (17) and (18) is confirmed. Appendix E shows that the bra evolved as in Eq. (18) is the instantaneous bra of the ket in Eq. (17). In Appendix C the situation for parallel transport is presented for clarity.

Crank Nicolson
Unlike the static-basis case, the Crank Nicolson approximate evolution of Eqs. (17) and (18) is not unitary regardless of dt. In order to see this let us proceed as follows. The Crank Nicolson algorithm for both equations is expressed as It is easy to see (using appropriate commutation) that, as defined, and, therefore, scalar products would seem to be preserved exactly. It is not so, however. The evolved bra ψ m• (t + dt) is only approximately the bra of the evolved ket ψ • m (t+dt), and therefore unitarity will only be approximate. This is a curvature effect; Appendix C shows it to happen for parallel transport as well. The actual bra of the ket in Eq. (19) is obtained by turning it around and applying metric tensors as needed, as where we use 1 for t and 2 for t + dt. Using Eq. (13) for S••(2) does not convert Eq. (21) into Eq. (20). The unitarity of the evolution is only approximate. It is tested below.

Tests: collision of two He atoms and H across graphite
The mentioned non-unitary evolution is explicitly shown numerically in Fig. 1. For the two occupied states, |ψ 1 and |ψ 2 in a collision between two He atoms, the quantities ψ 1 |ψ 1 − 1, ψ 2 |ψ 2 − 1, ψ 1 |ψ 2 , and ψ 2 |ψ 1 are plotted versus time for various values of dt, using the Crank Nicolson algorithm proposed in Eq. (19). One He atom is kept fixed in space while another moves past it on a fixed trajectory with an impact parameter of 0.5 Å, with a fixed velocity of 1 atomic unit.
The calculations are performed using the Siesta program, with a double-ζ polarised basis set. All the technical settings of the calculations are as in Ref. [24], for a box size of 10 Å. It is apparent how the overlap matrix for the two evolving states deviates from the starting unit matrix, depending on the size of dt as expected from the discussion above.
Using as reference the TS algorithm [6], which is unitary by construction and was extensively used in firstprinciples electronic stopping power calculations [21][22][23][24][25], the deterioration of orthonormality of the GP integrator is assessed by comparing both for exactly the same process and approximations (using the same basis). Fig. 1(c) clearly shows the creeping in of deviation from orthonormality of the evolving states for the GP algorithm of Eq. (19) as compared with the strictly unitary TS alternative, which is limited only by the accuracy in the diagonalization involved. The significant deviation between 0.1 and 0.35 fs is over the period where the two atoms are close enough to interact, i.e., when the basis states associated to the different atoms overlap. The magnitude of these larger deviations depends directly on the value of dt, although the shape is identical, as can be seen by comparing Figs. 1(a) and (b), but they return rapidly to close to zero once the atoms are further apart, with the final deviation from 0 at 0.5 fs also depending on dt at a much smaller scale.
That deviation from unitary evolution is behind the demand for smaller dt of the algorithm apparent in Fig. 2, where the convergence in energy transfer for the collision is shown (difference in electronic energy between times before the collision and after the collision).
The convergence, however, depends on system and property. Fig. 3 shows an analogous plot for the electronic stopping power S e , electronic energy uptake per unit length traversed by a proton projectile travelling across graphite along the (0001) direction (the details about this simulation can be found in Ref. [24]). ure shows that this property is quite similarly converged with both integrators. It should be noted, however, that the difference in computational effort for both integrators is substantial. The GP integrator following Eq. (19), requires the calculation of the connection D µνt , which represents two-centre integrals that are pre-calculated as tables at the beginning of a simulation, which then are interpolated and rotated as needed. It represents a very small part of the Siesta run, as the other two-centre integrals, such as the overlap and kinetic energy matrices. The TS integrator requires a diagonalization of the overlap matrix for the whole basis at every time step. The difference increases with system size, since the effort to calculation of D µνt scales linearly with system size (due to the sparsity of the matrix -same as for the overlap) while the said diagonalization scales with the cube power.
It is important to finish, however, revisiting the accuracy of the converged integration with both algorithms. It was already pointed out in Ref. [9] that the TS integrator, although perfectly unitary and displaying good convergence, does not produce the correct integration for high velocities. It worked well for electronic stopping power calculations for low velocity projectiles [21][22][23][24][25]. This is confirmed in Fig. 4, where the electronic stopping power for a proton across graphite is depicted, comparing both integrators with the empirical PSTAR data [32]. S e is indeed well reproduced by both algorithms for velocities below 1 atomic unit.
However, the deviation at higher velocities for the TS integrator is very apparent, very clearly confirming the formal results of Ref. [9], with a large overestimation of the electronic stopping power, by a factor of two around the Bragg peak, growing to a tenfold overestimation for velocities around ten atomic units. The simpler system of the two He atom collision is quite illustrative. Fig. 5 shows the electronic energy as a function of position of the projectile He atom, as it passes by an immobile one. The TS integrator converges better with dt, as shown in Fig. 2, but quite a lot of the physics is lost. In this case the key difference stems from the start of the evolution, which is an abrupt kick of the projectile nucleus. It takes a time for the electrons around that nucleus to respond, and there is a lag. Since there is nothing else until reaching the target He atom, the electronic cloud around the projectile oscillates back and forth, which becomes a more complex behaviour when colliding. This is completely missed by the TS algorithm, since it transposes the coefficients of the evolving states from the (orthogonalised) basis at t to the one at t + dt, which implies an instantaneous response to the initial kick, without oscillation. Very smooth, nicely converged, and quite unphysical. Admittedly, it is an example particularly ill-suited for the algorithm, but illustrative, nevertheless.

Orthonormalisation correction
Of course, the Crank-Nicolson propagation step of GP, in Eq. (19), can be made strictly unitary by force if adding a Löwdin orthonormalisation step using S −1/2 (t+ dt) as obtained from the diagonalization of S mn (t + dt) = ψ m (t + dt)|ψ n (t + dt) .
Notice that the S matrix is of N × N dimensions, N being the number of occupied propagating states, much smaller than the number of basis states. In practical terms, and given the satisfactory unitarity achieved directly by Eq. (19) for small dt, one can choose to evolve using the uncorrected algorithm for some time, then evaluating the S matrix once every number steps n s , and, whenever max mn {|S mn − δ mn |} exceeds some tolerance , an orthonormalisation step would be performed as described. The algorithm will be then optimised by choosing the best combination of dt, n s and , which will depend on the system under study.
We have shown above how different problems have different demands for convergence, depending, for instance, on their evolution being dominated by the basis motion and the connection, or by the Hamiltonian itself, either the evolution of the external potential, or the Hamiltonian effective spectral width for the evolving states. The interplay of those parameters can therefore vary quite substantially. The method would be nicely completed with a learning algorithm to adjust those parameters dynamically.

IV. CONCLUSIONS
The effect of the time evolution of the basis set and of the Hilbert space it spans is explored for the description of the time evolution of quantum states. The exploration is both formal and numerical, assessing the deterioration of the conservation of scalar products of evolving states (unitarity of the evolution) and its implication for convergence with time step in the integration of the dynamical equation by time discretisation.
Formally, a Crank-Nicolson algorithm using the connection of the manifold as gauge potential (GP) is shown to keep unitarity only approximately, unlike the integra-tor proposed in Ref. [9] for a moving basis within a static Hilbert space, and unlike the TS algorithm [6], which is perfectly unitary regardless of dt size by construction and works well at low velocities [21][22][23][24][25], but was suggested to describe unphysical evolution at higher velocities.
The most numerically convenient GP integrator [that of Eq. (19)] is tried on two systems, a collision between two He atoms and the passing of a constant-velocity H projectile across graphite. The unitarity and dt convergence for that algorithm is compared with the TS one, the latter displaying the expected better convergence with dt, requiring from two to ten times less time steps for a given simulation time. That advantage is however offset by the more efficient (and better basis-size scaling) GP algorithm of Eq. (19), which only demands, per time step, the calculation of a sparse matrix of two-centre integrals (in linear-scaling operations) instead of the overlap matrix diagonalization of TS.
The deviation from the physical evolution of the TS integration is confirmed for nuclear (and basis function) velocities comparable to or larger than the valence electron velocities. A very significant overestimation of the Bragg peak (a factor of two in both height and position) is observed for the electronic stopping power of the proton shooting through graphite along the centre of a (0001) channel. And a very clear modification of the expected physics is observed for the same integrator when describing the two-atom collision.
Although other routes for unitary integrators are proposed in this work, only the approximately unitary GP algorithm of Eq. (19) is found to be satisfactory. The work also provides a better perspective in the understanding of unitary evolution of quantum states with evolving basis sets, in the context of curved manifolds. It can be checked more conventionally, using for the bra coefficients the complex conjugate of the ones for the ket, ψ µ m = ψ m |e µ = ψ µ * m , as follows: where a filled bullet has been introduced for every upper (lower) index that contracts with a contiguous lower (upper) index, and where we have used the fact that ∂ t S νµ = D µtν + D µνt .
Appendix B: Non-orthogonal static basis 1. Exponential solution a. In the matrix representation A more traditional proof than the one in Eq. (10) is presented here for the unitarity of the evolution of states in a static Hilbert space Ω, for a non-orthogonal static basis. It is well known, but it serves to set up the scene for later manipulations. The ψ µ m elements are actually the conventional expansion coefficients of |ψ m in the {|e µ } basis, as |ψ m = |e µ e µ |ψ m = |e µ ψ µ m . Scalar products remain constant: where filled bullets indicate contracted indices that are contracted as in Appendix A. In the fourth line, 1 = S•• S •• was introduced, and in the fifth line we made use of the following relationship The first equality is easily checked by expanding the exponential, and using the fact that where we have just complex-conjugated the whole equation, and changed the order of ψ and the exponential to reflect the index contraction (ν) as the order in a matrix product. The same equation can be re-expressed as or, multiplying by the overlap on the right, nor D µνt are antihermitian, and it is hard to see how and why they would give unitary propagation. The parallel transport case can be illustrative, i.e., the evolution of states for zero covariant derivative. In an Euclidean space the states would not change in time, but they do in a curved manifold. Because of the Schrödinger equation, zero covariant derivative implies H µ ν (t) = 0. Parallel transport for ket and bra would then be described by For constant D µ νt the coefficients for the ket and the bra evolve as given by the solution of Eq. C1, namely, In the Ω(t) manifold, the scalar-product-preserving propagation ψ m |ψ n (t + dt) = ψ m |ψ n (t) becomes Using Eq. C2, Since e dt D ν µt e −dtD µ δt = δ ν δ , the unitarity of propagation in Eq. (C3) is demonstrated. The key is in Eq. (5).

Evolved bra
We have shown that the scalar product between the evolved representation of the ket ψ µ m (t) ∈ Ω(t) and the evolved representation of its bra, ψ n ν (t), is preserved at later times. We can also show the preservation of the scalar product directly for the bra ψ nν (t + dt), corresponding to ψ ν n (t + dt) ∈ Ω(t + dt). In other words, the evolved bra is the bra of the evolved ket, i.e., we want to check that for the ket represented by the bra would be For the purpose let us write down Eq. C4 as e µ , t + dt|ψ m , t + dt = e −dt e µ |∂teν e ν , t|ψ m , t and turn it around (and complex conjugate it) where we have changed the order of factors in the r.h.s. to reflect the summed index ν. Now, and Eq. C6 becomes ψ mσ (t + dt) S σµ (t + dt) = ψ mλ (t) S λν (t)e −dtD µ νt .
coinciding with the evolved bra in Eq. C5, as expected.

Crank Nicolson for parallel transport
Above we just saw that since parallel transport preserves scalar products, it also does it in the case of a constant connection, no matter for how long. The question is now what happens with approximate evolution over a finite dt.
Since the evolution of both ket and bra under parallel transport with constant connection [Eq. (C2)] is mathematically analogous to their evolution under Eq. (9) with fixed basis, the Crank Nicolson algorithm can be used to approximate their transport, As done in Section III C 2, it is apparent (using appropriate commutation) that ψ m• (t + dt)ψ • n (t + dt) = ψ m• (t)ψ • n (t) , and, therefore, scalar products would seem to be preserved exactly. But again, they are not. The evolved bra ψ m• (t + dt) as in Eq. (C8) is only approximately the bra of the evolved ket ψ • m (t + dt), and therefore unitarity will only be approximate. The actual bra of the ket in Eq. (C7) is rather Using Eq. (13) for S••(t + dt) does not convert Eq. (C9) into Eq. (C8).
Appendix D: Alternative connection and integrator for finite dt When travelling between t and t + dt, S µν (t) and S µν (t + dt) can be calculated directly in ambient space using Eq. (15). Assuming a constant connection we could then use the D µνt arising from the solution of Eq. 13, instead of calculating it explicitly as in Eq. (16). However, the relation between overlaps alone should not be sufficient for the determination of the connection, since any rotation of the basis at t + dt should leave S µν (t + dt) unaltered while the transformation between t and t + dt would change.

Parallel transport transformation
Defining the transformation tensor, it performs the transformation from t to t + dt for any vector following parallel transport (see Appendix C), We could define A • • more generally as the one doing that operation regardless of D µ νt being constant or not. If not, Eq. (D1) would have to be replaced by a time-ordered integral, but once in possession of A µ ν , the rest of this section would be the same.
It is analogous to a basis set transformation, except that A µ ν = A µ ν ≡ e µ , t + dt | e ν , t , as defined in Ref. [9], which is the relevant transformation when Ω does not change with time.