Corrections to Wigner type phase space methods
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.
[section] [section] [section] [section] [section] [section]
Molecular quantum systems are described by an unbounded self-adjoint operator, the Schrödinger operator
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 . 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
with the appropriate -scaling of the time-derivative or – equivalently – by the action of the unitary evolution operator , since
Even though the Schrödinger equation 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
the associated Hamiltonian system
and the corresponding Hamiltonian flow
as the classical counterparts to the Schrödinger operator , the time-dependent Schrödinger equation , 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 Appendix 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  when developing the thermodynamics of quantum mechanical systems. Crucial properties of the Wigner function are the orthogonality relation
and the asymptotically classical time evolution
for the solution of the Schrödinger equation . The Wigner phase space method of E. Heller  and the statistical quasiclassical method of H. Lee and M. Scully , for example, use these properties for computing the transition probabilities from a given state to according to
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
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 . Often, the Weyl operator is also called the quantum observable assocaited with the classical observable . Besides the convenient relation , Weyl quantization enjoys the beautiful property that the trace of the product of Weyl operators can be expressed as the integral
From the Weyl point of view, the asymptotic time evolution of the Wigner function can equivalently be formulated as
This quantum-classical approximation of time evolved quantum observables is known as Egorov’s theorem and has first been formulated in  within the Hörmander theory of pseudo-differential operators, see also  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 , who approximate the time-dependent correlation function of and according to
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  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  have proposed to improve the second order approximation for the Wigner function of by with
where the derivative of the effective potential is given by
For the correcting factor they mention the three options
and they use the term Wigner trajectories for the solutions of the ordinary differential equation , , see also . Later, A. Donoso and C. Martens  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 .
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 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
where denotes a generalization of the usual Poisson bracket involving derivatives of the functions and up to order , see . In the context of response theory, M. Kryvohuz and J. Cao  use this expansion up to the fourth term for a systematic improvement of linear response computations over long times. The commutator expansion 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 .
For general Hamiltonian functions , the expansion implies a higher order version of Egorov’s theorem,
with leading order term and corrections
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 , when deriving higher order estimates for the Wigner function of the wave function .
Since the higher order Egorov expansion is built of even powers of the parameter , the first correction
provides the fourth order estimate
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
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
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 , which is then applied for numerical test calculations.
This paper is organized as follows. develops Egorov’s theorem to the next order with respect to the parameter and formulates this correction as a first order ordinary differential equation. discusses a discretization of this corrected approximation for the computation of expectation values. 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.
2Higher order Corrections
Let be a smooth function of subquadratic growth. That is, for all with there exists with
We look for an approximate classical observable such that
We require at time , rewrite the difference according to
and observe that the commutator
plays a crucial rule. Multiplying this commutator with , we obtain an asymptotic expansion in even powers of ,
denotes the th generalized Poisson bracket of two smooth functions , see Appendix A. Since the first generalized Poisson bracket coincides with the usual Poission bracket, , we have
Let be the flow associated with the classical Hamilton function , and set . Then,
and we obtain Egorov’s theorem,
2.1The first correction
The asymptotic commutator expansion allows to systematically derive higher order corrections to Egorov’s theorem. For constructing the first correction, we set . The computation
suggests to choose
Indeed, we compute the time derivative,
where the last equation uses that the classical flow as a symplectic transformation of phase space preserves the Poisson bracket. We therefore obtain
2.2Ordinary Differential Equations for the Correction
Our next aim is to reformulate the correction as
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.
We observe that we can write
where the 3-tensor is defined by
A closer look at this 3-tensor reveals that it is symmetric:
Let . For , we set , and otherwise , . Then,
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 ,
Since the 3-tensor is symmetric and , we obtain the reformulation
Next, we compute the time derivatives of integrals, which are of the form observed in the defining equations of the time-dependent tensors , , and .
We start with
For introducing the -derivative in the integral, we compute
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 . Our result shows nonetheless, that the computation is feasible.
The initial values are clear, since , and are defined in as an integral from to , over a smooth integrand. Next, we study the time derivative of . For this, we write the Hamiltonian system as
and apply the differential operator to obtain the Jacobi stability equation
Then, by Lemma ?,
So we can finish the proof for by using that . The proofs for and are analogous: We differentiate the Jacobi stability equation to obtain
Moreover, we use
as well as and .
2.3Vectorization: General Hamiltonians
For the numerical simulation of the ordinary differential system , we vectorize the tensors. Recall, that for matrices and , the Kronecker product is defined as
The vectorization of a k-tensor is defined as
That is, if
The following observation allows the vectorization of the products occuring on the right hand side of our differential equation .
We define the matrix . That is,
for and . Then,
Now we reformulate the results of Theorem ? in vectorized form.
Applying Lemma ? to the part of the ordinary differential equation, we obtain