Focus Topic: Propagators for the TDSE Part I: The Time-Evolution - - PowerPoint PPT Presentation

AARHUS UNIVERSITY

AU

Focus Topic: Propagators for the TDSE Part I: The Time-Evolution Operator and Crank-Nicolson

Kenneth Hansen

QUSCOPE Meeting, Aarhus University, Denmark

December 17, 2015

Kenneth Hansen December 17, 2015 1 / 44

AARHUS UNIVERSITY

AU

The Time-Evolution Operator

Time evolution in QM done through an operator: |Ψ(t) = ˆ U(t, t0) |Ψ(t0) This operator should fulfill: ˆ U(t, t0) = ˆ U(t, t1)ˆ U(t1, t0), ˆ U(t0, t0) = ˆ I, ˆ U†(t, t0) = ˆ U−1(t, t0) = ˆ U(t0, t) The Time-Dependent Schr¨

• dinger Equation (TDSE) gives us the

requirement that the time-evolution operator must satisfy: i ∂ ∂t ˆ U(t, t0) = ˆ H(t)ˆ U(t, t0) Which can be rewritten as: ˆ U(t, t0) = ˆ I − i t

t0

ˆ H(t1)ˆ U(t1, t0)dt1

Kenneth Hansen December 17, 2015 2 / 44

AARHUS UNIVERSITY

AU

Dividing the time-step in smaller steps so tk = t0 + k∆t for k = 0, 1, ..., N and by making a sequence by self-insertion in the integral form of the time-evolution operator we can expand the time-evolution operator: ˆ U(tk+1, tk) =

∞

• n=0

ˆ U

(n)

(tk+1, tk), With ˆ U

(0)

= ˆ I ˆ U

(n)

(tk+1, tk) = (−i)n tk+1

tk

dt1 t1

tk

dt2 . . . tn−1

tk

dtn ˆ H(t1)ˆ H(t2) . . . ˆ H(tn) n = 1, 2, ...

Kenneth Hansen December 17, 2015 3 / 44

AARHUS UNIVERSITY

AU

If the number of intervals is sufficiently large it will be a good approximation that, ˆ H(t) = ˆ H(tk) for tk ≤ t ≤ tk+1. The expansion of the time-evolution operator then collapses to, ˆ U(tk+1, tk) =

∞

• n=0

ˆ U

(n)

(tk+1, tk) =

∞

• n=0

(−i)(n) [(∆t)ˆ H(tk)]n n! . This is per definition an exponential function expansion and can be contracted as, ˆ U(tk+1, tk) = exp[−i∆t ˆ H(tk)]. Numerically evolving states is now performed using this exponential

• perator on a state. There are different ways of doing this and we will now

present some of them.

Kenneth Hansen December 17, 2015 4 / 44

AARHUS UNIVERSITY

AU

The Explicit Euler Method

Instead of working with ˆ H(tk) one can improve most methods by working with ˆ H(tk+1/2). (Working at mid-points/in the mean) The first order expansion of the exponential results in the explicit Euler scheme: |Ψ(tk+1/2) =

• ˆ

I − i∆t 2 ˆ H(tk+1/2)

• |Ψ(tk)

This scheme has an error of O((∆t)2) but is extremely unstable. One could go to a higher order in the expansion an gain an error in the norm of O((∆t)3), but since it still isn’t unitary and it still doesn’t suppress non physical states (these will grow exponentially) this is not a good path to go.

Kenneth Hansen December 17, 2015 5 / 44

AARHUS UNIVERSITY

AU

The Implicit Euler Method

Instead of increasing the order of the expansion we use that: ˆ U†(tk+1, tk) = ˆ U−1(tk+1, tk) = ˆ U(tk, tk+1) = exp[i∆t ˆ H(tk+1)], and get the implicit equation: ˆ U(tk, tk+1) |Ψ(tk+1) = |Ψ(tk) ⇒ |Ψ(tk+1) =

• ˆ

