%\documentclass[aps,prc,preprint,showpacs,groupedaddress]{revtex4-1}
\documentclass[aps,prc,reprint,noshowpacs,groupedaddress,onecolumn]{revtex4}
\usepackage{bm}% bold math
\usepackage{amsmath}
\usepackage{amsfonts}
\usepackage{amscd}
%\usepackage[caption=false]{subcaption}
\usepackage{graphicx}
%\usepackage{subcaption}
\usepackage{color}
\def\beq{\begin{equation}}
\def\eeq{\end{equation}}
\def\sumint{\sum\mspace{-25mu}\int}
\def\directsumint{\sum\mspace{-25mu}\int\,_{\oplus}}
\begin{document}

\title{Scattering using real-time path integrals}

\author{W.~N.~Polyzou}
\affiliation{Department of Physics and Astronomy, The University of
Iowa, Iowa City, IA 52242}

\author{Ekaterina Nathanson}
\affiliation{School of Science and Technology, Georgia Gwinnett
College, Lawerenceville, GA 30042}

\date{\today}

%\pacs{}

\begin{abstract}

  {\bf 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.  More importantly,
  quantum algorithms for calculating path integrals require real-time
  evolution.

  {\bf Purpose:} To test a recently introduced method for performing direct
  calculations of scattering observables using real-time path
  integrals in order to understand the challenges facing real-time
  path integral calculations of scattering observables.  

  {\bf 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 can in principle be applied
  to arbitrary short-range potentials.

  {\bf Results and Conclusions:} The method is applied to compute
  matrix elements of M{\o}ller wave operators applied to narrow wave
  packets.  These are used to calculate half-shell sharp-momentum
  transition matrix elements for one-dimensional potential scattering.
  The calculations for half-shell transition operator matrix elements,
  while not efficient, converge to the numerical solution of the
  Lippmann-Schwinger equation.  This provides a demonstration that 
  real-time path integrals can be used to compute scattering
  observables.

\end{abstract}

\maketitle

\section{Introduction}

The purpose of this paper is to explore the possibility of using
real-time path integral methods \cite{Feynman_1}\cite{Feynman_2} to
calculate scattering observables.  The proposed computational method
is based on a recent formulation of the path integral
\cite{Muldowney}\cite{Katya_1}\cite{Katya_2} that replaces the
integral over paths by the expectation value of the potential
functional, $F[\gamma] := e^{-i \int V(\gamma(t)) dt}$, with respect
to a complex probability distribution on a space of paths.  The space
of paths is cylinder sets of paths with fixed starting and end points.
The method is not limited to quadratic potentials; 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 \cite{Ulam}\cite{Metropolis}.
The most important experimental observables are scattering
observables.  While clever methods have been developed to extract
scattering observables from Euclidean path integrals
\cite{Luscher}\cite{Luscher2}\cite{Lee}\cite{Luu}\cite{Gordon},
imaginary time is not the most natural setting for computing
scattering observables.

In the future it is anticipated that path integral computations will
be performed by quantum computers.  These computations require real-time
evolution.  While the technology for quantum computation is not
mature, an investigation of real-time methods for the computation
of path integrals is warranted.  These investigations can be used to
study the efficiency of real-time computational methods in
anticipation of future developments in quantum computing.  While most
applications of real-time path integrals to scattering are limited to
generating perturbative expansions \cite{Faddeev}, there are related
works on the formulation of scattering theory based on real-time path
integrals \cite{Rosenfelder1} \cite{Rosenfelder2} \cite{Campbell}.

This paper investigates the application of the formulation of the path
integral in \cite{Muldowney}\cite{Katya_1}\cite{Katya_2} to directly
calculate scattering observables in real time.  While this test is
limited to potential scattering in one dimension, the computation uses
strongly convergent time-dependent methods that are applicable to
general scattering reactions in quantum mechanics and quantum field
theory.  While the calculations presented in this paper are not
applicable to these complex reactions, the long term goal of quantum
information science is be able to treat realistic reactions, so it
make sense to investigate methods that are not limited to the one-dimensional
case.

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 approximates a unitary transfer matrix that generates the
contribution due to multiple paths in parallel, one time step at a time.
%For fields or multiparticle systems this matrix is prohibitively large.

%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 is a path
%integral implementation of standard time-dependent scattering theory
%which can be directly generalized to construct multi-channel many-body
%scattering observables.  

%The theory of complex probabilities was first contemplated by R.
%Henstock \cite{Henstock}\cite{bartle}, based on a generalized theory of
%integration that he co-developed.  The Henstock-Kurzweil integral is
%similar to a Riemann integral, except for a given prescribed error, 
%uniform interval widths and random evaluation points are replaced by
%tagged partitions, consisting of intervals and evaluation points
%that are restricted by a so-called gauge function.  The existence of
%such a gauge function ensures that the sum over the product of the 
%interval widths with function values at the tagged evaluation points
%converges to a number(the integral), as the error is reduced.

The theory of complex probabilities was first contemplated by R.
Henstock \cite{Henstock}\cite{bartle}, 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 $\epsilon$-$\delta$ definition is replaced by: for any
prescribed error, $\epsilon$, there is a positive function,
$\delta (x)$ (called a gauge), where the intervals $I_n$ and
evaluation points $x_n$ in the Riemann sums satisfy
$I_n \in [x_n -\delta(x_n) , x_n+\delta (x_n)]$. 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, $\delta (x)$, 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 dimensional
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 \cite{Muldowney}\cite{Katya_1}\cite{Katya_2}.

%It appears reasonable to expect
%that this limit will be the same as the limit of Riemann sums as the
%size of cylinder sets get arbitrary small.

In classical probability theory the probability is a non-negative
function on a collection of sets.  It has the property that the
probability that a measurement will fall in a set or its complement is
1 and the probability that a measurements will fall in a countable
union of disjoint sets is the sum of the probabilities that the
measurement will fall in each set.  If the sets are generated by
intervals under complements and countable unions, the collection of
sets are Lebesgue measurable sets and probability is a Lebesgue
integrable function.  Because 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 {\it 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 explored by
P. Muldowney.  He developed a suggestion of Henstock that the Feynman
path integral could be understood using this framework.  Muldowney
\cite{Muldowney} proved that the method converges to a local solution
of the Schr\"odinger equation.  It was recently shown by
P. J{\o}rgensen and one of the authors \cite{Katya_1}\cite{Katya_2}
that the local solutions could be patched together to construct a
global solution of the Schr\"odinger 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 \cite{Katya_2}, leading to a rigorous
interpretation of the real-time 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 \cite{Gill} 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, although it is unlikely to be
competitive with standard methods using ordinary computers.

%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 \cite{Simon},
which is used to express the unitary time-evolution operator as a
strong limit
\beq
e^{-iHt} = \lim_{N \to \infty} (e^{-i {p^2\over2\mu}{t \over N} }e^{-iV{t \over N} })^N
\label{i.1}
\eeq
provided that ${{p}^2\over 2 \mu} +V$ 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 \cite{Katya_3}.

The next step is to insert complete sets of intermediate states that
separately diagonalize both the kinetic energy, $p^2/2\mu$, and
potential energy, $V$.  The integrals over the intermediate momenta
are Gaussian Fresnel integrals which can be performed analytically.
What remains is the large $N$ limit of $N$-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 \cite{Feynman_2}
and can be
found in any textbook that covers path integrals.  The new steps are:
\begin{itemize}

\item[1] Replace each of the $N$ integrals over the real line by a sum
  of integrals over small intervals.  These intervals represent
  windows that a 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.

\item[2] The next step is to perform the integrals over the different
  products of N intervals, one for each time slice, {\it assuming that
    the potential is zero}.  The results are complex quantities that
  are labeled by $N$ space intervals, $\{I_{n_1}, \dots I_{n_N}\}$;
  one interval for each time slice.  Each sequence of intervals
  defines a cylinder set.  A cylinder set is an ordered set of $N$
  windows at $N$ intermediate times.  Every path goes
  through a unique cylinder set.  Thus each cylinder set represents an
  equivalence class of paths.  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.

\item[3] 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
  $e^{-i \sum V(x_i) \Delta t}$, where $x_i$ is any sample point in
  the $i^{th}$ interval in the cylinder set, converges to the ``path
  integral'' defined by the Trotter product 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,

\end{itemize}

The new interpretation is that the ill-defined real-time 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 $N$-dimensional integrals and (2) even if they could be
computed either analytically or numerically, there are too many of
them.  If there are $M$ intervals in each of $N$ time slices, the number of
cylinder sets is $M^N$, where the final result is obtained in the
limit that both $M$ and $N$ 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,
%hey provide a different way of thinking about the problem.
%ore importantly,
they provide a laboratory for testing computational strategies that
could be implemented in the future on a quantum computer.

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 time-dependent method has the advantage that it
is in principle applicable to multi-particle scattering in three
dimensions, scattering with long-range potentials \cite{Dollard}
\cite{dollard:1978we} and scattering in quantum field theory
\cite{Haag:1958vt}\cite{Ruelle:1962}.  Using a method that can be
applied to a large class of scattering reactions is preferred for
developing algorithms that might be implemented on a quantum computer
in the future.  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.
 
\section{Scattering observables using time-dependent methods}

The application that will be considered in this work is scattering in
one-dimension using an attractive Gaussian potential, $V(x) = -v_0e^{-(x/r_0)^2}$.
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 multi-particle
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 $\vert \psi_i \rangle$ to a final state $\vert \psi_f \rangle$,
is
\beq
P = \vert \langle \psi_f (t) \vert \psi_i (t) \rangle \vert^2 =
\vert \langle \psi_f (0) \vert \psi_i (0) \rangle \vert^2 .
\label{s.1}
\eeq
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, $\vert \psi_i (t)
\rangle$, is simple long before the collision; it looks like a system
of free moving particles, $\vert \psi_{i0} (t) \rangle$.  Similarly,
the final state, $\vert \psi_f (t) \rangle$, has a simple form long
after the collision; it looks like a system of free moving particles,
$\vert \psi_{f0} (t) \rangle$.  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 $\vert \psi_{f0}(t) \rangle$ and
$\vert \psi_{i0}(t) \rangle$, is given by the scattering asymptotic
conditions:
\beq
\lim_{t \to \infty} \Vert \vert \psi_f(t) \rangle - \vert \psi_{f0}(t)
\rangle \Vert =0 
\qquad
\lim_{t \to -\infty} \Vert \vert \psi_i(t) \rangle - \vert \psi_{i0}(t)
\rangle \Vert =0 .
\label{s.2}
\eeq
In the general multichannel or field theoretic case these expressions
have the form
\beq
\lim_{t \to \infty} \Vert e^{-iHt} \vert \psi_f(0) \rangle - 
\Phi_\alpha e^{-iH_\alpha t}\vert \psi_{f\alpha }(0)
\rangle \Vert =0 
\qquad
\lim_{t \to -\infty} \Vert e^{-iH t} \vert \psi_i(0) \rangle - 
\Phi_\beta e^{-iH_\beta t}\vert \psi_{i\beta}(0)
\rangle \Vert =0 
\label{s.3a}
\eeq
where $H$ is the Hamiltonian of the system, $\alpha$ and $\beta$
represent the final and initial scattering channels, $\Phi_{\alpha}$
and $\Phi_{\beta}$ are injection operators that map spaces of
asymptotically free point particles into the physical Hilbert space by
including the internal structure of the asymptotic particles,
$H_\alpha$ and $H_\beta$ are the energies of the free asymptotic
particles and $\vert \psi_{i\beta}(0)\rangle$ and
$\vert \psi_{i\beta}(0)\rangle $ are products of sharp momentum wave
packets in the total momentum and spin of each asymptotic fragment.
In quantum field theories the $\Phi_{\alpha/ \beta}$ are products of
the quasilocal functions of fields that create one-body states out of
the vacuum in Haag-Ruelle theory \cite{Haag:1958vt}\cite{Ruelle:1962}.
There are also $\Phi_{\alpha/ \beta}$'s for Coulomb and other
long-range interactions \cite{Dollard} \cite{dollard:1978we}.

The unitarity of $e^{iHt}$ can be
used to write this as
\beq
\lim_{t \to \infty} \Vert\vert \psi_f(0) \rangle - 
e^{iHt}\Phi_\alpha e^{-iH_\alpha t}\vert \psi_{f\alpha }(0)
\rangle \Vert =0 
\qquad
\lim_{t \to -\infty} \Vert \vert \psi_i(0) \rangle - 
e^{iHt} \Phi_\beta e^{-iH_\beta t}\vert \psi_{i\beta}(0)
\rangle \Vert =0 
\label{s.3b}
\eeq
The important point is that in all cases the limits are strong and the
dynamics enters via $e^{-iHt}$ which can in principle be computed by a
path integral.  The only non-trivial input is the channel injection
operators $\Phi_\alpha$, which are products of bound states or
products of interpolating fields that create point-mass eigenstates
out of the vacuum.

For two-body scattering (\ref{s.3b})  
has the form
%\beq
%\lim_{t \to \infty} \Vert e^{-iHt} \vert \psi_f(0) \rangle - 
%e^{-iH_0t}\vert \psi_{f0}(0)
%\rangle \Vert =0 
%\qquad
%\lim_{t \to -\infty} \Vert e^{-iH t} \vert \psi_i(0) \rangle - 
%e^{-iH_0 t}\vert \psi_{i0}(0)
%\rangle \Vert =0 
%\label{s.3}
%\eeq
%where $H=H_0+V$ is the Hamiltonian of the system.
%The unitarity of the time-evolution operator, $e^{-iHt}$, 
%leads to the equivalent forms:
\beq
\lim_{t \to \infty} \Vert \vert \psi_f(0) \rangle - 
e^{iHt} e^{-iH_0t} \vert \psi_{f0}(0)
\rangle \Vert =0 
\qquad
\lim_{t \to -\infty} \Vert \vert \psi_i(0) \rangle - 
e^{iHt} e^{-iH_0t}  \vert \psi_{i0}(0)
\rangle \Vert =0 
\label{s.4}
\eeq
where $H=H_0+V$ is the Hamiltonian of the two-body system.
These formulas express the initial and final scattering states at the
common time, $t=0$, in terms of the corresponding non-interacting states at
the same time.

While these expressions formally involve the limits $t\to \pm \infty$,
if $t=0$ is taken as the time of the collision, the limit is realized
at the finite times when $\pm t$ 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 $t$ 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 $t$ or $-t$.  The minimum size of this $t$ 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 $t$ needs to be determined for each
calculation.

Because it appears in equation (\ref{s.4}) it is useful to define the
operator
\beq
\Omega (t) := e^{iHt} e^{-iH_0t} .
\label{s.5}
\eeq
The operator $\Omega (t)$ is a unitary operator that, for sufficiently
large $\pm t$, transforms the initial or final {\it non-interacting} wave
packet at time $t=0$ to the initial or final interacting wave packet at
$t=0$.  Equation (\ref{s.4}) shows
\beq
\vert \psi_i(0) \rangle \approx \Omega (-t) \vert \psi_{i0}(0) \rangle 
\label{s.5b}
\eeq
as $t$ gets sufficiently large.

The quantities of physical interest are the on-shell
transition matrix elements.   For the simplest case of
scattering by a short-range potential, $V(x)$, the transition matrix elements
are related to $\Omega (t)$ by 
\beq
\langle \mathbf{p}_f \vert T(E_i+i0^+) \vert \mathbf{p}_i \rangle
=
\lim_{t \to - \infty}
\langle \mathbf{p}_f \vert V \Omega (t) \vert \mathbf{p}_i \rangle
\label{s.6}
\eeq
where
\beq
T(z) = V +V (z-H)^{-1}V
\label{s.7}
\eeq
is the transition operator and
$E_i = {\mathbf{p}_i^2 \over 2m}= {\mathbf{p}_f^2 \over 2m}$.
Because the limit on the right is a strong limit, it only exists if 
the sharp-momentum generalized eigenstate, $\vert \mathbf{p}_i \rangle$, is
replaced by normalizable a wave packet.

A useful approximation is to use a narrow Gaussian wave packet, centered 
about the initial momentum, $\mathbf{p}_i$, with a delta-function normalization:
\beq
\int d\mathbf{p} \langle \mathbf{p} \vert \psi_{i0}\rangle =1 .
\label{s.8}
\eeq
With this choice \cite{Lippmann}, 
\beq
\langle \mathbf{p}_f \vert T(E_i) \vert \mathbf{p}_i \rangle
 \approx
\langle \mathbf{p}_f \vert V \Omega (t) \vert \mathbf{p}_i \rangle
\approx
\langle \mathbf{p}_f \vert V \Omega (t) \vert \psi_{i0}(0) \rangle
\label{s.9}
\eeq
in the large $t$, 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
$\langle \mathbf{p}_f \vert T(E_i) \vert \mathbf{p}_i \rangle$ 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 (\ref{s.9}).  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
\beq
\langle {p}_f \vert V \Omega (t) \vert \psi_{i0}(0) \rangle =
{1 \over \sqrt{2 \pi}}
\int dx e^{-ip_fx} \langle {x} \vert V \Omega (t) \vert \psi_{i0}(0) \rangle
\label{s.10}
\eeq
for a sufficiently large $t$ 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,
$\langle \psi_{f0}(0) \vert$, centered about ${p}_f$ with a
delta-function normalization:
\beq
\langle {p}_f \vert V \Omega (t) \vert \psi_{i0}(0) \rangle \approx 
\langle \psi_{f0}(0) \vert V \Omega (t) \vert \psi_{i0}(0) \rangle .
\label{s.11}
\eeq

The transition matrix elements in (\ref{s.9}) are directly related to 
scattering matrix elements,  cross sections,
phase shifts and, in this one-dimensional case,
to reflection and transmission coefficients.  In the 
one-dimensional case scattering matrix elements
are related to transition matrix elements by
\beq
\langle p \vert S \vert p' \rangle =
\delta (p-p') - 2 \pi i \delta (E-E') 
\langle p \vert T(E+i0^+) \vert p' \rangle  
\label{s.12}
\eeq
where $p=\pm p'$.
The phase shifts are defined by
\beq
\langle E(\pm) \vert S \vert E'\rangle = \delta (E-E') e^{2 i \delta_{\pm} (E)}
= \sqrt{{m \over p}} \langle p \vert S \vert p' \rangle  \sqrt{{m \over p'}}
.
\label{s.13}
\eeq
There are two values 
depending on whether the final state has momentum parallel or anti-parallel
to the initial state.

Using $\vert p \rangle = \vert E \rangle \sqrt{{p\over m}}$ in
(\ref{s.12}) gives
\beq
e^{2 i \delta_+ (E)}= 1 - 2 \pi i {m \over \vert p \vert }  
\langle p \vert T(E+i0^+) \vert p \rangle
\label{s.14}
\eeq
\beq
e^{2 i \delta_- (E)}= - 2 \pi i {m \over \vert p \vert }  
\langle -p \vert T(E+i0^+) \vert p \rangle .
\label{s.15}
\eeq
This leads to the following expressions for the phase shifts in terms of 
the real and imaginary parts of the on-shell transition matrix elements:
\beq
\delta_{+} (E) = {1 \over 2} \tan^{-1} ({-2 \pi m Re(\langle  p \vert
T(E+i0^+) \vert p
\rangle) 
\over \vert p\vert  + 2 \pi m 
Im (\langle  p \vert T(E+i0^+) \vert p \rangle)})
\label{s.16}
\eeq
\beq
\delta_{-} (E) = {1 \over 2} \tan^{-1} ({- Re(\langle - p \vert T (E+i0^+)
\vert p
\rangle) 
\over  
Im (\langle - p \vert T(E+i0^+) \vert p \rangle)}) .
\label{s.16}
\eeq
The on-shell transition matrix elements are related to the transmission
and reflection coefficients by
\beq
T=1-{2 \pi i p \over m} \langle p \vert T(E+i0^+) \vert p \rangle
\qquad
R= -{2 \pi i p \over m} \langle -p \vert T(E+i0^+) \vert p \rangle .
\label{s.16}
\eeq
The expression $\vert R\vert^2 + \vert T \vert^2=1$ expresses the unitarity
of the scattering operator
\beq
S_{++}^* S_{++} + S_{+-}^* S_{-+} =1
\label{s.16}
\eeq


\section{ Scattering using the Feynman path integral} 

The dynamical quantity needed as input to equation (\ref{s.9}) is 
\beq
\langle x \vert V \Omega (t) \vert \psi_{i0}(0) \rangle =
\langle x \vert V
e^{-iHt} \vert
\psi_{i0} (-t) \rangle .
\label{fs.1}
\eeq
for sufficiently large $t$.  The initial state at time zero is a 
Gaussian approximation to a delta function with the initial momentum,
$p_i$:
\beq
\langle p \vert \psi_{i0}  (0) \rangle =
{1 \over 2 \sqrt{\pi} \Delta p} e^{- {(p-p_i)^2 \over 4 (\Delta p)^2}} .
\label{fs.2}
\eeq
Here $\Delta p$ is the quantum mechanical uncertainty in $p$ for this
wave packet.  This wave packet needs to be evolved to $-t$ using the 
free time evolution.  The result is
\beq
\langle p \vert \psi_{i0} (-t) \rangle =
{1 \over 2 \sqrt{\pi} \Delta p} e^{ 
- {(p-p_i)^2 \over 4 (\Delta p)^2} + i {p^2 \over 2 \mu}t}.
\label{fs.3}
\eeq
These wave packets are needed in the coordinate basis in the path 
integral.  The Fourier transform of (\ref{fs.3}) can be computed 
analytically and expressed in terms of the momentum or coordinate uncertainty
of the wave packet:
\beq
\langle x \vert \psi_{i0} (t) \rangle =
(2 \pi)^{-1/2}
\sqrt{{1  \over 1  + i {2 (\Delta p)^2 t \over  \mu}}}
e^{-
{(\Delta p)^2 \over 1  + {4 (\Delta p)^4 t^2 \over  \mu^2}}
(x-{p_it\over \mu})^2}
e^{i{1 \over 1  + {4 (\Delta p)^4 t^2 \over  \mu^2}}
(x p_i   
+ 2 x^2 (\Delta p)^4 {t \over  \mu}
- { t \over  \mu}{p_i^2 \over 2 })} =
\label{fs.4}
\eeq
\beq
(2 \pi)^{-1/2}
\sqrt{{1  \over 1  + i { t \over 2 (\Delta x)^2  \mu}}}
e^{-{1 \over 4 \Delta x^2}
{1 \over 1  + { t^2 \over 4 \Delta x^4   \mu^2}}
(x-{p_it\over \mu})^2}
e^{i{1 \over 1  + {t^2 \over 4 \Delta x^4  \mu^2}}
(x p_i   
+ { x^2 t \over  8 \Delta x^4 \mu}
- { t \over  \mu}{p_i^2 \over 2 })} .
\label{fs.5}
\eeq
where $\Delta p \Delta x = {1 \over 2}$, since Gaussian wave functions 
represent minimal uncertainty states.  Equations (\ref{fs.4})
and (\ref{fs.5}) have the
form a Gaussian with a moving center multiplied by an oscillating
function.
These equations show that the center of this
initial wave packet moves with its classical velocity, $v=p_i/\mu$, so
the center of the wave packet is located at
${x}(t) = { {p}_i \over \mu}t$, where $\mu$ is the mass and $p_i$ 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 $t \approx {(R + \Delta x) \mu \over \vert \mathbf{p}_i \vert}$,
where $R$ is the range of the potential and $\Delta x$ is the width of
the wave packet.  This suggests that the asymptotic time for
scattering will be reached for
$t > {(R+ \Delta x)\mu \over \vert {p}_i\vert} $.  Because the
potential appears in (\ref{s.9}) in the expression for the transition
matrix elements, only the values of
$\langle {x} \vert e^{-iHt} \vert \psi_i(-t) \rangle$ 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 (\ref{fs.1}) can be expressed exactly as
\beq
\langle x \vert
V e^{-iHt} \vert 
\psi_{i0} (-t) \rangle =
\lim_{N \to \infty} 
\langle x \vert V
\left (e^{-iH{t/N}}\right )^N  \vert
\psi_{i0} (-t) \rangle . 
\label{fs.6}
\eeq
In the ``$x$'' representation
the wave function in (\ref{fs.6}) has the form
\beq
\int \langle x \vert e^{-i ({p^2 \over 2\mu} + V) t} \vert x_i \rangle
dx_i \langle x_i \vert \psi_{i0} (-t)\rangle  = 
\lim_{N\to \infty}
\int \langle x \vert (e^{-i ({p^2 \over 2\mu} + V)\Delta t})^N \vert x_i
\rangle dx_i \langle x_i \vert \psi_{i0} (-t)\rangle . 
\label{fs.7}
\eeq
where $\Delta t := t/N$. 

The only contributions to the large $N$ limit come from the
first-order terms in $\Delta t$.  This follows from the Trotter
product formula \cite{Simon}, which gives conditions for the operator
version of
\beq
e^x = \lim_{N \to \infty}(1+{x/N})^N 
\label{fs.8}
\eeq
to hold when $N\to \infty$ as a strong limit.

Using this property, the limit in (\ref{fs.7}) can be replaced by
\beq
\lim_{N\to \infty}
\langle x \vert 
\underbrace{e^{-i {p^2 \over 2\mu}\Delta t}
e^{-i V \Delta t}
\cdots 
e^{-i {p^2 \over 2\mu}\Delta t}
e^{-i V \Delta t}}_{N-\mbox{times}}
\vert x_i 
\rangle   dx_i \langle x_i \vert \psi_{i0} (-t)\rangle .
\label{fs.9}
\eeq
The following replacements were used in (\ref{fs.7}) 
to get (\ref{fs.9})
\beq
e^{-i ({p^2 \over 2\mu} + V)\Delta t} \to
1 -i ({p^2 \over 2\mu} + V)\Delta t  \to
(1 -i {p^2 \over 2\mu}\Delta t )(1-i V \Delta t )\to
e^{-i {p^2 \over 2\mu}\Delta t}
e^{-i V\Delta t} .
\label{fs.10}
\eeq
This representation has the advantage that unitary time evolution
is expressed as the limit of products of unitary operators.

The next step is to insert complete sets of intermediate position 
and momentum eigenstates
so ${p^2 \over 2\mu}$ and $V$ each become multiplication operators.  This 
leads to the expression 
\beq
\langle x_0 \vert
e^{-iHt} \vert
\psi_{i0} (-t) \rangle =
\lim_{N\to \infty}
\int \prod_{n=1}^N {dp_n dx_n \over 2\pi} 
e^{ i p_i (x_{n-1}-x_n) - i {p_n^2 \over 2\mu}\Delta t -i 
V(x_n) \Delta t} \langle x_N \vert \psi_{i0} (-t)\rangle .
\label{fs.11}
\eeq

The $p_n$ integrals are Gaussian Fresnel integrals and can be performed by
completing the square in the exponent
\beq
\langle x_0 \vert
e^{-iHt} \vert
\psi_{i0} (-t) \rangle =
\lim_{N\to \infty}
\int \prod_{n=1}^N {dp_n dx_n \over 2 \pi}
e^{
-i {\Delta t \over 2\mu} (p_n - {\mu \over \Delta t}(x_{n-1}-x_n))^2
+i {\mu \over 2\Delta t} (x_{n-1}-x_n)^2
-i 
V(x_n) \Delta t }
\langle x_N \vert \psi_{i0} (-t)\rangle. 
\label{fs.12}
\eeq
The general structure of the resulting momentum integrals is 
\beq
\int_{-\infty}^{\infty} e^{-i a p^2 + i bp } dp =
\sqrt{\pi \over i a} e^{i b^2/(4a) }.
\label{fs.13}
\eeq
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 $p$ axis and the other making a $45$ degree angle between the positive 
real and negative imaginary $p$ axis.

The resulting integral over the $N$ momentum variables is 
\[
\langle x_0 \vert e^{-iHt}  \vert \psi_{i0} (-t) \rangle =  
\]
\beq
\lim_{N\to \infty}
({\mu \over 2 \pi i \Delta t})^{N/2}
\int \prod_{i=n}^N dx_n
e^{
i {\mu \over 2\Delta t} (x_{n-1}-x_n)^2
-i V(x_n) \Delta t }
\langle x_N \vert \psi_{i0} (-t)\rangle.  
\label{fs.14}
\eeq
This is the standard form of the path integral derived by Feynman. 
The path integral interpretation is obtained by factoring 
a $\Delta t$ out of the sum in the exponent to get
\beq
\langle x_0 \vert e^{-itH}  \vert \psi_{i0} (-t) \rangle =  
\lim_{N\to \infty}
({\mu \over 2 \pi i \Delta t})^{N/2}
\int 
e^{
i \sum_{n=1}^N \left ({\mu \over 2 } ({x_{n-1}-x_n\over \Delta t})^2
- V(x_n)\right ) \Delta t }\prod_{m=1}^N dx_m
\langle x_N \vert \psi_{i0} (-t)\rangle . 
\label{fs.15}
\eeq
This looks like an integral over piece-wise linear paths between
points in the $N$ time slices, $(x_N \to x_{N-1} \to \cdots \to x_0)$,
weighted with the imaginary exponential of a finite difference 
``approximation'' of the action:
\beq
A \approx  \sum_{n=1}^N \left ({\mu \over 2 } ({x_{n-1}-x_n\over \Delta t})^2
- V(x_n)\right ) \Delta t .
\label{fs.16}
\eeq
This is in quotes because, due to the integrals, the numerator in the
finite difference does not get small as $\Delta t \to 0$, so the
interpretation of ${x_{n-1}-x_n\over \Delta t}$
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.

