Evaluation of Integrals with Highly Oscillatory Integrands:
Quantum Dynamics with Path Integrals
The solutions of many quantum dynamics problems can be formulated in terms of path integrals. Indeed, there is a formulation of quantum mechanics (less familiar than the Schrödinger wave function formulation, the Heisenberg matrix formulation, or the Dirac Hilbert space formulation) based entirely on the use of path integrals (Feynman and Hibbs, 1965). Path integrals are closely related to some of the mathematical approaches used to describe classical Brownian motion. A prototypical mathematical problem and a class of methods for solving it are described first; then a more general discussion of the problem is presented.
Prototypical Problem
This discussion of the prototypical problem follows closely the discussion of Doll et al. (1988). Consider the integral I(t) =
dx 
(x) exp(itf(x)) where t is a real parameter
corresponding to the time, x is the set of coordinates for a many-dimensional space, and (x) and f(x) are functions in that space. The function 
(x) is positive semidefinite and can without loss of generality be taken as normalized to unity. Thus, it can be regarded as a probability density in the space. The function f(x) is real (for this protot
ype problem).
For problems of interest, the dimensionality of the space can be very large (from 1 to a hundred to thousands), and the complexity of the function precludes analytic integration. For any particular x, (x) and f(x) can be evaluated numerically, and typically there is much more computation required for the evaluation of
than f.
The range of the x integration may be bounded, but more often it is unbounded. For typical problems of interest, however, integrals of the formdx (x) (f(x))
exist for all positive integers n.
The goal is usually to evaluate the integral for a specific value of t (or for a set of values of t) by making use of calculated values of and f at appropriately chosen points
and to estimate the accuracy of the answer. In some cases, the desired quantity is the area under the function
or perhaps the Fourier transform
Techniques for evaluation of such integrals of I would be extremely worthwhile, but this discussion focuses on the evaluation of I(t) for specific values of t.
Discussion of the Problem
A standard way of attacking this problem would be a Monte Carlo integration scheme in which points x are generated with a probability distribution of (x) and f is ev
aluated at these points. Then
where N points are generated andx
is the nth point. For small values of t, this can be a practical procedure. For large values of t, however, this procedure is extreme
ly inefficient because of the large amount of cancellation among the various terms.
If t is large and f is not a constant, it is clear that the integrand is highly oscillatory; thus, small regions of space where f is varying will tend to contribute small net amounts to the integrand for large t. Nev ertheless, individual points in such a region will contribute to the sum an exponential whose magnitude is unity. Cancellation among the various terms will give an accurate and small answer only when each such small region is sampled many times in the Mo nte Carlo evaluation.
The regions of the space in which df(x)/dx is zero or small will make the major contributions to the integral for large t. Thus, a scheme that concentrates the Monte Carlo sampling at and near such regions would be desira ble. Such schemes are usually called "stationary-phase" methods, since at such points the phase of the exponential is stationary.
Stationary-Phase Monte Carlo Methods
A variety of sampling methods have been developed that focus on the stationary-phase points (Doll, 1984; Filinov, 1986; Doll et al., 1988). Here the approach of Doll et al. (1988), which is similar to the others, is outlined.
The integral of interest can be rewritten identically as
where
and P(y) is an arbitrary normalized probability distribution. (This result is exact if the range of integration is over all positive and negative values of the integration variables.) The function P(y) is typically chosen to be a Gaussian distribution centered at y = 0. D(x,t) is called a "damping function," because for such a typical choice of P(y), D(x,t) is largest where f is stationary and hence it damps the integrand in Equation (29) away from the stationary points.
An evaluation of D(x,t) poses the same problems as evaluation of I(t). However, for specific choices of P, various approximations for D can be constructed and used as the basis for a more efficient sampling me
thod for the evaluation of I. For example, if P(y) is a narrow Gaussian function with a maximum at y = 0, then for large t an approximation can be constructed by replacing 
(x - y)/(x) by unity and performing a Taylor expansion of the exponent. One then obtains an approximation for D called the "first-order gradient approximation
":
Then
This integral, which is the major contribution to the correct answer, is generally more amenable to evaluation by sampling, in this case by using (x)D0(x,t) as the probabi
lity distribution. A complete method would require the development of procedures to estimate the correction terms that must be added on the right to make this an equality (see Doll et al., 1988, for additional details). The essential idea of these proce
dures is that use of a precise approximate formula for D automatically includes much of the cancellation responsible for decreasing the value of I.
The stationary-phase Monte Carlo methods provide some useful ideas for formulating tractable sampling methods for solving some problems of interest. It would be worthwhile to develop better versions of these methods both for I(t) and for its Fo urier transform.
Alternative Approaches to the Prototype Problem
In many cases of interest, the function I(t) is known to be an analytic function of t at the origin. This knowledge may be of use in developing approximate evaluation techniques based on analytic continuation. Most analytic continuat ion methods are based on a more detailed formulation of the problem of interest that is tied closely to the physical nature of the problem. But even for this simple formulation of the prototype problem, one might ask whether the fact that I(t ) is an analytic function of (complex) t for some region that includes the origin suggests any alternate ways of calculating I(t) accurately.
For example, one might consider the use of Padé approximants (i.e., ratios of polynomials; Baker and Gammel, 1970), which are often used for approximate analytical continuation of analytic functions, especially meromorphic functions. Expanding t he exponential in I(t) gives
where
The first few Taylor series coefficients can then be estimated by sampling on (x) in the usual way, but the difficulty of making accurate estimates grows as n increases. These coefficients
might then be used to fit I(t) with a Padé approximant.
This approach raises a number of mathematical questions.
n make it inappropriate to use Padé approximants to fit the function? Are there better alternatives?
n?
n or how to construct ap
propriate Padé approximants? Other Formulations and Solutions of the Basic Problem
v The nature and dimensionality of the x space
and the specific form of the (x) function are determined by the nature of the quantum mechanical problem to be solved. In genera
l, however, the coordinates of x are Cartesian (or other) coordinates for a mechanical system. Similarly, the form of the function f(x) is determined by the nature of the problem. Various analytic continuation techniques
have been applied to attack specific special cases of the prototype problem. The question then arises as to whether these are isolated solutions of specific problems or special cases of an underlying, not yet discovered, general method for solving proble
ms of this type. This is not the place to delve into the details of specific problems; instead, some of the relevant mathematical principles involved are highlighted. For additional discussion of related problems, the reader is referred to Doll (1984),
Makri (1991), and Wolynes (1987).
Analytic Continuation in Time. In some problems, the goal is to evaluate a function of time t of the form
More generally, there is a function of a complex variable z of the form
that is analytic in a region including the positive real axis and part of the positive imaginary axis, and for which h(t) = lim
->0+ H( + it). The function g(z,x) is an analytic function of all its arguments, and f(x) is real (for real x). In many cases, the formulation in Equ
ation (31) is required to make the problem well defined, because g(it,x) is purely imaginary and the integrals in Equation (30) are not absolutely convergent (for an infinite range of integration).
The calculation of H(z) for real z presents a tractable sampling problem because g(z,x) is real for real z. Then, sampling a set of points distributed with a probability density proportional to g(z,x) and calculating f at those points allow H to be calculated for real z. In fact, as long as z has a positive real part, the problem of evaluating H can be converted to a well-defined sampling problem, b
ecause 
g(z,x) -> -
for large positive and negative values of the integration variables. Then g(z,x) can be used as a probability distribution, and the right side of Equation (31) can be expressed as a ratio of two averages over this distribution:
A variety of analytic continuation methods have been used to solve problems of this type. One strategy is to perform tractable sampling calculations to obtain H(z) for real z and then use Padé approximant techniques to estimate H as z approaches the imaginary axis (Miller et al., 1983; Thirumalai and Berne, 1983; Jacquet and Miller, 1985; Yamashita and Miller, 1985; Makri and Miller, 1989). These techniques motivate the same sorts of mathematical questions as those listed above for a different approach to analytic continuation.
Deformation of the Integration Contours. Another technique has been used in certain situations in which the integrand is an analytic function of all the variables in x. The integration is from -
to along the real line for each of these variables. Deformation of these contours of integration and evaluation of the resulting integrals by sampling methods can be an effective techni
que if the deformed contours come near or pass through stationary-phase points of the integrand (Chang and Miller, 1987; Mak and Chandler, 1991). In one situation in which this was used (Mak and Chandler, 1991), the contour was described by a set of para
meters that were varied during the calculation to maximize some measure of the success of the sampling process being used to calculate the integral.
Most General Formulation. Perhaps the most general formulation of these problems is the following one. Given a set of complex variables z1, ..., zn and a function of these variables g(z1,..., zn) that is an analytic function of all these variables for a domain that includes the real axis for each of them, calculate the following complex function of m complex variables:
The method should work for highly oscillatory functions and should take into account the fact that each evaluation of g and the location of stationary-phase points are expensive. Here z1,. . . , zm are analogous to the t in previous examples, and zm+1, ..., zn are the integration variables.
References
Baker, Jr., G.A., and J.L. Gammel, 1970, The Padé Approximant in Theoretical Physics, Academic Press, New York.
Chang, J., and W.H. Miller, 1987, Monte-Carlo path integration in real-time via complex coordinates, J. Chem. Phys. 87:1648.
Doll, J.D., 1984, Monte Carlo Fourier path integral methods in chemical dynamics, J. Chem. Phys. 81:3536.
Doll, J.D., Thomas L. Beck, and David L. Freeman, 1988, Quantum Monte-Carlo dynamics--The stationary phase Monte-Carlo path integral calculation of finite temperature time correlation-functions, J. Chem. Phys. 89:5753-5763.
Feynman, R.P., and A.R. Hibbs, 1965, Quantum Mechanics and Path Integrals, McGraw- Hill, New York.
Filinov, V.S., 1986, Calculation of the Feynman integrals by means of the Monte-Carlo method, Nucl.Phys. B 271:717-725.
Jacquet, R., and W.H. Miller, 1985, J. Phys. Chem. 89:2139.
Mak, C.H., and D. Chandler, 1991, Coherent-incoherent transition and relaxation in condensed-phase tunneling systems, Phys. Rev. 44:2352-2369.
Makri, N., 1991, Feynman path integration in quantum dynamics, Computer Physics Communications 63:389.
Makri, N., and W.H. Miller, 1987, Monte-Carlo integration with oscillatory integrands: implications for Feynman path integration in real time, Chem. Phys. Lett. 139:10-14.
Makri, N., and W.H. Miller, 1989, Exponential power series expansion for the quantum time evolution operator, J. Chem. Phys. 90:904-911.
Miller, W.H., S.D. Schwartz, and J.D. Tromp, 1983, Quantum mechanical rate constants for bimolecular reactions, J. Chem. Phys. 79:4889-90.
Thirumalai, D., and B.J. Berne, 1983, On the calculation of time correlation functions in quantum systems: path integral techniques, J. Chem. Phys. 79:5029-33.
Wolynes, P.G., 1987, Imaginary time path integral Monte Carlo route to rate coefficients for nonadiabatic barrier crossing, J. Chem. Phys. 87:6559.
Yamashita, K., and W.H. Miller, 1985, "Direct" calculation of quantum mechanical rate constants via path integral methods: Application to the reaction path Hamiltonian, with numerical test for the H + H2 reaction in 3D, J. Chem. Phys. 82:5475-54 84.
NAS Home Page | NAP Home Page | Reading Room | Report Home Page