I + i∆t 2 ˆ H(tk+1/2) −1 |Ψ(tk) The implicit method is a lot more stable as unphysical components don’t grow exponentially like they do in the explicit method. It is still an error of O((∆t)2) method and still lacks unitarity.

Kenneth Hansen December 17, 2015 6 / 44

AARHUS UNIVERSITY

AU

The Crank-Nicolson Method

The implicit method is better than the explicit method but it still isn’t unitary. Combining the implicit and explicit method into one we get the Crank-Nicolson form: |Ψ(tk+1) =

• ˆ

I + i∆t 2 ˆ H(tk+1/2) −1 ˆ I − i∆t 2 ˆ H(tk+1/2)

• |Ψ(tk)

Which is unitary! (realize that ˆ U†(tk+1, tk) = ˆ U(tk+1, tk)) This is also known as Cayley’s form of the complex-exponential for time-evolution. The Crank-Nicolson scheme has an error of O((∆t)3) and is unconditionally stable.

Kenneth Hansen December 17, 2015 7 / 44

AARHUS UNIVERSITY

AU

t or n x or j FTCS (a) Fully Implicit (b) Crank-Nicolson (c)

0Numerical Recipes, 2. ed., page 850 Kenneth Hansen December 17, 2015 8 / 44

AARHUS UNIVERSITY

AU

The Crank-Nicolson Method - Numerically

The Crank-Nicolson method is used with a grid-based representation of the wave function. Remembering the Schr¨

• dinger Equation in a length gauge:

i ∂ ∂t Ψ(x, t) =

• −1

2 ∂2 ∂x2 + ˆ V (x) + ˆ E(t)ˆ x

• Ψ(x, t)

We then use the second-order central difference formula: ∂2 ∂x2 Ψ(xj, t) = Ψ(xj+1, t) − 2Ψ(xj, t) + Ψ(xj−1, t) (∆x)2 , and get evolution matrices U±(tk+1/2) whose elements are:

U±

j,j′(tk+1/2) =

     ±i∆t/[2(∆x)2], j = j′ + 1, j′ − 1, 1 ∓ i∆t[1/(∆x)2 + V (xj) + E(tk+1/2)xj], j = j′, 0, else.

Kenneth Hansen December 17, 2015 9 / 44

AARHUS UNIVERSITY

AU

Using the Crank-Nicolson scheme consists now in solving the linear equations with appropriate boundary conditions: U−(tk+1/2)Ψ Ψ Ψ(tk+1) = U+(tk+1/2)Ψ Ψ Ψ(tk). It is seen that U−(tk+1/2) is tridiagonal and an inverse can therefore be found without much computational effort.

Kenneth Hansen December 17, 2015 10 / 44

AARHUS UNIVERSITY

AU

Using the Crank-Nicolson method

Good things about this method: Accurate to the second order Unconditionally stable Unitary Can be computationally efficient (O(n2)) For this to be an effective method it has to be brought into a tridiagonal (or band diagonal/Toeplitz) form to simplify calculating the inverse (solving the implicit equation) This is almost done automatically in 1D, but will require effort in more dimensions!

Kenneth Hansen December 17, 2015 11 / 44

AARHUS UNIVERSITY

AU

Alternating Direction Implicit Method (ADI)

For multi-dimensional systems a split-step can be performed to simplify the problem.

i,j i,j-1 i,j+1 i,j i,j i,j-1 i,j+1 n+1 n+1/2 n

x y z

i-1,j i+1,j

Kenneth Hansen December 17, 2015 12 / 44

AARHUS UNIVERSITY

AU

References

A good introduction to numerical integration of wave functions: Atoms in Intense Laser Fields, C.J. Joachian, N.J. Kylstra, R.M. Potvliege, 2012 General introduction to Crank-Nicolson from computational standpoint: Numerical Recipes, W.H. Press et.al.

Kenneth Hansen December 17, 2015 13 / 44

AARHUS UNIVERSITY

AU

End of Part I Thank You!

Kenneth Hansen December 17, 2015 14 / 44