\section{The Muldowney-Nathanson-J{\o}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 \cite{Muldowney}, Nathanson and 
J{\o}rgensen \cite{Katya_1}\cite{Katya_2}.  This
provides a means for treating a large class of interactions, and
eliminates the questionable ``finite difference'' approximation of the
kinetic energy in (\ref{fs.15}).  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 
paths, and computing
the expectation value of the random variable $F[\gamma] = e^{-i \int
V(\gamma (t) ) dt}$ with respect to this complex probability
distribution.  In this expression $\gamma(t)$ is a path
between $x_N$ and $x_0$ , 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 $e^{-i \int V(\gamma
  (t) ) dt}$ differs from 1 only on the portion of the path, $\gamma
(t)$, 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 $N$-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. 

Finally, the use of the operator $\Omega (t)$ in (\ref{s.5}) means
that the quantity being computed in (\ref{s.6}) and (\ref{s.9}) 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 $t$ 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
some cancellations are expected.  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
$x_n$ in (\ref{fs.14}) into a sum of integrals over $M_n+1$ intervals,
$I_{mn}$,
at the $n^{th}$ time slice:
\beq
\int dx_n = \sum_{m=0}^{M_n} \int_{I_{mn}} dx_n .
\label{njm.1}
\eeq
The intervals $I_{mn}$ are chosen to be disjoint and cover the real line.  They
are taken to have the  general form
\beq
\underbrace{(-\infty, x_{1n})}_{I_{0n}}, 
\underbrace{[x_{1n},x_{2n})}_{I_{1n}} , \cdots , 
\underbrace{[x_{M-1,n},x_{M,n})}_{I_{M-1,n}},
\underbrace{[x_{M,n},\infty)}_{I_M,n} .
\label{njm.2}
\eeq
Using this decomposition the limit in (\ref{fs.15}) becomes
\beq
\langle x_0 \vert e^{-itH}  \vert \psi_{i0} (-t) \rangle =  
\lim_{N\to \infty} 
({\mu \over 2 \pi i \Delta t})^{N/2}
\sum_{m_1 \cdots m_N}
\prod_{n=1}^N  \int_{I_{m_n}} dx_n
e^{
i {\mu \over 2\Delta t} (x_{n-1}-x_n)^2
-i V(x_n) \Delta t }
\langle x_N \vert \psi_{i0} (-t)\rangle.  
\label{njm.3}
\eeq
The sum is over the $(M_1+1)\times \cdots \times (M_N+1)$ $N$-fold
Cartesian products of intervals (cylinder sets) for the $N$ time
slices.  Each continuous path from $x_N$ to $x_0$ goes through one
interval in each time slice and is thus an element of a unique
cylinder set.  The $m_n$ sums range over $0\leq m_n \leq M_n$, $1 \leq
n \leq N$.  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, $I_{mn}$.
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 $y_{mn} \in I_{mn}$ in the interval.  In the calculations exhibited in
section 6, $y_{mn}$ is taken to be the midpoint of the interval $I_{mn}$.
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
\beq
e^{-i \sum_{n=1}^N V(x_{n})\Delta t}\langle x_N \vert \psi_{i0} (-t)\rangle \to
e^{-i \sum_{n=1}^N V(y_{mn})\Delta t}\langle y_{mN} \vert \psi_{i0} (-t)\rangle .
\label{njm.4}
\eeq
in the integrand of equation (\ref{njm.3}) is expected to be 
a good approximation on
the cylinder set $I_{m_0}\times I_{m_1} \times \cdots \times I_{m_N}$.

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.
  
%Nathanson and
%J{\o}rgensen 
%showed that for suitable potentials the resulting
%expression is a global solution of the Schr\"odinger equation.

The replacement (\ref{njm.4}) leads to the following approximate expression
\beq
\langle x_0 \vert e^{-iHs} \vert \psi_{i0}(-t) \rangle \approx
\lim_{N\to \infty} \lim_{I_{mn}\to 0} 
({ \mu \over 2 \pi i \Delta t})^{N/2}
\sum_{m1 \cdots mN} \prod_{n=1}^N 
(
\int_{I_{mn}}
dx_n
e^{
i {\mu \over 2\Delta t} (x_{n-1}-x_n)^2}
)
e^{-i V(y_{mn}) \Delta t }
\psi_{i0} (y_{mN},-t).  
\label{njm.5}
\eeq

The integrals, 
\beq
%\boxed{
P(x_0,I_{m1} \cdots I_{mN}) := 
({\mu \over 2 \pi i \Delta t})^{N/2}
\prod_{n=1}^N  \int_{I_{mn}} dx_n
e^{
i {\mu \over 2\Delta t} (x_{n-1}-x_{n})^2},  
%}
\label{njm.6}
\eeq
are interpreted as complex probabilities to arrive at $x_0$ by passing
through the sequence of intervals $I_{mN} \cdots I_{m1}$. 
The probability interpretation follows because the sum of these quantities
over all intervals is 1, independent of $x_0$:
\beq
\sum_{{m1} \cdots {mN}} 
P(x_0,I_{m1}, \cdots ,I_{mN}) = 1 .
\label{njm.7}
\eeq
This is because, using a simple change of variables, 
the sum can be transformed to the product of $N$ Gaussian-Fresnel 
integrals that are normalized to unity.   

Specifically, $P(x_0,I_{m1}, \cdots ,I_{mN})$ is interpreted as the
``complex probability'' for a path to pass through $I_{mN}$ at time
$t_N$, $I_{mN-1}$ at time $t_{N-1}$, $\cdots$, $I_{m1}$ at time
$t_{m1}$, and end up at $x_0$ at time $t$.  The set of
paths that pass through $I_{mN}$ at time $t_N$, $ \cdots , I_{m1}$ at
time $t_1$ define a cylinder set of paths.  The right most
(initial time) interval only gets contributions from the sample points
$y_{mN}$ that are in the support of the initial wave packet.  Equation
(\ref{njm.7}) is consistent with the requirement that every
path goes through one and only one cylinder set with complex
probability 1.

In \cite{Nelson} 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
\beq
%\boxed{
\langle x_0 \vert e^{-iHt} \vert \psi_{i0}(-t) \rangle = 
\lim_{N\to \infty} \lim_{Vol(I_{mn})\to 0}
\sum_{m1 \cdots mN}  P(x_0, I_{m1}, \cdots ,I_{mN})
e^{
-i \sum_{n=1}^{N} V(y_{mn}) \Delta t }
\langle y_{mN} \vert \psi_{i0} (-t)\rangle . 
%}
\label{njm.8}
\eeq
For the half-infinite intervals, $I_{0n}$ and $I_{Mn}$ the upper and
lower boundaries increase (resp. decrease) in the limit.  
%These 
%^are not expected to contribute in this limit because of the 
%square integrability of the initial state and unitarity of $\Omega (t)$. 
Equation
(\ref{njm.8}) is like a Riemann integral with a complex volume element,
except it is interpreted as the expectation of the
random variable
$ e^{ -i \sum_{n=1}^{N} V(y_{mn}) \Delta t }
\langle y_{mN} \vert \psi_{i0} (-t)\rangle$ with
respect to the complex probability distribution
$P(x_0,I_{m1}, \cdots ,I_{mN})$.  Nathanson and J{\o}rgensen show that
the complex probability $P(x_0,I_{m1} \cdots I_{mN})$ is concentrated
on continuous paths and (\ref{njm.8}) converges to a global solution of the
Schr\"odinger 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 $e^{-i
%\int  V(\gamma(t)) dt}$.
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
\[
\int \langle p_f \vert x \rangle V(x) dx 
\sum_{m0,m1 \cdots mN} P(x,I_{m1}, \cdots , I_{mN})
e^{
-i \sum_{n=1}^{N} V(y_{mn}) \Delta t }
\langle y_{mN} \vert \psi_{i0} (-t)\rangle 
\approx
\]
\beq
{1 \over \sqrt{2 \pi}} \sum_{m0, m1, \cdots , mN} e^{-ip_f y_{m0}} \delta y
V(y_{m0}) P(y_{m0},I_{m_1} \cdots I_{m_N}) e^{ -i \sum_{n=1}^{N} V(y_{m_n}) \Delta t }
\langle y_{m_N} \vert \psi_{i0} (-t)\rangle 
\label{njm.9}
\eeq
where $\delta y$ is the width of the $I_{m_N}$ interval, and the 
sum is over the final sample points and the finite intervals, $I_{1_N} 
\cdots I_{M-1,N}$.

%which formally leads to 
%\[
%\langle \psi_{f0}(0) \vert V e^{-iHs} \vert \psi_{i0}(-s) \rangle = 
%\]
%\beq
%\lim_{N\to \infty} \lim_{Vol(I_{mn})\to 0}
%\sum_{m0,m1 \cdots mN}
%\chi^*( I_{m0})
%P(y_{m0}, I_{m1}, \cdots ,I_{mN})
%e^{
%-i \sum_{i=1}^{N} V(y_{mn}) \Delta t }
%\psi (y_{mN},-s) . 
%\label{njm.9}
%\eeq
%where
%\beq
%\chi^*(I_{m0}) :=
%\int_{I_{m0}}
%\langle \psi_{f0}(0) \vert x \rangle V(x) dx
%\label{njm.10}
%\eeq
%These integrals can be computed exactly for minimal uncertainty wave
%packets and gaussian interactions. 

%A direct application of this approach is not feasible due to the large
%number of cylinder sets that have to be summed ($\approx (M+1)^{N+1}$) over
%and the computational time needed to compute the complex probability
%associated with each of these sets.
  
\section{Factorization}

The input to the Muldowney-Nathanson-J{\o}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
$P(x_0,I_{m_11}, \cdots ,I_{m_NN})$ could be computed analytically,
there are $(M+1)^{N}$ cylinder sets in the limit that $M$ and $N$
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
\beq
\langle x_{N-1} \vert e^{-iH \Delta t} \vert
\psi_{i0} (-t) \rangle  \approx
\sum_m P(x_{N-1},I_{mN}) e^{
  -i V(y_{mN}) \Delta t }
\langle y_{mN} \vert \psi_{i0} (-t)\rangle . 
\label{f.1}
\eeq
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
\[
\langle x_{N-2} \vert e^{-iH 2\Delta t} \vert \psi_{i0} (-t) \rangle \approx
\sum_{m_{N-1}} P(x_{N-2},I_{m{(N-1)}}) e^{ -i V(y_{m{(N-1)}}) \Delta t
} \langle y_{m(N-1)} \vert e^{-iH \Delta t} \vert \psi_{i0} (-t) \rangle
\approx 
\]
\beq 
\sum_{m(N-1),{m{N}}} P(x_{N-2},I_{m({N-1}}) e^{-i
V(y_{m(N-1)}) \Delta t } P(y_{{(N-1)}} ,I_{mN}) e^{ -i
V(y_{mN}) \Delta t } \langle y_{mN} \vert \psi_{i0} (-t)\rangle . 
\label{f.2}
\eeq
Repeating this for all $N$ time steps gives the following approximation 
\beq
\langle x_0 \vert e^{-i H t} \vert \psi_{i0} (-t) \rangle \approx
\sum_{} P(x_0,I_{m1}) e^{ -i V(y_{m1}) \Delta t} 
\prod_{n=2}^N P(y_{m({n-1})},I_{nm}) e^{ -i V(y_{mn}) \Delta t}
\langle y_{mN} \vert \psi_{i0} (-t)\rangle 
\label{f.3}
\eeq
where the sum is over the cylinder sets.   The factorization
leads to the following approximation of the complex probability
on the cylinder set $\{I_{m1},\cdots , I_{mN}\}$:
\beq
P(x_0,I_{m1},\cdots ,I_{mN}) \approx 
P(x_0,I_{m1}) 
\prod_{n=2}^N P(y_{m({n-1})},I_{nm})
\label{f.4}
\eeq
With this approximation (\ref{s.9}) becomes 
\[
\langle x_0 \vert V   \vert e^{-i H t} \vert \psi_{i0} (-t) \rangle \approx
\]
\beq
\sum_{m1 \cdots mN}
V(x_0) 
P(x_0,I_{m1}) e^{ -i V(y_{m1}) \Delta t} 
\prod_{n=2}^N P(y_{m({n-1})},I_{nm}) e^{ -i V(y_{mn}) \Delta t}
\langle y_{mN} \vert \psi_{i0} (-t)\rangle . 
\label{f.5}
\eeq
This representation has a significant advantage over (\ref{njm.8}) 
because the matrix elements
\beq
K_{m,k} =
P(y_{m},I_{k}) e^{ -i V(y_{k}) \Delta t}
\label{f.6}
\eeq
where
\beq
P(y_{m},I_{k}) = 
({\mu \over 2 \pi i \Delta t})^{1/2}
\int_{x_k}^{x_{k+1}} dx
e^{
i {\mu \over 2\Delta t} (y_{m}-x)^2} 
=
\sqrt{{1  \over i \pi }}
\int_{\sqrt{{\mu \over 2\Delta t}}(x_k-y_m)}^{\sqrt{{\mu \over 2\Delta t}}(
x_{k+1}-y_m)}
e^{ 
i \beta^2} d\beta
\label{f.7}
\eeq
can  be computed analytically,  and  powers  of  this matrix  can  be
computed efficiently.  The integrals in (\ref{f.7}) for finite intervals 
are Fresnel integrals of the form
\beq
I[a,b] = \int_a^b e^{ix^2} dx =
%\sqrt{\pi \over 2} (C_{c}( x)+ i 
%S_{c}(x))\vert_a^b =
%\]
%\beq
\sqrt{\pi \over 2} (C_{c}(b)- C_{c}(a))
+i \sqrt{\pi \over 2} (S_{c}(b)- S_{c}(a)) .
\label{f.8}
\eeq
where
\beq
C_c(x) = \sqrt{{2 \over \pi}} \int_0^x \cos(t^2) dt
\qquad
S_c(x) = \sqrt{{2 \over \pi}} \int_0^x \sin(t^2) dt .
\label{f.9}
\eeq
Note that these definitions differ from the definitions of
Fresnel integrals given in \cite{Abramowitz} Abramowitz and Stegun.
They are related by
\beq
C_c (\sqrt{{\pi \over 2}}x) = C_{AS} (x)
\qquad
S_c (\sqrt{{\pi \over 2}}x) = S_{AS} (x).
\label{f.10}
\eeq

For the semi-infinite interval with  
$a=-\infty$
\beq
I[-\infty,b] = \int_{-\infty}^b e^{ix^2} dx =
{1 \over 2} \int_{-\infty}^{\infty} e^{ix^2} dx - \int_b^0 e^{ix^2} dx =
\sqrt{{\pi \over 2}} \left ({1+i \over 2} + C_c(b)+i S_c(b) \right ) .
\label{f.11}
\eeq
and for $b=\infty$
\beq
I[a,\infty] = \int_a^{\infty} e^{ix^2} dx =
{1 \over 2} \int_{-\infty}^{\infty} e^{ix^2} dx - \int_0^a e^{ix^2} dx =
\sqrt{{\pi \over 2}} \left ({1+i \over 2} - C_c(a)-i S_c(a) \right ) .
\label{f.11}
\eeq
Using these formulas leads to the following expressions for 
the one-step matrix $K_{mk}$ when $I_k$ is a finite interval:
\[
K_{mk} =
P(y_{m},I_{k})  e^{ -i V(y_{k}) \Delta t} =
\]
%\[ 
%\sqrt{{1  \over i \pi }}
%\left (
%\sqrt{\pi \over 2} (C_{c}(\sqrt{{\mu \over 2\Delta t}}(x_{k+1}-y_m))- 
%C_{c}(\sqrt{{m \over 2\Delta t}}(
%x_{k}-y_m)))
%\right . 
%\]
%\[
%\left .
%+i \sqrt{\pi \over 2} (S_{c}(\sqrt{{\mu \over 2\Delta t}}(x_{k+1}-y_m))- S_{c}(%\sqrt{{m \over 2\Delta t}}(
%x_{k}-y_m))) 
%\right ) e^{ -i V(y_{k}) \Delta t} =
%\]
\[ 
{1 \over 2}
\left (
(C_{c}(\sqrt{{\mu \over 2\Delta t}}(x_{k+1}-y_m))- 
C_{c}(\sqrt{{\mu \over 2\Delta t}}(x_{k}-y_m)) +
S_{c}(\sqrt{{\mu \over 2\Delta t}}(x_{k+1}-y_m))- 
S_{c}(\sqrt{{\mu \over 2\Delta t}}(x_{k}-y_m)))
\right . 
\]
\beq
\left .
+i (S_{c}(\sqrt{{\mu \over 2\Delta t}}(x_{k+1}-y_m))-
S_{c}(\sqrt{{\mu \over 2\Delta t}}(x_{k}-y_m))- 
C_{c}(\sqrt{{\mu \over 2\Delta t}}(x_{k+1}-y_m))+ 
C_{c}(\sqrt{{\mu \over 2\Delta t}}(x_{k}-y_m)))
\right ) e^{ -i V(y_{k}) \Delta t} . 
\label{f.12}
\eeq
For $x_{k}=x_0= -\infty $:
\[
K_{mk} =
P(y_{m},I_{k0})  e^{ -i V(y_{k}) \Delta t}
=
\]
%\[
%\sqrt{{1  \over i\pi}}
%\left ( 
%{1 \over 2} \int_{-\infty}^{\infty} e^{ix^2} dx - \int_b^0 e^{ix^2} dx
%\right ) e^{ -i V(y_{0}) \Delta t} =
%\]
%\[
% \sqrt{{1  \over i\pi}}
% \sqrt{{\pi \over 2}} \left ({1+i \over 2} + C_c(\sqrt{{\mu \over 2\Delta t}}
% (x_{1}-y_m))+
%i S_c(\sqrt{{\mu \over 2\Delta t}}(x_{1}-y_m)) \right )
%e^{ -i V(y_{k}) \Delta t} =
%\]
\beq
{1  \over 2}
\left ( 1 +
(C_{c}(\sqrt{{\mu \over 2\Delta t}}(x_1-y_m))+
S_{c}(\sqrt{{\mu \over 2\Delta t}}(x_1-y_m)))
+i (
S_{c}(\sqrt{{\mu \over 2\Delta t}}(x_1-y_m)) -
C_{c}(\sqrt{{\mu \over 2\Delta t}}(x_1-y_m)))
\right) 
e^{ -i V(y_{0}) \Delta t} 
\label{f.13}
\eeq
and for $x_{k+1}=x_{M+1} \infty $:
\[
K_{mM} =
P(y_{m},I_{M})  e^{ -i V(y_{M}) \Delta t} = 
\]
%\[
%\sqrt{{1  \over i \pi }}
%\left ( {1 \over 2} \int_{-\infty}^{\infty} e^{ix^2} dx - \int_0^a e^{ix^2} 
%dx \right ) e^{ -i V(y_{k}) \Delta t} =
%\]
%\[
% \sqrt{{1  \over i\pi}}
%\sqrt{{\pi \over 2}} \left ({1+i \over 2} - C_c(\sqrt{{m \over 2\Delta t}}(
%x_{M}-y_m)))-i S_c(\sqrt{{m \over 2\Delta t}}(
%x_{M}-y_m))) 
%\right )e^{ -i V(y_{M}) \Delta t} =
%\]
\beq
{1 \over 2} \left (
1 -( 
S_{c}(\sqrt{{m \over 2\Delta t}}(x_{M}-y_m))+ 
C_{c}(\sqrt{{m \over 2\Delta t}}(x_{M}-y_m))) 
-i(
S_{c}(\sqrt{{m \over 2\Delta t}}(x_{M}-y_m))- 
C_{c}(\sqrt{{m \over 2\Delta t}}(x_{M}-y_m)))
\right )
e^{ -i V(y_{M}) \Delta t} .
\label{f.14}
\eeq
Combining these approximations the expression for 
$\langle \psi_{f0} (0) \vert V \vert e^{-i H \Delta t} \vert \psi_{i0} (-t) \rangle$
is approximately given by
\beq
\langle x \vert  V \vert e^{-i H t} \vert \psi_{i0} (-t) \rangle \approx
\sum_{mk}
V(x)P(x,I_m)e^{-i V(y_m) t}K^{N-1}_{mk}
\langle y_{k} \vert \psi_{i0} (-t)\rangle 
\label{f.15}
\eeq
where $K_{mn}$ is the $(M+1)\times (M+1)$ matrix in equations
(\ref{f.12}-\ref{f.14}).  This matrix only requires $M+1$ values of the
potential.  The other elements needed for this computation are
the minimal uncertainty wave function at time $-t$ at the
same $M+1$ points and the potential.  

The above formulas are for one-dimensional scattering.  The
three-dimensional case involves products of these expressions.
Transition matrix elements can be extracted from (\ref{f.15}) using
(\ref{s.10}).  

\section{Computational considerations}

To test this method, approximate transition matrix elements are
computed for the example of a particle of mass $\mu$ scattering from
an attractive Gaussian potential in one dimension.  These computations
provide a laboratory that can be used to develop strategies for
formulating realistic computations on a quantum computer.

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) $N$ time slices
are chosen to be equally spaced and (2) the number, $M-1$, and width
$\delta x$ 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 $(p^2/2\mu)\Delta t$ and
$V\Delta t$ vanish.  In a computation these terms need to be small.
In the first term $p$ 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.  The potential needs to be
approximately constant on each spatial interval.  The Trotter condition
also requires that  $m \delta x^2 / 2 \Delta t $ is small.
%so the
%width of each interval $\delta x$  and the width of the time steps $\Delta t$
%are realted.  Specifically 
%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 $\Delta p$ small makes $\Delta x$ large
- which increases the time that the wave packet feels the potential.
Both $P(x,I_m)$ and $e^{i V(x) \Delta t}$ 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 ${\mu (\delta x)^2 \over 2 \Delta t}$ is small
requires $N/M^2$ to be small, where $N$ is the number of time steps
and $M$ is the number of intervals per time step.  This means that
shortening the time step may require 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.
\begin{table}
\caption{initial wave packet parameters}
\begin{tabular}{ll}
\hline
\hline
mass &$\mu = 1.0$\\
initial momentum &$p_i=   5.0$\\ 
momentum uncertainty &$\Delta p_i =0.25$\\ 
position uncertainty &$\Delta x= 2.0$\\
initial velocity &$v_i=p_i/\mu = 5.$\\
\hline
\end{tabular}
\end{table}
The particle scatters off of a Gaussian potential of range $r_0$ and strength
$-v_0$,
\beq
V(x) = -v_0 e^{-(x/r_0)^2}.
\label{cc.1}
\eeq
The values of the potential parameters used in the test
calculations are listed in table II.
\begin{table}
\caption{potential parameters}
\begin{tabular}{ll}
\hline
\hline 
strength &$ v0=5.0$\\ 
range &$r0=1.0$\\
\hline
\end{tabular}
\end{table}

The Trotter product formula is justified provided that the time steps
satisfy $(p^2/2m) \Delta t \approx 12.5 \Delta t$ and
$V \Delta t \approx 5 \Delta t$ are small.  The minimum total time is
$t= \mu (r_0+\Delta x)/p = .6$.  The calculations require a slightly
longer time for convergence, but convergence can be obtained with a
surprisingly large $\Delta t$.  The active volume is the sum of the
width of wave packet plus the range of the potential, which is about
$3.0$ units.  This must be decomposed into small intervals where the
potential and wave packet are approximately constant.

Graphical methods can be used to determine how long before the collision
the initial wave packet is out of the range of the interaction.
Figure 1 shows a plot of the potential and the initial wave packet at
$t=-3$, before it feels the effects of the potential.  The
solid curve represents the potential.  The dashed curve represents the
real part of the initial wave function at $t=-3$ and the dotted curve
represents the imaginary part of the initial wave function at $t=-3$.
The non-interacting wave packet is evolved to $t=-3$ from $t=0$ so it
includes the spreading of the wave packet.
%At $t=-3$ the wave packet is moving to the right and is not yet in the
%range of the potential.
This figure suggests that for this initial wave packet and potential that
$\Omega (-\infty) \vert \psi_{i0} \rangle \approx
\Omega (-3) \vert \psi_{i0} \rangle $.
%
%$t=-3$ is a sufficient time to start the
%computation 
%
%Figure 1: Shows the potential (\ref{cc.1}) in configuration space.  
%It gives a more precise measure of the full width of 
%the potential, which just over about $4.0$ units. 
%
%%figure 1:
%
%\begin{figure}
%\centering
%\includegraphics[scale=.5]{fig_1.pdf}
%\caption{V(x)}
%\label{figure 1}
%\end{figure}  
%
%Figure 2: Shows the real and imaginary parts of the initial
%free wave packet $\langle x \vert \psi_{0i}(0) \rangle $ and plots of the 
%potential on the same graph.
What is relevant
is the combined width of the initial wave packet and the potential
%In this figure the potential is represented by the solid line, the
%real part of the initial wave packet is represented by the dashed line
%and the imaginary part of the initial wave packet is represented by
%the dotted line.
which gives an estimate of the active volume
where the particle feels the effect of the potential. 
%needed for a calculation.
%An estimate of the minimum $t$ needed for
%convergence is the time it takes for the wave packet to no longer feel
%the effects of the potential.  Figure 2 shows that the width initial
%wave packet is about twice the width of the potential.
%\begin{figure}
%\centering
%\includegraphics[scale=.5]{fig_2.pdf}
%\caption{V(x),$\langle x \vert \psi_{0i}(0) \rangle$}
%\label{figure 2}
%\end{figure}  
Figure 1 indicates that the active volume is about 12
units.  For a wave packet moving with speed $v=p_i/m=-5$ the wave
packet will travel 15 units in a time $t=3.0$.  This should be
sufficiently long to move the potential out of the range of the
potential.
%Note that evolving the free wave packet to
%$t=-3$ accounts for the effects of the spreading of the wave packet. 
%
%if the spreading of the free wave packet is ignored.  The
%effect of the spreading can be included by comparing the free wave
%packet at time $t=-3$ with the potential. This is illustrated in
%figure 3 which plots the real and imaginary parts of the initial free
%wave packet at time $t=-3$ with the potential.  As in figure 2 the
%potential is represented by the solid line, the real part of the
%initial wave packet is represented by the dashed line and the
%imaginary part of the initial wave packet is represented by the dotted
%line.
Figure 1 shows that even with the effects of wave-packet
spreading, at $t=-3.0$ the wave packet has not reached the range of
the potential.  This suggests that $t=3.0$ is a good first guess at
the time $t$ sufficient for convergence of
$\Omega(-t)\vert \psi_{0i}(0) \rangle$.

%\begin{figure}
%\centering
%\includegraphics[scale=.5]{psi_minus_v_0_t2.pdf}
%\caption{V(x),$\psi(x,-2)$}
%\label{figure 3}
%\end{figure}  

\begin{figure}
\centering
%\includegraphics[scale=.5]{fig_3.pdf}
\includegraphics[scale=.5]{figure_1.pdf}
%\caption{Potential and initial wave packet at $t=-3$}
\caption{Plot of the potential and the incident initial free wave packet
at time $t=-3$.}
\label{figure 1}
\end{figure}  

Numerical calculations are possible because of the factorization of
the complex probabilities into products of matrix elements of one-step
probabilities.  Table 3 shows the sum of M=5000 one-step complex
probabilities computed at different points.  The computation shows
that these quantities behave like complex probabilities.  The real
part of the sum of 5000 complex probabilities is always 1 and the
imaginary part is always 0, independent of the final $x$ value.  The
table indicates the stability of the sum of these large numbers of
one-step complex probabilities since the cancellation of all of the
imaginary terms is accurate to 15-17 significant figures.
\begin{table}
\caption{Sums of the real and imaginary parts of the complex probabilities.} 
\begin{tabular}{lll}
\hline
\hline
$x$ & $\sum_n \Re (P_n(x))$ &$\sum_n \Im (P_n(x))$\\
\hline
-25.005001 & p=1.000000$\times 10^0$ & + i 5.273559$\times 10^{-16}$ \\
%-24.604921 p=1.000000e+00 + i(9.992007e-16)
%-24.204841 p=1.000000e+00 + i(5.828671e-16)
%-23.804761 p=1.000000e+00 + i(-1.415534e-15)
%-23.404681 p=1.000000e+00 + i(-1.193490e-15)
%-23.004601 p=1.000000e+00 + i(4.996004e-16)
%-22.604521 p=1.000000e+00 + i(6.383782e-16)
%-22.204441 p=1.000000e+00 + i(-1.387779e-16)
%-21.804361 p=1.000000e+00 + i(1.110223e-16)
%-21.404281 p=1.000000e+00 + i(1.942890e-16)
%-21.004201 p=1.000000e+00 + i(3.330669e-16)
%-20.604121 p=1.000000e+00 + i(-1.332268e-15) 
-20.204041 &p=1.000000$\times 10^0$ &- i 1.387779$\times 10^{-16}$ \\
%-19.803961 p=1.000000e+00 + i(-1.471046e-15)
%-19.403881 p=1.000000e+00 + i(-1.304512e-15)
%-19.003801 p=1.000000e+00 + i(-6.106227e-16)
%-18.603721 p=1.000000e+00 + i(-9.714451e-16)
%-18.203641 p=1.000000e+00 + i(-2.498002e-16)
%-17.803561 p=1.000000e+00 + i(-4.440892e-16)
%-17.403481 p=1.000000e+00 + i(4.440892e-16)
%-17.003401 p=1.000000e+00 + i(-2.775558e-16)
%-16.603321 p=1.000000e+00 + i(3.608225e-16)
%-16.203241 p=1.000000e+00 + i(1.221245e-15)
%-15.803161 p=1.000000e+00 + i(8.604228e-16)
%-15.403081 p=1.000000e+00 + i(2.220446e-16)
-15.003001 & p=1.000000$\times 10^0$ & + i 2.775558$\times 10^{-17}$ \\
%-14.602921 p=1.000000e+00 + i(-9.436896e-16)
%-14.202841 p=1.000000e+00 + i(1.249001e-15)
%-13.802761 p=1.000000e+00 + i(-8.326673e-17)
%-13.402681 p=1.000000e+00 + i(-2.775558e-16)
%-13.002601 p=1.000000e+00 + i(1.665335e-16)
%-12.602521 p=1.000000e+00 + i(2.498002e-16)
%-12.202440 p=1.000000e+00 + i(0.000000e+00)
%-11.802360 p=1.000000e+00 + i(-2.498002e-16)
%-11.402280 p=1.000000e+00 + i(-6.661338e-16)
%-11.002200 p=1.000000e+00 + i(-9.992007e-16)
%-10.602120 p=1.000000e+00 + i(-8.049117e-16)
-10.202040 & p=1.000000$\times 10^0$ &+ i 2.775558$\times 10^{-16}$ \\
%-9.801960  p=1.000000e+00 + i(2.775558e-17)
%-9.401880  p=1.000000e+00 + i(9.714451e-16)
%-9.001800  p=1.000000e+00 + i(1.110223e-15)
%-8.601720  p=1.000000e+00 + i(-6.383782e-16)
%-8.201640  p=1.000000e+00 + i(0.000000e+00)
%-7.801560  p=1.000000e+00 + i(3.608225e-16)
%-7.401480  p=1.000000e+00 + i(1.554312e-15)
%-7.001400  p=1.000000e+00 + i(5.551115e-17)
%-6.601320  p=1.000000e+00 + i(-4.996004e-16)
%-6.201240  p=1.000000e+00 + i(-3.053113e-16)
%-5.801160  p=1.000000e+00 + i(-1.942890e-16)
%-5.401080  p=1.000000e+00 + i(-6.383782e-16) 
-5.001000  & p=1.000000$\times 10^0$ &+ i 7.771561$\times 10^{-16}$ \\
%-4.600920  p=1.000000e+00 + i(8.049117e-16)
%-4.200840  p=1.000000e+00 + i(1.498801e-15)
%-3.800760  p=1.000000e+00 + i(1.110223e-15)
%-3.400680  p=1.000000e+00 + i(1.415534e-15)
%-3.000600  p=1.000000e+00 + i(1.471046e-15)
%-2.600520  p=1.000000e+00 + i(2.498002e-16)
%-2.200440  p=1.000000e+00 + i(-5.551115e-17)
%-1.800360  p=1.000000e+00 + i(5.273559e-16)
%-1.400280  p=1.000000e+00 + i(9.436896e-16)
%-1.000200  p=1.000000e+00 + i(-7.494005e-16)
%-0.600120  p=1.000000e+00 + i(-1.276756e-15)
%-0.200040  p=1.000000e+00 + i(-7.494005e-16)
0.200040   &p=1.000000$\times 10^0$ &+ i 3.608225$\times 10^{-16}$ \\
%0.600120   p=1.000000e+00 + i(-4.996004e-16)
%1.000200   p=1.000000e+00 + i(-3.330669e-16)
%1.400280   p=1.000000e+00 + i(6.938894e-16)
%1.800360   p=1.000000e+00 + i(1.137979e-15)
%2.200440   p=1.000000e+00 + i(1.748601e-15)
%2.600520   p=1.000000e+00 + i(9.992007e-16)
%3.000600   p=1.000000e+00 + i(1.471046e-15)
%3.400680   p=1.000000e+00 + i(-1.110223e-16)
%3.800760   p=1.000000e+00 + i(-2.220446e-16)
%4.200840   p=1.000000e+00 + i(-7.494005e-16)
%4.600920   p=1.000000e+00 + i(1.137979e-15)
5.001000   &p=1.000000$\times 10^0$ &+ i 1.665335$\times 10^{-16}$ \\
%5.401080   p=1.000000e+00 + i(-1.693090e-15)
%5.801160   p=1.000000e+00 + i(-1.443290e-15)
%6.201240   p=1.000000e+00 + i(-3.330669e-16)
%6.601320   p=1.000000e+00 + i(-8.881784e-16)
%7.001400   p=1.000000e+00 + i(2.775558e-16)
%7.401480   p=1.000000e+00 + i(1.415534e-15)
%7.801560   p=1.000000e+00 + i(1.859624e-15)
%8.201640   p=1.000000e+00 + i(6.383782e-16)
%8.601720   p=1.000000e+00 + i(1.165734e-15)
%9.001800   p=1.000000e+00 + i(9.159340e-16)
%9.401880   p=1.000000e+00 + i(2.498002e-16)
%9.801960   p=1.000000e+00 + i(-1.942890e-16)
10.202040  &p=1.000000$\times 10^0$ &+ i 5.273559$\times 10^{-16}$ \\
%10.602120  p=1.000000e+00 + i(3.053113e-16)
%11.002200  p=1.000000e+00 + i(2.220446e-16)
%11.402280  p=1.000000e+00 + i(-8.604228e-16)
%11.802360  p=1.000000e+00 + i(-3.608225e-16)
%12.202440  p=1.000000e+00 + i(-5.828671e-16)
%12.602521  p=1.000000e+00 + i(-1.554312e-15)
%13.002601  p=1.000000e+00 + i(5.551115e-17)
%13.402681  p=1.000000e+00 + i(4.163336e-16)
%13.802761  p=1.000000e+00 + i(6.661338e-16)
%14.202841  p=1.000000e+00 + i(-6.938894e-16)
%14.602921  p=1.000000e+00 + i(8.326673e-17)
15.003001  &p=1.000000$\times 10^0$ &+ i 8.049117$\times 10^{-16}$ \\
%15.403081  p=1.000000e+00 + i(1.415534e-15)
%15.803161  p=1.000000e+00 + i(7.494005e-16)
%16.203241  p=1.000000e+00 + i(0.000000e+00)
%16.603321  p=1.000000e+00 + i(1.498801e-15)
%17.003401  p=1.000000e+00 + i(-1.554312e-15)
%17.403481  p=1.000000e+00 + i(4.996004e-16)
%17.803561  p=1.000000e+00 + i(-3.608225e-16)
%18.203641  p=1.000000e+00 + i(2.498002e-16)
%18.603721  p=1.000000e+00 + i(1.526557e-15)
%19.003801  p=1.000000e+00 + i(-6.661338e-16)
%19.403881  p=1.000000e+00 + i(3.885781e-16)
%19.803961  p=1.000000e+00 + i(-7.494005e-16)
20.204041  &p=1.000000$\times 10^0$&+ i 1.137979$\times 10^{-15}$ \\
%20.604121  p=1.000000e+00 + i(-6.383782e-16)
%21.004201  p=1.000000e+00 + i(8.049117e-16)
%21.404281  p=1.000000e+00 + i(5.551115e-16)
%21.804361  p=1.000000e+00 + i(2.775558e-17)
%22.204441  p=1.000000e+00 + i(-3.053113e-16)
%22.604521  p=1.000000e+00 + i(2.775558e-17)
%23.004601  p=1.000000e+00 + i(-8.881784e-16)
%23.404681  p=1.000000e+00 + i(-9.714451e-16)
%23.804761  p=1.000000e+00 + i(-1.165734e-15)
%24.204841  p=1.000000e+00 + i(-1.165734e-15)
%24.604921  p=1.000000e+00 + i(5.551115e-16)
25.005001  &p=1.000000$\times 10^0$ &- i 2.164935$\times 10^{-15}$\\
\hline                                          
\end{tabular} 
\end{table} 

