Scattering using real-time path integrals

# Scattering using real-time path integrals

## Abstract

Background: Path integrals are a powerful tool for solving problems in quantum theory that are not amenable to a treatment by perturbation theory. Most path integral computations require an analytic continuation to imaginary time. While imaginary time treatments of scattering are possible, imaginary time is not a natural framework for treating scattering problems.

Purpose: To test a recently introduced method for performing direct calculations of scattering observables using real-time path integrals.

Methods: The computations are based on a new interpretation of the path integral as the expectation value of a potential functional on a space of continuous paths with respect to a complex probability distribution. The method has the advantage that it can be applied to arbitrary short-range potentials.

Results: The new method is tested by applying it to calculate half-shell sharp-momentum transition matrix elements for one-dimensional potential scattering. The calculations for half shell transition operator matrix elements are in agreement with a numerical solution of the Lippmann-Schwinger equation. The computational method has a straightforward generalization to more complicated systems.

## I Introduction

The purpose of this paper is to explore the possibility of using real-time path integral methods Feynman (1948)Feynman and Hibbs (1965) to calculate scattering observables. The proposed computational method is based on a recent formulation of the path integral Muldowney (2012)Nathanson and Jørgensen (2015)Nathanson (2015) that replaces the integral over paths by the expectation value of the potential functional, , with respect to a complex probability distribution on a space of paths. The space of paths is cylinder sets of continuous paths with fixed starting and end points. An advantage of the method is that it can be used to compute real-time path integrals for arbitrary short-range potentials.

The motivation for this work is that path integrals provide one of the most reliable methods to study strongly interacting quantum systems. In cases of physical interest the path integral is analytically continued to imaginary time, so the integrals can be approximated using Monte-Carlo methods Richtmyer et al. (1947)Metropolis and S. (1949). The most important experimental observables are scattering observables. While clever methods have been developed to extract scattering observables from Euclidean path integrals Luscher and Wolff (1990)Luscher (1986)Bing-Nan et al. (2016)Luu and Savage (2011)Aiello and Polyzou (2016), imaginary time is not the most natural setting for computing scattering observables. While most applications of real-time path integrals to scattering are limited to generating perturbative expansions Faddeev and Slavnov (1980), there are related works on the formulation of scattering theory based on real-time path integrals Rosenfelder (012701) Carron and Rosenfelder (2011) Campbell et al. (1975).

The purpose of this paper is to determine if the formulation of the path integral in Muldowney (2012)Nathanson and Jørgensen (2015)Nathanson (2015) can be used to directly calculate scattering observables in real time. While this initial test is limited to potential scattering in one dimension, the results suggest that the method can be extended to treat few-body scattering in three dimensions. The practical applicability of the method to more complex systems remains an open question. Compared to Euclidean methods, which use Monte Carlo integration, the real-time path integral in this work is approximated by matrix elements of powers of a large matrix, which can be computed efficiently. The matrix is related to a unitary transfer matrix that generates the contribution of multiple paths in parallel, one time step at a time.

In this work the real-time path integral is used to compute sharp-momentum transition matrix elements, since they are the dynamical input to multichannel scattering observables for any number of particles. The method used to construct transition matrix elements in this one-dimensional potential scattering example can be directly applied to construct multi-channel many-body scattering observables in three dimensions.

The theory of complex probabilities was first contemplated by R. Henstock Henstock (1963)Bartle (2001), based on a generalized theory of integration that he co-developed with J. Kurzweil. The Henstock-Kurzweil integral is similar to a Riemann integral, except the usual - definition is replaced by: for any prescribed error, , there is a positive function, (called a gauge), where the intervals and evaluation points in the Riemann sums satisfy . In the Henstock-Kurzweil case the points and intervals are correlated or tagged and the evaluation point is normally taken as one of the endpoints of the interval. Henstock-Kurzweil integral reduces to the Riemann integral when the gauge function, , is a constant.

One advantage of the Henstock-Kurzweil integral is that it can still be approximated to arbitrary accuracy by generalized Riemann sums. The class of Henstock-integrable functions contains the Lebesgue measurable functions as well as functions that are not absolutely integrable. These considerations generalize to infinite dimensionals integrals, which are inductive limits of finite dimensional integrals over cylinder sets that have a well-defined volume and correlated evaluation points and times. Precise definitions for the interested reader can be found in Muldowney (2012)Nathanson and Jørgensen (2015)Nathanson (2015).