AARHUS UNIVERSITY

AU

Focus Topic: Propagators for the TDSE Part II: Split-Operator Method

Haruhide Miyagi

QUSCOPE Meeting, Aarhus University, Denmark

December 17, 2015

Please have a look at the handouts!

Haruhide Miyagi December 17, 2015 15 / 44

AARHUS UNIVERSITY

AU

Focus Topic: Propagators for the TDSE Part III: Arnoldi-Lanczos Propagator

Chuan Yu

QUSCOPE Meeting, Aarhus University, Denmark

December 17, 2015

Chuan Yu December 17, 2015 16 / 44

AARHUS UNIVERSITY

AU

Arnoldi-Lanczos Propagator

Short-time propagator |Ψ(t + ∆t) = exp[−i ˆ H(t)∆t]|Ψ(t)

Exact calculation of matrix exponential scales cubically with dimension, feasible only for small systems.

Idea of Arnoldi-Lanczos method

1 Truncate Taylor expansion to some order L, and span Krylov subspace

by a small set of vectors ˆ Hk |Ψ , (k = 0, . . . , L) |Ψ(t + ∆t) ≈

L

• n=0

(n!)−1 −i ˆ H(t)∆t n |Ψ(t) .

2 Find an orthonormal set of vectors |Qk , (k = 0, . . . , L) of Krylov

subspace with modified Gram-Schmidt algorithm, starting from a normalized state |Q0 = |Ψ(t) /Ψ(t) = |Ψ(t) /

• Ψ(t)|Ψ(t).

Chuan Yu December 17, 2015 17 / 44

AARHUS UNIVERSITY

AU

Arnoldi-Lanczos Propagator (cont’d)

Idea of Arnoldi-Lanczos method (cont’d)

3 Approximate wavefunction and Hamiltonian as

|Ψ(t + ∆t) ≈

L

• j=0

|Qj Qj|Ψ(t + ∆t) , ˆ H ≈

L

• j=0

L

• k=0

|Qj Qj| ˆ H |Qk Qk| , where Qj| ˆ H |Qk is matrix element of reduced Hamiltonian HL

4 Diagonalize HL with a similarity transformation S−1HLS = λ 5 Evaluate coefficents cj = Qj|Ψ(t + ∆t) = Qj| exp(−i ˆ

H∆t) |Ψ(t)

Chuan Yu December 17, 2015 18 / 44

AARHUS UNIVERSITY

AU

Arnoldi-Lanczos Propagator (cont’d)

Final expression

|Ψ(t + ∆t) ≈

L

• j=0

|Qj Qj| exp(−i ˆ H∆t) |Ψ(t) ≈

L

• j=0

|Qj Qj|

L

• k=0

L

• l=0

|Qk [exp (−iHL∆t)]kl Ql|Q0 Ψ(t) =

L

• j=0

|Qj

L

• k=0

[S]jk exp [−iλk∆t]

• S−1

k0 Ψ(t).

Reduced Hamiltonian HL is of dimension (L + 1) × (L + 1), which is much smaller than Hamiltonian H of dimension N × N. Krylov dimension should be chosen to ensure convergence, i.e., |Ψ(t + ∆t) should be well-described by {|Q0 , · · · , |QL}.

Chuan Yu December 17, 2015 19 / 44

AARHUS UNIVERSITY

AU

Arnoldi-Lanczos Algorithm

Starting from a normalized vector |Q0 = |Ψ(t) /

• Ψ(t)|Ψ(t),

for k = 0, . . . , L do |Φ(0)

k+1 = ˆ

H |Qk for j = 0, . . . , k do βj,k = Qj|Φ(j)

k+1

|Φ(j+1)

k+1 = |Φ(j) k+1 − |Qj βj,k

end for βk+1,k = Φ(k+1)

k+1 =

• Φ(k+1)

k+1 |Φ(k+1) k+1

|Qk+1 = |Φ(k+1)

k+1 /βk+1,k