The accuracy of the numerical computation of the time evolution of the
initial wave packet depends on having sufficiently small $\delta x$
and $\Delta t$.  The limiting size depends on the initial wave packet.

Since the time evolution of the free wave packet can be computed
analytically, one test of accuracy of the free evolution based on
using products of one-step probabilities is to start with the exact
initial wave packet at $t=-3$ and transform it back to the initial
time, $t=0$, using multiple applications of the one-step probability
matrices.  This can then be compared to
the exact initial wave packet at $t=0$.  In this test $t=-3$ is chosen
because the wave packets should be in the asymptotic region
at that time.  This test uses 30 time slices separated by $\Delta_t=.1$ and
$5000$ spatial steps between $-25$ and $25$ corresponding to a spatial
resolution $\delta_x=.01$.  The result of this calculation is
shown in figure 2.  In these plots the dashed lines
represent the calculated $t=0$ wave function while the dotted line
represents the exact $t=0$ wave function given in (\ref{fs.4}).

These figures compare
\beq
\psi_{i0}(x,0) 
\qquad \mbox{to} \qquad
\sum P(x,I_{n1}) P(y_{n1},I_{n2}) \cdots
P(y_{n19 },I_{n20})\langle y_{n20} \vert \psi_{0i}(-3) \rangle .
\label{cc.2}
\eeq
In both panels the real and imaginary parts of both wave
functions fall on top of each other.  
%Table 4 shows the numerical values of some of the data plotted in these
%graphs

\begin{figure}
\begin{minipage}[t]{.45\linewidth}
%\begin{subfigure}[t]{.45\linewidth}
\centering
\includegraphics[scale=.5]{figure_2a.pdf}  
%\includegraphics[scale=.5]{fig_4.pdf}
%%\caption{Re( $\langle x \vert \psi_{0i}(0) \rangle$)
%and Re($\langle x \vert U_0 (3) \vert \psi_{0i}(-3) \rangle$)
%using 30 time steps}
%\caption{real part}
%\label{figure 2}
%\end{figure}  
\end{minipage}
%\end{subfigure}
%\begin{subfigure}
\begin{minipage}[t]{.45\linewidth}
% \begin{figure}\
\centering
\includegraphics[scale=.5]{figure_2b.pdf}  
%\includegraphics[scale=.5]{fig_5.pdf}
%\caption{Imaginary part}
%\caption{Im( $\langle x \vert \psi_{0i}(0) \rangle$)
%  and Im$\langle x \vert U_0(3) \vert \psi_{0i}(-3) \rangle$)
%using 30 time steps}
%\label{figure 3}
\end{minipage}
\label{figure 2}
%\end{subfigure}
\caption{Plots of the real and imaginary parts of the free wave packet
at time $t=0$ and the real and imaginary parts of the free path integral
evolved wave packet  
from $t=-3$ to $t=0$.}
\end{figure}  
A more detailed examination of the difference appears in figure 3,
which shows the absolute value of the difference of the real and
imaginary parts of the exact and approximate evolved free wave
function divided by the real and imaginary parts of the exact wave
function.  The spikes in the relative error plots are located near the
zeros of the exact wave functions.  The relative error is smallest
between $-4.0$ and $4.0$, which is at the center of the wave packet.
In this region, with the exception of the spike in the imaginary part
at 0, the relative error for the real and imaginary parts of the wave
function in this region is bounded by .007.  Outside of this region
the relative error gets larger.  This is in part because the
path integral has a some fixed numerical uncertainty while   
exact wave function falls off exponentially, so the expression for
the relative error has a denominator that exponentially
approaches 0.    In order to calculate the transition matrix elements
using (\ref{s.10}) or (\ref{s.11}) the exact wave function (with
interaction) only needs to be accurate only in the range of the
potential.

\begin{figure}
\begin{minipage}[t]{.45\linewidth}  
\centering
%\includegraphics[scale=.5]{fig_4a_h9.pdf}
\includegraphics[scale=.5]{figure_3a.pdf}
%\caption{Relative error in the real part of the free path integral
%evolved wave packet
%from $t=-3$ to $t=0$.}
%\caption{Re($\vert \langle x \vert \psi_{0i}(0) \rangle -
%\langle x \vert U_0 (3) \vert \psi_{0i}(-3) \rangle\vert/.4$)
%using 30 time steps}
%\label{figure 4}
%\end{figure}  
\end{minipage}
%\begin{figure}
\begin{minipage}[t]{.45\linewidth}  
\centering
%\includegraphics[scale=.5]{fig_5a_h9.pdf}
\includegraphics[scale=.5]{figure_3b.pdf}
%caption{Relative error in the imaginary part of the free path integral evolved
%wave packet
%evolved from $t=-3$ to $t=0$.}
%\caption{Im($\vert \langle x \vert \psi_{0i}(0) \rangle -
% \langle x \vert U_0(3) \vert \psi_{0i}(-3) \rangle \vert /.4$)
% using 30 time steps}
\end{minipage}
\label{figure 3}
\caption{Relative error in the real and imaginary parts of the free path
integral
evolved wave packet
from $t=-3$ to $t=0$.}
\end{figure}  
These comparisons indicate that both the time steps and resolution are
sufficiently small to evolve the free wave packet for $t=3.0$

In these calculations there was no attempt at efficiency; however
because of the analytic expressions for the one-step
probabilities, the one-step probabilities were computed on the fly in
order to avoid storing large matrices.  The resolution was chosen
to be sufficiently fine to accurately represent the initial wave packet.

The next figure illustrates the effect of the potential on the
evolution of the wave function.  Figure 4 compares the real
and imaginary parts of
the scattering wave function
$\langle x \vert \psi_{i}(0)\rangle= \langle x \vert \Omega
(-3.0)\vert \psi_{i0}(0)\rangle$ (dotted line) at time $t=0$ to the
initial wave packet $\langle x \vert \psi_{i0} (t=0)\rangle$ at $t=0$.
The parameters used for this computation are $\Delta_t=.1$,
$\delta_x=.01$ and $x\in [-25,25]$, which are the same parameters used
to produce figures 2 and 3.  The functions are calculated at the
midpoint of each finite interval.  Figure 4
shows the change in phase as the interaction is turned on.  The dashed
lines in both panels represent the initial wave packet, while the
dotted lines represent the interacting wave packet.

These calculations approximate the incident scattered
wave function at the collision time $(t=0)$.
As discussed in section 2, the exact incident wave function at $t=0$
(or any common time) is needed to compute the differential cross section.  These
figures show that both the free and interacting $t=0$ wave
functions occupy approximately the same volume, (about 20 units) which
shows that using $\Omega(t)$ to calculate the scattering wave function
eliminates spreading of the wave function.
At time $t=0$ the phase of the incident wave is more
dominantly shifted in the forward (right) side of the wave packet, which
has spent the most time interacting with the interaction.
This why these wave
functions do not look like separated transmitted and reflected waves.
\begin{figure}
\begin{minipage}[t]{.45\linewidth}
\centering
%\includegraphics[scale=.5]{fig_6.pdf}
\includegraphics[scale=.5]{figure_4a.pdf}
%\caption{Shift in the phase of the real part of the incident wave packet at the% collision time $t=0$}
%\caption{Re($\langle x \vert \psi_{0i}(0) \rangle$)
%and Re($\langle x \vert \Omega(t=-3)\vert \psi_{i0}(0) \rangle$)}
%\label{figure 6}
%\end{figure}
\end{minipage}
\begin{minipage}[t]{.45\linewidth}
%\begin{figure}
\centering
%\includegraphics[scale=.5]{fig_7.pdf}
\includegraphics[scale=.5]{figure_4b.pdf}
%\caption{Shift in the phase of the imaginary part of the incident wave packet at the collision time $t=0$}
%\caption{Im($\langle x \vert \psi_{0i}(0) \rangle$)
%and Im($\langle x \vert \Omega(t=-3)\vert\psi_{i0}(0) \rangle$)}
%\label{figure 7}
\end{minipage}
\caption{Shift in the phase of the real and imaginary parts of the incident wave packet at the collision time $t=0$}
\label{figure 4}
\end{figure}  

One of the complications of performing scattering calculations using
real time path-integral methods is that the final expression for the
transition matrix elements involves several approximations that
need to be tested for convergence.
With the potential turned on it is important to check that (1) the
volume $[-25,25]$ is sufficiently large, the time step $\Delta t=.1$
is sufficiently small, the resolution $\delta_x=.01$ is sufficiently
small, the total time $t=3.0$ is sufficiently large,  and the calculation
is stable with respect to changing the sample point $y_i$ in the
interval $I_i$.  The plots in figures 5-10 investigate the sensitivity 
of the time $t=0$ interacting wave functions to variations of the
parameters used in the calculations shown in figure 4.

Note that all of the calculations assume that the evolved wave
functions vanish on the half infinite intervals.  This is justified
both graphically and because the wave packet remains square
integrable.

The Trotter product formula is justified in the small time step limit.
The calculations illustrated in figure 4 used a time step
$\Delta_t=.1$.  Figure 5 compares the real and imaginary parts of
the scattered wave function $\langle x \vert \psi_{i}(0)\rangle= 
\langle x \vert \Omega (-3.0)\vert \psi_{i0}(0)\rangle$
using a time step size of $\Delta_t=.1$ (dashed curve) with
corresponding calculations using a time step of half the size,
$\Delta_t=.05$ (dotted curve).  The plots of the real and imaginary
parts of the scattered $t=0$ wave functions for $\Delta_t=.1$ and
$\Delta_t=.05$ fall on top of each other.  This indicates that the
time resolution $\Delta_t$ is sufficiently fine for this calculation
\begin{figure}
\begin{minipage}[t]{.45\linewidth}
\centering
\includegraphics[scale=.5]{figure_5a.pdf}
%\includegraphics[scale=.5]{fig_8.pdf}
%\caption{Plots of the real part of the time $t=0$
%scattered wave packet
%calculated using Trotter time steps of $\Delta t =.1$ and
%$\Delta t=.05$.
%}
%\caption{Re($\langle x \vert \Omega (-3) \vert \psi_{i0} \rangle$) for time steps $\Delta t =.1,.05$}
%\label{figure 8}
%\end{figure}  
\end{minipage}
\begin{minipage}[t]{.45\linewidth}
%\begin{figure}
\centering
%\includegraphics[scale=.5]{fig_9.pdf}
\includegraphics[scale=.5]{figure_5b.pdf}
%\caption{Plots of the imaginary part of the time $t=0$
%scattered wave packet
%calculated using Trotter time steps of $\Delta t =.1$ and
%$\Delta t =.05.$
%}
%\caption{Im($\langle x \vert \Omega (-3) \vert \psi_{i0} \rangle$) for time steps $\Delta t =.1,.05$}
%\label{figure 9}
\end{minipage}
\caption{Plots of the real and imaginary parts of the time $t=0$
scattered wave packet
calculated using Trotter time steps of $\Delta t =.1$ and
$\Delta t =.05.$
}
\label{figure 5}
\end{figure}  

While the Trotter approximation requires a sufficiently small time
step, an accurate approximation of the wave operators requires a
sufficiently large time.  The initial choice of the approximating
$\Omega_- := \lim_{t \to -\infty}\Omega(t)$ by $\Omega (-3)$ was
determined by examining the range of the potential and width and speed
of the wave packet.  Graphical methods indicated that the incident wave
packed was still in the asymptotic region at $t=-3$.  
Figure 6 shows the
effect of increasing the time from $t=3.0$ to $t=6.0$, keeping the
size of the time step $\Delta_t=.1$ constant, on the real and
imaginary parts of the calculated scattering wave function.  The plots
of the wave functions for $t=-3.0$ (dashed curve) and $t=-6.0$ (dotted
curve) fall on top of each other.  This suggests that $t=3.0$ is
sufficient for convergence.

\begin{figure}
\begin{minipage}[t]{.45\linewidth}
\centering
%\includegraphics[scale=.5]{fig_10.pdf}
\includegraphics[scale=.5]{figure_6a.pdf}
% \caption{Plots of the real part of the scattered wave packet
%using approximate wave operators evaluated at $t=-3$ and $t=-6$.}  
%\caption{Re($\langle x \vert \Omega (-3) \vert \psi_{i0} \rangle$) and Re($\langle x \vert \Omega (-6) \vert \psi_{i0} \rangle$)}
%\label{figure 10}
%\end{figure}  
\end{minipage}
\begin{minipage}[t]{.45\linewidth}
%\begin{figure}
\centering
\includegraphics[scale=.5]{figure_6b.pdf}
%\includegraphics[scale=.5]{fig_11.pdf}
%\caption{Plots of the imaginary part of the scattered wave packet
%using approximate wave operators evaluated at $t=-3$ and $t=-6$.}  
%\caption{Im($\langle x \vert \Omega (-3) \vert \psi_{i0} \rangle$) and Im($\langle x \vert \Omega (-6) \vert \psi_{i0} \rangle$)}
%\label{figure 11}
\end{minipage}
\caption{Plots of the real and imaginary parts of the scattered wave packet
using approximate wave operators evaluated at $t=-3$ and $t=-6$.}  
\label{figure 6}
\end{figure}  

The other limit used in this formulation of the path integral is the
decomposition of the volume at each time slice into sufficiently small
intervals. The
spatial resolution of the intervals that define the cylinder sets
should be sufficiently small that the wave function is effectively
constant on them.  The required resolution must be small 
due to the oscillating nature of the wave function,
Figure 7 shows the effect of increasing the
spatial resolution from 5000 intervals ($\delta_x=.01$, dashed curves)
to 10000 intervals ($\delta_x=.005$, dotted curves) on the real and
imaginary parts of the scattering wave function.  In these
calculations $\Delta_t=.1$, $t=3.0$ and the interval is $[-25,25]$.
Again the curves for the real and imaginary part of the wave function
fall on top of each other.  This indicates that for this problem a
spatial resolution $\delta_x=.01$ (5000 intervals)) is sufficient for
convergence.

\begin{figure}
\begin{minipage}[t]{.45\linewidth}
\centering
%\includegraphics[scale=.5]{fig_12.pdf}
\includegraphics[scale=.5]{figure_7a.pdf}
%\caption{Plots of the real part of the scattered
%wave packet using cylinder sets with spacial resolutions of $.01$ and $.005$.}
%\caption{Re($\langle x \vert \Omega (-3) \vert \psi_{i0} \rangle$) and Re($\delta_x=.01, .005$).}
%\label{figure 12}
%\end{figure}  
\end{minipage}
\begin{minipage}[t]{.45\linewidth}
%\begin{figure}
\centering
%\includegraphics[scale=.5]{fig_13.pdf}
\includegraphics[scale=.5]{figure_7b.pdf}
%\caption{Plots of the imaginary part of the scattered
%wave packet using cylinder sets with spacial resolutions of $.01$ and $.005$.}
%%\caption{Im($\langle x \vert \Omega (-3) \vert \psi_{i0} \rangle$) and Im($\delta_x=.01, .015$)}
%\label{figure 13}
\end{minipage}
\caption{Plots of the real and imaginary parts of the scattered
wave packet using cylinder sets with spacial resolutions of $.01$ and $.005$.}
\label{figure 7}
\end{figure}  

While the previous plots suggest that wave functions vanish outside of
the volume $[-25,25]$, it is still important to check that the results
are stable with respect to increasing the active volume of the
calculation.  Figure 8 compares calculations of the real and
imaginary parts of the wave function where the volume is changed from
$[-25,25]$ (dotted curves) to $[-50,50]$ (dashed curves) keeping
$\delta_x=.01$, $\Delta_t=.1$, and $t=3.$ Again, the calculations
indicate the volume $[-25,25]$ is sufficient for convergence.


\begin{figure}
\begin{minipage}[t]{.45\linewidth}
\centering
%\includegraphics[scale=.5]{fig16.pdf}
%\includegraphics[scale=.5]{fig_14.pdf}
\includegraphics[scale=.5]{figure_8a.pdf}
%\caption{Plots of the real part of the scattered wave packet
%computed in the regions $[-25,25]$ and $[-50,50]$.
%}
%\caption{Re($\langle x \vert \Omega (-3) \vert \psi_{i0} \rangle$) for Vol$=\pm 25,\pm 50$}
%\label{figure 14}
%\end{figure}  
\end{minipage}
\begin{minipage}[t]{.45\linewidth}
% \begin{figure}
\centering
%\includegraphics[scale=.5]{fig17.pdf}
%\includegraphics[scale=.5]{fig_15.pdf}
\includegraphics[scale=.5]{figure_8b.pdf}
%\caption{Plots of the imaginary part of the scattered wave packet
%computed in the regions $[-25,25]$ and $[-50,50]$.
%}
%\caption{Im($\langle x \vert \Omega (-3) \vert \psi_{i0} \rangle$) for Vol$=\pm 25,\pm 50$}
%\label{figure 15}
\end{minipage}
\caption{Plots of the real and imaginary parts of the scattered wave packet
computed in the regions $[-25,25]$ and $[-50,50]$.
}
\label{figure 8}
\end{figure}

The size of the spatial intervals, $\delta x$, should be sufficiently
small that the potential and initial wave free wave packet are
approximately constant on the intervals.  Figure 9 compares
calculations where the potential is evaluated at the center or left
endpoint of each interval for $\delta x=.01$. Figure 9 shows the real
and imaginary parts of the wave function, where the sample points are at
the left
(dash-dot line) and the center (dashed line) of the interval.
This tests whether the potential is
locally constant on each interval.  These graphs show a small shift in
the overall phase of the wave function.  Figure 10 repeats the
calculations shown in Figure 9 by increasing the resolution
by a factor of 2.  These figures show a corresponding reduction in the
difference between the two curves.  The conclusion is that the
calculations in figure 4 seem to be most sensitive to the choice of
the position in each interval used to evaluate the potential.  This
sensitivity may be due to large slope
near the edges of the Gaussian potential shown in figure 1. 
\begin{figure}
\begin{minipage}[t]{.45\linewidth} 
\centering
%\includegraphics[scale=.5]{fig14a.pdf}
%\includegraphics[scale=.5]{fig_16.pdf}
\includegraphics[scale=.5]{figure_9a.pdf}
%\caption{Plots of the real part of the scattered wave packet
%with the potential evaluated the at the center or left end of each
%interval $\delta x=.01$.  
%}
%\caption{Re($\langle x \vert \Omega (-3) \vert \psi_{i0} \rangle$) for left vs center evaluation of $V(x)$}
%\label{figure 16}
%\end{figure}
\end{minipage}
\begin{minipage}[t]{.45\linewidth}
%\begin{figure}
\centering
%\includegraphics[scale=.5]{fig15a.pdf}
%\includegraphics[scale=.5]{fig_17.pdf}
\includegraphics[scale=.5]{figure_9b.pdf}
%\caption{Plots of the imaginary part of the scattered wave packet
%with the potential evaluated the at the center or left end of each
%interval for $\delta x=.01$. 
%}
%\caption{Im($\langle x \vert \Omega (-3) \vert \psi_{i0} \rangle$) for left vs center evaluation of $V(x)$}
%\label{figure 17}
\end{minipage}
\caption{Plots of the real and imaginary parts of the scattered wave packet
with the potential evaluated the at the center or left end of each
interval for $\delta x=.01$. 
}
\label{figure 9}
\end{figure}
 
\begin{figure}
\begin{minipage}[t]{.45\linewidth} 
\centering
%\includegraphics[scale=.5]{fig14.pdf}
%\includegraphics[scale=.5]{fig_18.pdf}
\includegraphics[scale=.5]{figure_10a.pdf}
%\caption{Plots of the real part of the scattered wave packet
%with the potential evaluated the at the center or left end of each
%interval $\delta x=.005$.
%}
%\caption{Re($\langle x \vert \Omega (-3) \vert \psi_{i0} \rangle$) for left vs center evaluation of $V(x)$}
%\label{figure 18}
% \end{figure}
\end{minipage}
\begin{minipage}[t]{.45\linewidth}
%\begin{figure}
\centering
%\includegraphics[scale=.5]{fig15.pdf}
%\includegraphics[scale=.5]{fig_19.pdf}
\includegraphics[scale=.5]{figure_10b.pdf}
%\caption{Plots of the imaginary part of the scattered wave packet
%with the potential evaluated the at the center or left end of each
%interval for $\delta x=.005$.
%}
%\caption{Im($\langle x \vert \Omega (-3) \vert \psi_{i0} \rangle$) for left vs center evaluation of $V(x)$}
%\label{figure 19}
\end{minipage}
\caption{Plots of the real and imaginary parts of the scattered wave packet
with the potential evaluated the at the center or left end of each
interval for $\delta x=.005$.
}
\label{figure 10}
\end{figure}

The goal of this work is to determine if these scattering wave
functions can be used to calculate sharp-momentum transition matrix
elements.  Because the limits in the Trotter product formula are
strong limits, the initial momentum was replaced by a narrow wave
packet.  To test the effect of the smearing, the sharp-momentum Born
approximation is compared to the Born approximation where the initial
sharp momentum state is replaced by a Gaussian delta function of width
.25, which was used in the calculations above.

This comparison is illustrated in figure 11.  The solid curve shows
the Gaussian approximation to a delta function with width $\Delta p=.25$.
The dotted curve shows the potential with initial momentum $p=5$,
as a function of the final momentum.  The dashed curve shows the 
same potential with the initial momentum replaced by the Gaussian 
delta function state centered at $p=5$ as a function of the final momentum. 
\begin{figure}
\centering
%\includegraphics[scale=.5]{fig19.pdf}
%\includegraphics[scale=.5]{fig_20.pdf}
\includegraphics[scale=.5]{figure_11.pdf}
\caption{Plots of the delta function normalized Gaussian,
the exact Born approximation to $T$, and the Born approximation
smeared with the initial Gaussian wave packet.}
%\label{figure 23}
%\label{figure 20}
\label{figure 11}
\end{figure}  
The figure shows a small decrease in the matrix elements due to 
the smearing near the on-shell value.
 
As another test of the numerical convergence,  the smeared Born 
approximation $\langle p \vert V \vert \Psi_{i0}(0)\rangle$ is compared to 
$\langle p \vert V \Omega_0 (V=0,t=3) \vert \Psi_{i0}(-3)\rangle$, where
the {\it potential is turned off} to compute $\Omega_0(3)$.  
In this calculation 
$\delta x= .01$ and $\Delta t =.1$.  This comparison is shown in figure
12.  The imaginary part of $\langle p \vert V \Omega_0 (V=0,t=-3) \vert \Psi_{i0}\rangle$, shown in the solid line should vanish, while the real parts of 
$\langle p \vert V \vert \Psi_{i0}(0)\rangle$ (dotted line) and 
$\langle p \vert V \Omega_0 (V=0,t=-3) \vert \Psi_{i0}\rangle$
(dashed line), should agree.  The figure shows that 
the calculation accurately approximates the smeared Born approximation.

\begin{figure}
\centering
%\includegraphics[scale=.5]{fig_21.pdf}
%\includegraphics[scale=.5]{fig_21.pdf}
\includegraphics[scale=.5]{figure_12.pdf}
\caption{Plots of the exact smeared Born approximation,
the real and imaginary parts of the Born approximation smeared
against the wave packet
computed using complex probabilities with V set to 0.} 
%\label{figure 21}
\label{figure 12}
%\label{figure 25}
\end{figure}  
The approximate calculations give
$V(x) \langle x \vert \Omega (-t) \vert \psi_{i0}\rangle$.  The
sharp-momentum transition matrix elements can be computed using a
direct Fourier transform (\ref{s.10}) or by integrating against
Gaussian delta function with the desired final momentum (\ref{s.11}).
These two methods of calculation are compared in figure 13.
In these calculations the $\Delta p$ of the final Gaussian is .25,
which is the same value used in the initial Gaussian.  The dashed
curves show the real and imaginary parts of the smeared transition
matrix elements computed using a direct numerical Fourier transform
(dashed curves) compared to the curves which show the corresponding
quantities that replace the final momentum by a Gaussian approximation
to a delta function (dotted curves).  For these calculations the time
step was taken to be $\Delta t =.025$, and $\delta_x= .003$ which is
smaller than the time step used in the calculations shown in figures 4-10.
This is because the $T$-matrix calculations, which involve
Fourier transforms or integration against oscillating wave packets,
are more sensitive to the accuracy of the scattering wave functions.
These figures show that both methods give results within a few percent
of each other.

\begin{figure}
\begin{minipage}[t]{.45\linewidth}
\centering
%\includegraphics[scale=.5]{fig21.pdf}
%\includegraphics[scale=.5]{fig_22.pdf}
\includegraphics[scale=.5]{figure_13a.pdf}
%\caption{Plots of the real part of the $p=5.0$ half-shell T matrix
%computed by numerical Fourier transform and
%by smearing with a delta-function normalized Gaussian wave packet.}
%\caption{Real(T) - Fourier transform vs wave packet} 
%\label{figure 27}
%\label{figure 22}
%\end{figure}  
\end{minipage}
\begin{minipage}[t]{.45\linewidth}
%\begin{figure}
\centering
%\includegraphics[scale=.5]{fig21.pdf}
%\includegraphics[scale=.5]{fig_23.pdf}
\includegraphics[scale=.5]{figure_13b.pdf}
%\caption{Plots of the imaginary part of the $p=5.0$ half-shell T matrix
%computed by numerical Fourier transform and
%by smearing with a delta-function normalized Gaussian wave packet.
%}
%\caption{Imaginary(T) - Fourier transform vs wave packet} 
%\label{figure 28}
%\label{figure 23}
\end{minipage}
\label{figure 13}
\caption{Plots of the real and imaginary parts of the $p=5.0$ half-shell T matrix
computed by numerical Fourier transform and
by smearing with a delta-function normalized Gaussian wave packet.
}
\end{figure}  

Figure 14 compares the real and imaginary parts of the smeared
half-shell transition matrix elements (using the direct Fourier
transform method (dotted curves) to sharp momentum half-shell
transition matrix elements computed by numerically solving the
Lippmann-Schwinger equation (dashed curves).
%(dashed curves) using the method in \cite{Rubtsova}.
As with the other calculation, errors
are a few percent at the on-shell point.  The comparison is between a
sharp-momentum matrix element and a matrix element where the initial
state is smeared with a narrow wave packet, so there will be some
residual difference due to the smearing.  The relative error
for the real and imaginary parts are illustrated in figure 15.
While the relative errors can large, there are a few percent near the
on shell point, $p=5$.  All of the errors can
be reduced by increasing the number of time steps and intervals.
 
\begin{figure}
\begin{minipage}[t]{.45\linewidth}
\centering
%\includegraphics[scale=.5]{2steps_lc_re.pdf}
%\includegraphics[scale=.45]{fig_24.pdf}
\includegraphics[scale=.45]{figure_14a.pdf}
%\caption{Plots of the real part of the
%half-shell transition
%matrix computed using the Fourier transform of the path integral
%calculation and by numerical solution of the Lippmann-Schwinger
%equation ($k=5.$).}
%\caption{Re(T) exact vs path integral}
%\label{figure 20}
%\label{figure 24}
%\end{figure}  
\end{minipage}
\begin{minipage}[t]{.45\linewidth}
%\begin{figure}
\centering
%\includegraphics[scale=.5]{2steps_lc_im.pdf}
%\includegraphics[scale=.45]{fig_25.pdf}
\includegraphics[scale=.45]{figure_14b.pdf}
%\caption{
%Plots of the imaginary part of the
%half-shell transition
%matrix computed using the Fourier transform of the path integral
%calculation and by numerical solution of the Lippmann-Schwinger
%equation ($k=5.$).}
%\caption{Im(T) exact vs. path integral}
%\label{figure 21}
%\label{figure 25}
\end{minipage}
\label{figure 14}
\caption{
Plots of the real and imaginary parts of the
half-shell transition
matrix computed using the Fourier transform of the path integral
calculation and by numerical solution of the Lippmann-Schwinger
equation ($k=5.$).}
\end{figure}

%The relative error of the real and imaginary parts of the
%half shell transition matrix elements are illustrated in figures 26 and 27.
%These plots show the absolute value of the difference between the two curves
%in figures 24 and 25 divided by the maximum magnitude of the half-shell
%transition matrix elements.
%The relative error in both cases peaks at about three percent.
\begin{figure}
\begin{minipage}[t]{.45\linewidth}  
\centering
%\includegraphics[scale=.5]{fig_21.pdf}
%\includegraphics[scale=.5]{plot_ex_r.pdf}
%\includegraphics[scale=.5]{fig_26_r_r.pdf}
%\includegraphics[scale=.5]{fig_26_9_r.pdf}
\includegraphics[scale=.5]{figure_15a.pdf}
% \caption{Plot of the absolute value of the relative error
%of the real part of the half-shell transition matrix
%path integral calculated using path integrals and
%the Lippmann-Schwinger equation where $k=5.0$  is the on-shell point.} 
%  (Re(T) calc - Re(T) exact)/( $\vert$T$\vert$ calc max) $p=5.0$ } 
%\label{figure 26}
%\label{figure 25}
%\end{figure}  
\end{minipage}
\begin{minipage}[t]{.45\linewidth}
%\begin{figure}
\centering
%\includegraphics[scale=.5]{fig_21.pdf}
%\includegraphics[scale=.5]{plot_ex_i.pdf}
%\includegraphics[scale=.5]{fig_27_r_i.pdf
\includegraphics[scale=.5]{figure_15b.pdf}
%\includegraphics[scale=.5]{fig_27_9_i.pdf}
%\caption{Plot of the absolute value of the relative error
%of the imaginary part of the half-shell transition matrix
%path integral calculated using path integrals and
%the Lippmann-Schwinger equation where $k=5.0$ is the on-shell point.}
%  (Im(T) calc - Im(T) exact)/( $\vert$T$\vert$ calc max) $p=5.0$ } 
%\label{figure 27}
% \label{figure 25}
\end{minipage}
\label{figure 15}
\caption{Plot of the absolute value of the relative error
of the real and imaginary parts of the half-shell transition matrix
path integral calculated using path integrals and
the Lippmann-Schwinger equation where $k=5.0$ is the on-shell point.}

\end{figure}  

In this one-dimensional problem there are two on-shell transition
matrix elements corresponding to the transmitted wave and reflected
wave.  The calculations exhibited in the series of figures above are
only for the transmitted wave.  The on-shell matrix element for the
reflected wave can be evaluated by changing the sign of the final
momentum either in the Fourier transform on in the Gaussian
approximation to the delta function.  At this energy the transition
matrix elements for the reflected wave are about five orders of
magnitude smaller that those of the transmitted wave and are not
reliable.

Phase shifts, transmission and reflection coefficients can be
computed from the transition matrix elements.  The results for
$\Delta t=.025$, $\delta_x= .003$ between $-15..0$ and $15.0$,
$\Delta p =.25$ are shown below. 
For the Fourier Transform method
the transition matrix elements are
\[
\langle p \vert T \vert p \rangle = -8.04\times 10^{-01} -i 8.17 \times 10^{-01}  \qquad
\langle -p \vert T \vert p \rangle = - 2.41 \times 10^{-06} + i3.16 \times 10^{-06} .
\]
The phase shift for forward scattering is 
\[
\delta (p) =23.0^{\deg} .
\]
A numerical calculation of the Lippmann-Schwinger equation
gives a phase shift of $\delta (p)= 25.4^{\deg}$
which is slightly larger than the computed phase shift.
The phase shift for the reflected wave involves a ratio of two small
numbers, both of which are smaller than computational errors,
so it was not computed. 

The transmission and reflection coefficient are
\[
T=-.027 + i 1.01  \qquad 
R=  3.03 \times 10^{-6} +i 3.99 \times 10^{-6}
\]
and the unitarity check gives
\[
\vert T\vert ^2+\vert R\vert ^2= 1.02
\]
For the Gaussian approximate delta function 
the transition matrix elements are
\[
\langle p \vert T \vert p \rangle = -7.85e^{-01} -i 7.96e^{-01} \qquad
\langle -p \vert T \vert p \rangle = -2.16e^{-06}  +i 3.35e^{-06} .
\]
and the phase shifts are
\[
\delta(p) = 23.0^{\deg} \qquad \delta(-p) = 0.0^{\deg} 
\]
which agree with the phase shifts calculated using
the Fourier transform method.
The reflection and transmission coefficients are
\[
T=-1.59\times 10^{-04} +i .986 \qquad
R =2.71\times 10^{-06} +i 4.2 \times 10^{-06}\\
\]
and the unitarity check gives
\[
\vert T\vert ^2+\vert R\vert ^2= .972
\]
Note that there are slight differences in the unitary checks for
the two methods of computation.