Becasue Henstock-Kurzweil integrable functions include Lebesgue integrable functions, Henstock proposed the replacement of a probability theory based on countably additive positive measures by one based on finitely additive generalized Riemann sums. The two formulations of probability coincide for real probabilities, but the Henstock formulation of probability extends to non-positive and complex probabilities. This was further developed by P. Muldowney. He developed a suggestion of Henstock that the Feynman path integral could be understood using this framework. Muldowney Muldowney (2012) proved that the method converges to a local solution of the Schrödinger equation. It was recently shown by P. Jørgensen and one of the authors Nathanson and Jørgensen (2015)Nathanson (2015) that the local solutions could be patched together to construct a global solution of the Schrödinger equation, and a unitary one-parameter time-evolution group. The unitary one-parameter time-evolution group is used to formulate the scattering problem. In addition, the support of the complex probability is on paths that are almost everywhere continuous Nathanson (2015), leading to the interpretation of the path integral as the expectation of a functional of continuous paths with respect to a complex probability distribution.

The approach taken in this paper is a numerical implementation of the formulation of the path integral as the expectation of a potential functional with respect to a complex probability distribution. Gill and Zachary Gill and Zachary (2008) investigated an alternative direct application of the Henstock-Kurzweil integral to path integrals.

While this is just a reinterpretation of the ordinary path integral, it leads to a computational method that is not limited to quadratic potentials. In this paper this method is applied to treat a one-dimensional scattering problem with a Gaussian potential. While this is a simple problem, from a path integral perspective it involves approximating a large-dimensional oscillating integral. The computational method used in this work computes the transition matrix elements in a straightforward, but inefficient way. The efficiency of the method can be significantly improved, however it is not known if it can scaled up to treat problems involving many particles or fields.

This purpose of this work is only to demonstrate that it is possible to perform non-perturbative scattering calculations using this method. An appealing property of this method is that the interpretation of the calculation as the expectation value of a potential functional with respect to a complex probability distribution on a space of continuous paths is not lost in the numerical implementation. In particular the computational algorithm can be deconstructed to find the approximate complex probability for any given cylinder set of paths.

Since the discussion that follows is somewhat detailed, a brief summary that outlines the essential elements of the computational method is given below. The new method starts like the ordinary path integral. It is based on the Trotter product formula Reed and Simon (1980), which is used to express the unitary time-evolution operator as a strong limit

 e−iHt=limN→∞(e−ip22μtNe−iVtN)N (1)

provided that is essentially self-adjoint. It has recently been shown that the limit exists for a class of Hamiltonians with singular potentials that have self-adjoint extensions Nathanson and Jørgensen (2017).

The next step is to insert complete sets of intermediate states that separately diagonalize both the kinetic energy, , and potential energy, . The integrals over the intermediate momenta are Gaussian Fresnel integrals which can be performed analytically. What remains is the large limit of -dimensional integrals over the real line. Each integral is interpreted as an integral over all space at a given time slice. These steps are standard Feynman and Hibbs (1965) and can be found in any textbook that covers path integrals. The new steps are:

• Replace each of the integrals over the real line by a sum of integrals over small intervals. These intervals represent windows that a continuous path passes through at each time slice. The limit that the width of the finite intervals vanish, and the finite endpoints of the semi-infinite intervals become infinite is eventually taken. For computational purposes the intervals should be sufficiently small so the potential is approximately constant on each interval and is approximately zero on the semi-infinite intervals.

• The next step is to perform the integrals over the different products of N intervals, one for each time slice, assuming that the potential is zero. The results are complex quantities that are labeled by space intervals, ; one interval for each time slice. Each sequence of intervals defines a cylinder set. A cylinder set is an ordered set of windows at intermediate times. Every continuous path goes through a unique cylinder set. It is elementary to show that the sum of these integrals over all possible N-fold products of intervals (cylinder sets) is 1. This allows these integrals to be interpreted as complex probabilities that a path is an element of the associated cylinder set.

• If the intervals are sufficiently small so the potential is approximately locally constant on each interval, then the sum of products of the complex probabilities with , where is any sample point in the interval in the cylinder set, converges to the “path integral” defined by the Trotter formula. Because the potential is only needed at a finite set of sample points, this last step makes the method applicable to any-short range potential,

The new interpretation is that the ill-defined path integral is replaced by a well-defined expectation value of the potential functional with respect to a complex probability distribution on cylinder sets of continuous paths. This brief summary is discussed in more detail in what follows.

