Corrections to Wigner type phase space methods

# Corrections to Wigner type phase space methods

Wolfgang Gaim and Caroline Lasser Fachbereich Mathematik, Eberhard Karls Universität Tübingen, 72076 Tübingen, Germany Zentrum Mathematik, Technische Universität München, 80290 München, Germany
###### Abstract

Over decades, the time evolution of Wigner functions along classical Hamiltonian flows has been used for approximating key signatures of molecular quantum systems. Such approximations are for example the Wigner phase space method, the linearized semiclassical initial value representation, or the statistical quasiclassical method. The mathematical backbone of these approximations is Egorov’s theorem. In this paper, we reformulate the well-known second order correction to Egorov’s theorem as a system of ordinary differential equations and derive an algorithm with improved asymptotic accuracy for the computation of expectation values. For models with easily evaluated higher order derivatives of the classical Hamiltonian, the new algorithm’s corrections are computationally less expensive than the leading order Wigner method. Numerical test calculations for a two-dimensional torsional system confirm the theoretical accuracy and efficiency of the new method.

###### ams:
81S30, 81Q20, 81-08, 65D30, 65Z05
###### pacs:
82.20.Ln, 03.65.Sq, 34.10.+x
: Nonlinearity

,

## 1 Introduction

Molecular quantum systems are described by an unbounded self-adjoint operator, the Schrödinger operator

 H=−ε22Δ+V,

acting on the Hilbert space of complex-valued square-integrable functions . Here, is a small positive parameter related to the inverse square root of the average nuclear mass, and is the nuclear potential function resulting from the Born–Oppenheimer approximation, see [20]. The relevant time scales of nuclear quantum motion are of the order . The time evolution of an initial wave function is governed by the time-dependent linear Schrödinger equation

 \rmiε∂tψt=Hψt (1)

with the appropriate -scaling of the time-derivative or – equivalently – by the action of the unitary evolution operator , since

 ψt=\rme−\rmiHt/εψ0(t∈R).

Even though the Schrödinger equation (1) is a linear partial differential equation, the numerical simulation of physical quantities derived from the wave function is notoriously difficult for two reasons: The dimension of the nuclear configuration space is large. If a molecule consists of nuclei, then . Just for a single water molecule, for example, we have . Moreover, nuclear quantum motion is highly oscillatory. Solutions typically oscillate with frequencies of the order in time and space, while ranges between and , depending on the molecular system under consideration. For the hydrogen molecule H, for example, one has , while for iodine monobromide IBr.

As an answer to these challenges, chemical physicists have developed approximations involving the nonlinear time evolution of classical mechanics. One uses the Hamiltonian function

 h:Rd×Rd→R,h(q,p)=12|p|2+V(q), (2)

the associated Hamiltonian system

 ˙qt=∂ph(qt,pt),˙pt=−∂qh(qt,pt), (3)

and the corresponding Hamiltonian flow

 Φt:R2d→R2d

as the classical counterparts to the Schrödinger operator , the time-dependent Schrödinger equation (1), and the unitary evolution operator . This quantum-classical correspondence is most elegantly elaborated by using Wigner functions and Weyl quantization.

The Wigner function of a square integrable function is a real-valued, square integrable, continuous function on phase space obtained by an inverse Fourier transform of the autocorrelation function of , see the short summary in A. The Wigner function can be thought of as a probability density on phase space , though in general it might attain negative values. It has been introduced by E. Wigner in [25] when developing the thermodynamics of quantum mechanical systems. Crucial properties of the Wigner function are the orthogonality relation

 |⟨ϕ,ψ⟩|2=(2πε)d∫R2dWϕ(z)Wψ(z)\rmdz(ϕ,ψ∈L2(Rd))

and the asymptotically classical time evolution

 Wψt=Wψ0∘Φ−t+O(ε2),ε→0, (4)

for the solution of the Schrödinger equation (1). The Wigner phase space method of E. Heller [8, 2] and the statistical quasiclassical method of H. Lee and M. Scully [13], for example, use these properties for computing the transition probabilities from a given state to according to

 |⟨ϕ,ψt⟩|2≈(2πε)d∫R2dWϕ(Φt(z))Wψ0(z)\rmdz.

Viewing the Wigner function of a square integrable function as a tempered distribution, that is, as a continuous linear mapping from the Schwartz functions to the real numbers, one arrives at the identity

 ∫R2da(z)Wψ(z)\rmdz=⟨ψ,op(a)ψ⟩, (5)