end for gives an orthonormal basis set {|Q0 , · · · , |QL} together with matrix elements Qj| ˆ H |Qk =

• βj,k

for j ≤ k + 1, else.

Chuan Yu December 17, 2015 20 / 44

AARHUS UNIVERSITY

AU

Properties

General Unitary Stable Krylov dimension is related to system Hamiltonian and time step size. For small time step size, usually Krylov dimension less than 20 can give converged results. For Hamiltonian matrix of dimension N × N, computational complexity is O(N2). If Hamiltonian matrix is sparse, e.g., in finite-element discrete variable representation (FEDVR), computational complexity is O(N).

Chuan Yu December 17, 2015 21 / 44

AARHUS UNIVERSITY

AU

References

For Arnoldi-Lanczos

Unitary quantum time evolution by iterative Lanczos reduction,

• T. J. Park and J. C. Light, J. Chem. Phys. 85, 5870 (1986)

Explicit schemes for time propagating many-body wave functions,

• A. L. Frapiccini, A. Hamido, S. Schr¨
• ter, D. Pyke, F. Mota-Furtado,
• P. F. O’Mahony, J. Madro˜

nero, J. Eiglsperger, and B. Piraux,

• Phys. Rev. A 89, 023418 (2014)

For FEDVR

Numerical grid methods for quantum-mechanical scattering problems,

• T. N. Rescigno and C. W. McCurdy, Phys. Rev. A 62, 032706 (2000)

Parallel solver for the time-dependent linear and nonlinear Schr¨

• dinger equation, B. I. Schneider, L. A. Collins, and S. X. Hu,
• Phys. Rev. E 73, 036708 (2006)

Chuan Yu December 17, 2015 22 / 44

AARHUS UNIVERSITY

AU

End of Part III Thank You!

Chuan Yu December 17, 2015 23 / 44

AARHUS UNIVERSITY

AU

Focus Topic: Propagators for the TDSE Part IV: Polynomial Propagators - Chebyshev (et al.)

Daniel Reich

QUSCOPE Meeting, Aarhus University, Denmark

December 17, 2015

Daniel Reich December 17, 2015 24 / 44

AARHUS UNIVERSITY

AU

Many Ways to Solve the Schr¨

• dinger Equation

solve TDSE ≡ compute unitary propagator ˆ U (t) = e−i ˆ

Ht

but we learned: matrix exponentials are expensive (O

• N3

)

1 Taylor series to first order, adjust for unitarity =

⇒ Crank-Nicolson

2 split exponential into two parts for which computation is simple

( ˆ T diagonal in k-space, ˆ V diagonal in r-space) = ⇒ Split-Step

3 compute Taylor series to higher order but employ projected

Hamiltonian on a subspace of full Hilbert space = ⇒ Arnoldi-Lanczos

4 ??? Daniel Reich December 17, 2015 25 / 44

AARHUS UNIVERSITY

AU

Polynomial Propagators

idea: expand exponential in polynomial e−i ˆ

Ht |Ψ = ∞

• n=0

cnpn

• −i ˆ

Ht

• |Ψ

where pn is a polynomial of degree n matrix exponentiation problem (O

• n3

) transformed to (repeated) matrix-vector multiplication (O

• n2

) important from a numerical perspective (O

• n3
• vs. O
• n2

) ˆ H2 |Ψ =

• ˆ

H · ˆ H

• |Ψ

ˆ H2 |Ψ = ˆ H

• ˆ

H |Ψ

• Daniel Reich

December 17, 2015 26 / 44

AARHUS UNIVERSITY

AU

Taylor Series

expand exponential in polynomial ... ... Taylor series! e−i ˆ

Ht |Ψ = ∞

• n=0
• −i ˆ

Ht n n! |Ψ simple form, convergence radius = ∞, trivial to compute corresponds to Arnoldi-Lanczos but without any Hilbert space reduction

Daniel Reich December 17, 2015 27 / 44

AARHUS UNIVERSITY

AU

We’re done!

e−i ˆ

Ht |Ψ = ∞