There are two problems that must be overcome to make this into a computational method. They can be summarized by noting that (1) it is not clear that the complex probabilities can be computed analytically; they involve -dimensional integrals and (2) even if they could be computed either analytically or numerically, there are too many of them. If there are intervals in each of time slices, the number of cylinder sets is , where the final result is obtained in the limit that both and become infinite.

The challenge of this work is to overcome these obstacles without giving up the interpretation of the path integral as the expectation value of a functional of paths with respect to a complex probability distribution of paths.

The virtue of this formulation of the “path integral” is that the non-existent path integral is replaced by the well-defined expectation value of a potential functional of continuous paths with respect to a complex probability distribution on cylinder sets of continuous paths.

While the test problem treated in this paper can be computed more efficiently by directly solving the Lippmann-Schwinger equation, path integrals are a powerful tool for solving problems in quantum field theory that are not limited by perturbation theory. While the methods discussed in this paper do not directly apply to the field theoretic case, they provide a different way of thinking about the problem that could lead to new computational strategies.

The paper is organized as follows. The next section includes a brief discussion of the scattering formalism that will be used in the rest of the paper. This formalism is based on standard time-dependent scattering theory. The method has the advantage that it also applicable to multi-particle scattering in three dimensions. Section three discusses the Feynman path integral formulation of the scattering problem. The fourth section introduces the reinterpretation of the path integral as the expectation value of the potential functional with respect to a complex probability on a space of continuous paths. Section five introduces a factorization of the complex probabilities that makes numerical computations possible. Section six analyzes the test calculation for scattering from a simple short-range potential in one dimension. A summary and concluding remarks appear in section seven.

## Ii Scattering observables using time-dependent methods

The application that will be considered in this work is scattering in one-dimension using a Gaussian potential, . The goal is to calculate sharp-momentum transition matrix elements using a path-integral treatment of time-dependent scattering. The sharp-momentum transition matrix elements are simply related to the scattering cross section. The method can be formally generalized to compute sharp-momentum transition matrix elements for multiparticle scattering in three-dimensions. Because this method is ultimately based on the Trotter product formula, which requires a strong limit, the desired matrix elements need to be extracted using narrow wave packets.

In quantum mechanics the probability for scattering from an initial state to a final state , is

 P=|⟨ψf(t)|ψi(t)⟩|2=|⟨ψf(0)|ψi(0)⟩|2. (2)

The time-independence of the scattering probability follows from the unitarity of the time-evolution operator. Since this probability is independent of time, it can be computed at any convenient common time. In a scattering experiment the initial state, , is simple long before the collision; it looks like a system of free moving particles, . Similarly, the final state, , has a simple form long after the collision; it looks like a system of free moving particles, . The difficulty is that there is no common time when both states have a simple form. Time-dependent scattering theory provides a means to express the initial and final scattering states at a common time in terms of states of asymptotically free particles at a common time. The free particle states are easily computed at any time.

The relation of the initial and final states to the asymptotic system of free particles, described by and , is given by the scattering asymptotic conditions

 limt→∞∥|ψf(t)⟩−|ψf0(t)⟩∥=0limt→−∞∥|ψi(t)⟩−|ψi0(t)⟩∥=0 (3)

which can be written as

 limt→∞∥e−iHt|ψf(0)⟩−e−iH0t|ψf0(0)⟩∥=0limt→−∞∥e−iHt|ψi(0)⟩−e−iH0t|ψi0(0)⟩∥=0 (4)

where is the Hamiltonian of the system. The unitarity of the time-evolution operator, , leads to the equivalent forms:

 limt→∞∥|ψf(0)⟩−eiHte−iH0t|ψf0(0)⟩∥=0limt→−∞∥|ψi(0)⟩−eiHte−iH0t|ψi0(0)⟩∥=0. (5)

These formulas express the initial and final scattering states at the common time, , in terms of the corresponding non-interacting states at the same time.

While these expressions formally involve the limits , if is taken as the time of the collision, the limit is realized at the finite times when are large enough for the interacting particles to be outside of the range of the interaction. In a real experiment the times when the beam and target are prepared and when the reaction products are seen in detectors are finite; the infinite-time limits are a simple way to ensure that is large enough to reach the limiting form. This means that in a realistic calculation the limit can be replaced by a direct evaluation at a sufficiently large finite or . The minimum size of this depends on the range of the interaction, the size of the wave packets, and momentum distribution in the wave packet. When finite times are used in calculations, the minimal size of needs to be determined for each calculation.

Because it appears in equation (5) it is useful to define the operator

 Ω(t):=eiHte−iH0t. (6)