For the calculations at $p=5$, the momentum was sufficiently high that
the amplitude of the reflected wave was negligible relative to the
transmitted wave.

Two additional calculations were performed at lower momenta ($p=2.5$
and $p=1.0$) that have larger amplitude reflected waves.  For these
calculations the width of the initial wave packet in momentum had to
be sufficiently small so the spreading in position would not outrun
the motion of the center of the wave packet.  Decreasing the momentum
width of the initial wave packet increases the spatial volume of the
wave packet.  This, along with the slower moving wave packet,
increases the time that the particle feels the potential.

The results of these calculations are summarized in table III.  The
calculations in the table used the direct Fourier transform method,
although it gives the same phase shifts as the Gaussian wave packet
method.  The calculations at $p=1.0$ used $15000$ spatial intervals
between $-35.0$ and $35.0$ with a spatial resolution of
$\delta_x=.0047$.  The total time was $t=45.0$ which was broken up
into and $1000$ time steps of size $\Delta t=.045$.  The width of the
initial wave packet was $\Delta p= .1$.

The calculations at $p=2.5$ used $16000$
spatial intervals between $-15$ and $15$ with a spatial resolution of
$\delta_x=.0019$.  The total time was $t=7.0$ which was broken up into
and $1000$ time steps of size $\Delta t= .007$.  The width of the
initial wave packet was $\Delta p= .2$.  The phase shifts and unitarity
checks for the three calculations are shown in table 2

\begin{table}
\caption{Phase shifts and unitarity checks}
\begin{tabular}{llllll}
\hline
\hline
$p$ & $\phi^t$ & $\phi^t_e$ & $\phi^r$ & $\phi^r_e$ & $R^2+T^2$ \\
\hline
1   & -49.1$^{\deg}$ & -44.9$^{\deg}$  & 0.4$^{\deg}$ & 0.1$^{\deg}$ & 0.7\\
2.5 & 82.7.$^{\deg}$ &82.9$^{\deg}$    & -74.8$^{\deg}$  & -52.1$^{\deg}$ & 0.9\\ 
5.  & 23.0$^{\deg}$ & 25.4$^{\deg}$ & 0.0$^{\deg}$ & 0.0$^{\deg}$ & 1.0 \\
\hline
\end{tabular}
\end{table}

The half-shell transition matrix elements for these two momenta are
shown in figures 16-19.  The real and imaginary parts of the
calculated $p=1.0$ half-shell transition matrix elements are compared
to the same quantities by numerically solving the Lippmann-Schwinger
equation in figure 16.  The path integral calculation reproduces the
qualitative features of the solution of the Lippmann Schwinger
equation.  The corresponding relative errors for these two
calculations are shown in figure 17.  The relative errors are between
.1 and .2 near the on-shell point.

These lower energy calculations required larger volumes and more time steps than
the calculation at $p=5.0$, resulting in a much slower convergence.
\begin{figure}
\begin{minipage}[t]{.45\linewidth} 
\centering
%\includegraphics[scale=.5]{2steps_lc_re.pdf}
%\includegraphics[scale=.5]{fig_32_1.pdf}
%\includegraphics[scale=.5]{fig_26.pdf}
\includegraphics[scale=.5]{figure_16a.pdf}
%\caption{Plots of the real part of the
%half-shell transition
%matrix computed using the Fourier transform of the path integral
%calculation and by numerical solution of the Lippmann-Schwinger
%equation ($p=1.0$).}
%\label{figure 30}
%\label{figure 28}
%\end{figure}  
\end{minipage}
\begin{minipage}[t]{.45\linewidth}
%\begin{figure}
\centering
%\includegraphics[scale=.5]{2steps_lc_im.pdf}
%\includegraphics[scale=.5]{fig_33_1.pdf}
%\includegraphics[scale=.5]{fig_27.pdf}
\includegraphics[scale=.5]{figure_16b.pdf}
%\caption{Plots of the imaginary part of the
%half-shell transition
%matrix computed using the Fourier transform of the path integral
%calculation and by numerical solution of the Lippmann-Schwinger
%equation ($k=1.$).}
%\label{figure 21}
%\label{figure 29}
\end{minipage}
\caption{Plots of the real and imaginary parts of the
half-shell transition
matrix computed using the Fourier transform of the path integral
calculation and by numerical solution of the Lippmann-Schwinger
equation ($k=1.$).}
%\label{figure 21}
\label{figure 16}
\end{figure}

\begin{figure}
\begin{minipage}[t]{.45\linewidth}  
\centering
%\includegraphics[scale=.5]{fig_21.pdf}
%\includegraphics[scale=.5]{fig_30_r_r_9.pdf}
\includegraphics[scale=.5]{figure_17a.pdf}
%\caption{Plot of the absolute value of the relative error
%of the real part of the half-shell transition matrix
%path integral calculated using path integrals and
%the Lippmann-Schwinger equation where $k=1.0$  is the on-shell point.
%} 
%\label{figure 30}
%\label{figure 25}
%\end{figure}
\end{minipage}
\begin{minipage}[t]{.45\linewidth}  
%\begin{figure}
\centering
%\includegraphics[scale=.5]{fig_21.pdf}
%\includegraphics[scale=.5]{fig_31_r_i_9.pdf}
\includegraphics[scale=.5]{figure_17b.pdf}
%\caption{Plot of the absolute value of the relative error
%of the imaginary part of the half-shell transition matrix
%path integral calculated using path integrals and
%the Lippmann-Schwinger equation where $k=1.0$ is the on-shell point.
%} 
%\label{figure 31}
% \label{figure 25}
\end{minipage}
\caption{Plot of the absolute value of the relative error
of the real and imaginary parts of the half-shell transition matrix
path integral calculated using path integrals and
the Lippmann-Schwinger equation where $k=1.0$  is the on-shell point.
} 
\label{figure 17}
\end{figure}  

The real and imaginary parts of the half-shell transition matrix
elements for $p=2.5$ computed with the path integral and Lippmann
Schwinger equation are compared in figures 18 and 19.  The real part
of the transition matrix element is almost a factor of 10
smaller than the imaginary part.  In addition, from the
graph, near the on shell point $p=-2.5$ the real and imaginary parts
of the transition operator are near 0 so the calculation of the phase
shift, which for the reflected wave is a function of the ratio of
these two small quantities, is dominated by computational errors.

\begin{figure}
\begin{minipage}[t]{.45\linewidth}  
\centering
%\includegraphics[scale=.5]{2steps_lc_re.pdf}
%\includegraphics[scale=.5]{fig_34_2_5_r.pdf}
%\includegraphics[scale=.5]{fig_28.pdf}
\includegraphics[scale=.5]{figure_18a.pdf}
%\caption{Plots of the real part of the
%half-shell transition
%matrix computed using the Fourier transform of the path integral
%calculation and by numerical solution of the Lippmann-Schwinger
%equation ($k=2.5$).}
%\label{figure 30}
%\label{figure 32}
%\end{figure}
\end{minipage}
\begin{minipage}[t]{.45\linewidth}  
%\begin{figure}
\centering
%\includegraphics[scale=.5]{2steps_lc_im.pdf}
%\includegraphics[scale=.5]{fig_34_2_5_i.pdf}
%\includegraphics[scale=.5]{fig_29.pdf}
\includegraphics[scale=.5]{figure_18b.pdf}
%\caption{Plots of the imaginary part of the
%half-shell transition
%matrix computed using the Fourier transform of the path integral
%calculation and by numerical solution of the Lippmann-Schwinger
%equation ($k=2.5$).}
%\label{figure 21}
%\label{figure 33}
\end{minipage}
\caption{Plots of the real and imaginary parts of the
half-shell transition
matrix computed using the Fourier transform of the path integral
calculation and by numerical solution of the Lippmann-Schwinger
equation ($k=2.5$).}
%\label{figure 21}
\label{figure 18}
\end{figure}

The corresponding relative error in the real and imaginary
parts of the transition matrix at $p=2.5$ are shown in figures 19.
They are about $.03$ at the on shell point.

%\begin{figure}
%\centering
%\includegraphics[scale=.5]{fig_21.pdf}
%\includegraphics[scale=.5]{err_t1r.pdf}
%\caption{(Re(T) calc - Re(T) exact)/( $\vert$T$\vert$ calc max) $p=1.0$ } 
%\label{figure 34}
%\label{figure 25}
%\end{figure}  
%
%\begin{figure}
%\centering
%\includegraphics[scale=.5]{fig_21.pdf}
%\includegraphics[scale=.5]{err_t1i.pdf}
%\caption{(Im(T) calc - Im(T) exact)/($\vert$T$\vert$ calc max) $p=1.0$ } 
%\label{figure 35}
%\label{figure 25}
%\end{figure}  
\begin{figure}
\begin{minipage}[t]{.45\linewidth}  
\centering
%\includegraphics[scale=.5]{fig_21.pdf}
\includegraphics[scale=.5]{figure_19a.pdf}
%\includegraphics[scale=.5]{fig_34_r_r_9.pdf}
%\caption{Plot of the absolute value of the relative error
%of the real part of the half-shell transition matrix
%path integral calculated using path integrals and
%the Lippmann-Schwinger equation where $k=2.5$  is the on-shell point
%} 
%\label{figure 34}
%\label{figure 25}
%\end{figure}
\end{minipage}
\begin{minipage}[t]{.45\linewidth}  
%\begin{figure}
\centering
% \includegraphics[scale=.5]{fig_21.pdf}
\includegraphics[scale=.5]{figure_19b.pdf}
%\includegraphics[scale=.5]{fig_35_r_i_9.pdf}
%\caption{Plot of the absolute value of the relative error
%of the imaginary part of the half-shell transition matrix
%path integral calculated using path integrals and
%the Lippmann-Schwinger equation where $k=5.0$ is the on-shell point.}
%\label{figure 35}
%\label{figure 25}
\end{minipage}
\caption{Plot of the absolute value of the relative error
of the real and imaginary parts of the half-shell transition matrix
path integral calculated using path integrals and
the Lippmann-Schwinger equation where $k=5.0$ is the on-shell point.}
\label{figure 19}
\end{figure}  


%The half shell transition matrix elements are physically relevant
%quantities.  They can be used to calculate the phase shifts,
%cross sections and relfection and transmission coefficients.
%The relative errors for the real and imaginary part of the half
%shell transition matrix elements are shown as functions of the half
%shell momentum for $p=1.0,2.5 and 5.0$.  The plots show the
%difference between the calculated value and exact value divided by the
%maximum of the calculated value.  


%\begin{figure}
%\centering
%\includegraphics[scale=.5]{2steps_lc_re.pdf}
%\includegraphics[scale=.5]{fig_30_5.pdf}
%\caption{Re(T) exact vs path integral(k=5.0)}
%\label{figure 30}
%\label{figure 30}
%\end{figure}  

%\begin{figure}
%\centering
%\includegraphics[scale=.5]{2steps_lc_im.pdf}
%\includegraphics[scale=.5]{fig_31_5.pdf}
%\caption{Im(T) exact vs. path integral(k=5.0)}
%\label{figure 21}
%\label{figure 31}
%\end{figure}

The results indicate that lower energy calculations require more
computational resources than the higher energy calculations. There are
a number of reasons for this.  The slow moving wave packet must outrun
the negative spreading of the initial wave packet.  In addition, the
higher momentum resolution needed for low-energy calculations
increases the relevant volume that needs to be decomposed into
products of intervals.  The exhibited calculations, especially the
lower energy calculations, already involved iterating large matrices
($15000 \times 15000$). Fortunately these matrices do not have to be
stored.  The calculations are approximations in the sense that they
mathematically converge to the exact result.  This means that in
principle they can be improved by using more computational resources
(larger matrices and more time steps).  We did not attempt to go
beyond calculations of this size for $p=1.0$ and $p=2.5$ in part
because with the brute force computational methods used in this work,
small increases in resolution result in much larger calculations.
For this reason convergence is slow and it makes sense to first
investigate methods to improve the efficiency of the calculations,
rather than to continue investigating the limits of these inefficient
methods.  It is clear from the above discussions, that calculations of
zero-energy observables, like scattering lengths, would be very difficult. 
This is because it would be necessary to use narrow wave packets
in momentum that would result require large volumes.


\section{Summary and Conclusions}
 
This work demonstrated both the feasibility and challenges
of performing scattering calculations using real-time path integrals.

The exhibited calculations are based on a reinterpretation of the
path integral as the expectation of a potential functional with
respect to a complex probability distribution on a space of paths.
The space of paths is constructed by dividing the total time into a
large number of small time slices.  At each time slice space is
divided up into a large number of small windows, including two
semi-infinite windows.  The space of paths can then be decomposed into
disjoint equivalence classes where equivalent paths pass thorough the
same intervals at each time slice.
 
A complex probability is assigned to each equivalence class.  It is
constructed by decomposing the propagation of a free quantum
mechanical system into a sum of parts associated with each equivalence
class. Interactions are introduced by considering the effects of the
potential at each window.  Mathematically the time evolution of an
initial wave packet is represented by an average of a path-dependent
potential functional on the space of paths with respect to the complex
probability distribution.  The virtue of this method is that
it is a real approximation in the sense that it converges to
a global solution of the Schr\"odinger equation.

For the convergence of the calculations the equivalence class of paths
must be sufficiently fine that the potential locally constant
on each spatial volume element.    In addition, the time
steps must be sufficiently small for the Trotter product formula to
converge.  The fact that the complex probabilities are largely
supported on continuous paths means that for smooth short-range
potentials the relevant volume is finite. 

One difficulty is the large number of equivalence classes.  Each
equivalence class can be represented as the product of $N$ 
one-step classes for each time step.  This can be used to make an
approximate factorization of the $N$-time step complex probability as
a product of $N$ one-time-step probabilities.  This approximation
replaces the sum over the large number of equivalence classes by
computing powers of a large matrix, which was used in the
calculations.  In this representation the contribution of
each equivalence class is approximated by a particular sequence of
products of matrix elements.  For example, if $M_{mn}$ represents
the one step probability matrix, the product of numbers
\beq
M_{n_0 n_1} M_{n_1 n_2} \cdots M_{n_{19} n_{20}} 
\label{con:1}
\eeq
is approximately the complex probability for passing through the
sequence of windows $I_{n_{20}}$ at time $t_1$, $I_{n_{19}}$ at time
$t_2, \cdots $.  This factorization represents a tremendous increase
in efficiency - by using matrix algebra to treat a large number of
equivalence classes of paths in parallel.  Computationally
this looks like replacing the
path integral by successive applications of a transfer matrix to an
initial state, except in this case it is possible to identify and
extract the contribution of each cylinder set of paths to the
free-particle time-evolution operator. 

The scattering calculations were performed by approximating M{\o}ller
wave operators \cite{Moller} applied to normalizable wave packets.
This has the advantage of removing wave packet spreading effects from
the scattering calculations.  The relevant volume is the region where
the particles interact.  It is determined by the range of the
interactions and size of the free wave packets.  One possible strategy
to improve the efficiency would be to reduce the size of the active
volume by building some of the asymptotic properties of scattering
wave function into the free wave functions. This is essential for
long-range potentials \cite{Dollard} \cite{Mulherin}\cite{Busalev},
but similar methods could be applied to short range potentials.
 
The $p=5.0$ calculations presented in this paper formally involve
applying the $30^{th}$ power of a $5000 \times 5000$ matrix to a fixed
vector (for the transition matrix element calculations the $100^{th}$
power of a $10000 \times 10000$ matrix). This corresponds to averaging
over $5000^{30}$ (resp $(10000)^{100}$) equivalence classes of paths.
Because the one-step probabilities could be computed analytically,
accurately and efficiently, matrix elements could be computed on the
fly, which means that the computer storage required for these
calculations amounted to storing few complex vectors of length
5000(10000).  One of the surprising aspects of these calculations is
the stability of the sums over the complex probabilities.

The $p=5.0$ calculations presented were performed in minutes on a
desktop computer.  The $p=1.0$ and $p=2.5$ calculations required
larger space-time volumes and more time steps and took more than a day
on the same computer.
In terms of the number of
cylinder sets the $p=1.0$ calculations presented used the equivalent of
$(15000)^{1000}$ cylinder sets while the $k=2.5$ calculations used the
equivalent of $(16000)^{1000}$ cylinder sets.  As a result they took
considerably longer.  No attempt was made to be efficient.  All
intervals and time steps were taken to be the same size.  This is the
analog of computing a Riemann integral with equally spaced intervals.
There is a great deal of freedom both in how to choose intervals and
time slices that was not exploited.  Some  strategies that could be
explored would be use more resolution of cylinder sets near classical
paths.  The structure of the potential could be used to determine 
and optimal decomposition of each time slice into intervals.
%Another strategy would be to utilize properties
%of the initial wave packet to choose identify the most important
%cylinder sets. 
The presented calculations simply applied the same matrix $N$ times to
the initial vector.  The application of the one-step probability
matrix to a localized vector could be made more efficient by
discarding small components of the resulting vector, reducing the size
of the vector that must be stored for each time step.  The one step
complex probabilities approximate a unitary transfer matrix for free
evolution, however exact unitarity is not forced on the resulting matrices.
Designing the one step probability so the matrix is exactly unitary
may reduce errors.

Beyond the numerical considerations, this framework is appealing in
that the input is a potential functional $F[\gamma]$; this picture 
is retained both exactly and in approximation.   This is in contrast to the 
usual path integral where the relevant weight functional formally 
looks like an action, but the terms that represent the time derivatives
have no legitimate interpretation as derivatives in the Trotter
product formula.  In the MNJ formulation, these terms do not appear 
explicitly; they are contained in the expression for the 
one-step probabilities

While the calculations in this paper are motivated by the 
complex probability interpretation, the computational strategy 
can be understood directly from Feynman's work.  His path integral 
results in the kernel $K(x,t;x',t')$ of the
time evolution operator (see equation (4.2) of \cite{Feynman_2})
\beq
\langle x \vert \psi (t) \rangle = \int K(x,t;x',t') dx' 
\langle x' \vert \psi (t') \rangle
\label{con:2}
\eeq
This can be expressed as the product of propagation over many 
time steps (see equation
(2.33) of \cite{Feynman_2}) :
\beq
\langle x \vert \psi (t) \rangle = \int K(x,t;x_1,t_1) dx_1
K(x_1,t_1;x_2,t_2) dx_2 \cdots K(x_{N-1},t_{N-1};x_N,t_N) dx_N 
\langle x_N \vert \psi (t_N) \rangle
\label{con:3}
\eeq
If the time intervals $t_{j+1}-t_j$, are sufficiently small then   
\beq
K(x_{j+1},t_{j+1};x_j,t_j) \approx
K_0(x_{j+1},t_{j+1};x_j,t_j)e^{-i V(x_j)(t_{j+1}-t_j)} 
\label{con:4}
\eeq 
where $K_0(x_{j+1},t_{j+1};x_j,t_j)$ is the free time-evolution 
kernel.  This is justified by the Trotter product formula.
Finally, if the integrals were replaced by numerical quadratures,
this would become
\[
\langle x \vert \psi (t) \rangle \approx \sum 
K_0(x,t ;x_{1n_1} ,t_1)e^{-i V(x_{1n_1})(t-t_1)} \Delta_{x_{1n_1}}
K_0(x_{1n_1},t_1 ;x_{2n_2} ,t_2)e^{-i V(x_{2n_2})(t_1-t_2)} \Delta_{x_{2n_2}}
\cdots \times
\]
\beq
K_0(x_{(N-1)n_{N-1}},t_{N-1} ;x_{Nn_N} ,t_N)e^{-i V(x_{Nn_N})(t_{N-1} -t_N)} \Delta_{x_{Nn_N}} \langle x_N \vert \psi (t_N) \rangle
\Delta_{x_{2n_2}}.
\label{con:5}
\eeq
The approximation
\beq
K_0(x_{1n_1},t_1 ;x_{2n_2} ,t_2)\Delta_{x_{2n_2}}
\approx 
\int_{x_{2n_2}-\Delta_{x_{2n_2}}/2}^{x_{2n_2}+\Delta_{x_{2n_2}}/2}
K_0(x_{1n_1},t_1 ;x ,t_2)dx
\label{con:6}
\eeq
corresponds to  the one-step probability  used in this work.   If this
replacement  is made  in  (\ref{con:5}) the  result  is equivalent  to
(\ref{njm.8}).  The important features to emphasize are (1) because of
the small time steps, the potential can be factored out and evaluated
at one of the quadrature points  and (2) the free  particle kernel is
known.  It is these two features that make  non-trivial calculations
possible.