• n=0
• −i ˆ

Ht n n! |Ψ

DONE!

Daniel Reich December 17, 2015 28 / 44

AARHUS UNIVERSITY

AU

We’re done?

e−i ˆ

Ht |Ψ = ∞

• n=0
• −i ˆ

Ht n n! |Ψ

DONE?

Daniel Reich December 17, 2015 29 / 44

AARHUS UNIVERSITY

AU

Uniformity Matters!

many polynomial approximations are troublesome at the domain boundary (cf. Runge’s phenomenon) issue with Taylor: each order has very unequal distribution along the interval (monomials xn have decreasing contribution in the center for increasing n) find a polynomial πn (x) = n

k=0 ckpk (x) that approximates a

function f uniformly, i.e. minimise for all n πn − f ∞ = max

x∈[−1,1] πn (x) − f (x)

make sure that increasing the order leads to uniform improvement make sure that pk∞ is small, i.e. added orders pk stay close around zero everywhere

Daniel Reich December 17, 2015 34 / 44

AARHUS UNIVERSITY

AU

Chebyshev Optimality

Chebyshev Optimality

Let Πn be the set of monic polynomials of order n in the interval [−1, 1], i.e. polynomials of the form pn (x) = xn + n−1

k=0 ckxk. Then, the

normalised Chebyshev polynomials, Tn (x) = 1 2n−1 cos (n arccos (x)) , minimise πn∞ among all members of Πn (with Tn (x) ∞ =

1 2n−1 ).

Daniel Reich December 17, 2015 35 / 44

AARHUS UNIVERSITY

AU

Chebyshev Propagator

e−i ˆ

Ht |Ψ = ∞

• n=0

cnpn

• −i ˆ

Ht

• |Ψ

domain of pn? what does pn(ˆ H) mean anyway? e−i ˆ

Ht |Ψ = ∞

• n=0
• l

cnpn (−iElt) |El El | Ψ for Hermitian ˆ H, domain of pn should be the imaginary axis map real axis to imaginary axis = ⇒ complex Chebyshev polynomials φn (x) = Tn (−ix) Tn ∈ [−i, i] = iR?, but in numerics the spectrum of ˆ H is compact! normalise Hamiltonian (spectral radius ∆, minimal eigenvalue Emin) ˆ Hnorm = 2 ˆ H − 1

2∆ + Emin

• ∆

Daniel Reich December 17, 2015 36 / 44

AARHUS UNIVERSITY

AU

Chebyshev Propagator

e−iHdt |Ψ =

• n

cnφn

• −i ˆ

Hnorm

• |Ψ

cn = (2 − δn0) e−i( ∆

2 +Emin)dtJn

∆ 2 dt

• Jn (x) are the Bessel functions of first kind

for n > ∆

2 dt : Jn

∆

2 dt

• −

→ 0 exponentially fast = ⇒ together with φn ∈ [−1, 1] a priori knowledge of truncation threshold (roughly nmax ≃ 2∆dt for standard applications) easy evaluation of φn with recursion relation φn(ˆ A) = 2ˆ Aφn−1(ˆ A) − φn−2(ˆ A), φ1(ˆ A) = ˆ A, φ0(ˆ A) = 1

• nly one additional appication of Hamiltonian per order, coefficients

can be reused when propagation is continued (but watch the spectral radius for TDSE solving!)

Daniel Reich December 17, 2015 37 / 44

AARHUS UNIVERSITY

AU

Algorithm

[from M. H. Goerz, “Optimizing Robust Quantum Gates in Open Quantum Systems”]

Daniel Reich December 17, 2015 38 / 44

AARHUS UNIVERSITY

AU

Advantages and Disadvantages

mathematically optimal polynomial propagator exact solution of the TDSE for time-independent Hamiltonians, no further information required (except for spectral range) independent on representation of wave function, just requires knowledge of application of ˆ H simple, easy-to-implement algorithm generalisable to imaginary time propagation (Tk (x) instead of φk (x)) very sensitive to wrongly estimated spectral range (evidence: loss of norm conservation, can be caught during algorithm with some

• verhead)
• nly usable for purely real (or purely imaginary) spectrum of ˆ