The operator is a unitary operator that, for sufficiently large , transforms the initial or final non-interacting wave packet at time to the initial or final interacting wave packet at . Equation (5) shows

 |ψi(0)⟩≈Ω(−t)|ψi0(0)⟩ (7)

as gets sufficiently large.

The primary quantities of interest are the on-shell transition matrix elements. For the simplest case of scattering by a short-range potential, , the transition matrix elements are related to by

 ⟨pf|T(Ei+i0+)|pi⟩=limt→−∞⟨pf|VΩ(t)|pi⟩ (8)

where

 T(z)=V+V(z−H)−1V (9)

is the transition operator and . Because the limit on the right is a strong limit, it only exists if the sharp-momentum generalized eigenstate, , is replaced by normalizable a wave packet.

A useful approximation is to use a narrow Gaussian wave packet, centered about the initial momentum, , with a delta-function normalization:

 ∫dpψi0(p)=1. (10)

With this choice Lippmann and Schwinger (1950),

 ⟨pf|T(Ei)|pi⟩≈⟨pf|VΩ(t)|pi⟩≈⟨pf|VΩ(t)|ψi0(0)⟩ (11)

in the large , narrow wave packet limit. Gaussian wave packets are minimal uncertainty states, which provide maximal control over both the momentum resolution and spatial width of the free wave packet. For computational purposes, the width of the wave packet in momentum should be chosen so varies slowly on the support of the wave packet.

In this work “path integral” methods are used to compute the right-hand side of equation (11). Because the path integral is formulated in terms of paths in the coordinate representation, the actual quantity that needs to be computed (in one dimension) is the Fourier transform

 ⟨pf|VΩ(t)|ψi0(0)⟩=1√2π∫dxe−ipfx⟨x|VΩ(t)|ψi0(0)⟩ (12)

for a sufficiently large and narrow wave packet. This can be computed by a direct computation of the Fourier transform or by integrating over a narrow final-state wave packet, , centered about with a delta-function normalization:

 ⟨pf|VΩ(t)|ψi0(0)⟩≈⟨ψf0(0)|VΩ(t)|ψi0(0)⟩ (13)

## Iii Scattering using the Feynman path integral

The dynamical quantity needed as input to equation (11) is

 ⟨x|VΩ(t)|ψi0(0)⟩=⟨x|Ve−iHt|ψi0(−t)⟩. (14)

for sufficiently large . The initial state at time zero is a Gaussian approximation to a delta function with the initial momentum, :

 ⟨p|ψi0(0)⟩=12√πΔpe−(p−pi)24(Δp)2. (15)

Here is the quantum mechanical uncertainty in for this wave packet. This wave packet needs to be evolved to using the free time evolution. The result is

 ⟨p|ψi0(−t)⟩=12√πΔpe−(p−pi)24(Δp)2+ip22μt. (16)

These wave packets are needed in the coordinate basis in the path integral. The Fourier transform of (16) can be computed analytically

 ⟨x|ψi0(t)⟩=(2π)−1/2 ⎷11+i2(Δp)2tμe−(Δp)21+4(Δp)4t2μ2(x−pitμ)2ei11+4(Δp)4t2μ2(xpi+2x2(Δp)4tμ−tμp2i2)= (17)
 (2π)−1/2 ⎷11+it2(Δx)2μe−14Δx211+t24Δx4μ2(x−pitμ)2ei11+t24Δx4μ2(xpi+x2t8Δx4μ−tμp2i2). (18)

where , since Gaussian wave functions represent minimal uncertainty states.

Equations (17) and (18) show that the center of this initial wave packet moves with its classical velocity, , so the center of the wave packet is located at , where is the mass and is the mean momentum of the initial wave packet. Ignoring the spreading of the wave packet, it will be in the range of the interaction for a time , where is the range of the potential and is the width of the wave packet. This suggests that the asymptotic time for scattering will be reached for . Because the potential appears in (11) in the expression for the transition matrix elements, only the values of inside the range of the potential are needed to calculate the transition matrix elements.

While the example of a single-particle interacting with a local potential in one dimension is used to test the proposed computational method, the computational method presented below formally generalizes to many-body reactions in three dimensions.

Equation (14) can be expressed exactly as

 (19)

In the “” representation equation (19) has the form

 ∫⟨x|e−i(p22μ+V)t|xi⟩dxi⟨xi|ψi0(−t)⟩=limN→∞∫⟨x|(e−i(p22μ+V)Δt)N|xi⟩dxi⟨xi|ψi0(−t)⟩. (20)

where .

