QQUANTUM DYNAMICS WITH THE PARALLEL TRANSPORTGAUGE
DONG AN ∗ AND
LIN LIN † Abstract.
The dynamics of a closed quantum system is often studied with the direct evolutionof the Schr¨odinger equation. In this paper, we propose that the gauge choice ( i.e. degrees of freedomirrelevant to physical observables) of the Schr¨odinger equation can be generally non-optimal fornumerical simulation. This can limit, and in some cases severely limit the time step size. We findthat the optimal gauge choice is given by a parallel transport formulation. This parallel transportdynamics can be simply interpreted as the dynamics driven by the residual vectors, analogous tothose defined in eigenvalue problems in the time-independent setup. The parallel transport dynamicscan be derived from a Hamiltonian structure, thus suitable to be solved using a symplectic andimplicit time discretization scheme, such as the implicit midpoint rule, which allows the usage ofa large time step and ensures the long time numerical stability. We analyze the parallel transportdynamics in the context of the singularly perturbed linear Schr¨odinger equation, and demonstrate itssuperior performance in the near adiabatic regime. We demonstrate the effectiveness of our methodusing numerical results for linear and nonlinear Schr¨odinger equations, as well as the time-dependentdensity functional theory (TDDFT) calculations for electrons in a benzene molecule driven by anultrashort laser pulse.
Key words.
Schr¨odinger equation; Quantum dynamics; Gauge; Parallel transport; Densitymatrix; von Neumann equation; Symplectic method; Singularly perturbed system; Time-dependentdensity functional theory; Adiabatic theorem
1. Introduction.
Consider the following set of coupled nonlinear Schr¨odingerequations i (cid:15)∂ t Ψ( t ) = H ( t, P )Ψ( t ) . (1.1)Here we assume 0 < (cid:15) (cid:28)
1. Ψ( t ) = [ ψ ( t ) , . . . , ψ N ( t )] are N time-dependent wavefunctions subject to suitable initial and boundary conditions. H ( t, P ) is a self-adjointtime-dependent Hamiltonian. P ( t ) is called the density matrix and defined as P ( t ) = Ψ( t )Ψ ∗ ( t ) = N (cid:88) j =1 ψ j ( t ) ψ ∗ j ( t ) . (1.2)Note that when the initial state Ψ(0) consists of N orthonormal functions, the func-tions in Ψ( t ) will remain orthonormal for all t , i.e. ( ψ i ( t ) , ψ j ( t )) = δ ij , where ( · , · )denotes a suitable inner product. Then P ( t ) = N (cid:88) j,k =1 ψ j ( t )( ψ j ( t ) , ψ k ( t )) ψ ∗ k ( t ) = N (cid:88) j =1 ψ j ( t ) ψ ∗ j ( t ) = P ( t ) , i.e. P ( t ) is a projector. The explicit dependence of the Hamiltonian on t is oftendue to the existence of an external field, and we assume the partial derivatives ∂ m H∂t m are of O (1) in some suitable norms for all m ≥
1. Hence when 0 < (cid:15) (cid:28)
1, the wave ∗ Department of Mathematics, University of California, Berkeley, Berkeley, CA 94720. Email: dong [email protected] † Department of Mathematics, University of California, Berkeley, Berkeley, CA 94720 and Com-putational Research Division, Lawrence Berkeley National Laboratory, Berkeley, CA 94720. Email: [email protected] a r X i v : . [ m a t h . NA ] F e b unctions can oscillate on a much smaller time scale than that of the external fields,and this is called the singularly perturbed regime [13].The equations (1.1) are rather general and appear in several fields of scientificcomputation. In the simplest setup when N = 1 and H ( t, P ) ≡ H ( t ), this is the linearSchr¨odinger equation. Another example is the nonlinear Schr¨odinger equation (NLSE)used for modeling nonlinear photonics and Bose-Einstein condensation process [10],i (cid:15)∂ t ψ ( t ) = H ( t ) ψ ( t ) + g | ψ ( t ) | ψ ( t ) , (1.3)where H ( t ) is a Hermitian matrix obtained by discretizing the linear operator − ∆+ V ( x, t ). Since N = 1, P ( t ) = ψ ( t ) ψ ∗ ( t ), and | ψ ( t ) | = diag[ P ( t )] is a nonlinear localpotential. When N >
1, the coupled set of Schr¨odinger equations must be solvedsimultaneously. This is the case in the time-dependent density functional theory(TDDFT) [34, 32].The simulation of Eq. (1.1) and in particular (1.3) has been studied via a widerange of numerical discretization methods, such as explicit Runge-Kutta methods [36],implicit Runge-Kutta methods [6], operator splitting methods [3, 27], Magnus expan-sion methods [6, 7], exponential time differencing methods [21], spectral deferredcorrection methods [19], dynamical low rank approximation [24], adiabatic state ex-pansion [18, 39], to name a few. What this paper focuses on is not to develop anothernumerical scheme to directly discretize (1.1), but to propose an alternative formula-tion that is equivalent to (1.1), and can be solved with improved numerical efficiencyusing existing discretization schemes.More specifically, note that if we multiply Ψ( t ) by a time-dependent unitarymatrix U ( t ) ∈ C N × N , the resulting set of rotated wave functions, denoted by Φ( t ) =Ψ( t ) U ( t ), yields the same density matrix as P ( t ) = Φ( t )Φ ∗ ( t ) = Ψ( t ) [ U ( t ) U ∗ ( t )] Ψ ∗ ( t ) = Ψ( t )Ψ ∗ ( t ) . (1.4)Since the unitary rotation matrix U ( t ) is irrelevant to the density matrix which isused to represent many physical observables, U ( t ) is called the gauge, and Eq. (1.4)indicates the density matrix is gauge-invariant . Furthermore, Eq. (1.1) can be directlywritten in terms of the density matrix asi (cid:15)∂ t P ( t ) = [ H ( t, P ) , P ( t )] , (1.5)where [ H, P ] := HP − P H is the commutator between H and P . Eq. (1.5) is calledthe von Neumann equation (or quantum Liouville equation), which can be viewedas a more intrinsic representation of quantum dynamics since the gauge degrees offreedom are eliminated completely.The simulation of the von Neumann equation can also be advantageous from theperspective of time discretization. Consider the simplified scenario that H ( t, P ) ≡ H ( P ) does not explicitly depend on t , and the initial state Ψ(0) consists of a set ofeigenfunctions of H , i.e. H [ P ] ψ j (0) = ψ j (0) λ j (0) , j = 1 , . . . , N, P = N (cid:88) j =1 ψ j (0) ψ ∗ j (0) . (1.6)Eq. (1.6) is a set of nonlinear eigenvalue equations. When solved self-consistently, thesolution to the Schr¨odinger equation (1.1) has an analytic form ψ j ( t ) = exp (cid:18) − i (cid:15) λ j (0) t (cid:19) ψ j (0) , j = 1 , . . . , N, (1.7) hich oscillates on the O ( (cid:15) ) time scale. Hence many numerical schemes still need toresolve the dynamics with a time step of O ( (cid:15) ). On the other hand, the right handside of the von Neumann equation vanishes for all t , and hence nominally can bediscretized with an arbitrarily large time step! Of course one can use techniques suchas integration factors [8] to make this simulation using the Schr¨odinger equation asefficient. However this example illustrates that the gap in terms of the size of the timestep generally exists between the Schr¨odinger representation and the von Neumannrepresentation.In this paper, we identify that such gap is solely due to the gauge degrees offreedom in the Schr¨odinger representation. By optimizing the gauge choice, one canpropagate the wave functions using a time step comparable to that of the von Neu-mann equation. We demonstrate that the optimized gauge is given by a paralleltransport (PT) formulation. We refer to this gauge as the parallel transport gauge,and the resulting dynamics as the parallel transport dynamics. Correspondingly thetrivial gauge U ( t ) ≡ I N in Eq. (1.1) is referred to as the Schr¨odinger gauge, and theresulting dynamics as the Schr¨odinger dynamics. We remark that the PT dynamicscan also be interpreted as an analytic and optimal way of performing the dynamicallow rank approximation [24] for Eq. (1.1). Note that the simulation of the von Neu-mann equation requires the explicit operation on the density matrix P ( t ). When alarge basis set such as finite elements or planewaves is used to discretize the partialdifferential equation, the storage cost of P ( t ) can be often prohibitively expensivecompared to that of the wave functions Ψ( t ). Hence the PT dynamics combines theadvantages of both approaches, namely to perform simulation using the time step sizeof the von Neumann equation, but with cost comparable to that of the Schr¨odingerequation.We analyze the effectiveness of the PT dynamics for the linear time-dependentSchr¨odinger equation in the near adiabatic regime. We remark that efficient numericalmethods have been recently developed in this regime based on the construction ofa set of instantaneous adiabatic states [18, 39]. The assumption is that the wavefunctions can be approximated by the subspace spanned by low energy eigenstates ofthe Hamiltonian at each t . The dimension of the subspace is often chosen to be cN ,where c is a relatively small constant. Compared to these methods, the PT dynamicsalways operates only on N wave functions, and therefore has reduced computationaland the storage cost. The PT dynamics is also applicable beyond the near adiabaticregime.By extending the quantum adiabatic theorem [29, 2] to the PT dynamics, weprove that the local truncation error of the PT dynamics gains an extra order ofaccuracy in terms of (cid:15) , when the time step is O ( (cid:15) ) or smaller. The PT dynamics,after a slight modification, can be derived from a Hamiltonian system similar to thatin the Schr¨odinger dynamics. Hence the gain of accuracy for the local truncation errorcan be directly translated to the global error as well for long time simulation.We demonstrate the effectiveness of the PT dynamics using numerical resultsof the model linear and nonlinear Schr¨odinger equations. We also perform time-dependent density functional theory (TDDFT) calculations for the electrons in abenzene molecule driven by an ultrashort laser pulse, near and beyond the adiabaticregime. When the spectral radius of the Hamiltonian is large, it is suitable to discretizethe PT dynamics using a symplectic and implicit time discretization scheme, such asthe implicit midpoint rule, and the resulting scheme can significantly outperform thesame scheme for the Schr¨odinger dynamics. We also find that other time-reversible nd implicit time discretization schemes, such as the Crank-Nicolson scheme, canyield similar performance as well. Numerical results confirm our analysis in the nearadiabatic regime, and indicate that the convergence of the PT dynamics can startwhen the time step size is much larger than O ( (cid:15) ). This is in contrast to the Schr¨odingerdynamics where the error stays flat until the time step reaches below O ( (cid:15) ). ForTDDFT calculations, we find that our discretized PT dynamics can achieve 31 . .
2. Parallel Transport Gauge.
Since the concept of the parallel transportgauge is associated with the time propagation instead of spatial discretization, forsimplicity of the presentation, unless otherwise specified, we assume that Eq. (1.1)represents a discrete, finite dimensional quantum system, i.e. for a given time t , ψ j ( t )is a finite dimensional vector, and H ( t, P ) is a finite dimensional matrix. If the quan-tum system is spatially continuous, we may first find a set of orthonormal bases func-tions { e j ( r ) } dj =1 satisfying (cid:82) e ∗ j ( r ) e j (cid:48) ( r ) d r = δ jj (cid:48) , and expand the continuous wave-function as (cid:101) ψ j ( r , t ) ≈ (cid:80) dj =1 ψ j ( t ) e j ( r ). Then after a Galerkin projection, Eq. (1.1)becomes a d -dimensional quantum system, and the inner product for the coefficients ψ j ( t ) becomes the standard (cid:96) -inner product as ( ψ j ( t ) , ψ k ( t )) := ψ ∗ j ( t ) ψ k ( t ) = δ jk .Hence we can use the linear algebra notation. The star notation is interpreted asthe complex conjugation when applied to a scalar, and Hermitian conjugation whenapplied to a vector or a matrix. For simplicity let us consider the case N = 1 first, where thegauge matrix U ( t ) simply becomes a phase factor c ( t ) ∈ C , | c ( t ) | = 1. Note that thegauge choice cannot affect physical observables such as the density matrix. Henceconceptually we may think that the time-dependent density matrix P ( t ) has alreadybeen obtained as the solution of the von Neumann equation (1.5) on some time interval[0 , T ]. Similarly the wave function ψ ( t ) satisfying the Schr¨odinger dynamics is alsoknown. Then the relation P ( t ) ϕ ( t ) = ϕ ( t ) , ϕ ( t ) = ψ ( t ) c ( t ) (2.1)is satisfied for any gauge choice. For simplicity we use the notation ˙ ϕ ( t ) = ∂ t ϕ ( t ), anddrop the explicit t -dependence in all quantities, as well as the P -dependence in theHamiltonian unless otherwise noted. Our goal is to find the time-dependent gaugefactor c ( t ) so that the rotated wave function ϕ ( t ) varies as slowly as possible . Thisgives rise to the following minimization problem,min c ( t ) (cid:107) ˙ ϕ ( t ) (cid:107) s.t. ϕ ( t ) = ψ ( t ) c ( t ) , | c ( t ) | = 1 . (2.2)In order to solve (2.2), note that P ( t ) is a projector, we split ˙ ϕ into two orthogonalcomponents, ˙ ϕ = P ˙ ϕ + ( I − P ) ˙ ϕ. (2.3) y taking the time derivative with respect to both sides of the first equation inEq. (2.1), we have ( I − P ) ˙ ϕ = ˙ P ϕ. (2.4)Then (cid:107) ˙ ϕ (cid:107) = (cid:107) P ˙ ϕ (cid:107) + (cid:107) ( I − P ) ˙ ϕ (cid:107) = (cid:107) P ˙ ϕ (cid:107) + (cid:107) ˙ P ϕ (cid:107) = (cid:107) P ˙ ϕ (cid:107) + (cid:107) ˙ P ψ (cid:107) . (2.5)In the last equality, we have used that | c ( t ) | = 1. Note that the term (cid:107) ˙ P ψ (cid:107) isindependent of the gauge choice, so (cid:107) ˙ ϕ (cid:107) is minimized when P ˙ ϕ = 0 . (2.6)Therefore instead of writing down the minimizer of Eq. (2.2) directly, we define thegauge implicitly through Eq. (2.6).Let us write down an equation for ϕ ( t ) directly. Combining equations (2.4), (2.6),(1.5) and (2.1), we have˙ ϕ = ˙ P ϕ = 1i (cid:15) [ H, P ] ϕ = 1i (cid:15) ( Hϕ − ϕ ( ϕ ∗ Hϕ )) , (2.7)or equivalently i (cid:15)∂ t ϕ = Hϕ − ϕ ( ϕ ∗ Hϕ ) . (2.8)For reasons that will become clear shortly, we refer to this gauge choice as the parallel transport gauge , and Eq. (2.8) as the parallel transport (PT) dynamics. Com-paring with the Schr¨odinger dynamics, we find that the PT dynamics only introducesone extra term ϕ ( ϕ ∗ Hϕ ). The right hand side of Eq. (2.8) takes the form of theresidual vector in the solution of eigenvalue problem of the form (1.6). Hence the PTdynamics can be simply interpreted as the dynamics driven by the residuals. There-fore we expect that the PT dynamics can be particularly advantageous in the nearadiabatic regime [18, 39], i.e. when ϕ is close to be the eigenstate of H , and all theresidual vectors are therefore small.Now we provide an alternative interpretation of the gauge choice using the paralleltransport formulation associated with a family of projectors. For simplicity let usassume H ( t ) is already discretized into a finite dimensional Hermitian matrix for each t and so is P ( t ). Given the single parameter family of projectors { P ( t ) } defined onsome interval [0 , T ], we define A ( t ) = i (cid:15) [ ∂ t P ( t ) , P ( t )] . (2.9)It can be directly verified that A ( t ) is a Hermitian matrix for each t , and induces adynamics i (cid:15)∂ t T ( t ) = A ( t ) T ( t ) , T (0) = I. (2.10) T ( t ) is a unitary matrix for each t . T ( t ) is called the parallel transport evolutionoperator (see e.g. [28, 9]). The connection between the parallel transport dynamicsand the parallel transport evolution operator is given in Proposition 1. Proposition 1.
Define ϕ ( t ) = T ( t ) ψ (0) where T ( t ) is the evolution operatorsatisfying (2.10) , and P ( t ) satisfies the von Neumann equation (1.5) . Then P ( t ) = ϕ ( t ) ϕ ∗ ( t ) , and ϕ ( t ) satisfies the parallel transport dynamics (2.8) . roof . First we prove the following relation P ( t ) T ( t ) = T ( t ) P (0) (2.11)by showing that both sides solve the same initial value problem. Note that T ( t ) P (0)satisfies a differential equation of the form (2.10) with the initial value T (0) P (0).We would like to derive the differential equation P ( t ) T ( t ) satisfies. Taking the timederivative on both sides of the identity P ( t ) = P ( t ), we yield two useful relations˙ P = ˙ P P + P ˙ P , P ˙ P P = 0 . (2.12)Then using Eq. (2.10),i (cid:15)∂ t ( P T ) = i (cid:15) ˙ P T + i (cid:15)P [ ˙ P , P ] T = i (cid:15) ˙ P P T . On the other hand, A ( P T ) = i (cid:15) ( ˙ P P P
T − P ˙ P P T ) = i (cid:15) ˙ P P T . Therefore i (cid:15)∂ t ( P T ) = A ( P T ) . (2.13)Hence P T also satisfies an equation of the form (2.10). This proves Eq. (2.11) bynoticing further the shared initial condition P (0) T (0) = T (0) P (0).Using Eq. (2.11), we have P ( t ) ϕ ( t ) = P ( t ) T ( t ) ψ (0) = T ( t ) P (0) ψ (0) = T ( t ) ψ (0) = ϕ ( t ) . (2.14)Since T ( t ) is unitary, we have (cid:107) ϕ ( t ) (cid:107) = 1 for all t . Hence P ( t ) = ϕ ( t ) ϕ ∗ ( t ) . (2.15)The only thing left is to show that the gauge choice in ϕ ( t ) is indeed the paralleltransport gauge. Using Eq. (2.11) and (2.13), we havei (cid:15)∂ t ϕ = i (cid:15)∂ t ( T ψ (0)) = i (cid:15)∂ t ( P T ) ψ (0) = i (cid:15) ˙ P P T ψ (0) = HP ϕ − P HP ϕ. (2.16)Here we have used the von Neumann equationi (cid:15) ˙ P = HP − P H.
Finally using Eq. (2.14) and (2.15), we havei (cid:15)∂ t ϕ = Hϕ − ϕ ( ϕ ∗ Hϕ ) , which is precisely the parallel transport dynamics.In order to see why the parallel transport gauge can be more advantageous, con-sider again the time-independent example (1.6) in the introduction for the case N = 1.We find that the right hand side of Eq. (2.8) vanishes, and the solution is simply ϕ ( t ) = ϕ (0) = ψ (0)for all t . This implies that the parallel transport gauge is c ( t ) = exp (cid:0) + i (cid:15) λ (0) t (cid:1) thatperfectly cancels with the rotating factor in (1.7). Hence the PT dynamics yields the a) wave functions (b) Evolution of centers Fig. 2.1: (a) Real parts of the wave functions at x = 25 with the Schr¨odinger gaugeand the PT gauge, respectively. (b) Centers of the wave functions. Parameters arechosen to be T = 1 , (cid:15) = 0 . h = 10 − .slowest possible dynamics by completely eliminating the time-dependent phase factor,and the time step for propagating the PT dynamics can be chosen to be arbitrarilylarge as in the case of the von Neumann equation.For a more complex example, consider a time-dependent nonlinear Schr¨odingerequation in one dimension to be further illustrated in Section 5. Fig. 2.1 (a) showsthe evolution of the real part of the solution ψ ( t ) from the Schr¨odinger dynamics, andthat of ϕ ( t ) from the PT dynamics, respectively. We find that the trajectory of ϕ ( t )varies considerably slower than that of ψ ( t ), which allows us to use a much largertime step for the simulation. Fig. 2.1 (b) measures the accuracy of the average of theorbital center (cid:104) x (cid:105) ( t ), using simulation with the implicit midpoint rule, also known asthe Gauss-Legendre method of order 2 (GL2) scheme. We compare the performanceof the GL2 scheme with the Schr¨odinger gauge (S-GL2) and that with the PT gauge(PT-GL2) with the same step size h = 0 . h = 10 − . We observe that the solution from PT-GL2agrees very well with the reference solution, while the phase error of the solution fromS-GL2 becomes noticeable already after t = 0 . For simplicity let us consider the linear Schr¨odingerequation, i.e. H ( t, P ) ≡ H ( t ), and assume H ( t ) is a real symmetric matrix for all t .It is well known that the Schr¨odinger dynamics is a Hamiltonian system [30, 31, 11].More specifically, we separate the solution ψ into its real and imaginary parts as ψ = q + i p. (2.17)The (cid:96) -inner product associated with real quantities such as p, q are denoted by( p, q ) := p T q . We also introduce the canonically conjugate pair of variables ( τ, E )to eliminate the explicit dependence of H ( t ) on time [5, 11]. This gives the followingenergy functional E ( τ, q, E, p ) = 12 (cid:15) (cid:2) q T H ( τ ) q + p T H ( τ ) p (cid:3) + E. (2.18) he Hamiltonian system corresponding to this energy functional is ∂ t τ = ∂ E ∂E = 1 ,∂ t q = ∂ E ∂p = 1 (cid:15) H ( τ ) p,∂ t E = − ∂ E ∂τ = − (cid:15) (cid:20) q T ∂H ( τ ) ∂τ q + p T ∂H ( τ ) ∂τ p (cid:21) ,∂ t p = − ∂ E ∂q = − (cid:15) H ( τ ) q. (2.19)Hence τ is simply the time variable, and − E is the usually defined energy of the systemup to a constant. By combining the equations for q, p we obtain the Schr¨odingerdynamics for ψ .Although the PT dynamics only differs from the Schr¨odinger dynamics by thechoice of the gauge, interestingly, the PT dynamics cannot be directly written as aHamiltonian system. To illustrate this, we first separate the real and imaginary partsof ϕ as in (2.17), and the PT dynamics can be written as ∂ t q = 1 (cid:15) ( Hp − ( q T Hq + p T Hp ) p ) ,∂ t p = 1 (cid:15) ( − Hq + ( q T Hq + p T Hp ) q ) . (2.20)If this dynamics can be derived from some energy functional E , then ∂ E ∂p = 1 (cid:15) ( Hp − ( q T Hq + p T Hp ) p ) ,∂ E ∂q = 1 (cid:15) ( Hq − ( q T Hq + p T Hp ) q ) . (2.21)Straightforward computation reveals that ∂ E ∂p∂q = ∂ E ∂q∂p is not true in general, andhence the PT dynamics (2.8) cannot be a Hamiltonian system.Fortunately, the PT dynamics can be slightly modified to become a Hamiltoniansystem. Consider the following modified energy functional E ( τ, q, E, p ) = 12 (cid:15) ( q T H ( τ ) q + p T H ( τ ) p )(2 − q T q − p T p ) + E. (2.22)The corresponding Hamiltonian equations are ∂ t τ = ∂ E ∂E = 1 ,∂ t q = ∂ E ∂p = 1 (cid:15) (cid:2) H ( τ ) p (2 − q T q − p T p ) − ( q T H ( τ ) q + p T H ( τ ) p ) p (cid:3) ,∂ t E = − ∂ E ∂τ ,∂ t p = − ∂ E ∂q = 1 (cid:15) (cid:2) − H ( τ ) q (2 − q T q − p T p ) + ( q T H ( τ ) q + p T H ( τ ) p ) q (cid:3) . (2.23)Again τ is the same as t , and the conjugate variable E ( t ) satisfies E ( t ) = − (cid:15) ( q T H ( t ) q + p T H ( t ) p )(2 − q T q − p T p ) + constant . ompared to the PT dynamics (2.20), we have an extra factor (2 − q T q − p T p ) in theequations and the energy. Proposition 2 states that the solution to the PT dynam-ics (2.20) is the same as the solution of the Hamiltonian system (2.23). Proposition 2. If ( τ, q, E, p ) solves the Hamiltonian system (2.23) with normal-ized initial value condition p T (0) p (0) + q T (0) q (0) = 1 , then ( q ( t ) , p ( t )) solves (2.20) with the same initial value condition, and ϕ ( t ) = q ( t ) + i p ( t ) solves the PT dynam-ics (2.8) .Proof . Comparing Eq. (2.23) with Eq. (2.20), we only need to show the identity p T p + q T q = 1holds for all t . By computing ddt ( p T p + q T q ) =2( p T ∂ t p + q T ∂ t q )= 1 (cid:15) ( − − q T q − p T p ) p T Hq + 2( q T Hq + p T Hp ) p T q + 2(2 − q T q − p T p ) q T Hp − q T Hq + p T Hp ) q T p ) = 0 , we find that p T p + q T q is invariant during the propagation. Together with the nor-malized initial condition, we complete the proof.Proposition 2 suggests that the Hamiltonian form of the PT dynamics isi (cid:15)∂ t ϕ = Hϕ (2 − ϕ ∗ ϕ ) − ϕ ( ϕ ∗ Hϕ ) , (2.24)which shares exactly the same solution with (2.8) using the condition ϕ ∗ ϕ = 1.At the end of this part, we briefly discuss the Hamiltonian structure of the non-linear Schr¨odinger equation and the associated PT dynamics. Let us consider thediscretized nonlinear Schr¨odinger equation (1.3), which can be reformulated as aHamiltonian system driven by the energy functional E ( τ, q, E, p ) = 12 (cid:15) (cid:104) q T H ( τ ) q + p T H ( τ ) p + g | q | + | p | ) ) (cid:105) + E. (2.25)The PT dynamics corresponding to Eq. (1.3) can be written asi (cid:15)∂ t ϕ = H ϕ + g | ϕ | ϕ − ϕ ( ϕ ∗ H ϕ ) − gϕ ( ϕ ∗ | ϕ | ϕ ) . (2.26)Similar to the linear case, the PT dynamics itself cannot be reformulated as a Hamilto-nian system in general, but can be slightly modified to become a Hamiltonian system.More precisely, define the energy functional E ( τ, q, E, p ) = 12 (cid:15) (cid:2) q T H ( τ ) q + p T H ( τ ) p + g Tr(( | q | + | p | ) ) (cid:3) (2 − q T q − p T p ) − g (cid:15) Tr(( | q | + | p | ) ) + E, (2.27)then the Hamiltonian system driven by this energy functional can be written asi (cid:15)∂ t ϕ = ( H ϕ + 2 g | ϕ | ϕ )(2 − ϕ ∗ ϕ ) − ϕ ( ϕ ∗ H ϕ ) − gϕ ( ϕ ∗ | ϕ | ϕ ) − g | ϕ | ϕ. (2.28)Again this equation shares the same solution with Eq. (2.26) using the condition ϕ ∗ ϕ = 1. .3. General case. The PT dynamics derived in the previous sections can bedirectly generalized to Eq. (1.1) with
N >
1. Define the transformed set of wave func-tions Φ( t ) = Ψ( t ) U ( t ) = [ ϕ ( t ) , . . . , ϕ N ( t )], where U ( t ) ∈ C N × N is a gauge matrix.Following the same derivation in Section 2.1, we find that the parallel transport gaugeis given by the condition P ˙Φ = 0 . (2.29)This gives rise to the following PT dynamicsi (cid:15)∂ t Φ( t ) = H ( t, P ( t ))Φ( t ) − Φ( t )[Φ ∗ ( t ) H ( t, P ( t ))Φ( t )] , P ( t ) = Φ( t )Φ ∗ ( t ) . (2.30)Again the PT dynamics is driven by the residual vectors as in eigenvalue problems.In addition, the Hamiltonian structure is also preserved for the PT dynamics.For simplicity let us consider the linear Hamiltonian H ( t ). We separate the set of PTwave functions Φ into real and imaginary parts asΦ( t ) = q ( t ) + i p ( t ) . Define the energy functional E ( τ, q , E, p ) = 12 (cid:15) Tr (cid:16) ( q T H ( τ ) q + p T H ( τ ) p )(2 I N − q T q − p T p ) (cid:17) + E. (2.31)The associated Hamiltonian system is ∂ t τ = ∂ E ∂E = 1 ,∂ t q = ∂ E ∂ p = 1 (cid:15) ( H ( τ ) p (2 I N − q T q − p T p ) − p ( q T H ( τ ) q + p T H ( τ ) p )) ,∂ t E = − ∂ E ∂τ ,∂ t p = − ∂ E ∂ q = 1 (cid:15) ( − H ( τ ) q (2 I N − q T q − p T p ) + q ( q T H ( τ ) q + p T H ( τ ) p )) . (2.32)Similar with the case when N = 1 (Proposition 2), we can show that p T p + q T q = I N provided the orthonormal initial value condition. Therefore the solution to the Hamil-tonian system (2.32) can exactly form a set of solutions to the PT dynamics.Due to the straightforward generalization as described above, unless otherwisenoted, we will focus on the case N = 1 for the rest of the paper.
3. Time discretization.
When the spectral radius of the Hamiltonian is rel-atively small and (cid:15) ∼ O (1), explicit time integrators such as the 4th order Runge-Kutta method (RK4) and the Strang splitting method can be very efficient, and canbe applied to both the Schr¨odinger dynamics and the PT dynamics. However, theadvantage of propagating the PT dynamics can become clearer when (cid:15) becomes smallor when the spectral radius of H becomes very large, which is typical in e.g. TDDFTcalculations. In this scenario, all explicit time integrators must take a very smalltime step, which may become very costly. It should be noted that in the Schr¨odingerdynamics, the solution often oscillates rapidly on the time scale of (cid:15) as indicated in q. (1.7). Standard implicit discretization schemes, such as the implicit midpoint ruleand the Crank-Nicolson scheme, aim at interpolating such rapidly moving curves bylow order polynomials. Therefore the time step must still be kept on the order of (cid:15) tomeet the accuracy requirement, even though the numerical scheme itself may have alarge stability region or even A-stable [12].On the other hand, as discussed in Section 2.1, the PT dynamics transforms thefast oscillating wave function ψ ( t ) into a potentially slowly oscillating wave function ϕ ( t ) (as in Fig. 2.1 (a)). This makes it feasible to approximate ϕ ( t ) using a low orderpolynomial approximation. This statement will be further quantified by numericalresults in Section 5. Combined with an implicit time discretization scheme with alarge stability region, we may expect that the PT dynamics can be discretized with amuch larger time step than that in the Schr¨odinger dynamics.The Hamiltonian structure of the PT dynamics further invites the usage of a sym-plectic scheme for achieving long time accuracy and stability. The simplest symplecticand implicit scheme is the implicit mid-point rule, also known as the Gauss-Legendremethod of order 2 (GL2). We use a uniform time discretization t n = nh , and h is thetime step size. With some abuse of notations, we denote by ϕ ( t n ) the exact solutionat t n , and ϕ n the numerical approximation to ϕ ( t n ). Correspondingly we define P n = ϕ n ϕ ∗ n , H n = H ( t n , P n ) . It would also be helpful to define the effective nonlinear Hamiltonian H e ( t, ϕ ) as H e = H (2 − ϕ ∗ ϕ ) − ( ϕ ∗ Hϕ ) I, for Eq. (2.24), H e = ( H + 2 g | ϕ | )(2 − ϕ ∗ ϕ ) − ( ϕ ∗ H ϕ ) I − g ( ϕ ∗ | ϕ | ϕ ) I − g | ϕ | , for Eq. (2.28).Then the Hamiltonian equations (2.24) and (2.28) can be written in a uniform formi (cid:15)∂ t ϕ = H e ϕ. (3.1)The PT-Ham-GL2 discretization for discretizing the Hamiltonian equation (2.24)and (2.28) therefore becomes ϕ n +1 = ϕ n + h i (cid:15) H en + (cid:101) ϕ, (cid:101) ϕ = 12 ( ϕ n + ϕ n +1 ) , (3.2)Here (cid:101) ϕ can be interpreted as the approximation to ϕ ( t n + ) at the half time step, and H en + := H e ( t n + , (cid:101) ϕ ) . Note that the normalization condition (cid:101) ϕ ∗ (cid:101) ϕ → h →
0, but (cid:101) ϕ ∗ (cid:101) ϕ (cid:54) = 1 in general. Eq. (3.2) is a set of nonlinear equations for ϕ n +1 , and need to besolved iteratively. This can be viewed as a fixed point problem of the form ϕ = F ( ϕ ) , where the mapping F is explicitly defined as F ( ϕ ) = ϕ n + h i (cid:15) H en + (cid:101) ϕ, (cid:101) ϕ = 12 ( ϕ n + ϕ ) . (3.3) ssuming the fixed point exists and is unique, we may associate ϕ n +1 with the fixedpoint, and then move to the next time step. We may use any nonlinear equationsolving technique to solve such fixed point problem [23]. In this work, we use theAnderson mixing [1] method, which is a simplified Broyden-type method widely usedin electronic structure calculations [26].The PT-Ham-GL2 scheme can be simplified by directly applying the GL2 dis-cretization to the PT dynamics (2.8) and (2.26), with the efficient Hamiltonians tobe defined as H e = H − ( ϕ ∗ Hϕ ) I, for Eq. (2.8), H e = H + g | ϕ | − ( ϕ ∗ H ϕ ) I − g ( ϕ ∗ | ϕ | ϕ ) I, for Eq. (2.26).Again note that, unlike the continuous case, PT-GL2 is not equivalent to PT-Ham-GL2 since (cid:101) ϕ ∗ (cid:101) ϕ (cid:54) = 1 in general. Nevertheless, the norm of the numerical solutionsobtained by GL2 at the discretized time points t n are indeed conserved, which issummarized in the following proposition. Proposition 3.
Suppose ϕ n is the numerical solution obtained by applying GL2to one of the following PT dynamics, (2.24) , (2.28) , (2.8) and (2.26) . Assume that I − h (cid:15) H en + is always invertible in each step, then (cid:107) ϕ n (cid:107) = (cid:107) ϕ (cid:107) .Proof . We consider the GL2 scheme (3.2) for the uniform form (3.1). It sufficesto prove that (cid:107) ϕ n +1 (cid:107) = (cid:107) ϕ n (cid:107) for any n . We first substitute (cid:101) ϕ by ( ϕ n + ϕ n +1 ) andrewrite GL2 as (cid:18) I − h (cid:15) H en + (cid:19) ϕ n +1 = (cid:18) I + h (cid:15) H en + (cid:19) ϕ n . Note that for all defined H e , H e ∗ = H e , then ϕ ∗ n +1 ϕ n +1 = ϕ ∗ n (cid:18) I + h (cid:15) H en + (cid:19) ∗ (cid:18) I − h (cid:15) H en + (cid:19) ∗− (cid:18) I − h (cid:15) H en + (cid:19) − (cid:18) I + h (cid:15) H en + (cid:19) ϕ n = ϕ ∗ n (cid:18) I − h (cid:15) H en + (cid:19) (cid:18) I + h (cid:15) H en + (cid:19) − (cid:18) I − h (cid:15) H en + (cid:19) − (cid:18) I + h (cid:15) H en + (cid:19) ϕ n = ϕ ∗ n (cid:18) I − h (cid:15) H en + (cid:19) (cid:18) I − h (cid:15) H en + (cid:19) − (cid:18) I + h (cid:15) H en + (cid:19) − (cid:18) I + h (cid:15) H en + (cid:19) ϕ n = ϕ ∗ n ϕ n , where the second to the last line uses the fact that I − h (cid:15) H en + and I + h (cid:15) H en + commute.Similarly we may use other time-reversible (but not symplectic) schemes [11], suchas the trapezoidal rule discretization (known in this context as the Crank-Nicolsonmethod). So the PT-CN scheme becomes ϕ n +1 = ϕ n + h (cid:15) H en ϕ n + h (cid:15) H en +1 ϕ n +1 , (3.4)Here H en = H e ( t n , ϕ n ) , H en +1 = H e ( t n +1 , ϕ n +1 ). In both PT-GL2 and PT-CNschemes, we need to solve ϕ n +1 with nonlinear equation solvers as before. Althoughthese schemes are not symplectic schemes and the 2-norm of the numerical solution y PT-CN is not strictly conserved as in PT-Ham-GL2, numerical results in Section 5indicate that the performance of all the three schemes can be very comparable inpractice.Following the discussion above, we may readily obtain the corresponding schemefor N >
4. Analysis in the near adiabatic regime.
In this section, we demonstratethe advantage of the PT dynamics by analyzing the accuracy of the discretized PTdynamics in the near adiabatic regime. Our main result is that for h ≤ O ( (cid:15) ), a properdiscretization of the PT dynamics gains one extra order of accuracy in (cid:15) compared tothat of the Schr¨odinger dynamics.We extend the quantum adiabatic theorem [22, 2, 38] to the PT dynamics, whichshows that the PT wave function ϕ ( t ) can be decomposed into a component of whichthe oscillation is independent of (cid:15) and the magnitude is O (1), and a component thatis highly oscillatory with O ( (cid:15) ) magnitude. This leads to the desired result in terms ofthe local truncation error. We then obtain the global error estimate from the standardresults of symplectic integrators due to the Hamiltonian structure of the dynamics.Again, we restrict the scope of the theoretical analysis to the time-dependentlinear system with N = 1. While the generalization to the case N > H : [0 , T ] → C d × d is a Hermitian-valued and smooth map. The norms (cid:107) H ( t ) (cid:107) and (cid:107) H ( k ) ( t ) (cid:107) for all the time derivatives are bounded independentlyof (cid:15) and t ∈ [0 , T ].2. There exists a continuous function λ ( t ) ∈ spec( H ( t )) which is a simple eigen-value of H ( t ) and stays separated from the rest of the spectrum, i.e. thereexists a positive constant ∆ such thatdist( λ ( t ) , spec( H ( t )) \{ λ ( t ) } ) ≥ ∆ , ∀ t ∈ [0 , T ] . (4.1)3. The initial state ϕ (0) is the normalized eigenvector of H (0) associated withthe eigenvalue λ (0).The assumption 1 ensures that the solutions of both the Schr¨odinger dynamics andthe PT dynamics are smooth with respect to t . The assumption 2 is called the gapcondition [38].Before we continue, we would like to investigate a useful conclusion which can bedirectly derived from the above assumptions. Let Q ( t ) denote the projector on theeigenspace corresponding to λ ( t ). Q ( t ) can be expressed by the Riesz representationof the projector as Q ( t ) = − π i (cid:90) Γ( t ) R ( z, t ) dz (4.2) n which R ( z, t ) = ( H ( t ) − z ) − is the resolvent at time t and the complex contourcan be chosen as Γ( t ) = { z ∈ C : | z − λ ( t ) | = ∆ / } . Note that the assumption 2assures that such representation is well-defined and, together with assumption 1, Q ( t )is actually also a smooth bounded map, which is summarized in the following lemma. Lemma 4.
The norms of all time derivatives (cid:107) Q ( k ) ( t ) (cid:107) are bounded independentlyof (cid:15) .Proof . We follow the technique in [38]. The boundedness of Q ( t ) directly followsfrom the Riesz representation (4.2) and the boundedness of R ( z, t ) over the contourΓ( t ). The contour Γ( t ) depends on t . To avoid taking time derivatives over thecontour, note that the continuity of λ ( t ) implies that for any s ∈ [0 , T ], there exists aneighborhood B ( s, δ s ) such that | z − λ ( t ) | ≥ ∆ / , ∀ t ∈ B ( s, δ s ) ∩ [0 , T ] , z ∈ Γ( s ) . By finding a finite cover (cid:83) nj =1 B ( s j , δ s j ) ⊃ [0 , T ], for each t ∈ [0 , T ], there exists a s j such that t ∈ B ( s j , δ s j ) and we can rewrite Q ( t ) as Q ( t ) = − π i (cid:90) Γ( s j ) R ( z, t ) dz. (4.3)Such s j remains unchanged locally, hence Q ( k ) ( t ) = − π i (cid:90) Γ( s j ) R ( k ) ( z, t ) dz. The boundedness of Q ( k ) ( t ) can be directly assured by the boundedness of H ( k ) ( t ). First let us define the adiabatic evolution ϕ A ( t ) asthe solution to the following initial value problemi (cid:15)∂ t ϕ A = i (cid:15) [ ˙ Q, Q ] ϕ A , ϕ A (0) = ϕ (0) . (4.4)Since the matrix i (cid:15) [ ˙ Q, Q ] is Hermitian, (cid:107) ϕ A (cid:107) = 1 holds for all t . Following thesame proof of Eq. (2.14) in Proposition 1, we find that ϕ A is an eigenvector of H ( t )corresponding to λ ( t ), i.e. Q ( t ) ϕ A ( t ) = ϕ A ( t ) holds for all t ∈ [0 , T ].In the near adiabatic regime, we may separate ϕ ( t ) into the smooth component ϕ A and a remainder term. This is called the adiabatic theorem and is given in Theorem 5. Theorem 5.
Let ϕ ( t ) follow the PT dynamics (2.8) , and let ϕ A ( t ) follow theadiabatic evolution as defined in Eq. (4.4) . Then the following decomposition ϕ ( t ) = ϕ A ( t ) + (cid:15)ϕ R ( t ) (4.5) holds up to time T = O (1) . Furthermore, ϕ R ( t ) is infinitely differentiable, and (cid:107) ϕ R ( t ) (cid:107) is bounded independently of (cid:15) .Proof . The proof is organized according to the following three steps.1. Define another adiabatic evolution ϕ B , which satisfies an equation that re-sembles the PT dynamics.2. Prove the adiabatic decomposition with respect to ϕ B , i.e. there exists aninfinitely differentiable function η ( t ) such that ϕ ( t ) = ϕ B ( t ) + (cid:15)η ( t ) , ∀ t ∈ [0 , T ] , where (cid:107) η ( t ) (cid:107) is bounded independently of (cid:15) . . Prove that the difference between ϕ B and ϕ A is of O ( (cid:15) ).1. Define T B as the solution to the initial value problemi (cid:15)∂ t T B = ( H − ϕ ∗ Hϕ + i (cid:15) [ ˙ Q, Q ]) T B , T B (0) = I, (4.6)We define ϕ B according to ϕ B ( t ) := T B ( t ) ϕ (0) , which solves the initial value problemi (cid:15)∂ t ϕ B = ( H − ϕ ∗ Hϕ + i (cid:15) [ ˙ Q, Q ]) ϕ B , ϕ B (0) = ϕ (0) . (4.7)Since the matrix ( H − ϕ ∗ Hϕ + i (cid:15) [ ˙ Q, Q ]) is Hermitian, T B is a unitary evolution, and ϕ B is a normalized vector.Next we show that ϕ B ( t ) is an eigenvector of H ( t ) corresponding to λ ( t ), i.e. Q ( t ) ϕ B ( t ) = ϕ B ( t ) . (4.8)This can be done by showing that Qϕ B and ϕ B solve the same initial value problem.By the Leibniz rule and Eq. (4.7), we have ∂ t ( Qϕ B ) = ˙ Qϕ B + Q ˙ ϕ B = ˙ Qϕ B + Q [ ˙ Q, Q ] ϕ B − i (cid:15) Q ( H − ϕ ∗ Hϕ ) ϕ B . Use the identities similar to (2.12),˙ Q = ˙ QQ + Q ˙ Q, Q ˙ QQ = 0 , Q = Q, we have ˙ Q + Q [ ˙ Q, Q ] = ˙ QQ + Q ˙ Q + Q ˙ QQ − Q ˙ Q = ˙ QQ = ( ˙ QQ − Q ˙ Q ) Q = [ ˙ Q, Q ] Q. Hence ∂ t ( Qϕ B ) = [ ˙ Q, Q ] Qϕ B − i (cid:15) Q ( H − ϕ ∗ Hϕ ) ϕ B . Together with the identity QH = HQ , we have ∂ t ( Qϕ B ) = [ ˙ Q, Q ] Qϕ B − i (cid:15) ( H − ϕ ∗ Hϕ ) Qϕ B = − i (cid:15) ( H − ϕ ∗ Hϕ + i(cid:15) [ ˙ Q, Q ])( Qϕ B ) . Furthermore, the initial condition satisfies Q (0) ϕ B (0) = ϕ B (0) = ϕ (0). Hence Qϕ B solves the same initial value problem (4.7) as ϕ B .In summary, in step 1 we define another adiabatic evolution ϕ B ( t ) which is alsoan eigenstate of H ( t ) corresponding to λ ( t ) (Eq. (4.8)). Therefore, ϕ A ( t ) and ϕ B ( t )are both eigenstates of H ( t ) differing at most by a choice of gauge. . Now we estimate the distance between ϕ ( t ) and ϕ B ( t ). This can be done bymimicking the standard proof of the adiabatic theorem [2] with some modifications.By the definition of ϕ B , (cid:107) ϕ ( t ) − ϕ B ( t ) (cid:107) = (cid:107) ϕ ( t ) − T B ( t ) ϕ (0) (cid:107) = (cid:107)T − B ( t ) ϕ ( t ) − ϕ (0) (cid:107) . Define w ( t ) = T − B ( t ) ϕ ( t ), then (cid:107) ϕ ( t ) − ϕ B ( t ) (cid:107) = (cid:107) w ( t ) − w (0) (cid:107) = (cid:13)(cid:13)(cid:13)(cid:13)(cid:90) t ˙ w ( s ) ds (cid:13)(cid:13)(cid:13)(cid:13) . (4.9)In order to estimate ˙ w , differentiate the equation T B w = ϕ and we get˙ w = −T − B [ ˙ Q, Q ] T B w. (4.10)Note that if we define X ( t ) = − π i (cid:90) Γ( s j ) R ( z, t ) ˙ Q ( t ) R ( z, t ) dz where Γ( s j ) and R ( z, t ) are defined in the proof of Lemma 4, then (cid:107) X (cid:107) and (cid:107) ˙ X (cid:107) are bounded independently of (cid:15) , and [2, 38][ ˙ Q, Q ] = [
H, X ] . Then ˙ w = −T − B [ H, X ] T B w = − ( T − B H ) X T B w + T − B X ( H T B ) w. (4.11)To compute the first part of Eq. (4.11), we first take the time derivative of the identity I = T − B T B and get T − B H = − i (cid:15)∂ t ( T − B ) + ( ϕ ∗ Hϕ ) T − B − i (cid:15) T − B [ ˙ Q, Q ] . (4.12)Then the first part of Eq. (4.11) can be rewritten as − ( T − B H ) X T B w = i (cid:15)∂ t ( T − B ) X T B w + i (cid:15) T − B [ ˙ Q, Q ] X T B w − ( ϕ ∗ Hϕ ) T − B X T B w. (4.13)To compute the second part of Eq. (4.11), rewrite Eq. (4.6) as H T B = i (cid:15) ˙ T B + ( ϕ ∗ Hϕ ) T B − i (cid:15) [ ˙ Q, Q ] T B , (4.14)and then T − B X ( H T B ) w = i (cid:15) T − B X ˙ T B w − i (cid:15) T − B X [ ˙ Q, Q ] T B w + ( ϕ ∗ Hϕ ) T − B X T B w. (4.15)Sum up Eq. (4.13) and (4.15), then Eq. (4.11) becomes˙ w = i (cid:15) ( ∂ t ( T − B ) X T B + T − B X ˙ T B ) w + i (cid:15) T − B [[ ˙ Q, Q ] , X ] T B w. (4.16)In Eq. (4.16), the second term of the right hand side is already of O ( (cid:15) ). Now we turnto the first term to treat the derivatives ∂ t ( T − B ) and ˙ T B . By repeated usage of theLeibniz rule, Eq. (4.16) becomes˙ w = i (cid:15)∂ t ( T − B X T B ) w − i (cid:15) T − B ˙ X T B w + i (cid:15) T − B [[ ˙ Q, Q ] , X ] T B w = i (cid:15)∂ t ( T − B X T B w ) − i (cid:15) T − B X T B ˙ w − i (cid:15) T − B ˙ X T B w + i (cid:15) T − B [[ ˙ Q, Q ] , X ] T B w = i (cid:15)∂ t ( T − B Xϕ ) + i (cid:15) T − B X [ H, X ] ϕ − i (cid:15) T − B ˙ Xϕ + i (cid:15) T − B [[ ˙ Q, Q ] , X ] ϕ. (4.17) n the last equation we use again Eq. (4.10). Substitute Eq. (4.17) back to Eq. (4.9),we get (cid:107) ϕ ( t ) − ϕ B ( t ) (cid:107) = (cid:107) (cid:90) t ˙ w ( s ) ds (cid:107) ≤ (cid:15) (cid:107) ( T − B Xϕ )( t ) − ( T − B Xϕ )(0) (cid:107) + (cid:15) (cid:13)(cid:13)(cid:13)(cid:13)(cid:90) t (cid:0) T − B X [ H, X ] ϕ − T − B ˙ Xϕ + T − B (cid:2) [ ˙ Q, Q ] , X (cid:3) ϕ (cid:1) ds (cid:13)(cid:13)(cid:13)(cid:13) = O ( (cid:15) ) . (4.18)Therefore there exists η ( t ) such that ϕ ( t ) = ϕ B ( t ) + (cid:15)η ( t ) , (4.19)where (cid:107) η ( t ) (cid:107) is bounded independently of (cid:15) . The differentiability of η ( t ) followsdirectly from that of ϕ ( t ) and ϕ B ( t ).3. Comparing Eq. (4.19) with our goal, the only thing that we need to proveis that the distance between ϕ B and ϕ A is also O ( (cid:15) ). Note that ϕ A can be writtenas [10] ϕ A ( t ) = T (cid:20) exp (cid:18)(cid:90) t [ ˙ Q ( s ) , Q ( s )] ds (cid:19)(cid:21) ϕ A (0) , (4.20)where T is the time ordering operator due to the explicit time dependence of Q . Usingthe power series representation, the time-ordered exponential is defined as T (cid:104) e (cid:82) t A ( s ) d s (cid:105) = I + (cid:90) t A ( s ) d s + 12! (cid:90) t (cid:90) t T [ A ( s ) A ( s )] d s d s + · · · , (4.21)where the time-ordered product of two matrices T [ A ( s ) A ( s )] is given by T [ A ( s ) A ( s )] = (cid:40) A ( s ) A ( s ) , s ≥ s ; A ( s ) A ( s ) , s < s . (4.22)Using Duhamel’s principle, we have from Eq. (4.4) and (4.7) ϕ B ( t ) = ϕ A ( t )+ (cid:90) t T (cid:20) exp (cid:18)(cid:90) ts [ ˙ Q ( s (cid:48) ) , Q ( s (cid:48) )] ds (cid:48) (cid:19)(cid:21) · (cid:15) ( H ( s ) − ϕ ∗ ( s ) H ( s ) ϕ ( s )) ϕ B ( s ) ds (4.23)By Eq. (4.8), (4.19), and the normalization condition of ϕ and ϕ B ,( H − ϕ ∗ Hϕ ) ϕ B = − λ ( (cid:15)η ∗ ϕ B + (cid:15)ϕ ∗ B η ) ϕ B − (cid:15) ( η ∗ Hη ) ϕ B = − λ [( ϕ B + (cid:15)η ) ∗ ( ϕ B + (cid:15)η ) − ϕ ∗ B ϕ B − (cid:15) η ∗ η ] ϕ B − (cid:15) ( η ∗ Hη ) ϕ B = (cid:15) λ ( η ∗ η ) ϕ B − (cid:15) ( η ∗ Hη ) ϕ B = O ( (cid:15) ) . (4.24)Hence Eq. (4.23) implies ϕ B − ϕ A = O ( (cid:15) ) . (4.25) herefore, ϕ R := η + ( ϕ B − ϕ A ) /(cid:15) is infinitely differentiable, and (cid:107) ϕ R ( t ) (cid:107) is boundedindependently of (cid:15) . This proves the decomposition of the solution to the PT dynamics ϕ = ϕ B + (cid:15)η = ϕ B + (cid:15)ϕ R − ( ϕ B − ϕ A ) = ϕ A + (cid:15)ϕ R . (4.26)Theorem 5 gives a decomposition near the adiabatic regime with respect to thePT wave function. As a corollary, we also have the adiabatic theorem with respect tothe projector. Corollary 6.
For the projector P ( t ) , there exists an infinitely differentiablematrix-valued function R ( t ) such that P ( t ) = Q ( t ) + (cid:15)R ( t ) (4.27) holds for all t up to T = O (1) , where (cid:107) R ( t ) (cid:107) is bounded independently of (cid:15) .Proof . This follows directly from theorem 5 P = ϕϕ ∗ = ( ϕ A + (cid:15)ϕ R )( ϕ A + (cid:15)ϕ R ) ∗ = Q + (cid:15) ( ϕ R ϕ ∗ A + ϕ A ϕ ∗ R + (cid:15)ϕ R ϕ ∗ R ) . (4.28) Remark 7.
The adiabatic theorem for the Schr¨odinger wave function ψ ( t ) hasbeen well established in the literature e.g. [22, 2, 38], where the decomposition takesthe form ψ = ψ A + (cid:15) ˜ ψ R , and the adiabatic evolution ψ A satisfies i (cid:15)∂ t ψ A = ( H + i (cid:15) [ ˙ Q, Q ]) ψ A . (4.29) We compare our result with previous well-established ones from two aspects. First,there is an important difference between the PT eigenfunction ϕ A , governed by Eq. (4.4) ,and the standard one ψ A , governed by Eq. (4.29) . Although both ϕ A and ψ A are eigen-functions of H ( t ) , their phase factors are different, resulting in different oscillatorybehavior. More specifically, the standard wavefunction ψ A oscillates on the scale of O ( (cid:15) − ) since (at least intuitively) Eq. (4.29) is just a small perturbation of the originalSchr¨odinger dynamics. The PT eigenfunction ϕ A does not depend on (cid:15) , and thus os-cillates on the scale of O (1) . When projected to the eigenspace, the PT dynamics leadsto the optimal phase factor, and this verifies the effectiveness of the definition of PT(to minimize unnecessary oscillations) and provides another theoretical explanationof the performance shown in Fig. 2.1a. Second, our proof largely follows the existingworks of the adiabatic theorem [22, 2, 38]. Our main modification is to address thespecial non-linear term in the PT dynamics, even though the original Schr¨odingerdynamics is linear. Remark 8.
As mentioned at the end of step 1, ϕ B is also an eigenstate, andEq. (4.19) indeed leads to another version of the adiabatic theorem, but with notabledifferences from the decomposition in Theorem 5. First, the definition of ϕ B stillrelies on the information of ϕ , and thus is not a self-contained equation. Second, thenorms of the derivatives of ϕ B still depend on (cid:15) (more precisely one can prove that (cid:107) ϕ ( k ) B (cid:107) ∼ O (1 /(cid:15) k − ) for k ≥ ), which indicates that the gauge choice of ϕ B is notoptimal either. .2. Local truncation error. In this section, we show that after time dis-cretization, the local truncation error of the discretized PT dynamics improves byone order in terms of (cid:15) compared to that of the discretized Schr¨odinger dynamics inthe near adiabatic regime. This is given in Lemma 9.For simplicity we will focus on the numerical integrators in the classes of Runge-Kutta methods and linear multistep methods, both of which are widely used forsimulating the Schr¨odinger equation. We will refer numerical integrator to eithera Runge-Kutta method or a linear multistep method in our context. Recall that anumerical integrator with a given time step h , denoted by I h , can be generally writtenas u n +1 = I h ( u n , · · · , u n − l ) , (4.30)for some integer l ≥
0, and u n is the numerical approximation to the true solution u ( t n ). If I h is of order k , then the local truncation error at step n + 1, defined as L n +1 = I h ( u ( t n ) , · · · , u ( t n − l )) − u n +1 , should satisfy (cid:107) L n +1 (cid:107) ≤ Ch k +1 (cid:107) u ( k +1) ( ξ n +1 ) (cid:107) , for some ξ n +1 ∈ [ t n , t n +1 ]. When applied to the Schr¨odinger dynamics, the PTdynamics, or the associated Hamiltonian form, we may identify u with ψ , ϕ , or theequivalent ( q, p ) representation. Lemma 9.
Apply a numerical integrator of order k to the Schr¨odinger dynam-ics or its Hamiltonian form (2.19) . Then the local truncation error is bounded by Ch k +1 /(cid:15) r up to the time T ∼ O (1) , with r = k + 1 and C is a constant independentof h and (cid:15) . The same result holds for the PT dynamics (2.8) or its correspondingHamiltonian form (2.23) with r = k .Proof . It is sufficient to show that the derivatives satisfy (cid:107) ψ ( k +1) (cid:107) ≤ O (1 /(cid:15) k +1 ),and (cid:107) ϕ ( k +1) (cid:107) ≤ O (1 /(cid:15) k ) for any k ≥
0. This can be proved by induction.1. For ψ , the case k = 0 directly follows from Eq. (1.1). Assume the estimateholds for all the integers smaller than k , differentiate the Schr¨odinger equation k timesand we get ψ ( k +1) = 1i (cid:15) k (cid:88) j =0 (cid:18) kj (cid:19) H ( k − j ) ψ ( j ) . (4.31)By the induction and the assumption 1, (cid:107) ψ ( k +1) (cid:107) ≤ C(cid:15) k (cid:88) j =0 (cid:18) kj (cid:19) (cid:15) j ∼ O ( (cid:15) − ( k +1) ) . (4.32)2. For ϕ , we first study the derivatives of P , and then use the PT condition (2.6)to obtain the derivatives of ϕ .By Corollary 6, the von Neumann equation (1.5) and the identity HQ = QH , thefirst order derivative of P satisfies (cid:107) ˙ P (cid:107) = 1 (cid:15) (cid:107) HP − P H (cid:107) = (cid:107) HR − RH (cid:107) ≤ O (1) . urthermore, differentiate the von Neumann equation (1.5) k times, we get P ( k +1) = 1i (cid:15) k (cid:88) j =0 (cid:18) kj (cid:19) [ H ( j ) , P ( k − j ) ] , (4.33)from which we can show by induction that (cid:107) P ( k +1) (cid:107) ≤ O ( (cid:15) − k ) . (4.34)Now use the PT condition P ˙ ϕ = 0, we find for k = 0,˙ ϕ = ∂ t ( P ϕ ) = ˙
P ϕ ≤ O (1) . (4.35)Furthermore, ϕ ( k +1) = k (cid:88) j =0 (cid:18) kj (cid:19) P ( j +1) ϕ ( k − j ) , (4.36)from which we can prove by induction and Eq. (4.34) that ϕ ( k +1) ≤ O ( (cid:15) − k ) . (4.37) The analysis of the local truncation error directly extends tothe global error up to T ∼ O ( (cid:15) ), following the classical stability analysis. However, theLipschitz constants corresponding to the right hand side of the Schr¨odinger dynamicsand the PT dynamics are generally O (1 /(cid:15) ), which leads to an exponentially growingfactor exp( T /(cid:15) ) in the global error bounds. Hence we cannot directly obtain the globalerror estimate up to O (1) time.However, if we adopt the Hamiltonian formulation of the dynamics and employ asymplectic integrator, we can indeed obtain long time error estimates. This is statedin Theorem 10, of which the proof directly follows from Lemma 9 and Theorem X.3.1in [11]. Theorem 10.
Apply a symplectic integrator of order k to the Hamiltonian system (2.19) and (2.23) , then there exist constants c, C , independent of h and (cid:15) , such thatfor the time step h ≤ c(cid:15) , the numerical solutions up to the time T ∼ O (1) satisfy (cid:107) ( q n , p n ) − ( q ( t ) , p ( t )) (cid:107) ≤ C h k (cid:15) r . (4.38) Here r = k + 1 for the Schr¨odinger dynamics (2.19) and r = k for the PT dynamics (2.23) . Remark 11.
In Theorem X.3.1 in [11], all terms are bounded by O (1) termsand there is no (cid:15) dependence. In order to adapt its proof to the current situation,we observe the key fact in Theorem X.3.1 in [11] that the global error of a symplecticintegrator accumulates linearly in time with no exponential growing factor. Thereforethe local truncation error which is O ( h k +1 /(cid:15) r ) directly sums up linearly to the globalerror of O ( h k /(cid:15) r ) . Remark 12.
The nontrivial restriction on the time step size h ≤ c(cid:15) is becauseTheorem X.3.1 in [11] holds only for sufficiently small time steps. In general, h mustbe no larger than c/L where L is the Lipschitz constant of the right hand side of the amiltonian system, and is O (1 /(cid:15) ) in the singularly perturbed regime. Nonetheless,numerical results in Section 5 indicate that the PT dynamics may admit a considerablylarger time step in practice. Remark 13.
When a symplectic integrator is used, Theorem 10 is directly appli-cable to the Schr¨odinger dynamics. However, the PT dynamics (2.8) and the Hamil-tonian system (2.23) share the same exact solution, but lead to different numericalschemes even when the same integrator is used. Despite such difference, numericalresults in Section 5 indicate that the symplectic integrators, and even certain non-symplectic schemes, can still perform very well in the PT dynamics (2.8) . Remark 14.
Theorem 10 also indicates that the PT dynamics is relatively moreeffective when combined with low order methods. For instance, if we would liketo achieve some desired accuracy δ (assuming δ is sufficiently small), then for theSchr¨odinger dynamics, we should choose the time step size to be h ∼ O ( δ k (cid:15) k ) . For the PT dynamics, we should choose h ∼ O ( δ k (cid:15) ) . From this perspective, the gain of the PT dynamics is less significant when k is large.
5. Numerical results.
In this section we study the effectiveness of the PTdynamics using three examples. The first one is a toy example, which is a linearSchr¨odinger equation in C . This example gives a clear illustration of the performanceof different numerical methods near and beyond the adiabatic regime. The secondexample is a nonlinear Schr¨odinger equation in a one-dimensional space, where wealso compare the computational cost between the propagation of the Schr¨odingerdynamics and the PT dynamics. In the end we study the electron dynamics of abenzene molecule driven by an ultrashort laser pulse described by the time-dependentdensity functional theory (TDDFT).The test programs in the first two examples are written in MATLAB. We imple-ment the PT dynamics for TDDFT in the PWDFT code, which performs planewavebased electronic structure calculations. PWDFT is a self-contained module in themassively parallel DGDFT (Discontinuous Galerkin Density Functional Theory) soft-ware package written in MPI and C++ [25, 15]. All calculations are carried out usingthe BRC High Performance Computing service. Each node consists of two Intel Xeon10-core Ivy Bridge processors (20 cores per node) and 64 gigabyte (GB) of memory.We use the Anderson mixing for solving all the nonlinear fixed point problems, includ-ing those in the PT dynamics, the nonlinear Schr¨odinger equation, and the TDDFTcalculations. Here no preconditioner is used for the first two examples. We use ashifted Laplace preconditioner for the TDDFT example, which can be implementedefficiently in the planewave basis set using the fast Fourier transform. First we present a linear example, in which H ( t ) is chosento be H ( t ) = (cid:18) t − t δδ − ( t − t ) (cid:19) . (5.1)Here H ( t ) has the eigenvalues λ , ( t ) = ∓ (cid:112) ( t − t ) + δ , where δ > δ is large, the dynamics stays a) δ = 1 (b) δ = 0 . Fig. 5.1: Eigenvalues of H ( t ) in the toy example with t = 0 . δ .closer to the adiabatic regime, while the dynamics can go beyond the adiabatic regimewith a smaller δ (see Fig. 5.1). The initial value is always chosen to be the normalizedeigenvector of H (0) corresponding to λ (0) = − (cid:112) t + δ . We propagate the wavefunctions up to T = 1. For the choices of the parameters in the Anderson Mixing inpropagating PT dynamics, the step length α = 1, the mixing dimension is 20, andthe tolerance is 10 − . First we consider the near adiabatic case with δ = 1. We compare the following numerical methods: • S-RK4: fourth order Runge-Kutta method (RK4) applied to the Schr¨odingerequation (1.1) • PT-RK4: fourth order Runge-Kutta method (RK4) applied to the PT dy-namics (2.8) • S-GL2: implicit midpoint rule (GL2) applied to the Schr¨odinger equation (1.1) • PT-Ham-GL2: implicit midpoint rule (GL2) applied to the PT Hamiltoniansystem (2.23) • PT-GL2: implicit midpoint rule (GL2) applied to the PT dynamics (2.8) • PT-CN: trapezoidal rule (or the Crank-Nicolson method, CN) applied to thePT dynamics (2.8)Fig. 5.2 compares the performance of different methods for this toy example. Thenumerical error is computed by e ( h, (cid:15) ) = max n s.t. nh ∈ [0 ,T ] (cid:107) u n − u ( t n ) (cid:107) where u denotes ψ for the Schr¨odinger dynamics, ϕ for the PT dynamics and ( q, p )for the Hamiltonian systems, respectively.We first consider the explicit numerical methods. Fig. 5.2a and 5.2b give a com-parison between S-RK4 and PT-RK4. Not surprisingly, as an explicit method, RK4 isnumerically unstable for large time steps under both cases, and achieves fourth orderconvergence for small time steps. Furthermore, when h is small enough, e ( h, (cid:15) ) of thePT dynamics is smaller than that of the Schr¨odinger dynamics. Fig. 5.3a presents astudy on how e ( h, (cid:15) ) depends on (cid:15) , which reveals that by propagating the PT dynam- a) (cid:15) = 0 .
01 (b) (cid:15) = 0 . (cid:15) = 0 .
01 (d) (cid:15) = 0 . Fig. 5.2: Numerical errors of different numerical methods in the near adiabatic regimeof the toy example. (a)(b) compare S-RK4 and PT-RK4 for (cid:15) = 0 . , . (cid:15) = 0 . , . (cid:15) . This agrees with the theoreticalresults in Section 4.Next we test GL2 as an example of implicit symplectic methods applied to theHamiltonian systems. Fig. 5.2c compares the numerical performances of S-GL2 andPT-Ham-GL2. For small h , we observe a smaller error using the PT formulation, i.e. e ( h, (cid:15) ) of S-GL2 is O ( h /(cid:15) ) and e ( h, (cid:15) ) of PT-Ham-GL2 is O ( h /(cid:15) ) (see Fig. 5.3bfor a study on the (cid:15) dependence). This verifies the estimate in Theorem 10. Despitethat GL2 is a numerically stable scheme with a large time step, the step size of S-GL2is constrained by the requirement of the accuracy, while the step size of PT-Ham-GL2can be chosen to be considerably larger.More specifically, let us define the “turning point” h T to be the largest time stepsize when a scheme starts to converge. Numerically for second order schemes theturning point can be computed as h T = arg max (cid:26) h ∈ [ h , h ] : ∂ (log e ) ∂ (log h ) > (cid:27) a) (b) Fig. 5.3: Relationship between the asymptotic errors and (cid:15) in the near adiabaticregime of the toy example. Here we fix the time step size to be h = 10 − for both(a)(b). (a) (b) Fig. 5.4: (a) Relationship between the turning points and (cid:15) in S-GL2 and PT-Ham-GL2 in the near adiabatic regime of the toy example. (b) Relationship between themagnitude of the plateau of the numerical error and (cid:15) in PT-Ham-GL2.where [ h , h ] is a suitable interval containing the convergence interval of interests.In Fig. 5.2c we mark the turning points in S-GL2 and PT-Ham-GL2, and study theirdependence on (cid:15) in Fig. 5.4a. For S-GL2, the convergence starts at h T = O ( (cid:15) / ).For PT-Ham-GL2, a two-stage convergence behavior is observed. As h decreases,the scheme first starts to converge with second order at a relatively large time step h T = O ( (cid:15) / ). This first stage ends at h = O ( (cid:15) ) when e ( h, (cid:15) ) reaches a plateau withits magnitude being O ( (cid:15) ) (see Fig. 5.4b). Then the second-stage convergence startsat h T = O ( (cid:15) / ).In the end we compare the schemes PT-Ham-GL2, PT-GL2 and PT-CN. Althoughwe only justified the behavior of the global error for PT-Ham-GL2, numerical resultsin Fig. 5.5a and 5.5b indicate that there is no essential difference among these methods a) (cid:15) = 0 .
01 (b) (cid:15) = 0 . Fig. 5.5: Performance of PT-Ham-GL2, PT-GL2 and PT-CN in the near adiabaticregime of the toy example.in practice.
As the value of δ is reduced, the secondeigenstate corresponding to λ may contribute significantly to the wave function,which leads to the violation of the adiabatic regime.Fig. 5.6 investigates the Schr¨odinger wave function and the PT wave functionwith (cid:15) = 0 . , δ = 0 .
05. Fig. 5.6a and 5.6b compare the real parts of the Schr¨odingerwave function and the PT wave function. When t < t = 0 .
5, the system stays closeto the adiabatic regime and the PT wave function is nearly flat. However, when t > t , the PT wave function starts to oscillate as well. Fig. 5.6c shows an orthogonaldecomposition of the PT wave function into two orthogonal eigenspaces. Fig. 5.6dshows the evolution of the probability that the eigenstate corresponding to λ ( t ) isoccupied, which can be computed as | c | = | ( ϕ ( t ) , e ( t )) | and e ( t ) is the normalizedeigenstate of H ( t ) corresponding to λ ( t ). These results confirm that the oscillatorybehavior originates from the excited state corresponding to λ .As discussed before, such oscillatory nature in the wave functions may increasethe computational difficulty and require a smaller time step even for the PT dy-namics. Fig. (5.7) compares e ( h, (cid:15) ) for S-GL2, PT-Ham-GL2, PT-GL2 and PT-CNrespectively. The results confirm that the PT dynamics is always more accurate thanthe Schr¨odinger dynamics using the same step size, but the gain becomes smaller as δ decreases. Next we studythe performance of the PT dynamics in a singularly perturbed nonlinear Schr¨odingerequation in one dimension.i (cid:15)∂ t ψ ( x, t ) = − ∂ x ψ ( x, t ) + V ( x, t ) ψ ( x, t ) + g | ψ ( x, t ) | ψ ( x, t ) , x ∈ [0 , L ] ψ ( x,
0) = ψ ( x ) ψ (0 , t ) = ψ ( L, t ) . (5.2) a) (b)(c) (d) Fig. 5.6: The Schr¨odinger and the PT wave functions beyond the adiabatic regime inthe toy example. In all sub-figures, parameters are chosen to be (cid:15) = 0 . , δ = 0 . h = 10 − . (a)(b) showthe first and second entry of the real part of the Schr¨odinger wave function and thePT wave function, respectively. (c) shows a decomposition of the PT wave functioninto the two orthogonal eigenspaces (in the sub-figure we only present the real partof the first entry). (d) shows the time evolution of the probability that the eigenstatecorresponding to λ is occupied.We set L = 50, and the external potential is chosen to be a time-dependent Gaussianfunction modeling a moving potential well (Fig. 5.8) V ( x, t ) = − exp( − . x − R ( t )) ) (5.3)with a time-dependent center R ( t ) = 25 + 1 . − t − . ) + exp( − t − . ) . (5.4)Note that R ( t ) varies on the O (1) time scale.We use equidistant nodes x k = kh x and the second-order finite difference schemefor spacial discretization, and we fix h x = 0 . g = 2 . , T = 1 , (cid:15) = 0 . a) δ = 0 .
07 (b) δ = 0 . δ = 0 .
03 (d) δ = 0 . Fig. 5.7: Numerical errors of different numerical methods beyond the adiabatic regimein the toy example. In all sub-figures (cid:15) = 0 . δ = 0 . , . , .
03, respectively.(d) gives a comparison of PT-Ham-GL2, PT-GL2 and PT-CN with δ = 0 . α = 1, the mixing dimension is 20, and thetolerance is 10 − . Fig. 5.9 compares e ( h, (cid:15) ) of S-GL2, PT-Ham-GL2, PT-GL2 andPT-CN, and confirms the same numerical behavior as in the toy example.Next we study the computational cost by comparing the total number of theAnderson mixing steps versus the numerical error e ( h, (cid:15) ) up to T = 1. Fig. 5.10 clearlydemonstrates that in order to achieve the same level of accuracy, all the methodspropagating the PT dynamics, including PT-Ham-GL2, PT-GL2 and PT-CN, aremuch more efficient than S-GL2. This is valid across the entire range of the step sizesunder study. Asthe last example, we demonstrate the performance of the PT dynamics for a ben-zene molecule driven by an ultrashort laser pulse using the time-dependent densityfunctional theory (TDDFT). The TDDFT equations arei ∂ t Ψ( t ) = H ( t, P )Ψ( t ) , P ( t ) = Ψ( t )Ψ ∗ ( t ) , (5.5) a) (b) Fig. 5.8: External potential and the time-dependent center for the nonlinearSchr¨odinger equation. (a) (b)
Fig. 5.9: Numerical errors of different numerical methods in the example of the non-linear Schr¨odinger equation. Parameters are chosen to be T = 1 , (cid:15) = 0 . ∂ t Φ( t ) = H ( t, P )Φ( t ) − Φ( t )Φ ∗ ( t ) H ( t, P )Φ( t ) , P ( t ) = Φ( t )Φ ∗ ( t ) . (5.6)The number of wavefunctions N is 15 for benzene example. Compared to the setupof singularly perturbed equations, here in the sense that the parameter (cid:15) is formallyset to 1 in TDDFT equations. However, as will be seen later, the PT dynamics canstill result in significant computational advantage. The Hamiltonian takes the form H ( t, P ) = −
12 ∆ + V ext ( r , t ) + V PP ( r ) + V Hxc [ ρ ( t )] . (5.7) ig. 5.10: Total numbers of the Anderson mixing versus the numerical error.Here V PP is the pseudopotential operator due to the electron-ion interaction, and weuse the Optimized Norm-Conserving Vanderbilt (ONCV) pseudopotential [14] witha kinetic energy cutoff E cut = 30 Hartree. After spatial discretization, V PP becomesa matrix independent of the time t and the density matrix P . V Hxc is the sum ofthe Hartree and exchange-correlation potentials. We use the Perdew-Burke-Ernzerhof(PBE) [33] exchange correlation potential that depends on the electron density ρ ( t ) =diag[ P ( t )]. The external potential V ext ( r , t ) = r · E ( t ) is given by a time-dependentelectric field E ( t ) = ˆ k E max exp (cid:104) − ( t − t ) a (cid:105) sin[ ω ( t − t )] , (5.8)where ˆ k is a unit vector defining the polarization of the electric field. The parame-ters a, t , E max , ω define the width, the initial position of the center, the maximumamplitude of the Gaussian envelope, and the frequency of the laser, respectively. Inpractice ω and a are often determined by the wavelength λ and the full width at halfmaximum (FWHM) pulse width [35], i.e. λω = 2 πc and FWHM = 2 a √ c is the speed of the light. In this example, the peak electric field E max is 1.0 eV/˚A,occurring at t = 15 . x axis (the benzene molecule is in x - y plane,see Fig. 5.11a). We consider one relatively slow laser with wavelength 800 nm, andanother faster laser with wavelength 250 nm, respectively (Fig. 5.11). The electrondynamics for the first laser is in the near adiabatic regime, where the system staysnear the ground state after the active time interval of the laser, while the second laserdrives electrons to excited states. We implement S-RK4 and PT-CN in the PWDFTpackage, and propagate TDDFT to T = 30 . α is 0 .
2, the mixing dimension is 10, and the tolerance is 10 − .We measure the accuracy using the dipole moment D ( t ) := Tr[ r P ( t )], as well as theenergy difference E ( t ) − E (0) along the trajectory.Figure 5.12 shows the numerical results for the 800 nm laser using S-RK4 with astep size 0.0005 fs and PT-CN with a step size 0.05 fs. In this case, the system staysnear the ground state after the active time interval of the laser. After 25.0 fs, the totalenergy for S-RK4 only increases by 2 . × − eV, and hence we may use the resultsfrom S-RK4 as our benchmark. We remark that S-RK4 becomes unstable at largetime step sizes. Even when increasing the time step to be 0.001 fs, S-RK4 blows up a)(b) (c) Fig. 5.11: (a) The benzene molecule. The direction of the external electric field isalong the x-axis. This figure is generated by VMD package [16]. (b)(c) The intensityof the electric field. The peak electric field E max is 1.0 eV/˚A, occurring at t = 15 . x direction, and the total energy difference. After25.0 fs, the total energy is nearly constant and only slightly increases by 2 . × − eV compared to that of the initial state.Since the computational cost of TDDFT calculations is mainly dominated by thecost of applying the Hamiltonian matrix to wave functions, we measure the numericalefficiency using the number of such matrix-vector multiplications. Although PT-CNrequires more matrix-vector multiplications in each time step, the total number ofmatrix-vector multiplications is still significantly reduced due to the larger time stepsize, and PT-CN usually achieves a significant speedup. More specifically, in thiscase, during the time interval for which the laser is active (from 5.5 fs to 24.5 fs),the average number of matrix-vector multiplications in each PT-CN time step is 12.6,and the total number of matrix-vector multiplications in the simulation is 4798. Onthe other hand, the number of matrix-vector multiplications in each S-RK4 time stepis 4, and the total number of matrix-vector multiplications during this period using a) (b) Fig. 5.12: (a) Dipole moment along the x-direction and (b) total energy differencewith the 800 nm laser.time step 0.0005 fs is 152000. Hence the overall speedup of PT-CN over RK4 is 31 . x directionoscillates more strongly due to the excitation. PT-CN needs to adopt a smaller timestep size 0.005 fs, and still gives a very good approximation to the electron dynamicscompared to S-RK4, For the dipole moment, PT-CN results match very well with S-RK4 benchmark during (Fig. 5.13b) and after (Fig. 5.13c and 5.13d) the active timeinterval of the laser. The total energy obtained by PT-CN matches very well withthat in S-RK4 benchmark during the active interval and stays at a constant level withan average increase of 0.5340 eV by the end of the simulation (Fig. 5.13e and 5.13f).In this case, PT-CN slightly overestimates the total energy after the laser’s action by7 . × − eV.For the computational costs within the period from 5.5 fs to 24.5 fs, the totalnumber of matrix-vector multiplications is still 152000 for S-RK4. The average num-ber of matrix-vector multiplications in each PT-CN time step is 7.5 due to the reducedstep size, and the total number of matrix-vector multiplications is 28610. Thereforein this case PT-CN achieves 5.3 times speedup over S-RK4.We remark that even the electron dynamics is beyond the adiabatic regime, PT-CN can still be stable with a larger time step. Table 5.1 measures the accuracy ofPT-CN with h = 0.005 fs, 0.0065 fs, 0.0075 fs, 0.01 fs and 0.02 fs, respectively. We findthat the number of matrix-vector multiplications systematically reduces as the stepsize increases. When the step size is 0.02 fs, the speed up over S-RK4 is 12 .
6, and thisis at the expense of overestimating the energy by 0 .
6. Conclusion.
Quantum dynamics can be equivalently written in terms of theSchr¨odinger equation for the wave function, and the von Neumann equation for thedensity matrix. However, the Schr¨odinger dynamics may require a very small time step ethod h (fs) AEI (eV) AOE (eV) MVM SpeedupS-RK4 0.0005 0.5260 / 152000 /PT-CN 0.005 0.5340 0.0080 28610 5.3PT-CN 0.0065 0.5347 0.0087 22649 6.7PT-CN 0.0075 0.5362 0.0102 21943 6.9PT-CN 0.01 0.5435 0.0175 15817 9.6PT-CN 0.02 0.5932 0.0672 12110 12.6Table 5.1: Accuracy and efficiency of PT-CN for the electron dynamics with the 250nm laser compared to S-RK4. The accuracy is measured using the average energyincrease (AEI) after 25.0 fs and the average overestimated energy (AOE) after 25.0fs. The efficiency is measured using the total number of matrix-vector multiplications(MVM) during the time interval from 5.5 fs to 24.5 fs, and the computational speedup.in numerical simulation due to the non-optimal gauge choice. In this paper, we proposeto close this gap by identifying the optimal gauge choice, which is obtained from theparallel transport formulation. The solution of the resulting parallel transport (PT)dynamics can be significantly less oscillatory to that of the Schr¨odinger dynamics,especially in the near adiabatic regime. The PT dynamics is suitable to be combinedwith implicit time integrators, which allows the usage of large time steps even whenthe spectral radius of the Hamiltonian is large, and/or when (cid:15) is small. Although ourglobal error analysis only applies to the Hamiltonian form of the PT dynamics withsymplectic integrators and a relatively small time step, our numerical results indicatethat the PT dynamics can be effectively discretized with more general numericalschemes and with much larger time steps. The mathematical understanding of thebehavior with a large time step is our future work. Combining the PT dynamics withnumerical schemes other than the Runge-Kutta methods and the linear multistepmethods, as well as more detailed numerical studies of the PT dynamics for the time-dependent density functional theory calculations are also under progress. Acknowledgments.
This work was partially supported by the National ScienceFoundation under Grant No. 1450372, No. DMS-1652330 (D. A. and L. L.), and bythe Department of Energy under Grant No. de-sc0017867, No. DE-AC02-05CH11231(L. L.). We thank the National Energy Research Scientific Computing (NERSC)center and the Berkeley Research Computing (BRC) program at the University ofCalifornia, Berkeley for making computational resources available. We thank StefanoBaroni, Roberto Car, Weile Jia, Christian Lubich, and Lin-Wang Wang for helpfuldiscussions.
REFERENCES[1]
D. G. Anderson , Iterative procedures for nonlinear integral equations , J. Assoc. Comput.Mach., 12 (1965), pp. 547–560.[2]
J. E. Avron and A. Elgart , Adiabatic theorem without a gap condition , Communications inMathematical Physics, 203 (1999), pp. 445–463.[3]
W. Bao, S. Jin, and P. A. Markowich , On time-splitting spectral approximations for theSchr¨odinger equation in the semiclassical regime , J. Comput. Phys., 175 (2002), pp. 487–524.[4]
F. A. Bornemann and C. Sch¨utte , On the singular limit of the quantum-classical moleculardynamics model , J. Appl. Math., 59 (1999), pp. 1208–1224.325]
H. Candy and W. Rozmus , A symplectic integration algorithm for seperate hamiltonian func-tions , J. Comput. Phys., 92 (1991), pp. 230–256.[6]
A. Castro, M. Marques, and A. Rubio , Propagators for the time-dependent Kohn-Shamequations , J. Chem. Phys., 121 (2004), pp. 3425–33.[7]
Z. Chen and E. Polizzi , Spectral-based propagation schemes for time-dependent quantumsystems with application to carbon nanotubes , Phys. Rev. B, 82 (2010), p. 205410.[8]
D. Cohen, T. Jahnke, K. Lorenz, and C. Lubich , Numerical integrators for highly oscil-latory hamiltonian systems: a review , in Analysis, modeling and simulation of multiscaleproblems, Springer, 2006, pp. 553–576.[9]
H. D. Cornean, D. Monaco, and S. Teufel , Wannier functions and z2 invariants in time-reversal symmetric topological insulators , Rev. Math. Phys., 29 (2017), p. 1730001.[10]
A. L. Fetter and J. D. Walecka , Quantum theory of many-particle systems , Courier Corp.,2003.[11]
E. Hairer, C. Lubich, and G. Wanner , Geometric numerical integration: structure-preserving algorithms for ordinary differential equations , Springer-Verlag Berlin Heidel-berg, second ed., 2006.[12]
E. Hairer, S. P. Nørsett, and G. Wanner , Solving ordinary differential equation I: nonstiffproblems , vol. 8, Springer, 1987.[13]
E. Hairer and G. Wanner , Solving ordinary differential equation II: stiff and differential-algebraic problems , vol. 8, Springer, 1991.[14]
D. R. Hamann , Optimized norm-conserving Vanderbilt pseudopotentials , Phys. Rev. B, 88(2013), p. 085117.[15]
W. Hu, L. Lin, and C. Yang , DGDFT: A massively parallel method for large scale densityfunctional theory calculations , J. Chem. Phys., 143 (2015), p. 124110.[16]
William Humphrey, Andrew Dalke, and Klaus Schulten , VMD – Visual Molecular Dy-namics , J. Molec. Graphics, 14 (1996), pp. 33–38.[17]
A. Iserles , A first course in the numerical analysis of differential equations , no. 44, CambridgeUniv. Pr., 2009.[18]
T. Jahnke and C. Lubich , Numerical integrators for quantum dynamics close to the adiabaticlimit , Numer. Math., 94 (2003), pp. 289–314.[19]
J. Jia and J. Huang , Krylov deferred correction accelerated method of lines transpose forparabolic problems , J. Comput. Phys., 227 (2008), pp. 1739–1753.[20]
C. F. Kammerer and A. Joye , Nonlinear quantum adiabatic approximation , arXiv:1906.11069,(2019).[21]
A.-K. Kassam and L. N. Trefethen , Fourth-order time-stepping for stiff PDEs , SIAM J. Sci.Comput., 26 (2005), pp. 1214–1233.[22]
T. Kato , On the adiabatic theorem of quantum mechanics , J. Phys. Soc. J. Jpn., 5 (1950),pp. 435–439.[23]
C. T. Kelley , Iterative methods for optimization , vol. 18, SIAM, 1999.[24]
O. Koch and C. Lubich , Dynamical low-rank approximation , SIAM J. Matrix Anal. Appl., 29(2007), pp. 434–454.[25]
L. Lin, J. Lu, L. Ying, and W. E , Adaptive local basis set for Kohn-Sham density functionaltheory in a discontinuous Galerkin framework I: Total energy calculation , J. Comput.Phys., 231 (2012), pp. 2140–2154.[26]
L. Lin and C. Yang , Elliptic preconditioner for accelerating self consistent field iteration inKohn-Sham density functional theory , SIAM J. Sci. Comp., 35 (2013), pp. S277–S298.[27]
C. Lubich , On splitting methods for Schrodinger-Poisson and cubic nonlinear Schrodingerequations , Math. Comp., 77 (2008), pp. 2141–2153.[28]
M. Nakahara , Geometry, topology and physics , CRC Press, 2003.[29]
G. Nenciu , Linear adiabatic theory. exponential estimates , Communications in MathematicalPhysics, 152 (1993), pp. 479–496.[30]
P. Nettesheim, F. A.Bornemann, B. Schmidt, and C. Sch¨utte , An explicit and symplec-tic integrator for quantum-classical molecular dynamics , Chem. Phys. Lett., 256 (1996),pp. 581–588.[31]
P. Nettesheim and C. Sch¨utte , Numerical integrators for quantum-classical molecular dy-namics , (1999), pp. 396–411.[32]
G. Onida, L. Reining, and A. Rubio , Electronic excitations: density-functional versus many-body Green’s-function approaches , Rev. Mod. Phys., 74 (2002), p. 601.[33]
J. P. Perdew, K. Burke, and M. Ernzerhof , Generalized gradient approximation madesimple , Phys. Rev. Lett., 77 (1996), pp. 3865–3868.[34]
E. Runge and E. K .U. Gross , Density-functional theory for time-dependent systems , Phys.Rev. Lett., 52 (1984), p. 997. 3335]
A. Russakoff, Y. Li, S. He, and K. Varga , Accuracy and computational efficiency of real-time subspace propagation schemes for the time-dependent density functional theory , J.Chem. Phys., 144 (2016), p. 204125.[36]
A. Schleife, E. W. Draeger, Y. Kanai, and A. A. Correa , Plane-wave pseudopotentialimplementation of explicit integrators for time-dependent Kohn-Sham equations in large-scale simulations , J. Chem. Phys., 137 (2012), p. 22A546.[37]
C. Sparber , Weakly nonlinear time-adiabatic theory , Ann. Henri Poinc´are, 17 (2016), pp. 913–936.[38]
S. Teufel , Adiabatic perturbation theory in quantum dynamics , Springer-Verlag Berlin Hei-delberg, first ed., 2003.[39]
Z. Wang, S.-S. Li, and L.-W. Wang , Efficient real-time time-dependent density functionaltheory method and its application to a collision of an ion with a 2D material , Phys. Rev.Lett., 114 (2015), pp. 1–5. 34 a) (b)(c) (d)(e) (f)