where denotes the bounded linear operator on obtained by the Weyl quantization of the function . The Weyl operator is a pseudo-differential operator, a generalized partial differential operator, which treats differentiation and multiplication by functions on the same footing, see the appendix. With an appropriate handling of operator domains, Weyl quantization also applies for unbounded linear operators, and the Schrödinger operator , for example, is the Weyl quantized classical Hamiltonian function defined in (2). Often, the Weyl operator is also called the quantum observable assocaited with the classical observable . Besides the convenient relation (5), Weyl quantization enjoys the beautiful property that the trace of the product of Weyl operators can be expressed as the integral

 \tr(op(a)op(b))=(2πε)−d∫R2da(z)b(z)\rmdz,

provided that the classical observables and satisfy appropriate regularity and growth conditions, see e.g. [4, Proposition 9.2] or [7, Proposition 284].

From the Weyl point of view, the asymptotic time evolution of the Wigner function (4) can equivalently be formulated as

 \rme\rmiHt/εop(a)\rme−\rmiHt/ε=op(a∘Φt)+O(ε2),ε→0. (6)

This quantum-classical approximation of time evolved quantum observables is known as Egorov’s theorem and has first been formulated in [5] within the Hörmander theory of pseudo-differential operators, see also [19, 1, 26] for refined error estimates in the context of semiclassical microlocal analysis. The computational power of classically propagating quantum observables has been recognized by W. Miller and H. Wang [17, 24], who approximate the time-dependent correlation function of and according to

 \tr(\rme\rmiHt/εop(a)\rme−\rmiHt/εop(b))≈(2πε)−d∫R2da(Φt(z))b(z)\rmdz.

In the chemical physics literature this approximation is referred to as the linearized semiclassical initial value representation (LSC-IVR), and it is derived from Fourier integral operator representations of the unitary evolution operator , see also [22] for a recent review. It seems that the link of LSC-IVR to Egorov’s theorem has not been noticed so far.

In a study of the one-dimensional Morse oscillator, H. Lee and M. Scully [14] have proposed to improve the second order approximation (4) for the Wigner function of by with

 ∂tW=−p∂qW+V′eff∂pW,

where the derivative of the effective potential is given by

 V′eff=V′−ε224γV′′′.

For the correcting factor they mention the three options

 γ1=∂3p(Wψ0∘Φ−t)/∂pW,γ2=∂3p(Wψ0∘Φ−t)/∂p(Wψ0∘Φ−t),γ3=∂3pW/∂pW,

and they use the term Wigner trajectories for the solutions of the ordinary differential equation , , see also [11, 12]. Later, A. Donoso and C. Martens [3] have defined their entangled classical trajectories by a fourth variant . Trajectories generated by and and generalisations thereof have also been incorporated for the numerical simulation of time-dependent correlation functions of one-dimensional systems by J. Liu and W. Willer [15].

Our aim here is also a higher order approximation of the dynamics. However, we will work with the rigorous mathematical framework of Egorov’s theorem (6) and construct phase space trajectories without any entanglement for arbitrary dimensions . The key element for proving Egorov’s theorem is an asympotic expansion of the commutator

 \rmiε[op(h),op(a)]∼∑k∈2N(ε2\rmi)kop({h,a}k+1). (7)

where denotes a generalization of the usual Poisson bracket involving derivatives of the functions and up to order , see Section 2. In the context of response theory, M. Kryvohuz and J. Cao [9] use this expansion up to the fourth term for a systematic improvement of linear response computations over long times. The commutator expansion (7) also reveals, that for Hamiltonians , wich are polynomials of degree less or equal than two, the remainder of Egorov’s theorem vanishes. This exact Egorov result is utilized by H. Waalkens, R. Schubert and S. Wiggins for their quantum normal form algorithm in dynamical transition state theory [23].

For general Hamiltonian functions , the expansion (7) implies a higher order version of Egorov’s theorem,

 \rme\rmiHt/εop(a)\rme−\rmiHt/ε∼∑k∈2Nεkop(ak(t)) (8)

with leading order term and corrections

 ak(t)=∑l∈{0,2,…,k−2}(\rmi2)k−l∫t0{h,al(τ)}k+1−l∘Φt−τ\rmdτ.

It is remarkable that the corrections are only built from derivatives of the observable , the Hamiltonian function  and the flow , which has also been emphasized by M. Pulvirenti in [18], when deriving higher order estimates for the Wigner function of the wave function .

Since the higher order Egorov expansion (8) is built of even powers of the parameter , the first correction

 a2(t)=−14∫t0{h,a∘Φτ}3∘Φt−τ\rmdτ (9)