The only contributions to the large limit come from the first-order terms in . This follows from the Trotter product formula Reed and Simon (1980), which gives conditions for the operator version of

 ex=limN→∞(1+x/N)N (21)

to hold when as a strong limit.

Using this property, the limit in (20) can be replaced by

 limN→∞⟨x|e−ip22μΔte−iVΔt⋯e−ip22μΔte−iVΔtN−times% |xi⟩dxi⟨xi|ψi0(−t)⟩. (22)

The following replacements were used in (20) to get (22)

 e−i(p22μ+V)Δt→1−i(p22μ+V)Δt→(1−ip22μΔt)(1−iVΔt)→e−ip22μΔte−iVΔt. (23)

The next step is to insert complete sets of intermediate position and momentum eigenstates so and each become multiplication operators. This leads to the expression

 ⟨x0|e−iHs|ψi0(−t)⟩=limN→∞∫N∏n=1dpndxn2πeipi(xn−1−xn)−ip2n2μΔt−iV(xn)Δt⟨xN|ψi0(−t)⟩. (24)

The integrals are Gaussian Fresnel integrals and can be performed by completing the square in the exponent

 ⟨x0|e−iHt|ψi0(−t)⟩=limN→∞∫N∏n=1dpndxn2πe−iΔt2μ(pn−μΔt(xn−1−xn))2+iμ2Δt(xn−1−xn)2−iV(xn)Δt⟨xN|ψi0(−t)⟩. (25)

The general structure of the resulting momentum integrals is

 ∫∞−∞e−iap2+ibpdp=√πiaeib2/(4a). (26)

These integrals are computed by completing the square in the exponent, shifting the origin, and evaluating the resulting integral by contour integration over a pie shaped path with one edge along the positive real axis and the other making a degree angle between the positive real and negative imaginary axis.

The resulting integral over the momentum variables is

 ⟨x0|e−iHt|ψi0(−t)⟩=
 limN→∞(μ2πiΔt)N/2∫N∏i=ndxneiμ2Δt(xn−1−xn)2−iV(xn)Δt⟨xN|ψi0(−t)⟩. (27)

This is the standard form of the path integral derived by Feynman. The path integral interpretation is obtained by factoring a out of the sum to get

 ⟨x0|e−itH|ψi0(−t)⟩=limN→∞(μ2πiΔt)N/2∫ei∑Nn=1(μ2(xn−1−xnΔt)2−V(xn))ΔtN∏m=1dxm⟨xN|ψi0(−t)⟩. (28)

This looks like an integral over piece-wise linear paths between points in the time slices, , weighted with the imaginary exponential of a finite difference “approximation” of the action:

 A≈N∑n=1(μ2(xn−1−xnΔt)2−V(xn))Δt. (29)

This is in quotes because, due to the integrals, the numerator in the finite difference does not get small as , so the interpretation of as an approximate derivative is not justified.

Irrespective of any concerns about the interpretation, this expression is mathematically well-defined as a limit of finite dimensional integrals. It converges as a result of the Trotter product formula, however it is not very useful for computational purposes because of the large dimensionality of the integrals needed for convergence, except in the case of quadratic interactions, where the integrals can be computed analytically.

## Iv The Muldowney-Nathanson-Jørgensen path integral

To compute the path integral for scattering it is necessary to overcome several obstacles. These include the large dimensionality of the integrals, the need to compute with general short-range interactions, the oscillatory nature of the integrals, and the spreading of the scattering wave packets. The purpose of this work is to investigate some methods that have the potential to overcome these obstacles.

The proposed computations are a consequence of the reformulation of the path integral due to Muldowney Muldowney (2012), Nathanson and Jørgensen Nathanson and Jørgensen (2015)Nathanson (2015). This provides a means for treating a large class of interactions, and eliminates the questionable “finite difference” approximation of the kinetic energy in (28). There still remain oscillations associated with the potential term; but they are only relevant inside of the finite range of the potential.

This method replaces the usual interpretation of the path integral by assigning a “complex probability” to subsets of continuous paths, and computing the expectation value of the random variable with respect to this complex probability distribution. In this expression is a continuous path between and , This differs from the standard interpretation in that the action functional is replaced by a potential functional, and the “measure” is replaced by a complex probability distribution. The random variable differs from 1 only on the portion of the path, , that is in the range of the potential. The potential functional does not suffer from the interpretational difficulties of the action functional in the standard path integral.