H requires approximation of explicit time-dependence in TDSE not applicable when inhomogeneities or non-linearities in the TDSE appear (e.g. Gross-Pitaevskii equation)

Daniel Reich December 17, 2015 39 / 44

AARHUS UNIVERSITY

AU

You might also find interesting ...

• R. Kosloff, “Time-Dependent Quantum-Mechanical Methods for Molecular

Dynamics”, J. Phys. Chem. 92, 2087 (1988) (nice, but slightly dated overwiew on propagation methods in the context of atomic and molecular physics)

• H. Tal-Ezer, R. Kosloff, “An accurate and efficient scheme for propagating the

time dependent Schr¨

• dinger equation”, JCP 81, 3967 (1984) (reference paper for

the Chebyshev propagator)

• M. H. Goerz, “Optimizing Robust Quantum Gates in Open Quantum Systems”,

Dissertation Universit¨ at Kassel (2015) (concise discussion of Chebyshev, Newton and Newton-Arnoldi propagator with pseudocodes for all algorithms)

• A. Gil et al., “Numerical Methods for Special Functions”, Society for Industrial

and Applied Mathematics (2007) (mathematical discussion on properties of Chebyshev polynomials and why they are useful)

• H. Tal-Ezer, “Polynomial approximation of functions of matrices and applications”,

Journal of Scientific Computing 4, 25 (1989) (mathematical discussion on polynomial approximation method for a variety of applications)

Daniel Reich December 17, 2015 40 / 44

AARHUS UNIVERSITY

AU

Extensions

Chebychev method can also be applied to other equations of motion, as long as an analytical solution can be written down and expanded in Chebychev polynomials expansion coefficients are not Bessel functions anymore, they must be derived using a cosine transform example: Chebychev propagator for the inhomogeneous TDSE [M. Ndong et al., JCP 130, 124108 (2009)] remarkable extension: use iterative time ordering to rewrite time-dependent TDSE as time-independent Hamiltonian TDSE with inhomogeneous source term = ⇒Chebyshev propagator for inhomogeneous Schr¨

• dinger equation

can eliminate requirement of approximating explicit time-dependence in ˆ H [M. Ndong et al., JCP 132, 064105 (2010)]

Daniel Reich December 17, 2015 41 / 44

AARHUS UNIVERSITY

AU

Newton Polynomials

what can we do for a truly complex spectrum, e.g. open quantum system evolutions following a Liouvillian instead of a Hamiltonian? ∂ ∂t ρ = Lρ solution still exponential (time-independent Liouvillian for simplicity) ρ (t) = eLtρ (0) approximate in complex plane by interpolating Newton polynomial f (z) =

• n

anRn (z) , Rn (z) =

n−1

• j=0

(z − zj) with sampling points {zj} fastest convergence when complex eigenvalue of L are used for zj but: diagonalisation problem for Liouvillians is even more complicated than for Hamiltonians!

Daniel Reich December 17, 2015 42 / 44

AARHUS UNIVERSITY

AU

Newton-Arnoldi

solution: estimate spectral domain, encircle it with a rectangle or ellipse, calculate expansion coefficients from sampling points on that boundary even better: use Arnoldi algorithm to obtain some approximate eigenvalues further reading:

• G. Ashkenazi et al., “Newtonian propagation methods applied to the

photodissociation dynamics of I −

3 ”, JCP 103, 10005 (1995)

(introduction to Newton propagator)

• H. Tal-Ezer, “On Restart and Error Estimation for Krylov

Approximation of w = f (A) v”, SIAM J. Sci. Comput. 29, 2426 (2007) (modification towards restarted Arnoldi-Newton propagator)

Daniel Reich December 17, 2015 43 / 44

AARHUS UNIVERSITY

AU

Thank You!

Daniel Reich December 17, 2015 44 / 44