The short time kernel
\[
K_0(x,t;x',t+\Delta t) 
\]
can be computed analytically leading to the following
unitary kernel for the approximate transfer matrix:
\[ 
TM(x,t;x',t+\Delta t) = 
K_0(x,t;x',t+\Delta t)e^{-i V(x')}
\]
\[
P(x'\to x, dx, \Delta t) := dx TM(x,t;x',t+\Delta t) 
\]
can be interpreted as a one-step complex probability
for a particle initially at $x'$ to end up within
$dx$ of $x$ in time $\Delta t$.  A direct numerical calculation
with this kernel may be more efficient that the brute force
calculation illustrated above.

The interesting question is how this method scales to more
particles or fields?  While
it is difficult to answer this question
without performing a test calculation and taking advantage of efficiencies,
%one thing that helps is the trotter approximation.
For small time steps the Trotter formula justifies the replacements
\[
e^{-i (\sum H_i + \sum V_{ij})\Delta t} \to
\prod e^{-i H_i \Delta t} e^{-i \sum V_{ij} \Delta t} \to
\prod e^{-i H_i \Delta t} \prod e^{-iV_{ij} \Delta t} .
\]
In this form there are more products for a given time step, but the
non-trivial products only involve exponentials of single pairwise
potentials or one-body free propagators. Each of these operators,
$e^{-i H_i \Delta t}$ or $e^{-iV_{ij} \Delta t}$, only
act on the variables associated with particle $i$ or particles $i$ and $j$ in
the many-body wave function.  This suggest that the non-trivial part a
few-body calculation scales with the number of interacting pairs.
There are other factors, in that both the relevant volume and
interaction time will increase as more particles are added.  In
addition serious challenges with storing the vector
representing the evolved many-body wave function must be overcome..

The exhibited one-dimensional calculations are of comparable
difficulty to computing three-dimensional two-body scattering
calculations using partial waves in the center of mass frame.
Time-dependent methods based on numerical computation of wave
operators have been successfully applied to treat-three body $n$-$D$
breakup with strong $+$ Coulomb interactions \cite{kroger}.  These
calculation suggest that with some improvements in efficiency that
real path integral scattering calculations involving two or three
particles might be possible.

While this work demonstrates that scattering observables can be
directly computed using real-time path integrals, the calculations
used in this work cannot compete with conventional computational
methods.  These methods may provide a useful laboratory
for investigating the efficiency of different approximation
methods for treating real time evolution, which is needed
for scattering calculations on a quantum computer.

Conceptually, the one-step probability multiplied by the one-step
potential functional, $e^{-i V(y_i)\Delta t}$, which is the key to the
computational method, is essentially a transfer matrix, which is a
unitarized version of using the Hamiltonian to solve the Schr\"odinger
equation for finite time by taking many small time steps.  This has a
lot in common with evolving a product of smeared field operators with
the Heisenberg equations of motion.  Also, the one-step complex
probability is essentially free propagation over short time, which is
also well understood in the field-theory case.
%The connection with
%this interpretation of the real-time path integral with time evolution
%based on a transfer matrix may be valuable to explore in this context.

One of the authors (WP) would like to acknowledge the generous
support of the US Department of Energy, Office of Science, 
grant number DE-SC0016457 who supported this research effort.

\bibliography{henstock_paper.bib}
\end{document}

\begin{figure}
\begin{minipage}[t]{.45\linewidth}
\centering
\includegraphics[angle=000,scale=.25]{polyzou-fig1.eps}
\caption{\bf V(x),$\langle x \vert \psi_{0i}(0) \rangle$}
\label{fig:1}
\end{minipage}
\begin{minipage}[t]{.45\linewidth}
\centering
\includegraphics[angle=000,scale=.25]{polyzou-fig2.eps}
\caption{\bf V(x),$\langle x \vert \psi_{0i}(-3) \rangle$}
\label{fig:2}
\end{minipage}
\end{figure}


\begin{figure}
\centering
%\includegraphics[scale=.5]{fig18.pdf}
\includegraphics[scale=.5]{fig_22.pdf}
\caption{Initial wave packet}
\label{figure 22}
\end{figure}  

To get some indication of the effect the smearing with this Gaussian
wave packet has on the computation of the transition matrix elements, 
matrix elements of the potential are compared to the same quantity with the 
initial momentum smeared against this Gaussian wave packet in 
figure 20,
\beq
\langle p \vert V \vert p_i \rangle
\qquad \mbox{to} \qquad
\int \langle p \vert V \vert p' \rangle dp'
\langle p' \vert \psi_{i}(-2) \rangle
\eeq
The potential matrix elements are, up a multiplicative constant, the
Born approximation to the transition matrix elements.  The on-shell
point corresponds to the minimum.  An error of comparable size can be
expected when comparing to the exact transition operator.  This can be
improved by choosing a narrower wave packet.

%\begin{figure}
%\centering
%\includegraphics[scale=.5]{born_vs_exact.pdf}
%\caption{Born vs exact}
%\label{figure 17}
%\end{figure}  
%
%
%\begin{figure}
%\centering
%\includegraphics[scale=.5]{born_vs_exact_yuk.pdf}
%\caption{Born vs exact(yukawa)}
%\label{figure 18}
%\end{figure}  

%Figure 17 Compares sharp momentum potential matrix elements
%compared to integrating the initial momentum with the Gaussian
%approximate delta functions used in the displayed calculations.
%In the calculations the momentum undertainty in the gaussian
%wave packet was taken to be 1/10 of the on shell momentum.
%The graphs in the figure shows that the smearing leads to a
%slight reduction in the on shell value of the potential;
%A similar size error can be anticipated in the transition 
%matrix elements, however the can be systematically improved by 
%increasing the momentum resolution in the initial wave packet.
%
%

Figure 22 shows the real (dashed line) and imaginary (dotted line) 
parts of the calculated
transition matrix elements plotted against the Born approximation
(solid line). 

\begin{figure}
\centering
%\includegraphics[scale=.5]{fig20.pdf}
\includegraphics[scale=.5]{fig_24.pdf}
\caption{Calculated transition matrix elements vs Born approximation}
\label{figure 24}
\end{figure}  

Figure 23 compares the  exact smeared Born approximation (solid line)
compared to $\langle p \vert V  U_0(2)\vert \psi_{i0}(-2)
\rangle$, where $U_0(2)$ constructed using 
products of one step probabilities. The real part of the last 
expression agrees with the exact result.  The calculated term
has a small imaginary part that vanishes in the exact case.   



%\begin{figure}
%\centering
%\includegraphics[scale=.5]{fig21.pdf}
%\includegraphics[scale=.5]{fig_26.pdf}
%\caption{Comparison to exact calculations} 
%\label{figure 26}
%\end{figure}  





\section{Appendix: Complex probabilities}

This appendix provides a brief introduction to the elements of complex
probabilities as they apply to the formulation of the path integral used 
in this work.

The standard elements of probability theory are (1) a set (2) a set
of measurable subsets and (3) a real positive probability measure
normalized to unity.  Random variables are represented by measurable
functions, and their expectation value is given by the integral over
the set with respect to the probability measure.

A simple example consider the set $\Omega = [0,1]$, the unit interval..  
The subsets are the Borel sets 
generated from open sub intervals by countable unions and complements.
The measure is Lebesgue measure.  Random variables are measurable functions
on the interval and their expectation value is the integral with respect to
Lebesgue measure.

As a result of properties of the Lebesgue integral this system 
satisfies the Kolomogorov axioms of probability theory,   
which require
\beq
P(E) = \int_E  dP(x) > 0 \qquad \mbox{for} \qquad  E \subseteq \Omega
\eeq
\beq
P(\Omega) = \int_{\Omega} dP(x) =1 
\eeq
\beq
P(\cup_{n=0}^{\infty}E_n ) = \sum_{n=0}^\infty P(E_n) .
\qquad E_m \cap E_n = \{ \emptyset \} \qquad m \not=n
\eeq
which gives an abstract characterization of the essential feature of this
example.  

Henstock \cite{Henstock} contemplated an alternate formulation of the
theory of probability, based on a generalization of the Riemann
integral, referred to as the Henstock-Kurzweil integral.  This
replaces the usual $\epsilon-\delta$ formulattion of the Reimann
integral by one where for a given $\epsilon$,  the $\delta$'s depends
on the function and varies with the argument of the function.  It is
designed to integrate functions like $\sin (1/x)$, which one can
imagine integrating by choosing $\delta$'s smaller than the
width of each cycle of the function, and adding the contributions in
an order to account for cancellations.  The Henstock-Kurzweil integral
formalizes this process.  A systematic account can be found in
\cite{Muldowney}\cite{Muldowney2}.  From a computational point of view
it reduces the integral to a limit of Reimann sums, with possible
restrictions on how to choose the points and intervals.

In his formulation of probability theory, when applied to the same
example, the probability function $P(E)$ is not required to be
positive or real.  It still must satisfy the normalization condition
$P(\Omega)=1$.  The Lebesgue integral is replaced by the
Henstock-Kurzweil integral. The condition of countable additivity of
the probabilites of disjoint sets is replaced by the weaker condition
finite additivity.  The analog of the measurable sets are sets
generated from open sub-intervals by finite unions and complements.
The countable additivity in the conventional formulation of
probability theory is a property of the measure.  The important
property of the Henstock integral is that the properly chosen finite
Reimann sums converge to an integral, without needing this additional
property of countable additivity.  When the probability is real and
positive it is equivalent to the conventional formulation of
probability theory.




\section{Implementation}

\section{Test case}


This figure shows the potential and the position of the free
wave packet at $t=-2$.  The relevant observation is that
by $t=-2$ the wave packet is out of the range of the potential.

\begin{figure}
\centering
\includegraphics[scale=.5]{potential_psi2.pdf}
\caption{V(x)}
\label{fig:x}
\end{figure}


\begin{figure}
\centering
\includegraphics[scale=.5]{psi_0_p.pdf}
\caption{V(x)}
\label{fig:x}
\end{figure}

This figure shows the initial wave packet at $t=0$
compared to the scattered wave packet at $t=0$ taking
the time limit to $2$


This figure shows that the scattered wave packet $(t=-2)$
is unchanged in going from 20-40 time steps. 



This figure shows that the scattered wave packets $(t=-2)$
and $t=-2.5$. 

\begin{figure}
\centering
\includegraphics[scale=.5]{omega_20_v_25.pdf}
\caption{V(x)}
\label{fig:x}
\end{figure}  

This graph compares the calculated smeared transition
matrix elements to the smeared Born apprximation.

\begin{figure}
\centering
\includegraphics[scale=.5]{exact_v_ba.pdf}
\caption{V(x)}
\label{fig:x}
\end{figure}  

This figure shows the free wave packet at 0 compared
to evolving the free wave packet at -2 back to 0
using the product of 20 1 step probabilities with
$delta t=.1$.  The fact the curves are essentially
indistinguishible implies that thee product of 1 step
probabilities gives a good approximation to a sum
over $3000^{20}$ cylinder sets.


\begin{figure}
\centering
\includegraphics[scale=.5]{omega0_20.pdf}
\caption{V(x)}
\label{fig:x}
\end{figure}  



These figures show the effect of evaluating the potiential function at the
center or left end of the inteval.

\begin{figure}
\centering
\includegraphics[scale=.5]{relc.pdf}
\caption{V(x)}
\label{fig:x}
\end{figure}  

\begin{figure}
\centering
\includegraphics[scale=.5]{imlc.pdf}
\caption{V(x)}
\label{fig:x}
\end{figure}  





In this section we test this method by using it to compute transition 
matrix elements for a particle of mass $\mu$ scattering from a one-dimensional
Gaussian potential.
The potential is assumed to have the Gaussian form
\beq
V(x) = \lambda e^{-x^2 /r_0^2}.
\label{tc.1}
\eeq
with strength $\lambda$ and range $r_0$.  The potential with strength
$\lambda =1$ and range $r_0=1$ is plotted in figure 1.
figure 6: potential

The one-dimensional problem with this potential can be solved
numerically.  Figure 2 shows plots of the potential,
$\langle k \vert V \vert k'\rangle$ , half-shell K matrix,
$\langle k \vert K(E') \vert k'\rangle$ and the real and imaginary
parts of the half-shell transition matrix
$\langle k \vert T(E') \vert k'\rangle$ as functions of $k$ for 
$k'=5.0$.


We choose a minimal uncertainty wave packet with width $\Delta_p=.1$.
It is plotted on the same scale as the off shell $T$ and potential.




The initial and final wave free particle states are minimal uncertainty states
with mean momenta $p_{0i}$ and $p_{0f}$ and uncertainties $\Delta p_i$
and $\Delta p_f$.  The form of these functions are given in the appendix.


We consider the integral (\ref{f.6}) with $N$ time slices.  In order
to utilize matrix multiplication the same partition is used on each
time slice.  For simplicity the partiton is taken to have $M$ equally
spaced points that end at the left and right endpoints $x=x_l$ and
$x=x_r$
\beq
x[m] = x_l + {m-1 \over M-1}(x_r-x_l) \qquad 1 < m \leq M . 
\label{tc.2}
\eeq
While this could certainly be improved upon, but it is a pretty robust method
provided we have a sufficiently fine partition.

Quadrature points are chosen to be at the mid point of each of these
intervals.  Quadrature points for the seminfinite are intervals are
chosen so all of quadrature points are also equally spaced: 
intervals  
\beq
y[0]= x_l - {1 \over 2(M-1)}  (x_r-x_l)
\label{tc.3}
\eeq
\beq
y[m] = x[m]+ {1 \over 2(M-1)}  (x_r-x_l)   \qquad 1 \leq m \leq M
\label{tc.4}
\eeq
\beq
dy[m] ={1 \over 2(M-1)}  (x_r-x_l) 
\label{tc.5}
\eeq
This choice is crude but robust.  Since the center of the 
initial wave packet
$\langle x \vert \psi_{0i} (t) \rangle$ moves with speed $p_{0i}/\mu$,
it will be outside of the domain for this quadrature at $t$ gets
large.  This does not happen if instead we compute
$\langle x \vert \Omega(t) \vert\psi_{0i} (-0) \rangle$

Figures 1 - 4 illustrate the movement of the center of initial free
minimal uncertainty wave packet, 
$\langle x \vert \psi_{0i}(-t) \rangle$, with time for a particle of
mass $\mu=1$, momentum $k=5.0$ and uncertainty in $k$, $\Delta k =.5$
for $-t=0,-.2,-.4,-.6$ and $-1.0$

Given the parameters  
We choose the following parameters
\begin{center}
\begin{tabular}{|l|l|}
\hline
$r_0$ &1.0\\
$\lambda$& 1.0\\
mass& 1.0 \\
$p_0$ & 5.\\
$\Delta p$ &.5\\
$\Delta x$ &1.0\\
$t$ &.6\\
\hline
\end{tabular}
\end{center}
For these parameters we can estimate that wave functions will be out 
of the scattering volume for $t > 1(1+1)/5 = 4.$  


It is useful to keep the calculations centered near the origin
so we can work with a fixed momentum grip.  To do this  
we define the one step probability matrix
\beq
K_{mn} = P(y_m,I_n) 
\label{tc.5}
\eeq
\beq
U_{n}= e^{-i V(y[n] )\Delta t}
\label{tc.6}
\eeq
and the 
\beq
P(y_m,I_n)^* = \int_{I_{n}} \langle y_m \vert e^{iH_0\Delta t} \vert
y \rangle dy . 
\label{tc.7}
\eeq 

The calculation is performed recursively 
\beq
M^1_{mn}  = \sum_{n}  K_{mn} U_n P^*_{nk}
\label{tc.8}
\eeq
\beq
M^l_{mn}  = \sum_{n}  K_{mr} U_r M^{l-1}_{rs}P^*_{sk}.
\label{tc.9}
\eeq
The final result is applied to the initial time zero 
minimal uncertainty state
\beq
\langle y_m \Vert \Omega (t) \Vert \psi_{i0} \rangle   
\approx \sum_n M^N_{mn} \psi (y[n])
\label{tc.9}
\eeq
The advantage of this method of computing is that at each 
stage of the iteration the wave function remains centered at
the origin, which allows us to work with a single quadrature.
Becasue there will be some spreading the the volume of the 
spatial quadrature needs to be sufficiently large.

The results of the first few iterations with time step $=$ 
show that the wave function remains centered near the origin:


This has to be integrated against the product of the final potential 
and final wave function in order to extract the transition operator.  


We have performed calculations with different numbers of time
steps using the above parameters in (\ref{}) to calculate (\ref{}). 

t=.2 .1 -10 10
tmatttt= 2.284787e-01 + i -2.512350e-02 
     t = 2.281e-01 +i  -2.508e-02 

t=.2 .05 -10 10
tmatttt= 2.172345e-01 + i -2.554324e-02 
     t = 2.169e-01 +    i -2.550e-02 

t=.2 .05 -12 12
tmatttt= 2.197880e-01 + i -2.583910e-02 
     t = 2.194e-01 +   i  -2.580e-02 

t=.2 .05 -14 14
tmatttt= 2.233609e-01 + i -2.634134e-02 
     t = 2.230e-01 +i     -2.630e-02 

t=.2 .05 -16 16
tmatttt= 2.243734e-01 + i -2.641181e-02 
     t = 2.240e-01 +i  -2.637e-02 

t=.2 .025  -10 10
tmatttt= 2.203772e-01 + i -2.582029e-02 
     t = 2.200e-01    +  i-2.578e-02 

t=.3 .1 -10 10
tmatttt= 2.283689e-01 + i -3.018381e-02 
t = 2.280e-01 +i  -3.014e-02 

t=.3 .05 -10 10
tmatttt= 2.081177e-01 + i -2.967450e-02 
     t = 2.078e-01 +    i -2.963e-02 

t=.3 .05 -16 16

t=.3 .025 -10 10
tmatttt= 2.083399e-01 + i -2.883398e-02 
     t = 2.080e-01 +i     -2.879e-02 

t=.4 .1  -10 10
tmatttt= 2.298457e-01 + i -7.060607e-03 
      t= 2.294e-01    +i  -7.048e-03 

t=.4 .05  -10 10
tmatttt= 2.001207e-01 + i -3.034409e-02 
     t = 1.998e-01 +i     -3.029e-02 

t=.4 .05 -16 16
tmatttt= 2.168274e-01 + i -3.286948e-02 
     t = 2.165e-01    +i  -3.282e-02 

t=.4 .025 -10 10
tmatttt= 1.961018e-01 + i -2.824561e-02 
     t = 1.958e-01 +i     -2.820e-02 

t.=5 .05  -10 10
tmatttt= 1.931390e-01 + i -2.879154e-02 
t = 1.928e-01 +i  -2.874e-02 

t=.5 .05 -16 16
tmatttt= 2.126872e-01 + i -3.279319e-02 
     t = 2.123e-01 +i  -3.274e-02 

t= .6 .05  -10 10
tmatttt= 1.870646e-01 + i -2.663089e-02 
t = 1.868e-01 +i  -2.658e-02 

t=.6 .05 -16 16
tmatttt= 2.078539e-01 + i -3.203504e-02 
t = 2.075e-01 +i  -3.199e-02 

t=.7 .05  -10 10
tmatttt= 1.817903e-01 + i -2.506831e-02 
t = 1.815e-01 +i  -2.502e-02 

t=.7 .05 -16 16
tmatttt= 2.024542e-01 + i -3.104170e-02 
t = 2.021e-01 +i  -3.099e-02 

t=.8 .05  -10 10
tmatttt= 1.768226e-01 + i -2.417508e-02 
t = 1.765e-01 +i  -2.413e-02 

t=.8 .05 -16 16

t=.9 .05  -10 10
tmatttt= 1.707988e-01 + i -2.332959e-02 
t = 1.705e-01 +i  -2.329e-02 

t=.9 .05 -16 16

t=1. .05  -10 10
tmatttt= 1.599348e-01 + i -2.190350e-02 
t = 1.597e-01 +i  -2.186e-02 

t=1.0 .05 -16 16

t= 1.1

t=1.2

t=.0 .05 -16 16  .25 width  
2.370117e-01 + i -3.693586e-02 

\begin{center}
\begin{tabular}{|l|l|}
\hline
800& points between -16: 16\\
\hline
%2&  3.617352e-01 + i -3.977637e-02\\ 
2& 2.24e-01 + i -2.64e-02 \\
%3&  3.615612e-01 + i -4.778803e-02\\
3& 2.20e-01 + i -3.13e-02 \\
%4&  3.622517e-01 + i -5.119849e-02\\
4& 2.17e-01 + i -3.29e-02 \\
%5&  3.629555e-01 + i -5.235822e-02\\
5&2.13e-01 + i -3.28e-02 \\
%6&  3.633625e-01 + i -5.262710e-02\\
6&2.08e-01 + i -3.20e-02 \\
%7&  3.635104e-01 + i -5.262668e-02\\
7& 2.02e-01 + i -3.10e-02 \\
%8&  3.635182e-01 + i -5.256206e-02\\
8& 1.97e-01 + i -3.00e-02 \\
%9&  3.634785e-01 + i -5.248025e-02\\
9& 1.91e-01 + i -2.89e-02 \\
\hline
%600 & points between -10: 10\\
%\hline
%9&  3.635550e-01 + i -5.256142e-02\\
\hline
\end{tabular}
\end{center}
While this looks like it is converging, it does not
agree with the exact result.

Decreasing the width of the gaussian gives numbers closer to the 
expected result (800 points on -16:16, with delta t=.05 and t=.9  
delta - = .1) 
\beq
t= 2.51e-01 + i -5.98e-02
\eeq
for $t=.8$: 
\beq
t= 2.61e-01 + i -5.74e-02
\eeq
for $1=1.0$:
\beq
t= 2.48e-01 + i -6.11e-02 
\eeq
A direct numerical calculation gives for 
$k'=4.99$,  
$\langle k' \Vert T(E'+i0^+) \Vert k' \rangle =2.80\times 10^{-01} - 
i 5.10\times 10^{-02}$ and 
and for $k'= 5.01$  
$\langle k' \Vert T(E'+i0^+) \Vert k' \rangle =2.80\times 10^{-01} - 
i 5.08\times 10^{-02}$

The potential, $K$-matrix and the real and imaginary parts of the 
half shell transition matrix are plotted for $k'=5.0$ in figure 8.



%
%exact 4.99 -2:10 .025 480   2.800891e-01 -5.102973e-02
%      5.01 -2:10 .025 480   2.801041e-01 -5.081778e-02
%      5.01 -10:10 .025 800  2.801041e-01 -5.081778e-02


\section{appendix- minimal uncertainty wave functions} 

The unit normalized minimal uncertainty state with mean momentum $p_0
:= \langle p \rangle$ 
and momentum uncertainty $\Delta p := \langle (p-p_0)^2 \rangle^{1/2}$ 
has the momentum-space wave function 
\beq
\langle p \vert \psi_0 (0) \rangle = (2 \pi)^{-1/4} (\Delta p)^{-1/2}
e^{- {(p-p_0)^2 \over 4 (\Delta p)^2}}.
\label{ap.1}
\eeq
The uncertainties in position and momentum are related by  
\beq
\Delta p \Delta x = {1 \over 2}
\label{ap.2}
\eeq
The Fourier transform of this state is 
\beq  
\langle x \vert \psi_0 (0) \rangle = (2 \pi)^{-1/4} (\Delta x)^{-1/2}
e^{ - {x^2 \over 4 (\Delta x)^2} + i p_0 x }
\label{ap.3}
\eeq
which has mean position $\langle x \rangle =0$.  

The corresponding time evolved states are
\beq
\langle p \vert \psi (t) \rangle = (2 \pi)^{-1/4} (\Delta p)^{-1/2}
e^{- {(p-p_0)^2 \over 4 (\Delta p)^2}- i {p^2 \over 2 \mu} t}
\label{ap.4}
\eeq
and 
\[
\langle x \vert \psi_0 (t) \rangle =
(2 \pi)^{-1/4}
( 2\Delta p)^{1/2}\sqrt{{1  \over 1  + i {2 (\Delta p)^2 t \over  \mu}}}
\times
\]
\beq
e^{-
{(\Delta p)^2 \over 1  + {4 (\Delta p)^4 t^2 \over  \mu^2}}
(x-{p_0t\over \mu})^2}
e^{i{1 \over 1  + {4 (\Delta p)^4 t^2 \over  \mu^2}}
(x p_0   
+ 2 x^2 (\Delta p)^4 {t \over  \mu}
- { t \over  \mu}{p_0^2 \over 2 })} =
\label{ap.5}
\eeq
\[
(2 \pi)^{-1/4}
( \Delta x)^{-1/2}\sqrt{{1  \over 1  + i { t \over 2 (\Delta x)^2  \mu}}}
\times
\]
\beq
e^{-{1 \over 4 \Delta x^2}
{1 \over 1  + { t^2 \over 4 \Delta x^4   \mu^2}}
(x-{p_0t\over \mu})^2}
e^{i{1 \over 1  + {t^2 \over 4 \Delta x^4  \mu^2}}
(x p_0   
+ { x^2 t \over  8 \Delta x^4 \mu}
- { t \over  \mu}{p_0^2 \over 2 })}
\label{ap.6}
\eeq

The integral of the product of the final time $0$ wave function with the
Gaussian potential used in section X 
over the finite interval, $I_{m}$, that appears
in equation (\ref{}) is 
\beq
\chi_m = \int_{I_m}  V(x) \langle x \vert \psi_{f0}(0)\rangle dx =
\label{ap.7}
\eeq
\beq
(2 \pi)^{-1/4} (\Delta x)^{-1/2} \int_{I_m} 
e^{ - {x^2 \over 4 (\Delta x)^2} + i p_0 x }
\lambda e^{- (x/r_0)^2 } dx =
\label{ap.8}
\eeq
\beq
(2 \pi)^{-1/4} (2 \Delta p)^{1/2} \int_{I_m} 
e^{ - (\Delta p)^2 x^2  + i p_0 x }
\lambda e^{- (x/r_0)^2 } dx =
\label{ap.9}
\eeq
\beq
{1 \over 2} \sqrt{{\pi \over \alpha}} e^{-p_0^2 \over 4 \alpha}
\left (
\mbox{erf}(\sqrt{\alpha}(x_{m+1}- i {p_0 \over 2 \alpha})-
\mbox{erf}(\sqrt{\alpha}(x_{m}- i {p_0 \over 2 \alpha})
\right )
\label{ap.10}
\eeq
where
\beq
\alpha = {1 \over 4 \Delta x^2 } + {1 \over r_0^2}=
\Delta p^2  + {1 \over r_0^2}
\label{ap.11}
\eeq
The integrals needed to extract the sharp momentum transition matrix
elements from $\langle \psi_{f0}(0) \Vert V \Vert \Omega(t) \vert
\psi_{i0}(0) \rangle$ in (\ref{s.19}) are 
\beq
\int \langle p \vert \psi_0 (0) \rangle dp = (2 \pi)^{-1/4} (\Delta p)^{-1/2}
e^{- {(p-p_0)^2 \over 4 (\Delta p)^2}} dp =
(2\pi)^{1/4} (2 \Delta p)^{1/2}   
\label{ap.12}
\eeq

\end{document}

Thus for the purpos
Equation (\ref{tr.16}) can be expressed in matrix form using
\beq
U_n= e^{-i V(v[n] )\Delta t}
\eeq
\beq
\psi_n= \psi (v[n])
\label{tr.30}
\eeq
\[
\langle v[n]  \vert e^{-i H_0 \Delta t} \vert \psi \rangle  \approx
\]
\beq
\sum_{k=0}^M K_{nk}U_k \psi_k  
\label{tr.31}
\eeq
For N time steps
\beq
\langle v[n]  \vert e^{-i H_0 \Delta t} \vert \psi \rangle  \approx
\sum_{k=0}^M (KU)^N_{nk} \psi_k  
\label{tr.32}
\eeq
or
\beq
\langle x \vert e^{-i H_0 N\Delta t} \vert \psi \rangle  \approx 
\sum_{k,n=0}^M (KU)_{\alpha(x)n}(KU)^{N-1}_{nk} \psi_k  
\label{tr.33}
\eeq
The advanatage of this is that if $N=2^L$ then we can compute
\beq
(KU)^2 = (KU)(KU)
\cdots
(KU)^{2^m} = (KU)^{2^{m-1}}(KU)^{2^{m-1}}
\label{tr.34}
\eeq
which allow us to treat multiple time steps geometrically.

Because we still have to deal with an initial time evolved
wave packet.  Since it can be expressed in terms of the 
inverse free propagation
\beq
\psi_k (-t) \approx (K^{-1})N psi_k (0)
\eeq
we can make calcellations at each step:
\beq
(KU)^M \psi_k(-t) = 
KU ( KU ( \cdots KUK^{_1})K^{-1} \cdots K^{-1}) \psi_k(0)
\eeq
It is useful to numerically test the unitarity of $K$


The multistep evolution can be represented using the complex probabilities
\[
\langle x \vert e^{-iH t} \vert \Psi \rangle
\approx
\sum_{k_1=1}^{M-1} \cdots \sum_{k_N=1}^{M-1}
P[\sqrt{{m \over 2\Delta t}}(u[k_N+1]-x), \sqrt{{m \over 2\Delta t}}(u[k_N]-x)]
e^{-i V(v[k_N]) \Delta t} \times
\cdots
\]
\[
P[\sqrt{{m \over 2\Delta t}}(u[k_2+1]-v[k_3]), 
\sqrt{{m \over 2\Delta t}}(u[k_2]-v[k_3])]
e^{-i V(v[k_2]) \Delta t} \times
\]
\beq
P[\sqrt{{m \over 2\Delta t}}(u[k_1+1]-v[k_2]), 
\sqrt{{m \over 2\Delta t}}(u[k_1]-v[k_2])]
e^{-i V(v[k_1]) \Delta t}\Psi (v[k_1]) 
\label{tr.34.1}
\eeq



moves, while  



The results of the previous section show that if we are given
a function $f(x)$, we can express 
\beq
\langle x \vert e^{-iH_0 \Delta t}e^{-iV \Delta t} \vert \psi \rangle 
\label{tr.1}
\eeq
in terms of the expectation values of a random variable with the 
complex probability distribution.  The input is the interaction 
and the initial wave packet
\beq
V(x)  \qquad \langle x \vert \psi \rangle .
\label{tr.2}
\eeq
The treatment of a single time step can be done using the Fresnel
integrals.  The data needed to complete this calculation are the
partition of the real line and the values of $V(x)$ and $\langle x
\vert \psi \rangle$ at a finite collection of points.  The result is a
new function of $x$.  We can repeat this process by applying $e^{-iH_0
\Delta t}e^{-iV \Delta t}$ again to the evolved function.

This is much simpler than summing over cylinder sets.  This is
because the evolution through $k$ time steps is a $k$-fold
matrix multiplication.  In the case that all of the matrices
are identical, high powers can be computed geometrically.

 
The starting point is to choose a set of $M+1$ equally 
spaced intervals on
$[-1,1]$.  The finite boundaries of each interval are 
\beq
x[m] = -1 + {2m \over M+1}  \qquad 1 \leq m \leq M
\label{tr.3}
\eeq
and the mid-points of each interval are
\beq
y[m] = -1 + {2m+1 \over M+1} \qquad 0 \leq m\leq M
\label{tr.4}
\eeq
The next step is to construct a mapping to the 
real line.  This could be in principle be different for 
each time step.  We leave them the same and use
a tangent map for the purpose of illustration:
\beq
u[m] = s \tan ({\pi \over 2}x[m]) \qquad 1 \leq m \leq M
\label{tr.5}
\eeq
\beq
v[m] = s \tan ({\pi \over 2}y[m]) \qquad 0 \leq m \leq M
\label{tr.6}
\eeq
where $s$ is a scale factor.  
The $u[m]$ are the transformed boundary points and the
$v[m]$ are the transformed quadrature points.
We assume $N$ time slices so a single time step
is given by
\beq
\Delta t = (t_f - t_i)/N .
\label{tr.7}
\eeq
The points $v[m]$ play the same role a quadrature points.

For a single time step the unitary evolution is approximately
\beq
\langle x \vert e^{-i H \Delta t} \vert \psi \rangle \approx 
\label{tr.8}
\eeq
\beq
\int \langle x \vert e^{-i {p^2 \over 2m}\Delta t} \vert
y \rangle dy e^{-i V(y)\Delta t} \psi (y) =
\label{tr.9}
\eeq
\beq
\int {1 \over 2\pi}e^{i p (x-y) -i {p^2 \over 2m}\Delta t}
dp dy e^{-i V(y)\Delta t} \psi (y) .
\label{tr.10}
\eeq
Next we complete the square in the exponent
\beq
\int {1 \over 2\pi} dp dy 
e^{
- (i {\Delta t \over 2m})(p- {m \over \Delta t} (x-y))^2 +
i ({\Delta t \over 2m}){m^2 \over \Delta t^2} (x-y)^2}
 e^{-i V(y)\Delta t} \psi (y) =
\label{tr.11}
\eeq
\beq
\int {1 \over 2\pi} dp dy
e^{
- (i {\Delta t \over 2m})(p- {m \over \Delta t} (x-y))^2 +
i {m \over 2 \Delta t} (x-y)^2}
e^{-i V(y)\Delta t} \psi (y)
\label{tr.12}
\eeq
We do the $p$ integral exactly to get 
\beq
\int dy
\sqrt{{m  \over 2 \pi i \Delta t}} 
e^{ 
i {m \over 2\Delta t} (x-y)^2}
e^{-i V(y)\Delta t} \psi (y) dy
\label{tr.13}
\eeq

The next step is to replace the $y$ integral by a sum of integrals over
intervals
\[
\int dy
\sqrt{{m  \over 2 \pi i \Delta t}} 
e^{ 
i {m \over 2\Delta t} (x-y)^2}
e^{-i V(y)\Delta t} \psi (y)  =
\]
\beq
\sum_{k=0}^M \int_{u[k]}^{u[k+1]} dy
\int dy
\sqrt{{m  \over 2 \pi i \Delta t}} 
e^{ 
i {m \over 2\Delta t} (x-y)^2}
e^{-i V(y)\Delta t} \psi (y) dy 
\label{tr.14}
\eeq
For sufficiently small intervals we can approximate
$e^{-i V(y)\Delta t} \psi (y)$ by 
$e^{-i V(v[k])\Delta t} \psi (v[k])$
on the
$k$-th interval
\beq
\sum_{k=0}^M \int_{u[k]}^{u[k+1]} dy
\int dy
\sqrt{{m  \over 2 \pi i \Delta t}} 
e^{ 
i {m \over 2\Delta t} (x-y)^2}
e^{-i V(y[k])\Delta t} \psi (y[k])  .
\label{tr.15}
\eeq
What remains are Fresnel integrals with a shifted limits of integration:
\beq
\sum_{k=0}^M \int_{u[k]-x}^{u[k+1]-x} dy
\sqrt{{m  \over 2 \pi i \Delta t}} 
e^{ 
i {m \over 2\Delta t} y^2} dy
e^{-i V(y[k])\Delta t} \psi (y[k])  
\label{tr.16}
\eeq
The Fresnel integrals can be done analytically.
To compute these integrals we change variables
\beq
\beta = \sqrt{{m \over 2\Delta t}}(y-x)
\qquad 
dy = \sqrt{{2\Delta t \over m }}d\beta 
\label{tr.17}
\eeq
With these changes
\[
\langle x \vert e^{-i H_0 \Delta t} \vert \psi \rangle =
\]
\[
\sum_{k=0}^M \int_{\sqrt{{m \over 2\Delta t}}(u[k]-x)}^{\sqrt{{m \over 2\Delta t}}(u[k+1]-x)}
\sqrt{{2\Delta t \over m }}d\beta 
\sqrt{{m  \over 2 \pi i \Delta t}} 
e^{ 
i \beta^2}
e^{-i V( v[k] )\Delta t} \psi (v[k])
=
\]
\beq
\sqrt{{1  \over i \pi }}
\sum_{k=0}^M \int_{\sqrt{{m \over 2\Delta t}}(u[k]-x)}^{\sqrt{{m \over 2\Delta t}}(u[k+1]-x)}
e^{ 
i \beta^2} d\beta
e^{-i V( v[k] )\Delta t} \psi (v[k])
\label{tr.18}
\eeq

In the language of complex probabilities the one step transfer matrix 
evolution can be expressed as
\beq
\boxed{
\langle x \vert e^{-iH_0 \Delta t} \vert \Psi \rangle
\approx
\sum_{k=0}^M
P[\sqrt{{m \over 2\Delta t}}(u[k+1]-x), \sqrt{{m \over 2\Delta t}}(u[k]-x)]
\Psi (v[k])
} 
\label{tr.18.1}
\eeq
and
\beq
\boxed{
\langle x \vert e^{-iH \Delta t} \vert \Psi \rangle
\approx
\sum_{k=0}^M
P[\sqrt{{m \over 2\Delta t}}(u[k+1]-x), \sqrt{{m \over 2\Delta t}}(u[k]-x)]
e^{-i V(v[k]) \Delta t}\Psi (v[k]) 
}
\label{tr.18.2}
\eeq
If we replace $v[0]=-\infty$ and $v[M]=+\infty$ then we can remove the 
first and last term from the sum - provided the inital wave packet 
vanishes at $\pm \infty$. 


If we define 
\beq
w[k+1](x):=\sqrt{{m \over 2\Delta t}}(u[k+1]-x)
\label{tr.19}
\eeq
and
\beq
w[k](x):=\sqrt{{m \over 2\Delta t}}(u[k]-x)
\label{tr.20}
\eeq
the integral becomes
\[
\sqrt{{1 \over 2i}}
[C_c(w[k+1](x)) - C_c(w[k](x)) 
+i (S_c (w[k+1](x)) - S_c(w[k](x))] =
\]
\[
{1 \over 2}[ (C_c(w[k+1](x))-C_c(w[k](x))+
 S_c(w[k+1](x))-S_c(w[k](x)))
\]
\beq
+i (S_c(w[k+1](x))-S_c(w[k](x)) -C_c(w[k+1](x))+C_c(w[k](x)))]
\label{tr.21}
\eeq
for $0<k<M$.  For $k=0$ it
is replaced by
\[
\sqrt{{1  \over  2 i }}
[{ 1+i \over 2} +  C_c ( w[1](x) ))
+ i S_c (w[1](x)))
] =
\]
\beq
{1 \over 2}[1+C_c(w[1](x))+
S_c(w[1](x)))
+i (S_c(w[1](x)) -C_c(w[k+1](x))]
\label{tr.22}
\eeq
and for $k=M$ by
\[
\sqrt{{1  \over  2 i }}
[{ 1+i \over 2} 
- C_c (w[M](x))
- i S_c  (w[M](x))     
] =
\]
\beq
{1 \over 2}[1-C_c(w[M](x))-
S_c(w[M](x)))
-i (S_c(w[M](x)) -C_c(w[M](x))]
\label{tr.23}
\eeq
where we have used 
\beq
\sqrt{{1  \over  2 i }} = {1 \over 2}( 1 -i)  
\label{tr.24}
\eeq
which leads to 
\[
\langle x \vert e^{-i H_0 \Delta t} \vert \psi \rangle  \approx
\]
\[ 
{1 \over 2}[1+C_c(w[1](x))+
S_c(w[1](x)))
\]
\[
+i (S_c(w[1](x)) -C_c(w[k+1](x))] e^{-i V(v[0] )\Delta t} \psi 
(
v[0]) +  
\]
\[
\sum_{k=1}^{M-1} 
{1 \over 2}[ (C_c(w[k+1](x))-C_c(w[k](x))+
S_c(w[k+1](x))-S_c(w[k](x)))
\]
\[
+i (S_c(w[k+1](x))-S_c(w[k](x)) -C_c(w[k+1](x))+C_c(w[k](x)))]
e^{-i V(v[k] )\Delta t} \psi 
(v[k]) +  
\]
\[
{1 \over 2}[1-C_c(w[M](x))-
S_c(w[M](x)))
\]
\beq
-i (S_c(w[M](x)) -C_c(w[M](x))]
e^{-i V(v[M] )\Delta t} \psi 
(v[M])
\label{tr.25}
\eeq
This gives the value of the time evolved state in terms of the
wave function and interaction on the grid 
$v[k]$.  The result is a
new function that we can time evolve to the next step.  This requires
evaluating the above at
\beq
x=  v[k]
\qquad
\label{tr.26}
\eeq
This leads to the matrix
\[
K_{nk} = 
{1 \over 2}[ (C_c(\sqrt{{m \over 2\Delta t}}(u[k+1]-v[n]))-C_c(
\sqrt{{m \over 2\Delta t}}(u[k]-v[n]))+
\]
\[
S_c(\sqrt{{m \over 2\Delta t}}(u[k+1]-v[n]))-S_c(\sqrt{{m \over 2\Delta t}}(u[k]-v[n])))
\]
\[
+i (S_c(\sqrt{{m \over 2\Delta t}}(u[k+1]-v[n]))-S_c(\sqrt{{m \over 2\Delta t}}(u[k]-v[n]))
\]
\beq
 -C_c(\sqrt{{m \over 2\Delta t}}(u[k+1]-v[n]))
+C_c(\sqrt{{m \over 2\Delta t}}(u[k]-v[n])))]
\label{tr.27}
\eeq
for $0 < k < M$
\[
K_{n0} =
{1 \over 2}[1+C_c(\sqrt{{m \over 2\Delta t}}(u[1]-v[n]))+
S_c(\sqrt{{m \over 2\Delta t}}(u[1]-v[n])))
\]
\beq
+i (S_c(\sqrt{{m \over 2\Delta t}}(u[1]-v[n])) -C_c(\sqrt{{m \over 2\Delta t}}
(u[1]-v[n]))]
\label{tr.28}
\eeq
for $k=0$ and 
\[
K_{nM}={1 \over 2}[1-C_c(\sqrt{{m \over 2\Delta t}}(u[M]-v[n]))-
S_c(\sqrt{{m \over 2\Delta t}}(u[M]-v[n])))
\]
\beq
-i (S_c(\sqrt{{m \over 2\Delta t}}(u[M]-v[n])) -C_c(\sqrt{{m \over 2\Delta t}}(u[M]-v[n]))]
\label{tr.29}
\eeq
for $k=M$
where we have replaced $x$ by $v[n]$.

Equation (\ref{tr.16}) can be expressed in matrix form using
\beq
U_n= e^{-i V(v[n] )\Delta t}
\eeq
\beq
\psi_n= \psi (v[n])
\label{tr.30}
\eeq
\[
\langle v[n]  \vert e^{-i H_0 \Delta t} \vert \psi \rangle  \approx
\]
\beq
\sum_{k=0}^M K_{nk}U_k \psi_k  
\label{tr.31}
\eeq
For N time steps
\beq
\langle v[n]  \vert e^{-i H_0 \Delta t} \vert \psi \rangle  \approx
\sum_{k=0}^M (KU)^N_{nk} \psi_k  
\label{tr.32}
\eeq
or
\beq
\langle x \vert e^{-i H_0 N\Delta t} \vert \psi \rangle  \approx 
\sum_{k,n=0}^M (KU)_{\alpha(x)n}(KU)^{N-1}_{nk} \psi_k  
\label{tr.33}
\eeq
The advanatage of this is that if $N=2^L$ then we can compute
\beq
(KU)^2 = (KU)(KU)
\cdots
(KU)^{2^m} = (KU)^{2^{m-1}}(KU)^{2^{m-1}}
\label{tr.34}
\eeq
which allow us to treat multiple time steps geometrically.

Because we still have to deal with an initial time evolved
wave packet.  Since it can be expressed in terms of the 
inverse free propagation
\beq
\psi_k (-t) \approx (K^{-1})N psi_k (0)
\eeq
we can make calcellations at each step:
\beq
(KU)^M \psi_k(-t) = 
KU ( KU ( \cdots KUK^{_1})K^{-1} \cdots K^{-1}) \psi_k(0)
\eeq
It is useful to numerically test the unitarity of $K$


The multistep evolution can be represented using the complex probabilities
\[
\langle x \vert e^{-iH t} \vert \Psi \rangle
\approx
\sum_{k_1=1}^{M-1} \cdots \sum_{k_N=1}^{M-1}
P[\sqrt{{m \over 2\Delta t}}(u[k_N+1]-x), \sqrt{{m \over 2\Delta t}}(u[k_N]-x)]
e^{-i V(v[k_N]) \Delta t} \times
\cdots
\]
\[
P[\sqrt{{m \over 2\Delta t}}(u[k_2+1]-v[k_3]), 
\sqrt{{m \over 2\Delta t}}(u[k_2]-v[k_3])]
e^{-i V(v[k_2]) \Delta t} \times
\]
\beq
P[\sqrt{{m \over 2\Delta t}}(u[k_1+1]-v[k_2]), 
\sqrt{{m \over 2\Delta t}}(u[k_1]-v[k_2])]
e^{-i V(v[k_1]) \Delta t}\Psi (v[k_1]) 
\label{tr.34.1}
\eeq


The difficulty with implementation is that one scale is set by the
initial wave packet - which means that the distribution of 
$v[k]$ points should be chosen base on the support of the
initial wave function.  This fixes distributions of the 
points $u[k]$ and $v[k]$.  The arguments of the 
Fresnel integrals are  $\sqrt{{m \over 2\Delta t}}(u[M]-v[n])$.
These variable get spread apart as the time step vanishes.  
For the integrals to be accurate we need both both distributions 
to have the property that the kernel and wave function at
adjacent points do not vary too quickly.  It gets harder to 
meet both of these conditions as $\Delta t \to 0$. 

To test this we compare to the exact evolution for different numbers of
time steps.  We turn off the potential and use a Gaussian
wave function for this test.

If we turn off the potential the exact evolution of a gaussian wave packet 
is
\beq
\langle x \vert \psi(t) \rangle = 
\int \langle x \vert e^{-i(p^2/2m)t} \vert p \rangle \pi^{-1/4} e^{-p^2/2} dp =
\eeq
\beq
\pi^{-1/4} {1 \over \sqrt{1+it/m}}e^{- {x^2 \over 2(1+it/m)}} 
\label{tr.35}
\eeq
where
\beq
{1 \over \sqrt{1+it/m}} =
\sqrt{{\sqrt{1+(t/m)^2} +1 \over 2(1+(t/m)^2)}}
-i\sqrt{{\sqrt{1+(t/m)^2} -1 \over 2(1+(t/m)^2)}}
\label{tr.36}
\eeq
Putting everything together gives
\[
\pi^{-1/4} (\sqrt{{\sqrt{1+(t/m)^2} +1 \over 2(1+(t/m)^2)}}
-i\sqrt{{\sqrt{1+(t/m)^2} -1 \over 2(1+(t/m)^2)}})
%{1 \over \sqrt{1+it/m}}
e^{- {x^2 \over 2(1+it/m)}} 
\]
\[
\pi^{-1/4} (\sqrt{{\sqrt{1+(t/m)^2} +1 \over 2(1+(t/m)^2)}}
-i\sqrt{{\sqrt{1+(t/m)^2} -1 \over 2(1+(t/m)^2)}})
%{1 \over \sqrt{1+it/m}}
e^{- {x^2 \over 2(1+(t/m)^2)}(1-it/m)} 
\]
\[
\pi^{-1/4} (\sqrt{{\sqrt{1+(t/m)^2} +1 \over 2(1+(t/m)^2)}}
-i\sqrt{{\sqrt{1+(t/m)^2} -1 \over 2(1+(t/m)^2)}})
%{1 \over \sqrt{1+it/m}}
e^{- {x^2 \over 2(1+(t/m)^2)}} \times
\]
\beq
(\cos ({t \over m}{x^2 \over 2(1+(t/m)^2)}  )
+ i \sin ({t \over m} {x^2 \over 2(1+(t/m)^2)}))=
\label{tr.39}
\eeq
\[
\pi^{-1/4}e^{- {x^2 \over 2(1+(t/m)^2)}} \times 
\]
\[
([\sqrt{{\sqrt{1+(t/m)^2} +1 \over 2(1+(t/m)^2)}}
\cos ({t \over m}{x^2 \over 2(1+(t/m)^2)}  )
+
\sqrt{{\sqrt{1+(t/m)^2} -1 \over 2(1+(t/m)^2)}}
\sin ({t \over m} {x^2 \over 2(1+(t/m)^2)})]
\]
\beq
-i
[
\sqrt{{\sqrt{1+(t/m)^2} -1 \over 2(1+(t/m)^2)}}
\cos ({t \over m}{x^2 \over 2(1+(t/m)^2)}  )
-\sqrt{{\sqrt{1+(t/m)^2} +1 \over 2(1+(t/m)^2)}}
\sin ({t \over m} {x^2 \over 2(1+(t/m)^2)})])
\label{tr.40}
\eeq
We compare the values at
\beq 
x[k]  = \sqrt{2 \Delta t \over m} v[k] 
\label{tr.41}
\eeq
for different numbers of time steps.

Thus we need to compare
\beq
\langle \sqrt{2 \Delta t \over m} v[k] \vert
\psi (n\Delta t) \rangle 
\label{tr.42}
\eeq
with
\beq
\sum_{l=0}^M (K)^n_{k,l} \psi_l
\label{tr.43}
\eeq


\section{plots of relevant quantities}

We choose the following parameters
\begin{center}
\begin{tabular}{|l|l|}
\hline
$r_0$ &1.0\\
$\lambda$& 1.0\\
mass& 1.0 \\
$p_0$ & 5.\\
$\Delta p$ &.5\\
$Delta x$ &1.0\\
$t$ &.6\\
\hline
\end{tabular}
\end{center}

ex code  4.924210e-01 -1.424500e+00 

deltap =.5  p = 5

\begin{center}
\begin{tabular}{|l|l|}
\hline
500& points between -10: 10\\
\hline
2&  3.617352e-01 + i -3.977637e-02\\ 
3&  3.615612e-01 + i -4.778803e-02\\
4&  3.622517e-01 + i -5.119849e-02\\
5&  3.629555e-01 + i -5.235822e-02\\
6&  3.633625e-01 + i -5.262710e-02\\
7&  3.635104e-01 + i -5.262668e-02\\
8&  3.635182e-01 + i -5.256206e-02\\
9&  3.634785e-01 + i -5.248025e-02\\
\hline
600 & points between -10: 10\\
\hline
9&  3.635550e-01 + i -5.256142e-02\\
\hline
\end{tabular}
\end{center}


exact 4.99 -2:10 .025 480  2.800891e-01 -5.102973e-02
      5.01 -2:10 .025 480  2.801041e-01 -5.081778e-02

figure 8 exact transition matrix elements

\begin{figure}
\centering
\includegraphics[scale=.7]{figure_8.pdf}
\caption{V(k,k'),K(k,k') and T(k,k',k')}
\label{fig:8}
\end{figure}  



figure 7 final chi (I)

\begin{figure}
\centering
\includegraphics[scale=.7]{figure_7.pdf}
\caption{$\int_{I_n} \psi(x) V(x) dx$}
\label{fig:7}
\end{figure}  



In figure 8 we take the on shell point to be $k'=5.1$.  The 
value f the transition matrix at this point is 
$T(k,k,k) = .4924210 - i 1.424500$ 


t = 1.781e-01 +i  -1.427e-02  -5 5 N=10 dt .1
t = 1.362e-01 +i   3.289e-02       N=12
t = 1.899e-01 +i  -2.200e-02       N=20 dt .05
xxxxxxxxxxx

\[
\langle x \vert \psi (t)  \rangle =
\]
\[
\int (2 \pi)^{-1/2} e^{ipx - i{p^2 \over 2 \mu}t}
(2 \pi)^{-1/4} (\Delta p)^{-1/2}
e^{- {(p-p_0)^2 \over 4 (\Delta p)^2}} dp =
\]
\[
\int (2 \pi)^{-1/2} (2 \pi)^{-1/4} (\Delta p)^{-1/2}
e^{ipx - i{p^2 \over 2 \mu}t - {(p-p_0)^2 \over 4 (\Delta p)^2}} dp =
\]
\[
(2 \pi)^{-1/2} (2 \pi)^{-1/4} (\Delta p)^{-1/2}
\int
e^{
-p^2 {1 \over 4 (\Delta p)^2} (1  + i {2 (\Delta p)^2 t \over  \mu} )+
p ({p_0 \over 2  (\Delta p)^2} +ix )
- {p_0^2 \over 4 (\Delta p)^2}} = 
\]
This is gaussian
\[
\int e^{-ap^2 +bp} = \sqrt{{\pi \over a}}e^{b^2 \over 4a}
\]
\[
a={1 \over 4 (\Delta p)^2} (1  + i {2 (\Delta p)^2 t \over  \mu} )
\]
\[
b=({p_0 \over 2  (\Delta p)^2} +ix ) 
\]
\[
{1 \over a} = {4 (\Delta p)^2 \over 1  + i {2 (\Delta p)^2 t \over  \mu}}
\]
\[
b^2 = -x^2 + ix {p_0 \over  (\Delta p)^2} + {p_0^2 \over 4  (\Delta p)^4}
\]
\[
{b^2 \over 4a} =
\]
\[
{(\Delta p)^2 \over 1  + i {2 (\Delta p)^2 t \over  \mu}}
\left (-x^2 + ix {p_0 \over   (\Delta p)^2} + {p_0^2 \over 4  (\Delta p)^4}
\right ) =  
\]
\[
{(\Delta p)^2 \over 1  + {4 (\Delta p)^4 t^2 \over  \mu^2}}
\left (-x^2 + { t \over  \mu}2x p_0
+ {p_0^2 \over 4  (\Delta p)^4}
+ i(x {p_0 \over  (\Delta p)^2} 
+  x^2 {2 (\Delta p)^2 t \over  \mu}
- { t \over  \mu}{p_0^2 \over 2  (\Delta p)^2})
\right ) =
\]
=  
\[
{(\Delta p)^2 \over 1  + {4 (\Delta p)^4 t^2 \over  \mu^2}}
\left (-(x-{p_0t\over \mu})^2 + {p_0^2 t^2 \over \mu^2} 
+ {p_0^2 \over 4  (\Delta p)^4}
+ i(2x {p_0 \over 2  (\Delta p)^2} 
+  x^2 {2 (\Delta p)^2 t \over  \mu}
- { t \over  \mu}{p_0^2 \over 2  (\Delta p)^2})
\right ) =
\]
\[
{(\Delta p)^2 \over 1  + {4 (\Delta p)^4 t^2 \over  \mu^2}}
\left (-(x-{p_0t\over \mu})^2 
+ i(2x {p_0 \over 2  (\Delta p)^2} 
+  x^2 {2 (\Delta p)^2 t \over  \mu}
- { t \over  \mu}{p_0^2 \over 2  (\Delta p)^2})
\right )+
\]
\[
{(\Delta p)^2 \over 1  + {4 (\Delta p)^4 t^2 \over  \mu^2}}
{p_0^2 \over 4  (\Delta p)^4} (1+ {4 (\Delta p)^4 t^2 \over \mu^2}) =
\]
\[
{(\Delta p)^2 \over 1  + {4 (\Delta p)^4 t^2 \over  \mu^2}}
\left (-(x-{p_0t\over \mu})^2 
+ i(2x {p_0 \over 2  (\Delta p)^2} 
+  x^2 {2 (\Delta p)^2 t \over  \mu}
- { t \over  \mu}{p_0^2 \over 2  (\Delta p)^2})
\right )
\]
\[
+ {p_0^2 \over 4  (\Delta p)^2} 
\]
With this evaluation
\[
\langle x \vert \psi (t) \rangle =
\]
\[
(2 \pi)^{-1/2} (2 \pi)^{-1/4} (\Delta p)^{-1/2}
e^{- {p_0^2 \over 4 (\Delta p)^2}}
e^{p_0^2 \over 4  (\Delta p)^2}
2(\Delta p)\sqrt{{\pi  \over 1  + i {2 (\Delta p)^2 t \over  \mu}}}
\times
\]
\[
e^{
{(\Delta p)^2 \over 1  + {4 (\Delta p)^4 t^2 \over  \mu^2}}
\left (-(x-{p_0t\over \mu})^2
+ i(2x {p_0 \over 2  (\Delta p)^2} 
+  x^2 {2 (\Delta p)^2 t \over  \mu}
- { t \over  \mu}{p_0^2 \over 2  (\Delta p)^2})
\right )} =
\]
\[
(2 \pi)^{-3/4}
2(\Delta p)^{1/2}\sqrt{{\pi  \over 1  + i {2 (\Delta p)^2 t \over  \mu}}}
\times
\]
\[
e^{
{(\Delta p)^2 \over 1  + {4 (\Delta p)^4 t^2 \over  \mu^2}}
\left (-(x-{p_0t\over \mu})^2
+ i(2x {p_0 \over 2  (\Delta p)^2} 
+  x^2 {2 (\Delta p)^2 t \over  \mu}
- { t \over  \mu}{p_0^2 \over 2  (\Delta p)^2})
\right )} =
\]
\[
(2 \pi)^{-1/4}
( 2\Delta p)^{1/2}\sqrt{{1  \over 1  + i {2 (\Delta p)^2 t \over  \mu}}}
\times
\]
\[
e^{-
{(\Delta p)^2 \over 1  + {4 (\Delta p)^4 t^2 \over  \mu^2}}
(x-{p_0t\over \mu})^2}
e^{i{1 \over 1  + {4 (\Delta p)^4 t^2 \over  \mu^2}}
(x p_0   
+ 2 x^2 (\Delta p)^4 {t \over  \mu}
- { t \over  \mu}{p_0^2 \over 2 })}
\]
The phase terms are unimportant - If we ignore them
\[
(2 \pi)^{-1/4}
( 2\Delta p)^{1/2}\sqrt{{1  \over 1  + i {2 (\Delta p)^2 t \over  \mu}}}
\times
\]
\[
e^{-
{(\Delta p)^2 \over 1  + {4 (\Delta p)^4 t^2 \over  \mu^2}}
(x-{p_0t\over \mu})^2}
\]
We also note that $(2 \Delta p)^{1/2} = (\Delta x)^{-1/2}$ 

Classically
\[
{\mu \over p_0}(2 \Delta x+ R) =
{\mu \over p_0}({1 \over \Delta p}+ R) =s
\]
is a lower bound on the scattering time, $s$, where $R$ is the range of 
the potential.

In three dimensions the above becomes
\[
\langle \mathbf{x} \vert \psi (t)  \rangle =
\]
\[
(2 \pi)^{-3/4}
( 2\Delta p)^{3/2}(\sqrt{{1  \over 1  + i {2 (\Delta p)^2 t \over  \mu}}})^3
\times
\]
\[
e^{-
{(\Delta p)^2 \over 1  + {4 (\Delta p)^4 t^2 \over  \mu^2}}
(\mathbf{x}-{\mathbf{p}_0t\over \mu})^2}
\]
These state have mean momentum $\mathbf{p}_0$ mean coordinate
at time $t=0$:
$\langle \mathbf{x}(t=0) \rangle =0$, and undertainty in
momentum $\Delta p$ and undertainty in coordinate
$\Delta x = {1 \over 2 \Delta p}$
 
And for multiparticle scattering we use symmetrized products of the
above with different $\mathbf{p}_0$
  



\section{grids}

In order to choose how to partition the real line we start with a gaussian 
potential of the form 

In this section we discuss the choice of ``quadrature points''.
we consider the simplest case of a gaussian potential
\beq
V(x) = \lambda e^{-x^2/r_0^2}. 
\label{g.1}
\eeq
The intervals $I_k$ are chosen so the potential and wave function
are slowly varying on $I_k$.  

We divide the range of the potential, $[0,V_{max}]$, into
equally spaced intervals.  This gives a maximal variation in the
potential of $V_{max}/(M+1)$ where $M+1$ is the number of intervals.
The inverse image of each of these provides a decomposition of the
real line into intervals.

As an example consider the one-dimensional Gaussian potential
with strength $V_0$ and range $R$: 
\beq
V(x) = V_0 e^{-x^2/R^2}
\qquad 
x= R \sqrt{-\ln (V/V_0)}
\label{g.3}
\eeq
This suggests a symmetric distribution of points about the
origin given by 
\beq
x[n] = \pm R \sqrt{-\ln ({n \over 2N+1}) }
\qquad
-N \leq n \leq N .
\label{g.4}
\eeq
Including the half infinite intevals we have
$M+1=2N +2$ intervals, with finite boundaries.
\beq
m= N+n+1 :1 \to M = 2N+1 \qquad n = m -N -1
\label{g.5}
\eeq
We relabel these using 
\beq
x[m] = \pm R \sqrt{-\ln ({m-N-1 \over 2N+1})}
\qquad
1 \leq m \leq M=2N+1
\label{g.6}
\eeq
We choose the quadrature points
\beq
y[m] = \pm R \sqrt{-\ln ({m-N-1/2 \over 2N+1})}
\qquad
0 \leq m \leq M=2N+1 .
\label{g.7}
\eeq
This grid is chosen so the potential varies slowly on the grid.
It is also importnant that initial and
final wave functions also vary slowly on
this grid.  This can be checked and additional points can 
included if necessary.

\section{Numerical tests}

\begin{figure}
\centering
\includegraphics[scale=.7]{potential.pdf}
\caption{potential: lambda=1,V0=1, r0=1}
\label{fig:1}
\end{figure}  

Here $r0=1$, $V0=1$, $delta_{p_i}=1.0$, $p_{0_i}=3.0$:
\begin{figure}
\centering
\includegraphics[scale=.7]{wave_function.pdf}
\caption{initial wave function:t=0}
\label{fig:1}
\end{figure}

does not appear to converge

delta t =.1 

t =  9.063e-02 +i  -1.408e-02   2
t =  8.774e-02 +i  -1.878e-02   3
t =  8.422e-02 +i  -2.244e-02   4
t =  8.014e-02 +i  -2.544e-02   5
t =  7.557e-02 +i  -2.795e-02   6
t =  7.055e-02 +i  -3.008e-02   7
t =  6.513e-02 +i  -3.186e-02   8
t = -2.234e-02 +i  -4.714e-03  30
t = -2.187e-02 +i  -4.086e-03  31 
t = -8.515e-03 +i  1.293e-03   40
Given grid (\ref{}) we wat to see if it looks adequate for 
computation.  We also want to use semiclassical analysis 
to estimate the minimal scatering time $2s$.  Formally we 
expect that $p_0t\over \mu \gtrsim 2r 2 \Delta x$.  
For a first guess we try 
\[
s > {4 \mu r \over (p_0 \Delta x)}
\]
If we take $\Delta p = \Delta x = 1/2$, $p_0=1$,$r=1$ and $\mu=1$ 
this suggests that we might expect to see convergence for $s>8$.


\section{Appendix: Asymptotic  wave packets}

A natural choice for the free wave packets $\vert \psi_{i,0}(t)
\rangle$ that appear in (\ref{s.2}) are minimal uncertainty states
with narrow momentum distributions.  These are Gaussian wave packets
or coherent states.  For a given small momentum spread these have the
smallest spatial distribution consistent with the constraints imposed
by the uncertainty principle.  It is desirable to make both of these
dimensions small; the small momentum spread is needed to get a
wave-packet-independent phase shift, while a narrow wave packet in
coordinate space minimizes the ``time'' $s$ needed for convergence.
The uncertainty principle implies that decreasing one of these
increases the other, which justifies our preference for minimal
uncertainty states.  The other important advantage of minimal
uncertainty states is because they are Gaussian, they can be computed
analytically.

Free minimal uncertainty wave packets in $n$ dimensions are products of 
$n$ one-dimensional minimal uncertainty wave packets.  Multiparticle 
asymptotic states can be taken to be products of single-particle minimal 
uncertainty states. 

The one-dimensional a minimal uncertainty wave packet with mean
momentum $p_0$, mean position $0$ at time $t=0$ and uncertainty
$\Delta p$ in $p$ and $\Delta x= {1 \over 2 \Delta p}$ in $x$, has the
form
\beq
\langle x\vert \psi_{0,p_0} (t) \rangle = 
\langle x \vert e^{-i H_0 t} \vert \psi_{0,p_0} (0) \rangle =
\int \langle x \vert p \rangle e^{-i p^2 t/ 2\mu} dp \langle p 
\vert \psi_{0,p_0} (0) \rangle
\label{aw.1}
\eeq
where the initial unit normalized wave function is 
\beq  
\langle p \vert \psi_{0,p_0} (0) \rangle = (2 \pi)^{-1/4} (\Delta p)^{-1/2}
e^{- {(p-p_0)^2 \over 4 (\Delta p)^2}}.
\label{aw.2}
\eeq
Performing the momentum integral gives
\[
\langle x\vert \psi_{0,p_0} (t) \rangle =
\]
\beq
(2 \pi)^{-1/4}
\sqrt{{2\Delta p  \over 1  + i 2 (\Delta p)^2 t/ \mu}}
\exp ({-
{(\Delta p)^2 \over 1  + 4 (\Delta p)^4 t^2 / \mu^2}
(x-{p_0t\over \mu})^2})
\exp (i{x p_0   
+ 2 x^2 (\Delta p)^4 t/  \mu
- p_0^2 t /(2 \mu) \over 1  + 4 (\Delta p)^4 t^2 / \mu^2})
\label{aw.3}
\eeq
The corresponding minimal uncertainty wave packet in three dimensions
is the product of three such factors,
\[
\langle \mathbf{x} \vert \psi (t)  \rangle =
(2 \pi)^{-3/4}
\left ({2\Delta p  \over 1  + i 2 (\Delta p)^2 t/ \mu}\right )^{3/2}
\exp (-
{(\Delta p)^2 \over 1  + 4 (\Delta p)^4 t^2 /  \mu^2}
(\mathbf{x}-{\mathbf{p}_0t\over \mu})^2)
\times
\]
\beq
\exp ({i{
\mathbf{x}\cdot \mathbf{p}_0   
+ 2 \mathbf{x}^2 (\Delta p)^4 t /\mu
-  t \mathbf{p}_0^2 / (2\mu) 
\over 1  + 4 (\Delta p)^4 t^2/  \mu^2}
} )
\label{aw.4}
\eeq
Note that while the minimal uncertainty states are computationally
useful, the phase shifts extracted from (\ref{s.7}) are independent of
the wave packet provided it is sufficiently narrow.  We can also see
that these free wave packets do not exhibit strong oscillations.  
The $1/t^2$ in the denominator suppresses the effects of the 
phase for large $t$. 


We also need integrals of these wave packets over finite or
semi infinite intervals to compute $\psi_{f0}^* (s,I_{m})$.

In order to compute these integrals it is useful to simplify
(\ref{aw.3}).  Formally this integral has the form
\beq
\langle x \vert \psi(t) \rangle =
A e^{-ax^2 +bx +c} =
A e^{c+{b^2 \over 4a}} e^{-a (x-{b \over 2a})^2}
\label{aw.4}
\eeq
To compute this first define
\beq
(1 + i {2 (\delta p)^2 t \over \mu}) = \lambda e^{i \phi}
\label{aw.5}
\eeq
where
\beq
\lambda = \sqrt{1 + {4 (\delta p)^4t^2 \over \mu^2}}
\qquad
\tan (\phi ) = {2 (\delta p)^2 t \over \mu}
\label{aw.6}
\eeq
In terms of these quantities the parameters in (\ref{aw.4}) are
\beq
A = (2 \pi)^{-1/4}\sqrt{2 (\delta p) \over \lambda} e^{-i \phi/2}
\label{aw.7}
\eeq
\beq
a= {(\delta p)^2 \over \lambda^2} - i {2 t (\delta p)^4 \over \mu \lambda^2}=
{(\delta p)^2 \over \lambda^2}(1 -i {2 (\delta p)^2 t \over \mu})
= {(\delta p)^2 \over \lambda}e^{-i \phi}
\label{aw.8}
\eeq
\beq
b= {2 (\delta p)^2 p_0 t \over \mu \lambda^2} + i {p_0 \over \lambda^2}=
i {p_0 \over \lambda^2}(1 -i {2 (\delta p)^2t \over \mu}) =
i {p_0 \over \lambda^2} \lambda e^{-i \phi} =
i {p_0 \over \lambda}e^{-i \phi}
\label{aw.9}
\eeq
\beq
c= - {1 \over \lambda^2} ({(\delta p)^2 p_0^2 t^2 \over \mu^2}
+i {p_0^2 t \over 2\mu}) =
-i {p_0^2 t \over 2 \mu \lambda^2}(1 - i {2 (\delta p)^2 t \over \mu})=
-i {p_0^2 t \over 2 \mu \lambda^2}\lambda e^{-i \phi} =
-i {p_0^2 t \over 2 \mu \lambda} e^{-i \phi}
\label{aw.10}
\eeq
\beq
{b^2 \over 4a} +c =
- {p_0^2 \over \lambda^2}e^{-2i\phi}{\lambda \over 4 (\delta p)^2}e^{i \phi}
-i {p_0^2 t \over 2 \mu \lambda} e^{-i \phi} =
- {e^{-i \phi}p_0^2 \over 4 \lambda (\delta p)^2 }(1 +i {2(\delta p)^2t \over \mu)}=
- {e^{-i \phi}p_0^2 \over 4 \lambda (\delta p)^2 }\lambda e^{i\phi} =
-{p_0^2 \over 4 (\delta p)^2 }
\label{aw.11}
\eeq
\beq
{b \over 2a} = i {p_0\over \lambda}e^{-i \phi} {\lambda \over 2 (\delta p)^2}
 e^{i\phi} = i {p_0 \over 2 (\delta p)^2}
\label{aw.12}
\eeq
This leads to a simple form for the wave function
\beq
\langle x \vert \psi (t) \rangle =
(2 \pi)^{-1/4} \sqrt{{2 (\delta p) \over \lambda }}e^{-i \phi/2}
e^{- {(\delta p)^2 \over \lambda} e^{-i \phi}
(x - i {p_0 \over 2 (\delta p)^2})^2 - {p_0^2 \over 4 (\delta p)^2}}
\label{aw.13}
\eeq
Note that the phase is always less than ${\pi \over 2}$ so the
function falls off for large $x$.  In three dimensions this becomes
\beq
\langle \mathbf{x} \vert \psi (t) \rangle =
(2 \pi)^{-3/4} ({{2 (\delta p) \over \lambda }})^{3/2}e^{-i 3\phi/2}
e^{- {(\delta p)^2 \over \lambda} e^{-i \phi}
(\mathbf{x} - i {\mathbf{p}_0 \over 2 (\delta p)^2})^2 - {\mathbf{p}_0^2 \over 4 (\delta p)^2}}
\label{aw.14}
\eeq
We are now in a position to perform the integral
\[
\int_{x_l}^{x_h} \langle x \vert \psi (t) \rangle dx =
\]
\[
(2 \pi)^{-1/4} \sqrt{2 (\delta p) \over \lambda }e^{-i \phi/2}
e^{- {p_0^2 \over 4 (\delta p)^2}}
\int_{x_l}^{x^h} e^{- {(\delta p)^2 \over \lambda} e^{-i \phi}
(x - i {p_0 \over 2 (\delta p)^2})^2}=
\]
\beq
(2 \pi)^{-1/4} \sqrt{2 (\delta p) \over \lambda }e^{-i \phi/2}
e^{- {p_0^2 \over 4 (\delta p)^2}}
{\sqrt{\lambda \over (\delta p)^2}}e^{i \phi/2}
\int_{y_l}^{y_h} e^{- {y^2} }dy
\label{aw.15}
\eeq
where
\beq
y = {(\delta p)\over \sqrt{\lambda}}e^{-i \phi/2}(x-i{p_0 \over 2 (\delta p)^2})
\label{aw.16}
\eeq
Using
\beq
\int_0^z e^{-y^2}dy = {\sqrt{\pi} \over 2} \mbox{erf}(z)
\label{aw.17}
\eeq
which gives
\[
\int_{x_l}^{x^h} \langle x \vert \psi (t) \rangle dx =
\]
\[
(2 \pi)^{-1/4} \sqrt{{2 (\delta p) \over \lambda }}
{\sqrt{\pi} \over 2} e^{-i \phi/2}
e^{- {p_0^2 \over 4 (\delta p)^2}}
{\sqrt{\lambda \over (\delta p)^2}}e^{i \phi/2} \times
\]
\[
\left (\mbox{erf}({(\delta p)\over \sqrt{\lambda}}e^{-i \phi/2}
(x_h-i{p_0 \over 2 (\delta p)^2})
- \mbox{erf}({(\delta p)\over \sqrt{\lambda}}e^{-i \phi /2}
(x_l-i{p_0 \over 2 (\delta p)^2}))
\right )=
\] 
\beq
(2 \pi)^{-1/4} \sqrt{{\pi \over 2 (\delta p)}} 
e^{- {p_0^2 \over 4 (\delta p)^2}}
\left  (\mbox{erf}({(\delta p)\over \sqrt{\lambda}}e^{-i \phi/2}
(x_h-i{p_0 \over 2 (\delta p)^2}))
- \mbox{erf}({(\delta p)\over \sqrt{\lambda}}e^{-i \phi /2}
(x_l-i{p_0 \over 2 (\delta p)^2}))
\right )
\label{aw.18}
\eeq
Here the only complex quantites are the argument of the
error functions.
\[
{(\delta p)\over \sqrt{\lambda}}(\cos (\phi/2) - i \sin (\phi/2))
(x -i{p_0 \over 2 (\delta p)^2}) =
\]
\beq
{(\delta p)\over \sqrt{\lambda}}\cos (\phi/2) x - 
{p_0 \over 2 \delta p \sqrt{\lambda}}\sin (\phi/2)
-i {(\delta p)\over \sqrt{\lambda}} x \sin (\phi/2)
-i {p_0 \over 2(\delta p) \sqrt{\lambda}}\cos (\phi/2))
\eeq
These numbers appear in the error function

xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
These integrals have the general form
\[
A \int_{x_1}^{x2}  e^{-ax^2 + bx +c} dx =
\]
\beq
{A\over 2} {\sqrt{\pi \over a}} e^{{b^2 \over 4a}+c}\left (
\mbox{erf}(\sqrt{a} (x_2 - {b \over 2a}) - 
\mbox{erf}(\sqrt{a} (x_1 - {b \over 2a}) \right ) 
\label{aw.5}
\eeq
where
\beq
A= (2 \pi)^{-1/4}
\sqrt{{2\Delta p  \over 1  + i 2 (\Delta p)^2 t /  \mu}}
\label{aw.6}
\eeq
\beq
a_0= {(\Delta p)^2 \over 1  + 4 (\Delta p)^4 t^2 / \mu^2}
\label{aw.7}
\eeq
\beq
a= a_0(1 +i 2  (\Delta p)^2 t/\mu)
\label{aw.8}
\eeq
\beq
b = 2a_0p_0 t/\mu + i a_0 p_0/ (\Delta p)^2
\label{aw.9}
\eeq
\beq
c=-a_0 p_0^2 t^2/\mu^2 - i {a_0 tp_0^2 \over 2 \mu (\Delta p)^2}  
\label{aw.10}
\eeq
This gives analytic expressions for the integrals over the
finite intervals.  For the half-infinite intervals the
relevant quantities can be computed by subtracting the
finite interval result from the Gaussian integral on the
line:
\[
A \int_{x_1}^{\infty}  e^{-ax^2 + bx +c} dx =
\]
\beq
{A\over 2} {\sqrt{\pi \over a}} e^{{b^2 \over 4a}+c}\left (
1 - 
\mbox{erf}(\sqrt{a} (x_1 - {b \over 2a}) \right ) 
\label{aw.11}
\eeq
\[
A \int_{-\infty}^{x2}  e^{-ax^2 + bx +c} dx =
\]
\beq
{A\over 2} {\sqrt{\pi \over a}} e^{{b^2 \over 4a}+c}\left (
\mbox{erf}(\sqrt{a} (x_2 - {b \over 2a}) -1 \right ) 
\label{aw.12}
\eeq

The three dimensional form of this equation is the product of
three 1-dimensional expressions of this form.  In three dimensions
there are $(M+1)^3$ distinct windows per time slice.

The other factor needed is the denominator of (\ref{s.7}).
For these states this can be expressed as
\[
\langle \psi_{0,p_f} (0)
\vert \delta ({p^2 - p_0^2\over 2 \mu}) 
\vert \psi_{0,p_i} (0) \rangle =
\]
\[
\int dp
(2 \pi)^{-1/4} (\Delta p_f)^{-1/2}
e^{- {(p-p_{0f})^2 \over 4 (\Delta p_f)^2}}
{\mu \over p_0} (\delta (p-p_0) + \delta (p+p_0))   
(2 \pi)^{-1/4} (\Delta p_i)^{-1/2}
e^{- {(p-p_{0i})^2 \over 4 (\Delta p_i)^2}} =
\]
\beq
{\mu \over p_0} (2 \pi)^{-1/2} (\Delta p_f \Delta p_i)^{-1/2}
\left (e^{- {(p_0-p_{0f})^2 \over 4 (\Delta p_f)^2} -
{(p_0-p_{0i})^2 \over 4 (\Delta p_i)^2}} +
e^{- {(p_0+p_{0f})^2 \over 4 (\Delta p_f)^2} -
{(p_0+p_{0i})^2 \over 4 (\Delta p_i)^2}}
\right ).
\label{aw.13}
\eeq
This is the one-dimensional form which can be written as
\[
\langle \psi_{0,p_f} (0)
\vert \delta ({p^2 - p_0^2\over 2 \mu}) 
\vert \psi_{0,p_i} (0) \rangle =
\]
\beq
{\mu \over p_0} (2 \pi)^{-1/2} (\Delta p_f \Delta p_i)^{-3/2}  
e^{-p_0^2 ({1 \over 4 (\Delta p_f)^2} + {1\over 4 (\Delta p_i)^2})}
e^{- ({p_{0f}^2 \over 4 (\Delta p_f)^2} + {p_{0i}^2\over 4 (\Delta p_i)^2})}
\cosh ({p_0 w \over 2})
\label{aw.14}
\eeq
where
\beq
w = \vert ({\mathbf{p}_{0f} \over (\Delta p_f)^2} + 
{{\mathbf{p}_{0i}\over (\Delta p_i)^2})} \vert .
\label{aw.12}
\eeq
In three dimensions the integral can also be performed - in this case
there are additional phase space integrals.
\[
\int 
(2 \pi)^{-3/4} (\Delta q)^{-3/2}
e^{- {(\mathbf{p}-\mathbf{p}_{0})^2 \over 4 (\Delta p)^2}}
d\mathbf{p} 
\delta ({\mathbf{p}^2 \over 2 \mu} - {\mathbf{q}^2 \over 2 \mu} )
d\mathbf{q}
(2 \pi)^{-3/4} (\Delta q)^{-3/2}
e^{- {(\mathbf{q}-\mathbf{q}_{0})^2 \over 4 (\Delta q)^2}}= 
\]
\[
\int 
(2 \pi)^{-3/4} (\Delta q)^{-3/2}
e^{- {(p^2 - 2 p p_0 u_1 + p_0^2) \over 4 (\Delta p)^2}}
du_1 d\phi_1  du_2 d \phi_2 p^4 dp dq {\mu \over p} \delta (p-q)
(2 \pi)^{-3/4} (\Delta q)^{-3/2}
e^{- {(p^2 - 2 {p}{q}_{0}u_2 + q_0^2 ) \over 4 (\Delta q)^2}}=
\]
\[
\int 
\mu (2 \pi)^{2-3/2} (\Delta q \Delta q )^{-3/2}  
e^{- {(p^2 + p_0^2) \over 4 (\Delta p)^2}}
e^{- {(p^2 + q_0^2) \over 4 (\Delta q)^2}}
p^3 dp 
{2 (\Delta p)^2 \over p p_0} 
(e^{pp_0 \over 2 (\Delta p)^2} - 
 e^{-pp_0 \over 2 (\Delta p)^2})  
{2 (\Delta q)^2 \over p q_0} 
(e^{pq_0 \over 2 (\Delta q)^2} - 
 e^{-pq_0 \over 2 (\Delta q)^2})  =
\]
\beq
\mu (2 \pi)^{2-3/2} (\Delta p \Delta q )^{-3/2}  
{2 (\Delta p)^2 \over  p_0} 
{2 (\Delta q)^2 \over  q_0}  
e^{- {p_0^2 \over 4 (\Delta p)^2}}
e^{- {p_0^2 \over 4 (\Delta q)^2}}
\int 
e^{- p^2 ({1 \over 4 (\Delta p)^2} + {1 \over 4 (\Delta q)^2})}
p dp 
(e^{pp_0 \over 2 (\Delta p)^2} - 
e^{-pp_0 \over 2 (\Delta p)^2})  
(e^{pq_0 \over 2 (\Delta q)^2} - 
e^{-pq_0 \over 2 (\Delta q)^2})  
\label{aw.13}
\eeq
The remaining inegrals are all Gaussian of the general form 
\beq
\int e^{-ap^2 +bp} p dp = \sqrt{\pi \over a} {b \over 2a} e^{b^2 \over 4a} 
\label{aw.14}
\eeq
where
\beq
a= ({1 \over 4 (\Delta p)^2} + {1 \over 4 (\Delta q)^2})
\label{aw.15}
\eeq
\beq
b_1= {p_0 \over 2 (\Delta p)^2} + {q_0 \over 2 (\Delta q)^2}
\label{aw.16}
\eeq
\beq
b_2= -1({p_0 \over 2 (\Delta p)^2} + {q_0 \over 2 (\Delta q)^2})=- b1
\label{aw.17}
\eeq
\beq
b_3= {p_0 \over 2 (\Delta p)^2} - {q_0 \over 2 (\Delta q)^2}
\label{aw.18}
\eeq
\beq
b_4= -{p_0 \over 2 (\Delta p)^2} + {q_0 \over 2 (\Delta q)^2}= - b4
\label{aw.19}
\eeq
So the full integral becomes
\[
\int 
(2 \pi)^{-3/4} (\Delta q)^{-3/2}
e^{- {(\mathbf{p}-\mathbf{p}_{0})^2 \over 4 (\Delta p)^2}}
d\mathbf{p} 
\delta ({\mathbf{p}^2 \over 2 \mu} - {\mathbf{q}^2 \over 2 \mu} )
d\mathbf{q}
(2 \pi)^{-3/4} (\Delta q)^{-3/2}
e^{- {(\mathbf{q}-\mathbf{q}_{0})^2 \over 4 (\Delta q)^2}}= 
\]
\[
{2 \mu (2 \pi)^{1/2} (\Delta p \Delta q )^{1/2}  
\over  p_0 q_0} \sqrt{\pi \over a^3} 
e^{- {p_0^2 \over 4 (\Delta p)^2}}
e^{- {p_0^2 \over 4 (\Delta q)^2}}
((b_1 e^{b_1^2 \over 4a} -e^{-b_1^2 \over 4a}) -
b_3  (e^{b_3^2 \over 4a}- e^{-b_3^2 \over 4a}))
\]
\beq
{4 \mu (2 \pi)^{1/2} (\Delta p \Delta q )^{1/2}  
\over  p_0 q_0} \sqrt{\pi \over a^3} 
e^{- {p_0^2 \over 4 (\Delta p)^2}}
e^{- {p_0^2 \over 4 (\Delta q)^2}}
((b_1 \sinh ({b_1^2 \over 4a}) -
b_3 \sinh ({b_3^2 \over 4a}))
\label{aw.20}
\eeq



{\color{blue}
\bigskip  
\hrule
\begin{center}
{\bf Computation of the integrals - not for publication:}
\end{center}
\bigskip
\[  
\langle p \vert \psi \rangle = (2 \pi)^{-1/4} (\Delta p)^{-1/2}
e^{- {(p-p_0)^2 \over 4 (\Delta p)^2}}
\]  
\[
\Delta p \Delta x = {1 \over 2}
\]
\[  
\langle x \vert \psi \rangle = (2 \pi)^{-1/4} (\Delta x)^{-1/2}
e^{i p_0 x - {x^2 \over 4 (\Delta x)^2}}
\]  
\[
\langle x \vert \psi (t)  \rangle =
\]
\[
\int (2 \pi)^{-1/2} e^{ipx - i{p^2 \over 2 \mu}t}
(2 \pi)^{-1/4} (\Delta p)^{-1/2}
e^{- {(p-p_0)^2 \over 4 (\Delta p)^2}} dp =
\]
\[
\int (2 \pi)^{-1/2} (2 \pi)^{-1/4} (\Delta p)^{-1/2}
e^{ipx - i{p^2 \over 2 \mu}t - {(p-p_0)^2 \over 4 (\Delta p)^2}} dp =
\]
\[
(2 \pi)^{-1/2} (2 \pi)^{-1/4} (\Delta p)^{-1/2}
\int
e^{
-p^2 {1 \over 4 (\Delta p)^2} (1  + i {2 (\Delta p)^2 t \over  \mu} )+
p ({p_0 \over 2  (\Delta p)^2} +ix )
- {p_0^2 \over 4 (\Delta p)^2}} = 
\]
This is gaussian
\[
\int e^{-ap^2 +bp} = \sqrt{{\pi \over a}}e^{b^2 \over 4a}
\]
\[
a={1 \over 4 (\Delta p)^2} (1  + i {2 (\Delta p)^2 t \over  \mu} )
\]
\[
b=({p_0 \over 2  (\Delta p)^2} +ix ) 
\]
\[
{1 \over a} = {4 (\Delta p)^2 \over 1  + i {2 (\Delta p)^2 t \over  \mu}}
\]
\[
b^2 = -x^2 + ix {p_0 \over  (\Delta p)^2} + {p_0^2 \over 4  (\Delta p)^4}
\]
\[
{b^2 \over 4a} =
\]
\[
{(\Delta p)^2 \over 1  + i {2 (\Delta p)^2 t \over  \mu}}
\left (-x^2 + ix {p_0 \over   (\Delta p)^2} + {p_0^2 \over 4  (\Delta p)^4}
\right ) =  
\]
\[
{(\Delta p)^2 \over 1  + {4 (\Delta p)^4 t^2 \over  \mu^2}}
\left (-x^2 + { t \over  \mu}2x p_0
+ {p_0^2 \over 4  (\Delta p)^4}
+ i(x {p_0 \over  (\Delta p)^2} 
+  x^2 {2 (\Delta p)^2 t \over  \mu}
- { t \over  \mu}{p_0^2 \over 2  (\Delta p)^2})
\right ) =
\]
=  
\[
{(\Delta p)^2 \over 1  + {4 (\Delta p)^4 t^2 \over  \mu^2}}
\left (-(x-{p_0t\over \mu})^2 + {p_0^2 t^2 \over \mu^2} 
+ {p_0^2 \over 4  (\Delta p)^4}
+ i(2x {p_0 \over 2  (\Delta p)^2} 
+  x^2 {2 (\Delta p)^2 t \over  \mu}
- { t \over  \mu}{p_0^2 \over 2  (\Delta p)^2})
\right ) =
\]
\[
{(\Delta p)^2 \over 1  + {4 (\Delta p)^4 t^2 \over  \mu^2}}
\left (-(x-{p_0t\over \mu})^2 
+ i(2x {p_0 \over 2  (\Delta p)^2} 
+  x^2 {2 (\Delta p)^2 t \over  \mu}
- { t \over  \mu}{p_0^2 \over 2  (\Delta p)^2})
\right )+
\]
\[
{(\Delta p)^2 \over 1  + {4 (\Delta p)^4 t^2 \over  \mu^2}}
{p_0^2 \over 4  (\Delta p)^4} (1+ {4 (\Delta p)^4 t^2 \over \mu^2}) =
\]
\[
{(\Delta p)^2 \over 1  + {4 (\Delta p)^4 t^2 \over  \mu^2}}
\left (-(x-{p_0t\over \mu})^2 
+ i(2x {p_0 \over 2  (\Delta p)^2} 
+  x^2 {2 (\Delta p)^2 t \over  \mu}
- { t \over  \mu}{p_0^2 \over 2  (\Delta p)^2})
\right )
\]
\[
+ {p_0^2 \over 4  (\Delta p)^2} 
\]
With this evaluation
\[
\langle x \vert \psi (t) \rangle =
\]
\[
(2 \pi)^{-1/2} (2 \pi)^{-1/4} (\Delta p)^{-1/2}
e^{- {p_0^2 \over 4 (\Delta p)^2}}
e^{p_0^2 \over 4  (\Delta p)^2}
2(\Delta p)\sqrt{{\pi  \over 1  + i {2 (\Delta p)^2 t \over  \mu}}}
\times
\]
\[
e^{
{(\Delta p)^2 \over 1  + {4 (\Delta p)^4 t^2 \over  \mu^2}}
\left (-(x-{p_0t\over \mu})^2
+ i(2x {p_0 \over 2  (\Delta p)^2} 
+  x^2 {2 (\Delta p)^2 t \over  \mu}
- { t \over  \mu}{p_0^2 \over 2  (\Delta p)^2})
\right )} =
\]
\[
(2 \pi)^{-3/4}
2(\Delta p)^{1/2}\sqrt{{\pi  \over 1  + i {2 (\Delta p)^2 t \over  \mu}}}
\times
\]
\[
e^{
{(\Delta p)^2 \over 1  + {4 (\Delta p)^4 t^2 \over  \mu^2}}
\left (-(x-{p_0t\over \mu})^2
+ i(2x {p_0 \over 2  (\Delta p)^2} 
+  x^2 {2 (\Delta p)^2 t \over  \mu}
- { t \over  \mu}{p_0^2 \over 2  (\Delta p)^2})
\right )} =
\]
\[
(2 \pi)^{-1/4}
( 2\Delta p)^{1/2}\sqrt{{1  \over 1  + i {2 (\Delta p)^2 t \over  \mu}}}
\times
\]
\[
e^{-
{(\Delta p)^2 \over 1  + {4 (\Delta p)^4 t^2 \over  \mu^2}}
(x-{p_0t\over \mu})^2}
e^{i{1 \over 1  + {4 (\Delta p)^4 t^2 \over  \mu^2}}
(x p_0   
+ 2 x^2 (\Delta p)^4 {t \over  \mu}
- { t \over  \mu}{p_0^2 \over 2 })}
\]
The phase terms are unimportant - If we ignore them
\[
(2 \pi)^{-1/4}
( 2\Delta p)^{1/2}\sqrt{{1  \over 1  + i {2 (\Delta p)^2 t \over  \mu}}}
\times
\]
\[
e^{-
{(\Delta p)^2 \over 1  + {4 (\Delta p)^4 t^2 \over  \mu^2}}
(x-{p_0t\over \mu})^2}
\]
We also note that $(2 \Delta p)^{1/2} = (\Delta x)^{-1/2}$ 

Classically
\[
{\mu \over p_0}(2 \Delta x+ R) =
{\mu \over p_0}({1 \over \Delta p}+ R) =s
\]
is a lower bound on the scattering time, $s$, where $R$ is the range of 
the potential.

In three dimensions the above becomes
\[
\langle \mathbf{x} \vert \psi (t)  \rangle =
\]
\[
(2 \pi)^{-3/4}
( 2\Delta p)^{3/2}(\sqrt{{1  \over 1  + i {2 (\Delta p)^2 t \over  \mu}}})^3
\times
\]
\[
e^{-
{(\Delta p)^2 \over 1  + {4 (\Delta p)^4 t^2 \over  \mu^2}}
(\mathbf{x}-{\mathbf{p}_0t\over \mu})^2}
\]
These state have mean momentum $\mathbf{p}_0$ mean coordinate
at time $t=0$:
$\langle \mathbf{x}(t=0) \rangle =0$, and undertainty in
momentum $\Delta p$ and undertainty in coordinate
$\Delta x = {1 \over 2 \Delta p}$
 
And for multiparticle scattering we use symmetrized products of the
above with different $\mathbf{p}_0$
\hrule
}  


\section{appendix 1}



Steven G. Johnson, A. Cervellino, J. Wuttke: libcerf, numeric library
for complex error functions, version […],
http://apps.jcns.fz-juelich.de/libcerf


\section{Fresnel integrals:}
\bigskip
One of the most important pieces of input to this method are analytic
expressions for the complex one-step probabilities (\ref{njm.6}).
These quantities are expressible in terms of the Fresnel integrals of
the general form:
\beq
\int_{I_k[n_k]}  ds_k e^{i s_k^2}
\label{fi.1}
\eeq
The basic Fresnel integrals are defined as special functions.
They are given in Abramowitz and Stegun as:
\beq
C_{AS}(x) := \int_0^x \cos ({\pi \over 2} t^2) dt
\qquad
S_{AS}(x) := \int_0^x \sin ({\pi \over 2} t^2) dt .
\label{fi.2}
\eeq
We use programs that compute a related
special function:
\beq
C_{c}(x) := \sqrt{{2\over \pi}} \int_0^x \cos (t^2) dt
\qquad
S_{c}(x) := \sqrt{{2\over \pi}} \int_0^x \sin ( t^2) dt
\label{fi.3}
\eeq
{\color{red} 
(Source code for the programs we used are open source codes from
ALGLIB - Proprietor - Bochkanov Sergey Anatolyevich).  The Fresnel
integral can be expressed in terms of these functions
\beq
\sqrt{{2 \over \pi}} \int_0^x e^{i t^2} dt =
\sqrt{{\pi \over 2}}(C_{c}(x) +i S_{c}(x)
\label{fi.4}
\eeq
}
These functions related to the A\&S functions with the same name by 
\beq
C_{AS}(x) = C_{c}(\sqrt{{\pi \over 2}}x)
\qquad
S_{AS}(x) = S_{c}(\sqrt{{\pi \over 2}}x)
\label{fi.5}
\eeq
This indicates that the integrals that appear in the ``complex
probabilities'' can be accurately computed.  Complex probabilities are
computed from integrals over finite or half-infinite intervals.  For
$a$ and $b$ finite the relevant integral can be expressed in terms of
the Fresnel sine and cosine integrals (\ref{fi.2}\ref{fi.3}):
\[
I[a,b] = \int_a^b e^{ix^2} dx =
\sqrt{\pi \over 2} (C_{AS}(\sqrt{{ 2 \over \pi}} x)+ i 
S_{AS}(\sqrt{{ 2 \over \pi}}x))\vert_a^b =
\]
\[
\sqrt{\pi \over 2} (C_{AS}(\sqrt{{ 2 \over \pi}}b)- 
C_{AS}(\sqrt{{ 2 \over \pi}}a))
+i \sqrt{\pi \over 2} (S_{AS}(\sqrt{{ 2 \over \pi}}b)- S_{AS}
(\sqrt{{ 2 \over \pi}}a)) =
\]
\beq
\sqrt{\pi \over 2} (C_{c}(b)- C_{c}(a))
+i \sqrt{\pi \over 2} (S_{c}(b)- S_{c}(a)) .
\label{fi.6}
\eeq
For $a=-\infty$
\beq
I[-\infty,b] = \int_{-\infty}^b e^{ix^2} dx =
{1 \over 2} \int_{-\infty}^{\infty} e^{ix^2} dx - \int_b^0 e^{ix^2} dx =
\sqrt{{\pi \over 2}} \left ({1+i \over 2} + C_c(b)+i S_c(b) \right ) .
\label{fi.7}
\eeq
For $b=\infty$
\beq
I[a,\infty] = \int_a^{\infty} e^{ix^2} dx =
{1 \over 2} \int_{-\infty}^{\infty} e^{ix^2} dx - \int_0^a e^{ix^2} dx =
\sqrt{{\pi \over 2}} \left ({1+i \over 2} - C_c(a)-i S_c(a) \right ) .
\label{fi.8}
\eeq
Even though all of the Fresnel integrals in the intervals $I_{[a,b]}$
can be performed exactly, there are very large number of cylinder sets
and any practical evaluation needs to make use of alternative methods.

\section{test calculation}

In this section we preform numerical tests of the computational
methods proposed above.  For the initial test we work in one
dimension.  We use a gaussian potential.

\section{Complex probabilities:}

The first step in the construction of a the complex probability
distribution is to select partitions of the real line that can 
be refined to make arbitrarily small partitions.  A simple 
choice is to start by dividing the interval $[-1,1]$ into $M+1$ equally 
spaced intervals with $M$ interior boundary points:
\beq
x[n] = -1 + {2n \over M+1} \qquad 1 \leq n \leq M .
\label{cp.1}
\eeq
The width of each interval vanishes as $M$ gets large.

Also define $M+1$ points at the midpoint of each of these
intervals that will be transformed to quadrature
points
\beq
y[n] = -1 + {2n \over M+1} + {1 \over M+1}  \qquad 0 \leq n \leq M .
\label{cp.2}
\eeq
The labels of these points are in a one-to-one correspondence 
with the labels of the intervals of the partition defined by
the boundary points (\ref{cp.1}).
Note that there is one more $y$ than $x$.  Our convention is that 
the $y$ indices start at $0$ while the $x$ indices start at $1$.
This convention is also used for labeling the intervals.
 
Next the points on $[-1,1]$ are transformed to $[-\infty,\infty]$.
The choice of the mapping is important.  Most of the points need to be
where the integrand is largest.  We used the combination of equally
spaced points on the interval and the tangent map (below) to get the
results in the table II.  In general the results of using different
maps should be tested and compared.

The tangent mapping transforms the 
equally spaced points on $[-1,1]$ 
to  transformed boundary points 
on $\mathbb{R}=[-\infty,\infty]$: 
\beq
u[n] = s \tan (x[n])  \qquad 1 \leq n \leq N
\label{cp.3}
\eeq
and transformed interior points
\beq
v[n] = s \tan (y[n])  \qquad 0 \leq n \leq N .
\label{cp.4}
\eeq
Here $s$ is a scale factor that is chosen so most of the
quadrature points lie where the integrand is large or has
the most structure.

The tangent map transforms the equally space intervals on $[-1,1]$
into the 
$M+1$ intervals:
\beq
\underbrace{[-\infty,u[1]]}_{I[0]}
\underbrace{[u[1],u[2]]}_{I[1]}
\cdots
\underbrace{[u[M-1],u[M]]}_{I[N-1]}
\underbrace{[u[M],\infty]}_{I[N]}
\label{cp.5}
\eeq
on the real line with a transformed quadrature point in the
interior of each interval
\beq
v[n] \in I[n] \qquad 0 \leq n \leq M .
\label{cp.6}
\eeq
These partition get smaller in the limit that $M$ gets large.

The complex probability associated with the $n^{th}$
interval is
\beq
P[n] = {\int_{I[n]} e^{ix^2 }dx
\over \int_{\mathbb{R}} e^{ix^2 }dx} =
%{1 \over \sqrt{2 \pi}}(1-i) \int_{I[n]} e^{ix^2 }dx
{1 \over \sqrt{i\pi}} \int_{I[n]} e^{ix^2 }dx .
\label{cp.7}
\eeq
By construction 
\beq
\sum_{n=0}^M P[n]=1 .
\label{cp.8}
\eeq
These probabilities can be computed analytically in terms of
the Fresnel sine and cosine integrals
\beq
P[0] =
{1 \over 2} \left (1 + C_c(u[1])+  S_c(u[1])
+i (S_c(u[1]) - C_c(u[1])) \right )  
\label{cp.9}
\eeq
\[
P[n] =
{1 \over 2} (C_{c}(u[n+1])- C_{c}(u[n]) 
+ S_{c}(u[n+1])- S_{c}(u[n]))
\]
\beq
+i {1\over 2} (S_{c}(u[n+1])- S_{c}(u[n])
-C_{c}(u[n+1])+ C_{c}(u[n]))
\label{cp.10}
\eeq
for $\qquad 0 < n < M$ and for $n=M$:
\beq
P[M] =
{1 \over 2} \left (1 - C_c(u[M])- S_c(u[M]) 
- i (S_c(u[M] - C_c(u[M]) \right ).
\label{cp.11}
\eeq
The expressions for the half infinite intervals used the fact that the
integrand in the Fresnel integral is even.  Then the integrals on
$[0,\infty]$ and $[-\infty,0]$ are half of the full integral.  The
boundary terms are computed using
\beq
\int_a^\infty = \int_0^\infty - \int_0^a = {1\over 2}\int_{-\infty}^\infty - \int_0^a
\qquad
\int_{-\infty}^b = \int_{-\infty}^0 - \int_b^0
= {1\over 2}\int_{-\infty}^\infty-\int_b^0
\label{cp.12}
\eeq
The definitions imply the boundary values $Ci(0) = Si (0) =0$. 

The complex probabilities are associated with partitions of the
real line into disjoint intervals.  Formally the probability associated
with the interval $[a,b]$ is defined by
\beq
P[a,b] = {\int_a^b e^{ix^2} dx \over \int_{-\infty}^{\infty} e^{iy^2} dy }
\label{cp.13}
\eeq
Note that it is enough to work with (\ref{cp.13}) 
since rescaling the exponent simply shifts the endpoints:
\beq
{\int_a^b e^{ix^2/r^2} dx \over \int_{-\infty}^{\infty} e^{iy^2/r^2} dy }=
P[a/r,b/r]
\label{cp.14}
\eeq

\section{Expectation values of random variables:}

The Riemann integral
\beq
I = \int_{\mathbf{R}} e^{ix^2} f(x) dx
\label{ev.1}
\eeq
can be approximated as the expectation value of the random variable
$f(x)$ with respect to the complex probability distribution
$P_{I_n}$:
\[
I = \int_{\mathbf{R}} e^{ix^2} f(x) dx =
\]
\beq
(\int_{\mathbf{R}} e^{ix^2} dx) 
\lim_{M\to \infty} \sum_{n=0}^M P[n] f(v[n])
\label{ev.2}
\eeq
where the normalization integral is 
\beq
\int_{\mathbf{R}} e^{ix^2} dx =
%\sqrt{\pi \over 2}(1+i)
\sqrt{i \pi}
\label{ev.3}
\eeq
which is needed to correctly normalize this expression when $f(x)=1$.
Here the limit is the limit that $M\to \infty$, when the partition is
constructed from equally spaced intervals on $[-1,1]$.

For a first test consider integrating the function:
\beq
f(x) = e^{-x^2 +x}.
\label{ev.4}
\eeq
This choice is probably too simple because it does not oscillate,
however it can be computed exactly, providing a means to test the
convergence of the complex expectation values to the integral.
 
The exact value of the integral is  
\beq
\int_{\mathbb{R}} e^{ix^2 -x^2 +x} dx = 
\sqrt{{\pi \over (1-i)}}e^{{1 \over 4(1-i)}}
=
\sqrt{{\pi (1+i) \over 2}}e^{{1+i \over 8}}=
\sqrt{{\pi \over \sqrt{2}}}e^{{1 \over 8}} \left (\cos ({1 + \pi \over 8})
+ i \sin ({1 +\pi \over 8})\right ) 
\label{ev.5}
\eeq
%{1+i \over 8}}=
%\]
%
%\sqrt{{\pi \over 2}} e^{1/8} (\cos ({1 \over 8})-
%\sin ({1 \over 8})) 
%+i \sqrt{{\pi \over 2}} e^{1/8} (\cos ({1 \over 8})+
%\sin ({1 \over 8})) 
%\label{a.22}
%\eeq


A practical evaluation of the expectation value of a random variable
over products of a large number of complex probability distributions,
like (\ref{d.30}), requires summing over a prohibitively
large number of cylinder sets.

This suggests the possibility of using statistical methods that
selectively choose the most important cylinder sets for a given random
variable.  While this will be challenging for a random variable like
$F(M,\mathbf{y}(\mathbf{s[\mathbf{n}]}))$, an important first test is
to see if statistical methods can be used for functions of
1 variable.

The most trivial way to proceed is to select a random point in the
interval $[-1,1]$ and use a tangent mapping to 
transform to the real line.  This has probability $1/(M+1)$ of being
transformed to the $i^{th}$ interval on the real line.  If we
sum over a very large number $(K>> M)$ random evaluations then
\[
I \approx  \sqrt{{\pi \over 2}}(1+i){\sum_{k=1}^K P[n(k)] f(y[k])
\over \sum_{k=1}^K P[n(k)]} =
\]
\beq
\sqrt{{\pi \over 2}}(1+i){\sum_{k=1}^K P[n(k)] f(y[k])
\over K/(M+1)} =
\sqrt{{\pi \over 2}}(1+i)\sum_{k=1}^K {M+1\over K} P[n(k)] f(y[k])
\label{ev.7}
\eeq
where $n(k)$ is the interval containing the $k^{th}$ point and
and $y[k]$ is the quadrature point in the $n(k)$ interval.




\section{phase shift calculation}

\[
\langle \psi_{f0} (t) \vert e^{-2iHt} \vert \psi_{i0}(-t) \rangle =
\]
\[
\sum \int_{I_{m0}} \psi^* (x,t) dx 
\langle y_{m_00} \vert e^{-i2 H t} \vert \psi_0 (-t) \rangle \approx
\sum_{} P(x_0,I_{m_{m_11}}) e^{ -i V(y_{m_11}) \Delta t} 
\prod_{n=2}^N P(y_{m_{n-1} n-1},I_{n_{n_1n}}) e^{ -i V(y_{m_nn}) \Delta t}
\psi (y_{m_NN},-t)
\]

In terms of the matrix 37 the formula for phase shift becomes
\[
e^{2i \delta (E)} = 
\]
\beq
{\sum_{mk} \psi_{0f}^* (I_m,t) 
M^N_{m,k} \psi_{0i} (y_{k},-t)
\over
\langle \psi_{0f}(0) \vert \delta (E-H_0) \vert
\psi_{0i}(0) \rangle }
\eeq
where 
\beq
\psi_{0f}^* (I_m,t) = \int_{I_m} \psi_{0f}^* (x,t) dx  
\eeq
and each one of these terms can be computed analytically for
any potential.


The two integrals that need to be performed can be computed analytically. 

A\&s 7.4.32:
\beq
\int_{x_1}^{x_2} 
e^{-a x^2 + 2bx +c}dx = 
{1 \over 2} \sqrt{\pi \over a} 
e^{{b^2-ac \over a}}
\left (
\mbox{erf} (\sqrt{a} x_2 + {b \over \sqrt{a}})
-
\mbox{erf} (\sqrt{a} x_1 + {b \over \sqrt{a}})
\right )
\eeq
can be used to integrate the final wave function over the 
intervals that we used to partition the time slices.


This allows us to calculate the full integral using 
\[
\int \phi_{f0}(x_0,t) \langle x_0 \vert e^{-2iHt}
\vert \psi_{i0}(t) \rangle dx_0 =
\]
\beq
\sum_{m} \int_{I_m} \phi_{f0}(x_0,t) dx_0
\langle y_m \vert e^{-2iHt}
\vert \psi_{i0}(t) \rangle dx_0
\eeq
In order to extract the phase shift we need to comptue the 
overlap of the asymptotic states at time 0:
\[
\int NN'\int dp  
e^{-(p-p_{0f})^2 -(p-p_{0i})^2}
\delta (E-p^2 \over 2 \mu )  =
\]
\[
NN'\sqrt{{\mu \over 2E}}  
e^{-(\sqrt{2 \mu E}-p_{0f})^2 -(\sqrt{2 \mu E}-p_{0i})^2}  =
\]
\beq
{1 \over 
\sqrt{\pi}\sqrt{\Delta p_f \Delta p_i}} 
\sqrt{{\mu \over 2E}}  
e^{-(\sqrt{2 \mu E}-p_{0f})^2 -(\sqrt{2 \mu E}-p_{0i})^2} 
\eeq
In three dimensions this becomes



\section{Quantum fields}

The interesting question is whether these ideas can be generalized to
treat local quantum fields.
The starting point is a pair of canonical fields satisfying
\beq
[ \phi (\mathbf{x},t) ,\pi (\mathbf{y},t)] = i \delta (\mathbf{x}-\mathbf{y})
\label{qf:1}
\eeq
If $\{ \xi_n (\mathbf{x}) \}$ is any orthonormal basis of real valued functions
of three-vector variables we can define the algebra of smeared fields 
\beq
\phi_n (t) := \int \xi_n (\mathbf{x}) \phi (\mathbf{x} ,t ) d\mathbf{x}
\eeq 
\beq
\pi_n (t) := \int \xi_n (\mathbf{x}) \pi (\mathbf{x} ,t ) d\mathbf{x}
\label{qf:2}
\eeq 
It follows from (\ref{qf:1}) and (\ref{qf:2}) that 
\beq
[ \phi_n (t) ,\pi_m (t)] = i \delta_{nm}  
\label{qf:3}
\eeq
Note that 
\beq
\pi (\mathbf{x},t) = \dot{\phi} (\mathbf{y},t)
\quad
\to 
\quad
\pi_n (t) = \dot{\phi}_n (t) .
\label{qf:4}
\eeq
If $\phi_n (t) = \phi_n (t)^{\dagger}$ and 
$\pi_n (t) = \pi_n (t)^{\dagger}$ then they both have real eigenvalues
\beq
\phi_n (t) \vert \lambda_\phi \rangle =
\lambda \vert \lambda_\phi \rangle
\qquad
\pi_n (t) \vert \lambda_\pi \rangle =
\lambda \vert \lambda_\pi \rangle
\label{qf:5}
\eeq
To determine the spectrum of $\phi_n(t)$ note that the commutation relations 
(\ref{qf:3}) imply
\beq
\phi_n (t) e^{-i \pi_n (t) \alpha} \vert \lambda \rangle =
(\lambda + \alpha) e^{-i \pi_n (t) \alpha} \phi_n (t) \vert \lambda \rangle
\label{qf:6}
\eeq
This shows that the spectrum of $\phi_n(t)$ is the real line,
since $\alpha$ can be any real number.

We introduce a basis where the $\phi_n (t)$ are multiplication operators
\beq
\pmb{\phi} := (\phi_1(t),\phi_2(t),\cdots )
\label{qf:7}
\eeq
\beq
\pmb{\phi} (t)\vert \mathbf{f} \rangle =
\mathbf{f} \vert \mathbf{f} \rangle
\label{qf:8}
\eeq 
where in principle the vector $\mathbf{f}$ has an infinite number of 
components.  For the moment we assume that the number of degrees of 
freedom are cut off.  We can also do the same for the $\pi_n(t)$:
\beq
\pmb{\pi} := (\pi_1(t),\pi_2(t),\cdots )
\label{qf:9}
\eeq
\beq
\pmb{\pi} (t)\vert \mathbf{g} \rangle =
\mathbf{g} \vert \mathbf{g} \rangle
\label{qf:10}
\eeq 
The Hamiltonian in this discrete representation has the form 
\beq
H = \sum_n {1 \over 2} \pi_n (t) \pi_n (t) +
\sum_{nk} {1 \over 2} \phi_n (t) D_{nk}\phi_k (t)+
\sum_{nk} {m^2 \over 2} \phi_n (t) \phi_n (t) +
\sum_{n_1 \cdots n_4} \Gamma_{n_1 \cdots n_4} 
\phi_{n_1 }(t) \cdots \phi_{n4} (t)
\label{qf:8}
\eeq 
where the coefficients are
\beq
D_{nk} = \int \pmb{\nabla}\chi_n (\mathbf{x}) \cdot
\pmb{\nabla}\chi_k (\mathbf{x})  d\mathbf{x}
\qquad
\Gamma_{n_1 \cdots n_4} = 
\int \chi_{n_1 }(\mathbf{x}) \cdots \chi_{n4} (\mathbf{x}) d\mathbf{x} 
\label{qf:9}
\eeq
Note
\beq
\langle \mathbf{f} \vert \mathbf{g} \rangle
= ({1 \over \sqrt{2 \pi}})^M
e^{i \mathbf{f}\cdot \mathbf{g} }
\label{qf:10}
\eeq
where in principle $M$ can be infinite.  Here w imagine that there is 
always cutoff in $M$. The problem is to learn how to take the limit as
$M\to \infty$.  We can expresse the Hamiltonian in vector form a
\beq
H = {1 \over 2} \pmb{\pi} (t) \cdot \pmb{\pi} (t) +
{1 \over 2} \pmb{\phi}(t)\cdot \mathbb{D}\cdot \pmb{\phi} (t)+
{m^2 \over 2} \pmb{\phi} (t) \cdot \pmb{\phi}(t) +
\Gamma (\pmb{\phi}(t)) 
\label{qf:11}
\eeq 
The time evoluton operator is 
\[
\langle \mathbf{f},t \vert
e^{-iH t} \vert \psi \rangle =
\]
\beq
\langle \mathbf{f},t \vert
(e^{-iH t/N})^N \vert \psi \rangle 
\label{qf:12}
\eeq
Let $\Delta t = t/N$.  Using the explicit form of the Hamiltonian 
\[  
\langle \mathbf{f},t \vert
(e^{-i(
{1 \over 2}\sum_n  \pi_n (t) \pi_n (t) +
{1 \over 2} \sum_{nk}  \phi_n (t) D_{nk}\phi_k (t)+
{m^2 \over 2} \sum_{nk} \phi_n (t) \phi_n (t) +
\sum_{n_1 \cdots n_4} \Gamma_{n_1 \cdots n_4} 
\phi_{n_1 }(t) \cdots \phi_{n4} (t))
\Delta t })^N \vert \psi \rangle .
\]
Next use the trotter prodduct formula to write 
\[  
\lim_{N\to \infty}
\langle \mathbf{f},t \vert [e^{-i \sum_n {1 \over 2} \pi_n (t) \pi_n
(t) \Delta t} e^{-i (\sum_{nk} {1 \over 2} \phi_n (t) D_{nk}\phi_k
(t)+ \sum_{n} {m^2 \over 2} \phi_n (t) \phi_n (t))\Delta t} 
\times 
\]
\[
e^{-i
\sum_{n_1 \cdots n_4} \Gamma_{n_1 \cdots n_4} \phi_{n_1 }(t) \cdots
\phi_{n_4} (t)) \Delta t }]^N \vert \psi \rangle =
\]
Next insert complete sets of intermediate states
\[
\int \langle \mathbf{f}_0 \vert \mathbf{g}_1 \rangle
d \mathbf{g}_1  d \mathbf{f}_1
e^{- {i \over 2} \mathbf{g}_1 \cdot 
\mathbf{g}_1 \Delta t} 
\langle \mathbf{g}_1 \vert \mathbf{f}_1 \rangle
e^{- {i \over 2} (\mathbf{f}_1\cdot \mathbb{D} \cdot
\mathbf{f}_1 + m^2 \mathbf{f}_1 \cdot \mathbf{f}_1 )\Delta t}
e^{-i
\Gamma (\mathbf{f}_1)  \Delta t } \cdots 
\]
\beq
\langle \mathbf{f}_{N-1} \vert \mathbf{g}_N \rangle
d \mathbf{g}_N  d \mathbf{f}_N
e^{- {i \over 2} \mathbf{g}_N \cdot 
\mathbf{g}_N \Delta t} 
\langle \mathbf{g}_N\vert \mathbf{f}_N \rangle
e^{-i {1 \over 2} (\mathbf{f}_N\cdot \mathbb{D} \cdot
\mathbf{f}_N + m^2 \mathbf{f}_N \cdots \mathbf{f}_N \Delta t}
e^{-i
\sum_{n_1 \cdots n_4} \Gamma ( \mathbf{f}_N) \Delta t } 
\psi (\mathbf{f}_N =
\label{qf:x}
\eeq
\[
({ 1 \over \sqrt{2 \pi}})^M \int 
e^{i (\mathbf{f}_0 - \mathbf{f}_1)\cdot \mathbf{g}_1}
d \mathbf{g}_1  d \mathbf{f}_1
e^{- {i \over 2} \mathbf{g}_1 \cdot 
\mathbf{g}_1 \Delta t} \times
\]
\[
e^{- {i \over 2} (\mathbf{f}_1\cdot \mathbb{D} \cdot
\mathbf{f}_1 + m^2 \mathbf{f}_1 \cdot \mathbf{f}_1 )\Delta t}
e^{-i
\Gamma (\mathbf{f}_1) \Delta t } \cdots 
\]
\beq
e^{i (\mathbf{f}_{N-1} - \mathbf{f}_N)\cdot \mathbf{g}_N}
d \mathbf{g}_N  d \mathbf{f}_N
e^{- {i \over 2} \mathbf{g}_N \cdot 
\mathbf{g}_N \Delta t} 
e^{- {i \over 2} (\mathbf{f}_{N}\cdot \mathbb{D} \cdot
\mathbf{f}_N + m^2 \mathbf{f}_{N} \cdots \mathbf{f}_N )\Delta t}
e^{-i
\Gamma (\mathbf{f}_{N }\Delta t}
\psi (\mathbf{f}_N)
\label{qf:x}
\eeq
As in the one diemnsional we can perform the integrals over 
$\mathbf{g}_i$ to get
\[
({ 1 \over \sqrt{i \Delta t}})^{MN} \int 
d \mathbf{f}_1
e^{ {i \over 2\Delta t } (\mathbf{f}_- - \mathbf{f}_1)^2 } 
e^{- {i \over 2} (\mathbf{f}_1\cdot \mathbb{D} \cdot
\mathbf{f}_1 + m^2 \mathbf{f}_1 \cdot \mathbf{f}_1 )\Delta t}
e^{-i
\Gamma (\mathbf{f}_1) \Delta t } \cdots 
\]
\beq
d \mathbf{f}_N
e^{ {i \over 2\Delta t} (\mathbf{f}_N - \mathbf{f}_{N-1})^2  } 
e^{- {i \over 2} (\mathbf{f}_{N}\cdot \mathbb{D} \cdot
\mathbf{f}_N + m^2 \mathbf{f}_{N} \cdots \mathbf{f}_N )\Delta t}
e^{-i
\Gamma (\mathbf{f}_{N }\Delta t}
\psi (\mathbf{f}_N)
\label{qf:x}
\eeq
which we write as
\beq
\lim_{n \to \infty}
({ 1 \over \sqrt{i \Delta t}})^{MN} \int 
\prod_{j=1}^N d\mathbf{f}_j
e^{ {i \over 2\Delta t }\sum_n  (\mathbf{f}_{n-1} - \mathbf{f}_n)^2 } 
e^{- {i \over 2} \sum_n (\mathbf{f}_n\cdot \mathbb{D} \cdot
\mathbf{f}_n + m^2 \mathbf{f}_n \cdot \mathbf{f}_n )\Delta t}
e^{-i \sum_n
\Gamma (\mathbf{f}_n) \Delta t }
\psi (\mathbf{f}_N)
\label{qf:x}
\eeq
The remaining integral involves a quadradic form in the $\mathbf{f}_n$ 
also with an interaction term.

The generalization of the DNP integral is to replace the 
$NM$ integrals over the $M$ field amplitudes at each of the $N$ 
times by integrals over $K$ sub-intervals.  The means that there are
intergrals over $NMK$ boxes.

Following what we did before we only integrate over the quadratic form.
These integrals can be normalized to represent complex probabilities.
The normalization integral can be done analytically.  What remains is 
a normalization factor, a sum over complex probabilities that 
each of the $M$ field degrees of freedom lies in a cylinder set
and this is multiplies by complex weight function
\beq
e^{-i \sum_n
\Gamma (\mathbf{f}^*_n) \Delta t }
\psi (\mathbf{f}^*_N)
\eeq
where the $*$ indicates an evaluation at a point in each the cylinder set.

The next step is to consider the corresponding computation of the
transfer matrix.  This be made to work like the one dimensional
case if we include the non-diagonal part of $\mathbb{D}$ in the 
interaction term.  This means that we used the generalized momenta
and mass terms to define the 1 step probility. While the 
one step matrix is $MK \times MK$ in the limit that both numbers
become infinite,  multiple steps involve matrix multiplication. 

This will be developed later.  The point is that with this choice 
we are not limited to trivial interactions, however there are 
serious issues with how to renormalize the system.    
  

\section{Appendix: Asymptotic  wave packets}

A natural choice for the free wave packets $\vert \psi_{i,0}(t)
\rangle$ that appear in (\ref{s.2}) are minimal uncertainty states
with narrow momentum distributions.  These are Gaussian wave packets
or coherent states.  For a given small momentum spread these have the
smallest spatial distribution consistent with the constraints imposed
by the uncertainty principle.  It is desirable to make both of these
dimensions small; the small momentum spread is needed to get a
wave-packet-independent phase shift, while a narrow wave packet in
coordinate space minimizes the ``time'' $s$ needed for convergence.
The uncertainty principle implies that decreasing one of these
increases the other, which justifies our preference for minimal
uncertainty states.  The other important advantage of minimal
uncertainty states is because they are Gaussian, they can be computed
analytically.

Free minimal uncertainty wave packets in $n$ dimensions are products of 
$n$ one-dimensional minimal uncertainty wave packets.  Multiparticle 
asymptotic states can be taken to be products of single-particle minimal 
uncertainty states. 

The one-dimensional a minimal uncertainty wave packet with mean
momentum $p_0$, mean position $0$ at time $t=0$ and uncertainty
$\Delta p$ in $p$ and $\Delta x= {1 \over 2 \Delta p}$ in $x$, has the
form
\beq
\langle x\vert \psi_{0,p_0} (t) \rangle = 
\langle x \vert e^{-i H_0 t} \vert \psi_{0,p_0} (0) \rangle =
\int \langle x \vert p \rangle e^{-i p^2 t/ 2\mu} dp \langle p 
\vert \psi_{0,p_0} (0) \rangle
\label{aw.1}
\eeq
where 
\beq  
\langle p \vert \psi_{0,p_0} (0) \rangle = (2 \pi)^{-1/4} (\Delta p)^{-1/2}
e^{- {(p-p_0)^2 \over 4 (\Delta p)^2}}.
\label{aw.2}
\eeq
Performing the momentum integral gives (up to phase)
\beq
\langle x\vert \psi_{0,p_0} (t) \rangle = 
(2 \pi)^{-1/4}
( 2\Delta p)^{1/2}(1  + 2i { (\Delta p)^2 t \over  \mu})^{-1/2}
\exp (-
{(\Delta p)^2 \over 1  + {4 (\Delta p)^4 t^2 \over  \mu^2}}
(x-{p_0t\over \mu})^2 )
\label{aw.3}
\eeq
This expression
shows that the center of this wave packet is moving 
with speed $p_0/\mu$ to the right.  It is centered at
the origin at $t=0$.

In three dimensions the corresponding minimal uncertainty
wave packet is the product of three such factors,
\beq
\langle \mathbf{x} \vert \psi (t)  \rangle =
(2 \pi)^{-3/4}
( 2\Delta p)^{3/2}(1  + 2i { (\Delta p)^2 t \over  \mu})^{-3/2}
\exp (-
{(\Delta p)^2 \over 1  + {4 (\Delta p)^4 t^2 \over  \mu^2}}
(\mathbf{x}-{\mathbf{p}_0t\over \mu})^2)
\label{aw.4}
\eeq
Note that while the minimal uncertainty states are computationally
useful, the phase shifts extracted from (\ref{s.7}) are independent of
the wave packet provided it is sufficiently narrow.

Two quantities are needed in order to calculate the
phase shifts using (\ref{s.7}).  One of these quantities are the
integrals of these free wave packets over finite or semi-infinite
space intervals.  These integrals can be expressed in terms of error
functions
\beq
\int_{x_1}^{x_2} N e^{-a(x-vt)^2}dx =
N{\sqrt{a\pi}\over 2}
\left (\mbox{erf} (\sqrt{a} (x_2-vt) )- \mbox{erf} (\sqrt{a} (x_1-vt) )
\right )
\label{aw.5}
\eeq

\beq
\int_{x_1}^{\infty} N e^{-a(x-vt)^2}dx =
N{\sqrt{a\pi}\over 2}
\left (1- \mbox{erf} (\sqrt{a} (x_1-vt) ) =
\right )
= {N\sqrt{a}}{\sqrt{\pi}\over 2}
\mbox{erfc} (\sqrt{a} (x_1-vt) ) 
\label{aw.6}
\eeq

\beq
\int_{-\infty}^{x_2} N e^{-a(x-vt)^2}dx =
N{\sqrt{a\pi}\over 2}
\left (\mbox{erf} (\sqrt{a} (x_2-vt) )- 1 )
\right )
=- N{\sqrt{a\pi}\over 2}
\mbox{erfc} (\sqrt{a} (x_2-vt) )
\label{aw.7}
\eeq
where
\beq
N=(2 \pi)^{-1/4}
( 2\Delta p)^{1/2}(1  + 2i { (\Delta p)^2 t \over  \mu})^{-1/2}
\label{aw.8}
\eeq
\beq
a={(\Delta p)^2 \over 1  + {4 (\Delta p)^4 t^2 \over  \mu^2}}
\qquad
\mathbf{v}={\mathbf{p}_0t\over \mu}
\label{aw.9}
\eeq

The other factor needed is the denominator of (\ref{s.7}).
For these states this can be expressed as
\[
\langle \psi_{0,p_f} (0)
\vert \delta ({p^2 - p_0^2\over 2 \mu}) 
\vert \psi_{0,p_i} (0) \rangle =
\]
\[
\int dp
(2 \pi)^{-1/4} (\Delta p_f)^{-1/2}
e^{- {(p-p_{0f})^2 \over 4 (\Delta p_f)^2}}
{\mu \over p_0} (\delta (p-p_0) + \delta (p+p_0))   
(2 \pi)^{-1/4} (\Delta p_i)^{-1/2}
e^{- {(p-p_{0i})^2 \over 4 (\Delta p_i)^2}} =
\]
\beq
{\mu \over p_0} (2 \pi)^{-1/2} (\Delta p_f \Delta p_i)^{-1/2}
\left (e^{- {(p_0-p_{0f})^2 \over 4 (\Delta p_f)^2} -
{(p_0-p_{0i})^2 \over 4 (\Delta p_i)^2}} +
e^{- {(p_0+p_{0f})^2 \over 4 (\Delta p_f)^2} -
{(p_0+p_{0i})^2 \over 4 (\Delta p_i)^2}}
\right )
\label{aw.10}
\eeq
This is the one dimensional form.  This can also e computed analytically 
in three dimensions.  In this case the integral 
becomes
\[
\langle \psi_{0,p_f} (0)
\vert \delta ({p^2 - p_0^2\over 2 \mu}) 
\vert \psi_{0,p_i} (0) \rangle =
\]
\beq
{\mu \over w} (2 \pi)^{-1/2} (\Delta p_f \Delta p_i)^{-3/2}  
e^{-p_0^2 ({1 \over 4 (\Delta p_f)^2} + {1\over 4 (\Delta p_i)^2})}
e^{- ({p_{0f}^2 \over 4 (\Delta p_f)^2} + {p_{0i}^2\over 4 (\Delta p_i)^2})}
\sin ({p_0v \over 2})
\label{aw.11}
\eeq
where
\beq
w = \vert ({\mathbf{p}_{0f} \over (\Delta p_f)^2} + 
{{\mathbf{p}_{0i}\over (\Delta p_i)^2})} \vert
\label{aw.12}
\eeq
In three dimensions the integral is 
\[
\int 
(2 \pi)^{-3/4} (\Delta q)^{-3/2}
e^{- {(\mathbf{p}-\mathbf{p}_{0})^2 \over 4 (\Delta p)^2}}
d\mathbf{p} 
\delta ({\mathbf{p}^2 \over 2 \mu} - {\mathbf{q}^2 \over 2 \mu} )
d\mathbf{q}
(2 \pi)^{-3/4} (\Delta q)^{-3/2}
e^{- {(\mathbf{q}-\mathbf{q}_{0})^2 \over 4 (\Delta q)^2}}= 
\]
\[
\int 
(2 \pi)^{-3/4} (\Delta q)^{-3/2}
e^{- {(p^2 - 2 p p_0 u_1 + p_0^2) \over 4 (\Delta p)^2}}
du_1 d\phi_1  du_2 d \phi_2 p^4 dp dq {\mu \over p} \delta (p-q)
(2 \pi)^{-3/4} (\Delta q)^{-3/2}
e^{- {(p^2 - 2 {p}{q}_{0}u_2 + q_0^2 ) \over 4 (\Delta q)^2}}=
\]
\[
\int 
\mu (2 \pi)^{2-3/2} (\Delta q \Delta q )^{-3/2}  
e^{- {(p^2 + p_0^2) \over 4 (\Delta p)^2}}
e^{- {(p^2 + q_0^2) \over 4 (\Delta q)^2}}
p^3 dp 
{2 (\Delta p)^2 \over p p_0} 
(e^{pp_0 \over 2 (\Delta p)^2} - 
 e^{-pp_0 \over 2 (\Delta p)^2})  
{2 (\Delta q)^2 \over p q_0} 
(e^{pq_0 \over 2 (\Delta q)^2} - 
 e^{-pq_0 \over 2 (\Delta q)^2})  =
\]
\[
\int 
\mu (2 \pi)^{2-3/2} (\Delta p \Delta q )^{-3/2}  
{2 (\Delta p)^2 \over  p_0} 
{2 (\Delta q)^2 \over  q_0}  
e^{- {p_0^2 \over 4 (\Delta p)^2}}
e^{- {p_0^2 \over 4 (\Delta q)^2}}
\int 
e^{- p^2 ({1 \over 4 (\Delta p)^2} + {1 \over 4 (\Delta q)^2})}
p dp 
(e^{pp_0 \over 2 (\Delta p)^2} - 
e^{-pp_0 \over 2 (\Delta p)^2})  
(e^{pq_0 \over 2 (\Delta q)^2} - 
e^{-pq_0 \over 2 (\Delta q)^2})  
\]
The remaining inegrals are all Gaussian of the general form 
\[
\int e^{-ap^2 +bp} d dp = \sqrt{\pi \over a} {b \over 2a} e^{b^2 \over 4a} 
\]
where
\[
a= ({1 \over 4 (\Delta p)^2} + {1 \over 4 (\Delta q)^2})
\]
\[
b_1= {p_0 \over 2 (\Delta p)^2} + {q_0 \over 2 (\Delta q)^2}
\]
\[
b_2= -1({p_0 \over 2 (\Delta p)^2} + {q_0 \over 2 (\Delta q)^2})=- b1
\]
\[
b_3= {p_0 \over 2 (\Delta p)^2} - {q_0 \over 2 (\Delta q)^2}
\]
\[
b_4= -{p_0 \over 2 (\Delta p)^2} + {q_0 \over 2 (\Delta q)^2}= - b4
\]
So the full integral becomes
\[
\int 
(2 \pi)^{-3/4} (\Delta q)^{-3/2}
e^{- {(\mathbf{p}-\mathbf{p}_{0})^2 \over 4 (\Delta p)^2}}
d\mathbf{p} 
\delta ({\mathbf{p}^2 \over 2 \mu} - {\mathbf{q}^2 \over 2 \mu} )
d\mathbf{q}
(2 \pi)^{-3/4} (\Delta q)^{-3/2}
e^{- {(\mathbf{q}-\mathbf{q}_{0})^2 \over 4 (\Delta q)^2}}= 
\]
\[
{2 \mu (2 \pi)^{1/2} (\Delta p \Delta q )^{1/2}  
\over  p_0 q_0} \sqrt{\pi \over a^3} 
e^{- {p_0^2 \over 4 (\Delta p)^2}}
e^{- {p_0^2 \over 4 (\Delta q)^2}}
((b_1 e^{b_1^2 \over 4a} -e^{-b_1^2 \over 4a}) -
b_3  (e^{b_3^2 \over 4a}- e^{-b_3^2 \over 4a}))
\]
\[
{4 \mu (2 \pi)^{1/2} (\Delta p \Delta q )^{1/2}  
\over  p_0 q_0} \sqrt{\pi \over a^3} 
e^{- {p_0^2 \over 4 (\Delta p)^2}}
e^{- {p_0^2 \over 4 (\Delta q)^2}}
((b_1 \sinh ({b_1^2 \over 4a}) -
b_3 \sinh ({b_3^2 \over 4a}))
\]






{\color{blue}
\bigskip  
\hrule
\begin{center}
{\bf Computation of the integrals - not for publication:}
\end{center}
\bigskip
\[  
\langle p \vert \psi \rangle = (2 \pi)^{-1/4} (\Delta p)^{-1/2}
e^{- {(p-p_0)^2 \over 4 (\Delta p)^2}}
\]  
\[
\Delta p \Delta x = {1 \over 2}
\]
\[  
\langle x \vert \psi \rangle = (2 \pi)^{-1/4} (\Delta x)^{-1/2}
e^{i p_0 x - {x^2 \over 4 (\Delta x)^2}}
\]  
\[
\langle x \vert \psi (t)  \rangle =
\]
\[
\int (2 \pi)^{-1/2} e^{ipx - i{p^2 \over 2 \mu}t}
(2 \pi)^{-1/4} (\Delta p)^{-1/2}
e^{- {(p-p_0)^2 \over 4 (\Delta p)^2}} dp =
\]
\[
\int (2 \pi)^{-1/2} (2 \pi)^{-1/4} (\Delta p)^{-1/2}
e^{ipx - i{p^2 \over 2 \mu}t - {(p-p_0)^2 \over 4 (\Delta p)^2}} dp =
\]
\[
(2 \pi)^{-1/2} (2 \pi)^{-1/4} (\Delta p)^{-1/2}
\int
e^{
-p^2 {1 \over 4 (\Delta p)^2} (1  + i {2 (\Delta p)^2 t \over  \mu} )+
p ({p_0 \over 2  (\Delta p)^2} +ix )
- {p_0^2 \over 4 (\Delta p)^2}} = 
\]
This is gaussian
\[
\int e^{-ap^2 +bp} = \sqrt{{\pi \over a}}e^{b^2 \over 4a}
\]
\[
a={1 \over 4 (\Delta p)^2} (1  + i {2 (\Delta p)^2 t \over  \mu} )
\]
\[
b=({p_0 \over 2  (\Delta p)^2} +ix ) 
\]
\[
{1 \over a} = {4 (\Delta p)^2 \over 1  + i {2 (\Delta p)^2 t \over  \mu}}
\]
\[
b^2 = -x^2 + ix {p_0 \over  (\Delta p)^2} + {p_0^2 \over 4  (\Delta p)^4}
\]
\[
{b^2 \over 4a} =
\]
\[
{(\Delta p)^2 \over 1  + i {2 (\Delta p)^2 t \over  \mu}}
\left (-x^2 + ix {p_0 \over   (\Delta p)^2} + {p_0^2 \over 4  (\Delta p)^4}
\right ) =  
\]
\[
{(\Delta p)^2 \over 1  + {4 (\Delta p)^4 t^2 \over  \mu^2}}
\left (-x^2 + { t \over  \mu}2x p_0
+ {p_0^2 \over 4  (\Delta p)^4}
+ i(x {p_0 \over  (\Delta p)^2} 
+  x^2 {2 (\Delta p)^2 t \over  \mu}
- { t \over  \mu}{p_0^2 \over 2  (\Delta p)^2})
\right ) =
\]
=  
\[
{(\Delta p)^2 \over 1  + {4 (\Delta p)^4 t^2 \over  \mu^2}}
\left (-(x-{p_0t\over \mu})^2 + {p_0^2 t^2 \over \mu^2} 
+ {p_0^2 \over 4  (\Delta p)^4}
+ i(2x {p_0 \over 2  (\Delta p)^2} 
+  x^2 {2 (\Delta p)^2 t \over  \mu}
- { t \over  \mu}{p_0^2 \over 2  (\Delta p)^2})
\right ) =
\]
\[
{(\Delta p)^2 \over 1  + {4 (\Delta p)^4 t^2 \over  \mu^2}}
\left (-(x-{p_0t\over \mu})^2 
+ i(2x {p_0 \over 2  (\Delta p)^2} 
+  x^2 {2 (\Delta p)^2 t \over  \mu}
- { t \over  \mu}{p_0^2 \over 2  (\Delta p)^2})
\right )+
\]
\[
{(\Delta p)^2 \over 1  + {4 (\Delta p)^4 t^2 \over  \mu^2}}
{p_0^2 \over 4  (\Delta p)^4} (1+ {4 (\Delta p)^4 t^2 \over \mu^2}) =
\]
\[
{(\Delta p)^2 \over 1  + {4 (\Delta p)^4 t^2 \over  \mu^2}}
\left (-(x-{p_0t\over \mu})^2 
+ i(2x {p_0 \over 2  (\Delta p)^2} 
+  x^2 {2 (\Delta p)^2 t \over  \mu}
- { t \over  \mu}{p_0^2 \over 2  (\Delta p)^2})
\right )
\]
\[
+ {p_0^2 \over 4  (\Delta p)^2} 
\]
With this evaluation
\[
\langle x \vert \psi (t) \rangle =
\]
\[
(2 \pi)^{-1/2} (2 \pi)^{-1/4} (\Delta p)^{-1/2}
e^{- {p_0^2 \over 4 (\Delta p)^2}}
e^{p_0^2 \over 4  (\Delta p)^2}
2(\Delta p)\sqrt{{\pi  \over 1  + i {2 (\Delta p)^2 t \over  \mu}}}
\times
\]
\[
e^{
{(\Delta p)^2 \over 1  + {4 (\Delta p)^4 t^2 \over  \mu^2}}
\left (-(x-{p_0t\over \mu})^2
+ i(2x {p_0 \over 2  (\Delta p)^2} 
+  x^2 {2 (\Delta p)^2 t \over  \mu}
- { t \over  \mu}{p_0^2 \over 2  (\Delta p)^2})
\right )} =
\]
\[
(2 \pi)^{-3/4}
2(\Delta p)^{1/2}\sqrt{{\pi  \over 1  + i {2 (\Delta p)^2 t \over  \mu}}}
\times
\]
\[
e^{
{(\Delta p)^2 \over 1  + {4 (\Delta p)^4 t^2 \over  \mu^2}}
\left (-(x-{p_0t\over \mu})^2
+ i(2x {p_0 \over 2  (\Delta p)^2} 
+  x^2 {2 (\Delta p)^2 t \over  \mu}
- { t \over  \mu}{p_0^2 \over 2  (\Delta p)^2})
\right )} =
\]
\[
(2 \pi)^{-1/4}
( 2\Delta p)^{1/2}\sqrt{{1  \over 1  + i {2 (\Delta p)^2 t \over  \mu}}}
\times
\]
\[
e^{-
{(\Delta p)^2 \over 1  + {4 (\Delta p)^4 t^2 \over  \mu^2}}
(x-{p_0t\over \mu})^2}
e^{i{1 \over 1  + {4 (\Delta p)^4 t^2 \over  \mu^2}}
(x p_0   
+ 2 x^2 (\Delta p)^4 {t \over  \mu}
- { t \over  \mu}{p_0^2 \over 2 })}
\]
The phase terms are unimportant - If we ignore them
\[
(2 \pi)^{-1/4}
( 2\Delta p)^{1/2}\sqrt{{1  \over 1  + i {2 (\Delta p)^2 t \over  \mu}}}
\times
\]
\[
e^{-
{(\Delta p)^2 \over 1  + {4 (\Delta p)^4 t^2 \over  \mu^2}}
(x-{p_0t\over \mu})^2}
\]
We also note that $(2 \Delta p)^{1/2} = (\Delta x)^{-1/2}$ 

Classically
\[
{\mu \over p_0}(2 \Delta x+ R) =
{\mu \over p_0}({1 \over \Delta p}+ R) =s
\]
is a lower bound on the scattering time, $s$, where $R$ is the range of 
the potential.

In three dimensions the above becomes
\[
\langle \mathbf{x} \vert \psi (t)  \rangle =
\]
\[
(2 \pi)^{-3/4}
( 2\Delta p)^{3/2}(\sqrt{{1  \over 1  + i {2 (\Delta p)^2 t \over  \mu}}})^3
\times
\]
\[
e^{-
{(\Delta p)^2 \over 1  + {4 (\Delta p)^4 t^2 \over  \mu^2}}
(\mathbf{x}-{\mathbf{p}_0t\over \mu})^2}
\]
These state have mean momentum $\mathbf{p}_0$ mean coordinate
at time $t=0$:
$\langle \mathbf{x}(t=0) \rangle =0$, and undertainty in
momentum $\Delta p$ and undertainty in coordinate
$\Delta x = {1 \over 2 \Delta p}$
 
And for multiparticle scattering we use symmetrized products of the
above with different $\mathbf{p}_0$
\hrule
}  

\bigskip

\hrule

\bigskip




This avoids  both complex numerical quadratures
and the  need to directly  sum over all  cylinder sets.  In  this case
each matrix  multiplication treats all  of the intervals at  each time
step  in parallel.   In this  way  the final  result is  given by  the
$N-{th}$ power of  an $M+1 \times M+1$ matrix.  Thus (\ref{}) becomes
\beq
\langle \psi_{f0} (s)  \vert e^{-i2 H s} \vert \psi_{i0} (-s) \rangle \approx
\sum_{mk}
\int_{I_{m}}
\psi_{i0} (x,-s)dx (M^N)_{mk}
\psi_{i0} (y_{mN},-s)
\label{f.8}
\eeq
Which reduces the computation of the path integral to the
computation powers of a matrix.  It is instructive to summarize the
input needed to compute the phase shift.  It includes
\begin{itemize}
\item [1.] A finite partition on the real line into M+1 sub intervals
where the potential is approximately constant on each one.

\item [2.] Values of the potential at M+1 points, one in each of
the intervals in [1.]

\item [3.] A fixed analytically computable $(M+1)\times (M+1)$ matrix.  

\item [4.] The ability to compute large powers of the matrix.

\item [5.] The choice of partition depens on the structure of the 
potential.  On simple treatment is to divide the range of the 
potential into a finite number of evenlly space intervals.
Then the break points $x_n$ can be chosen to satisfy $y_n= V(x_n)$
where the $y_n$ are breakpoints in the range of $V$.   
\end{itemize} 

Computation of powers of a matrix can be computed efficiently by many
different methods.  This still retains the complex probability
interpretation, except now these are approximate complex
probabilities.

The matrix $P(y_{m},I_{k})$ has one of the following three analytic forms
depending on whether the intervals are finite or half infinite: 
\[
P(y_{m},I_{k}) = 
\sqrt{{1  \over i \pi }}
\left (
\sqrt{\pi \over 2} (C_{c}(\sqrt{{\mu \over 2\Delta t}}(x_k-y_m))- 
C_{c}(\sqrt{{m \over 2\Delta t}}(
x_{k+1}-y_m)))
\right . 
\]
\beq
\left .
+i \sqrt{\pi \over 2} (S_{c}(\sqrt{{\mu \over 2\Delta t}}(x_k-y_m))- S_{c}(\sqrt{{m \over 2\Delta t}}(
x_{k+1}-y_m))) .
\right ) 
\label{f.9}
\eeq
For $x_{k}=-\infty (k=0)$:
\beq
P(y_{m},I_{k}) =
\sqrt{{1  \over i \pi }}
{1 \over 2} \int_{-\infty}^{\infty} e^{ix^2} dx - \int_b^0 e^{ix^2} dx =
\sqrt{{\pi \over 2}} \left ({1+i \over 2} + C_c(\sqrt{{\mu \over 2\Delta t}}(x_k-y_m))+
i S_c(\sqrt{{\mu \over 2\Delta t}}(x_k-y_m)) \right ) .
\label{f.10}
\eeq
For $x_{k+1}=\infty (k=M)$:
\beq
\sqrt{{1  \over i \pi }}
{1 \over 2} \int_{-\infty}^{\infty} e^{ix^2} dx - \int_0^a e^{ix^2} dx =
\sqrt{{\pi \over 2}} \left ({1+i \over 2} - C_c(\sqrt{{m \over 2\Delta t}}(
x_{k+1}-y_m)))-i S_c(\sqrt{{m \over 2\Delta t}}(
x_{k+1}-y_m))) 
\right )
\label{f.11}
\eeq
In the next section we compute analytic expressions for (\ref{f.6}). 


\end{document}