This is still a computationally intractable problem. In order to make this computationally tractable, the complex probability is approximately factored into a product of complex probabilities for each “time step”. These one-step “complex probabilities” have the advantage that they can be computed analytically and that the analytic calculation treats the free propagation exactly. The important simplification is that the one-step complex probability can be approximated by a matrix and the multi-step probability is approximately the -fold product of the same matrix. This reduces the calculation of the transition matrix elements to the computation of powers of a matrix applied to a vector. In this case the usual Monte Carlo integration is replaced by matrix multiplication, which can be performed efficiently even for large matrices.

Finally, the use of the operator in (6) means that the quantity being computed in (8) and (11) is a deformation of the initial free wave packet, at the time of collision, to the corresponding interacting packet at the same time. In this case both the free and interacting wave packets remain localized near the origin and the parameter interpolates between the free and interacting localized states. The spreading of the wave packet is only relevant during the time of collision, and even during that time one expects some cancellations. Thus, the relevant parts of the calculation take place in a finite space-time volume.

The fundamental new idea that is the key to the computational strategy, proposed by Muldowney, is to decompose the integral over each in (27) into a sum of integrals over intervals, , at the time slice:

 ∫dxn=Mn∑m=0∫Imndxn. (30)

The intervals are chosen to be disjoint and cover the real line. They are taken to have the general form

 (−∞,x1n)I0n,[x1n,x2n)I1n,⋯,[xM−1,n,xM,n)IM−1,n,[xM,n,∞)IM,n. (31)

Using this decomposition the limit in (28) becomes

 ⟨x0|e−itH|ψi0(−t)⟩=limN→∞(μ2πiΔt)N/2∑m1⋯mNN∏n=1∫Imndxneiμ2Δt(xn−1−xn)2−iV(xn)Δt⟨xN|ψi0(−t)⟩. (32)

The sum is over the -fold Cartesian products of intervals (cylinder sets) for the time slices. Each continuous path from to goes through one interval in each time slice and is thus an element of a unique cylinder set. The sums range over , . Up to this point everything is independent of how the intervals are chosen. For potentials and initial wave packets that are smooth, it is enough to choose the intervals sufficiently small so that the interaction and initial free wave packet are approximately constant on each interval, . Then the contribution from the potential can be factored out of the integral over the interval, and be replaced by evaluating the potential at any point in the interval. In the calculations exhibited in section 6, is taken to be the midpoint of the interval . Because of this, the potential no longer explicitly appears in the integrand, opening up the possibility to treat a large class of potentials. In the limit of small intervals this becomes exact. Thus, the replacement

 e−i∑Nn=1V(xn)Δt⟨xN|ψi0(−t)⟩→e−i∑Nn=1V(ymn)Δt⟨ymN|ψi0(−t)⟩. (33)

in the integrand of equation (32) is expected to be a good approximation on the cylinder set .

Formally the Henstock theory of integration, which is the basis of the probabilistic interpretation, restricts the choice of intervals, evaluation points and time slices needed for convergence. However, for smooth short-ranged potentials and wave packets, the Henstock integrals reduce to ordinary Riemann integrals. Motivated by this, it is assumed that convergence can be achieved using uniformly spaced time slices and intervals of fixed size. Numerical convergence provides an indication of the validity of this assumption.

The replacement (33) leads to the following approximate expression

 ⟨x0|e−iHs|ψi0(−t)⟩≈limN→∞limImn→0(μ2πiΔt)N/2∑m1⋯mNN∏n=1(∫Imndxneiμ2Δt(xn−1−xn)2)e−iV(ymn)Δtψi0(ymN,−t). (34)

The integrals,

 P(x0,Im1⋯ImN):=(μ2πiΔt)N/2N∏n=1∫Imndxneiμ2Δt(xn−1−xn)2, (35)

are interpreted as complex probabilities to arrive at by passing through the sequence of intervals . The probability interpretation follows because the sum of these quantities over all intervals is 1, independent of :

 ∑m1⋯mNP(x0,Im1,⋯,ImN)=1. (36)

This is because, using a simple change of variables, the sum can be transformed to the product of Gaussian-Fresnel integrals that are normalized to unity.

Specifically, is interpreted as the “complex probability” for a path to pass through at time , at time , , at time , and end up at at time . The set of continuous paths that pass through at time , at time define a cylinder set of continuous paths. The right most (initial time) interval only gets contributions from the sample points that are in the support of the initial wave packet. Equation (36) is consistent with the requirement that every continuous path goes through one and only one cylinder set with complex probability 1.

In Nelson (1964) Nelson defines a path integral by analytically continuing the mass in the kinetic energy term. His probability is related to the analytic continuation in the mass of Muldowney’s complex probability.