provides the fourth order estimate

 \rme\rmiHt/εop(a)\rme−\rmiHt/ε=op(a∘Φt+ε2a2(t))+O(ε4).

The aim of this paper is the exemplary analysis of this correction and its discretization for the numerical computation of expectation values. In a first step, we reformulate it as

 a2(t)=Λt1(Da∘Φt)+Λt2(D2a∘Φt)+Λt3(D3a∘Φt),

where

 Λtk:R2d×⋯×2d→R(k=1,2,3),

is an explicitly defined linear mapping from the space of -tensors to the real numbers, which depends on derivatives of the Hamiltonian function and the flow , but not on the observable . In spirit, this formulation of is close to the Wigner trajectories generated by the first correction factor . In the next step, we derive a first order ordinary differential system for the components of , , and . The vectorization of this equation results in

 ddtΛt=NtΛt+Ct, (10)

where the building blocks of the matrix and the vector are components of the tensors , , and evaluated along the flow . In the final step we derive a fourth order splitting scheme for the discretization of the ordinary differential equation (10), which is then applied for numerical test calculations.

This paper is organized as follows. Section 2 develops Egorov’s theorem to the next order with respect to the parameter and formulates this correction as a first order ordinary differential equation. Section 3 discusses a discretization of this corrected approximation for the computation of expectation values. Section 4 provides numerical experiments for a two-dimensional torsional quantum system, which confirm the theoretical considerations. The appendix summarizes basic properties of Wigner functions and Weyl operators.

## 2 Higher order Corrections

Let be a smooth function of subquadratic growth. That is, for all with there exists with

 ∥Dγh∥≤Cγ.

Let be a Schwartz function. Then, the following formal considerations can be turned into a proof according to [19, Théorème IV.10] or [1, Theorem 1.2].

We look for an approximate classical observable such that

 \rme\rmiHt/εop(a)\rme−\rmiHt/ε≈op(aappr(t)).

We require at time , rewrite the difference according to

 \rme\rmiHt/εop(a)\rme−\rmiHt/ε−op(aappr(t))=∫t0\rmd\rmdτ(\rme\rmiHτ/εop(aappr(t−τ))\rme−\rmiHτ/ε)\rmdτ = ∫t0\rme\rmiHτ/ε(\rmiε[op(h),op(aappr(t−τ))]−\rmd\rmdtop(aappr(t−τ)))\rme−\rmiHτ/ε\rmdτ,

and observe that the commutator

 [op(h),op(aappr(t))]=op(h)op(aappr(t))−op(aappr(t))op(h)

plays a crucial rule. Multiplying this commutator with , we obtain an asymptotic expansion in even powers of ,

 \rmiε[op(h),op(aappr(t))]∼∑k∈2N(ε2i)kop({h,aappr(t)}k+1), (11)

where

 {f,g}k=∑|α+β|=k(−1)|β|α!β!∂αq∂βpg∂βq∂αpf.