In this notation the “path integral” becomes

 ⟨x0|e−iHt|ψi0(−t)⟩=limN→∞limVol(Imn)→0∑m1⋯mNP(x0,Im1,⋯,ImN)e−i∑Nn=1V(ymn)Δt⟨ymN|ψi0(−t)⟩. (37)

For the half-infinite intervals, and the upper and lower boundaries increase (resp. decrease) in the limit. Equation (37) is like a Riemann integral with a complex measure, except it is interpreted as the expectation of the random variable with respect to the complex probability distribution . Nathanson and Jørgensen show that the complex probability is concentrated on continuous paths and (37) converges to a global solution of the Schrödinger equation in the limit of finer partitions and more time slices.

This reformulation of Feynman’s original path integral provides a justification to represent time evolution in quantum mechanics as an average over paths with complex probabilities. Scattering wave functions are expectation values of the potential functional . Transition matrix elements require an additional multiplication by the potential followed by the Fourier transform of the resulting quantity. For the case of equally spaced sample points this becomes

 ∫⟨pf|x⟩V(x)dx∑m0,m1⋯mNP(x,Im1,⋯,ImN)e−i∑Nn=1V(ymn)Δt⟨ymN|ψi0(−t)⟩≈
 1√2π∑m0,m1,⋯,mNe−ipfym0δyV(ym0)P(ym0,Im1⋯ImN)e−i∑Nn=1V(ymn)Δt⟨ymN|ψi0(−t)⟩ (38)

where is the width of the interval, and the sum is over the final sample points and the finite intervals, .

## V Factorization

The input to the Muldowney-Nathanson-Jørgensen formulation of the path integral is the complex probabilities that a path will be in a particular cylinder set of paths. Even if the probabilities could be computed analytically, there are cylinder sets in the limit that and become infinite. Summing over all of these configurations is not computationally feasible.

On the other hand, for the case of a single time step, the same approximations that were made for multiple time steps lead to the following

 ⟨xN−1|e−iHΔt|ψi0(−t)⟩≈∑mP(xN−1,ImN)e−iV(ymN)Δt⟨ymN|ψi0(−t)⟩. (39)