denotes the th generalized Poisson bracket of two smooth functions , see A. Since the first generalized Poisson bracket coincides with the usual Poission bracket, , we have

 \rmiε[op(h),op(aappr(t))]=op({h,op(aappr(t)})+O(ε2).

Let be the flow associated with the classical Hamilton function , and set . Then,

 \rmd\rmdtaappr(t)={h,aappr(t)},

and we obtain Egorov’s theorem,

 \rme\rmiHt/εop(a)\rme−\rmiHt/ε=op(a∘Φt)+O(ε2).

### 2.1 The first correction

The asymptotic commutator expansion (11) allows to systematically derive higher order corrections to Egorov’s theorem. For constructing the first correction, we set . The computation

 \rmiε[op(h),op(aappr(t))]=op({h,aappr(t)})−ε24op({h,a0(t)}3)+O(ε4)

suggests to choose

 a0(t)=a∘Φt,a2(t)=−14∫t0{h,a0(τ)}3∘Φt−τ\rmdτ.

Indeed, we compute the time derivative,

 \rmd\rmdta2(t) = −14({h,a0(t)}3+∫t0\rmd\rmdt({h,a0(τ)}3∘Φt−τ)\rmdτ) = −14({h,a0(t)}3+∫t0{h,{h,a0(τ)}3}∘Φt−τ\rmdτ) = −14({h,a0(t)}3+{h,∫t0{h,a0(τ)}3∘Φt−τ\rmdτ}),

where the last equation uses that the classical flow as a symplectic transformation of phase space preserves the Poisson bracket.We therefore obtain

 \rmd\rmdtaappr(t)={h,aappr(t)}−ε24{h,a0(t)}3.

Consequently,

 \rme\rmiHt/εop(a)\rme−\rmiHt/ε=op(aappr(t))+O(ε4).
###### Remark 2.1.

If is an observable invariant along the Hamiltonian flow, that is, for all , and additionally , then the corrected approximation is also invariant,

 aappr(t)=a∘Φt−ε24∫t0{h,a∘Φτ}3∘Φt−τ\rmdτ=a

for all . In particular, mass () and energy () are conserved.

### 2.2 Ordinary Differential Equations for the Correction

Our next aim is to reformulate the correction as

 a2(t)=−14((D3a∘Φt)ijkΛtkji+3(D2a∘Φt)ijΓtji+(Da∘Φt)iΞti),

where the components of the time-dependent tensors , , and satisfy a first order system of coupled ordinary differential equations, that is independent of the observable  and can be efficiently solved alongside the Hamiltonian flow . Here and in the following, we use Einstein’s summation convention for notational brevity.

Let

 J=(0Id−Id0)∈R2d×2d.

We observe that we can write

 a2(t) = −14∫t0{h,a∘Φτ}3∘Φt−τ\rmdτ = −14∫t0((D3(a∘Φτ))lmn(J˜D3h)nml)∘Φt−τ\rmdτ,

where the 3-tensor is defined by

 (˜D3h)ijk=⎧⎪ ⎪⎨⎪ ⎪⎩16(D3h)ijk,i=j=k12(D3h)ijk,i=j≠k or i=k≠j or k=j≠i(D3h)ijk,else

and

 (J˜D3h)ijk=JilJjmJkn(˜D3h)lmn. (12)

A closer look at this 3-tensor reveals that it is symmetric:

###### Lemma 2.1.

Let and be a -tensor. Then,

 Ji1l1⋯JiklkAl1⋯lkis symmetric⇔Ais symmetric.

In particular, the 3-tensor defined in (12) is symmetric.

###### Proof.

For every

 JimlAi1⋯im−1lim+1…ik={Ai1⋯im−1(im+d)im+1…ik,im≤d−Ai1⋯im−1(im−d)im+1…ik,else.

Let . For , we set , and otherwise , . Then,

 Ji1l1⋯JiklkAl1⋯lk=(−1)α1⋯(−1)αkAβ1…βk.

Hence, a permutation of results in a permutation of and . Thus, is symmetric if and only if is symmetric. ∎

We compute the third derivative of by the chain rule, where we denote the derivatives of by ,

 (D3(a∘Φτ))lmn= (D3a∘Φτ)ijk(DΦτ)il(DΦτ)jm(DΦτ)kn +(D2a∘Φτ)ij(D2Φτ)ilm(DΦτ)jn+(D2a)ij(D2Φτ)iln(DΦτ)jm +(D2a∘Φτ)ij(D2Φτ)imn(DΦτ)jl+(Da∘Φτ)i(D3Φτ)ilmn.

Since the 3-tensor is symmetric and , we obtain the reformulation

 a2(t)=−14((D3a∘Φt)ijkΛtkji+3(D2a∘Φt)ijΓtji+(Da∘Φt)iΞti)

with

 Λtijk=∫t0[(DΦτ)il(DΦτ)jm(DΦτ)kn(J˜D3h)nml]∘Φt−τ\rmdτ,Γtij=∫t0[(D2Φτ)ikl(DΦτ)jm(J˜D3h)mlk]∘Φt−τ\rmdτ,Ξti=∫t0[(D3Φτ)ijkl(J˜D3h)lkj]∘Φt−τ\rmdτ. (13)

Next, we compute the time derivatives of integrals, which are of the form observed in the defining equations of the time-dependent tensors , , and .

###### Lemma 2.2.

Let and smooth functions. Then,

 \rmd\rmdt∫t0(bf(τ))∘Φt−τ\rmdτ=∫t0(b\rmd\rmdτf(τ))∘Φt−τ\rmdτ+(bf(0))∘Φt.
###### Proof.

 \rmd\rmdt∫t0(bf(τ))∘Φt−τ\rmdτ=bf(t)+∫t0b∘Φt−τ(Df(τ))T∘Φt−τ\rmd\rmdtΦt−τ\rmdτ +∫t0(Db)T∘Φt−τ\rmd\rmdtΦt−τf(τ)∘Φt−τ\rmdτ.

For introducing the -derivative in the integral, we compute

 \rmd\rmdτ(bf(τ)∘Φt−τ)=b∘Φt−τ((\rmd\rmdτf(τ))∘Φt−τ−(Df(τ))T∘Φt−τ\rmd\rmdtΦt−τ) −(Db)T∘Φt−τ\rmd\rmdtΦt−τf(τ)∘Φt−τ.

Therefore,

 \rmd\rmdt∫t0(bf(τ))∘Φt−τ\rmdτ =bf(t)+∫t0(b\rmd\rmdτf(τ))∘Φt−τ\rmdτ−∫t0\rmd\rmdτ((bf(τ))∘Φt−τ)\rmdτ =∫t0(b\rmd\rmdτf(τ))∘Φt−τ\rmdτ+(bf(0))∘Φt.

Now we are ready to formulate and prove our first main result, the explicit system of ordinary differential equations describing the second order correction to Egorov’s theorem. In M. Zworski’s recent monograph, the higher order terms are referred to “as difficult to compute”, see [26, §11.1]. Our result shows nonetheless, that the computation is feasible.

###### Theorem 2.1 (Second Correction).

Let be a smooth function of subquadratic growth, , and the Hamiltonian flow associated with . Let be a Schwartz function. Then, there exists a constant such that for all

 ∥∥\rme\rmiop(h)t/εop(a)\rme−\rmiop(h)t/ε−op(a0(t)+ε2a2(t))∥∥≤Cε4

where and

 a2(t)=−14((D3a∘Φt)ijkΛtkji+3(D2a∘Φt)ijΓtji+(Da∘Φt)iΞti)

and the functions , , , , solve the ordinary differential system

 \rmd\rmdtΛtijk=MtilΛtljk+MtjlΛtilk+MtklΛtijl+C1ijk(t),\rmd\rmdtΓtij=(C2i(t))klΛtlkj+MtilΓtlj+MtjlΓtil,\rmd\rmdtΞti=(C3i(t))jklΛtlkj+3(C2i(t))jkΓtkj+MtilΞtl,Λ0ijk=Γ0ij=Ξ0i=0, (14)

with and

 C1ijk(t)=(J˜D3h∘Φt)ijk,(C2i(t))jk=(J⋅D3h∘Φt)ijk, (C3i(t))jkl=(J⋅D4h∘Φt)ijkl,

where for a matrix and -tensor , is given by and is defined by (12).

###### Proof.

The initial values are clear, since , and are defined in (13) as an integral from to , over a smooth integrand. Next, we study the time derivative of . For this, we write the Hamiltonian system as

 \rmd\rmdtΦt=JDh∘Φt, (15)

and apply the differential operator to obtain the Jacobi stability equation

 \rmd\rmdt(DΦt)ij=(J⋅D2h∘Φt)im(DΦt)mj.

Then, by Lemma 2.2,

 \rmd\rmdtΛtijk = +[(DΦ0)il(DΦ0)jm(DΦ0)kn(J˜D3h)lmn]∘Φt = ∫t0[(J⋅D2h∘Φτ)iν(DΦτ)νl(DΦτ)jm(DΦτ)kn(J˜D3h)lmn]∘Φt−τ\rmdτ +∫t0[(DΦτ)il(J⋅D2h∘Φτ)jν(DΦτ)νm(DΦτ)kn(J˜D3h)lmn]∘Φt−τ\rmdτ +∫t0[(DΦτ)il(DΦτ)jm(J⋅D2h∘Φτ)kν(DΦτ)νn(J˜D3h)lmn]∘Φt−τ\rmdτ +(J˜D3h)ijk∘Φt.

So we can finish the proof for by using that . The proofs for and are analogous: We differentiate the Jacobi stability equation to obtain

 \rmd\rmdt(D2Φt)ijk = (J⋅D3h∘Φt)imn(DΦt)nk(DΦt)mj+(J⋅D2h∘Φt)im(D2Φt)mjk, \rmd\rmdt(D3Φt)ijkl = (J⋅D4h∘Φt)imnμ(DΦt)nk(DΦt)mj(DΦt)μl +(J⋅D3h∘Φt)imn(DΦt)mj(D2Φt)nkl +(J⋅D3h∘Φt)imn(DΦt)nk(D2Φt)mjl +(J⋅D3h∘Φt)imn(DΦt)nl(D2Φt)mjk+(J⋅D2h∘Φt)im(D3Φt)mjkl.

Moreover, we use

 \rmd\rmdτ((D2Φτ)ikl(