This approximates the transformed wave function after one time step. The factorization follows if this wave function is used as the initial state in the transformation to the next time step

 ⟨xN−2|e−iH2Δt|ψi0(−t)⟩≈∑mN−1P(xN−2,Im(N−1))e−iV(ym(N−1))Δt⟨ym(N−1)|e−iHΔt|ψi0(−t)⟩≈
 ∑m(N−1),mNP(xN−2,Im(N−1)e−iV(ym(N−1))ΔtP(y(N−1),ImN)e−iV(ymN)Δt⟨ymN|ψi0(−t)⟩. (40)

Repeating this for all time steps gives the following approximation

 ⟨x0|e−iHt|ψi0(−t)⟩≈∑P(x0,Im1)e−iV(ym1)ΔtN∏n=2P(ym(n−1),Inm)e−iV(ymn)Δt⟨ymN|ψi0(−t)⟩ (41)

where the sum is over the cylinder sets. The factorization leads to the following approximation of the complex probability on the cylinder set :

 P(x0,Im1,⋯,ImN)≈P(x0,Im1)N∏n=2P(ym(n−1),Inm) (42)

With this approximation (11) becomes

 ⟨x0|V|e−iHt|ψi0(−t)⟩≈
 ∑m1⋯mNV(x0)P(x0,Im1)e−iV(ym1)ΔtN∏n=2P(ym(n−1),Inm)e−iV(ymn)Δt⟨ymN|ψi0(−t)⟩. (43)

This representation has a significant advantage over (37) because the matrix elements

 Km,k=P(ym,Ik)e−iV(yk)Δt (44)

where

 P(ym,Ik)=(μ2πiΔt)1/2∫xk+1xkdxeiμ2Δt(ym−x)2=√1iπ∫√μ2Δt(xk+1−ym)√μ2Δt(xk−ym)eiβ2dβ (45)

can be computed analytically, and powers of this matrix can be computed efficiently. The integrals in (45) for finite intervals are Fresnel integrals of the form

 I[a,b]=∫baeix2dx=√π2(Cc(b)−Cc(a))+i√π2(Sc(b)−Sc(a)). (46)

where

 Cc(x)=√2π∫x0cos(t2)dtSc(x)=√2π∫x0sin(t2)dt. (47)

Note that these definitions differ from the definitions of Fresnel integrals given in Abramowitz and Stegun (1964) Abramowitz and Stegun. They are related by

 Cc(√π2x)=CAS(x)Sc(√π2x)=SAS(x). (48)

For the semi-infinite interval with

 I[−∞,b]=∫b−∞eix2dx=12∫∞−∞eix2dx−∫0beix2dx=√π2(1+i2+Cc(b)+iSc(b)). (49)

and for

 I[a,∞]=∫∞aeix2dx=12∫∞−∞eix2dx−∫a0eix2dx=√π2(1+i2−Cc(a)−iSc(a)). (50)

Using these formulas leads to the following expressions for the one-step matrix when is a finite interval:

 Kmk=P(ym,Ik)e−iV(yk)Δt=
 12((Cc(√μ2Δt(xk+1−ym))−Cc(√μ2Δt(xk−ym))+Sc(√μ2Δt(xk+1−ym))−Sc(√μ2Δt(xk−ym)))
 +i(Sc(√μ2Δt(xk+1−ym))−Sc(√μ2Δt(xk−ym))−Cc(√μ2Δt(xk+1−ym))+Cc(√μ2Δt(xk−ym))))e−iV(yk)Δt. (51)

For :

 Kmk=P(ym,Ik0)e−iV(yk)Δt=
 12(1+(Cc(√μ2Δt(x1−ym))+Sc(√μ2Δt(x1−ym)))+i(Sc(√μ2Δt(x1−ym))−Cc(√μ2Δt(x1−ym))))e−iV(y0)Δt (52)

and for :

 KmM=P(ym,IM)e−iV(yM)Δt=
 12(1−(Sc(√m2Δt(xM−ym))+Cc(√m2Δt(xM−ym)))−i(Sc(√m2Δt(xM−ym))−Cc(√m2Δt(xM−ym))))e−iV(yM)Δt. (53)

Combining these approximations the expression for is approximately given by

 ⟨x|V|e−iHt|ψi0(−t)⟩≈∑mkV(x)P(x,Im)e−iV(ym)tKN−1mk⟨yk|ψi0(−t)⟩ (54)

where is the matrix in equations (51-53). This matrix only requires values of the potential. The other elements needed for this computation are the minimal uncertainty wave function at time at the same points and the potential.

While the above formulas are for one-dimensional scattering, the generalization to three dimensions is straightforward. The three-dimensional case involves products of these formulas. Transition matrix elements can be extracted from (54) using (12).

## Vi Computational considerations

To test this method, approximate transition matrix elements are computed for the example of a particle of mass scattering from a repulsive Gaussian potential in one dimension.

From a mathematical perspective the complex probability interpretation assumes that all of the integrals are Henstock-Kurzweil integrals. This means that for given prescribed error, there are restrictions on how to choose the intervals and evaluation points. On the half-infinite intervals the potential is approximately zero and the Henstock-Kurzweil integral is a Fresnel integral that can be computed exactly, while for the finite intervals the Henstock-Kurzweil integrals are Riemann integrals, so convergence can be realized using sufficiently fine, uniformly spaced, space and time grids.

In order keep the analysis as simple as possible (1) time slices are chosen to be equally spaced and (2) the number, , and width of the finite intervals on each time slice are chosen to be identical.

There are a number of constraints that have to be satisfied in order to get a converged approximation. The Trotter product formula gives the exact result in the limit that and vanish. In a computation these terms need to be small. In the first term is an unbounded operator, but the limit is a strong limit, so most of the momentum will be centered near the mean momentum of the initial free particle state. This can also be controlled if . This suggests that both of these dimensionless quantities should be less than 1 for a converged calculation. A second constraint is that the uncertainty in the momentum of the initial state should be less than the momentum. This avoids having slow moving or backward moving parts of the wave packet that will feel the potential for long times. Because of the uncertainty principle, making small makes large - which increases the time that the wave packet feels the potential. Both and oscillate, so the widths on the intervals on each time slice need to be small enough so these quantities are approximately constant on each interval. These limits can be realized by choosing sufficiently small time steps and sufficiently narrow intervals. The cost is higher powers of larger matrices.

Both the range of the potential and spatial width of the initial wave packet determine the active volume that needs to be broken into small intervals. The velocity of the wave packet determines the elapsed time that the initial wave packet interacts with the potential.

The condition that is small requires to be small, where is the number of time steps and is the number of intervals per time step. This means that shortening the time step generally requires including more intervals at each time step.

For the test an initial wave packet with the dimensionless parameters used in the calculations are listed in table I.