%\documentclass[aps,prc,preprint,showpacs,groupedaddress]{revtex4-1}
\documentclass[aps,prc,reprint,showpacs,groupedaddress,onecolumn]{revtex4}
\usepackage{bm}% bold math
\usepackage{amsmath}
\usepackage{amsfonts}
\usepackage{amscd} 
\usepackage{graphicx}
\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 Henstock 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}

We investigate the possibility of using real-time path integrals
to calculate scattering observables.  The computational method is 
a numerical implementation of an 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. 
This investigation demonstrates the success of the method in a simple 
test model.  How far it can be extended remains an open question.

\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 continuous paths with fixed starting and
end points.

%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 when the
gauge function, $\delta (x)$, is a constant.
 
One advantage is that the integral 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 infinite dimensionals integrals, which are
inductive limits of finite dimensional integrals over cylinder
sets that have a well-defined volume and correlated evaluation points and
times.  Precise definitions for the interested reader can be found in
\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.

This motivated the replacement of a probability theory based on
countably additive positive measures by one based on finitely additive
generalized Riemann sums.  The two formulations of probability
coincide for real probabilities, but the Henstock formulation of
probability extends to non-positive and complex probabilities.  This
was further developed by P. Muldowney.  He developed a suggestion of
Henstock that the Feynman path integral could be understood using this
framework.  Muldowney \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 time-evolution
group is needed 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 the interpretation of
the path integral as the expectation of functionals of continuous
paths with respect to a complex probability distribution.

The approach taken in this paper is a numerical implementation of the
formulation based on computing 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 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
computing a many-dimensional oscillating integral.  The computational
method used in this work computes the transition matrix elements in a
straightforward, but inefficient way.  While the efficiency of the
method can be significantly improved, it remains to be seen whether
the method can be scaled up to problems involving many particles or
fields.

This purpose of this work is only to demonstrate that it is possible
to perform scattering calculations using this method.  An appealing
property of the 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 paths is not lost in
the numerical implementation.  In particular the computational
algorithm can be deconstructed to find the approximate 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.  
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 and can be performed
analytically.  What remains is the 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 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 continuous path passes through at each time slice.
  The limit that the width of the finite intervals vanish, and the finite
  endpoints of the semi-infinite intervals become infinite is
  eventually taken.  For computational purposes the intervals should
  be sufficiently small so that 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$ intervals, $\{I_{n_1}, \dots I_{n_N}\}$; one interval for
  each time slice.  Each of sequence of intervals defines a cylinder
  set.  A cylinder set is an ordered set of $N$ windows
  at $N$ intermediate times.  Every continuous path goes through a
  unique cylinder set.  It is elementary to show that the sum of these
  integrals over all possible N-fold products of intervals (cylinder
  sets) is 1.  This allows them 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 formula.

\end{itemize}

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

There are two problems that must be overcome to make this into a
computational method.  They can be summarized by noting that (1) it is
not clear that the complex probabilities can be computed analytically;
they involve $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 expectation value of a
potential functional of continuous paths with respect to a complex
probability distribution of cylinder sets of continuous paths.

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

The paper is organized as follows.  The next section includes a brief 
discussion of the scattering formalism that will be used in the 
rest of the paper.  Section three discusses the conventional
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. 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 a Gaussian potential.  The goal is to calculate
sharp-momentum transition matrix elements using a path-integral
treatment of time-dependent scattering.  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
which can be written as
\beq
\lim_{t \to \infty} \Vert e^{-iHt} \vert \psi_f(0) \rangle - 
e^{-iH_0t}\vert \psi_{f0}(t)
\rangle \Vert =0 
\qquad
\lim_{t \to -\infty} \Vert e^{-iH t} \vert \psi_i(t) \rangle - 
e^{-iH_0 t}\vert \psi_{i0}(t)
\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
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 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
problem.

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 primary quantities of 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 eigenstate, $\vert \mathbf{p}_i \rangle$, is
replaced by normalizable a wave packet.

A useful approximation is to use a Gaussian wave packet, centered 
about the initial momentum, $\mathbf{p}_i$, with a delta-function normalization:
\beq
\int d\mathbf{p} \psi_{i0}(\mathbf{p}) =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 wave packet.  For
computational purposes, the width of the wave packet 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 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

\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  
\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}) 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_f \vert V
\left (e^{-iH{t/N}}\right )^N  \vert
\psi_{i0} (-t) \rangle . 
\label{fs.6}
\eeq
In the ``$x$'' representation
equation (\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

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^{-iHs} \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
where $x_0=x$.

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
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 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_{i-1}-x_i\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, 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 interactions that are more
general than quadratic, 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 continuous 
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 continuous 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.  In this case the usual Monte Carlo integration
is replaced by matrix multiplication, which can be performed
efficiently even for large matrices.

Finally, the use of the operator $\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
one expects some cancellations.  Thus,
the relevant parts of the calculation take place in a finite
space-time volume.

The fundamental new idea that is the key to the computational strategy,
proposed by Muldowney, is to decompose the integral over each
$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 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 continuous
paths that pass through $I_{mN}$ at time $t_N$, $ \cdots , I_{m1}$ at
time $t_1$ define a cylinder set of continuous paths.  The right most
(initial) 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 continuous
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 measure,
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 MNJ formulation of the path integral are 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{{m \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{{m \over 2\Delta t}}(x_{k}-y_m)) +
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 . 
\]
\beq
\left .
+i (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))- 
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 ) 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.  

While the above formulas are for one-dimensional scattering, the
generalization to three dimensions is straightforward.  The
three-dimensional case involves products of these formulas.
Transition matrix elements can be extracted from (\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 a
repulsive Gaussian potential in one dimension.

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

In order keep the analysis as simple as possible (1) $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 good 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.  This can also be
controlled if $m \delta x^2 / 2 \Delta t <1 $.  This suggests that
both of these dimensionless quantities should be less than 1.  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 numerical 
complexity.

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 generally requires including more intervals
at each time step.

For the test an initial wave packet with the dimensionless
parameters used in the calculations are listed in table 1.
\begin{table}
\begin{tabular}{|l|l|}
\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}
\caption{initial wave packet parameters}
\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 2.
\begin{table}
\begin{tabular}{|l|l|}
\hline
strength &$ v0=5.0$\\ 
range &$r0=1.0$\\
\hline
\end{tabular}
\caption{potential parameters}
\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= m (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 vary slowly.
 
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.  This figure can be used to estimate the active volume
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}  
These graphs suggest that the active volume is actually 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 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 3 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}
\caption{V(x),$\langle x \vert \psi_{0i}(-3) \rangle$}
\label{figure 3}
\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}
\begin{tabular}{|l|l|l|}
\hline
$x$ & $\sum_n \Re (P_n(x))$ &$\sum_n \Im (P_n(x))$\\
\hline
-25.005001 & p=1.000000e+00 & + i(5.273559e-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.000000e+00 &+ i(-1.387779e-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.000000e+00 & + i(2.775558e-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.000000e+00 &+ i(2.775558e-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.000000e+00 &+ i(7.771561e-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.000000e+00 &+ i(3.608225e-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.000000e+00 &+ i(1.665335e-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.000000e+00 &+ i(5.273559e-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.000000e+00 &+ i(8.049117e-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.000000e+00 &+ i(1.137979e-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.000000e+00 &+ i(-2.164935e-15)\\
\hline                                          
\end{tabular}
\caption{sums of complex probabilities} 
\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 to the exact wave packet at $t=-3$.  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 have reached the asymptotic region by
then.  This test used 30 time steps 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 figures 4 and 5.  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
Figures 4 and 5 show the real and imaginary parts of the exact $t=0$
free wave function and the exact $t=-3.0$ free wave function evolved
back to $t=0$ with 30 time slices.  These calculations used
$\Delta t =.1$ and $\delta x=.01$.  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}).

%Table 4 shows the numerical values of some of the data plotted in these
%graphs

\begin{figure}
\centering
\includegraphics[scale=.5]{fig_4.pdf}
\caption{Real part of $\langle x \vert \psi_{0i}(0) \rangle$
and $\langle x \vert U_0(3) \vert \psi_{0i}(-3) \rangle$ }
\label{figure 4}
\end{figure}  

\begin{figure}
\centering
\includegraphics[scale=.5]{fig_5.pdf}
\caption{Imaginary part of $\langle x \vert \psi_{0i}(0) \rangle$
and $\langle x \vert U_0(3) \vert \psi_{0i}(-3) \rangle$ }
\label{figure 5}
\end{figure}  
 
In these figures both the real and imaginary parts of both wave
functions fall on top of each other.  This suggest 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 matrices.  The resolution was chosen
to be sufficiently fine to represent the initial wave packet.

The next two figures illustrate the effect of the potential on the
evolution of the wave function.  Figure 6 compares the real part 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 figure 4.  The functions are calculated at the midpoint of 
each finite interval.  
Figure 7 shows the corresponding plot for the
imaginary part of the wave function.
Figures 6 and 7 show the change in phase as the interaction is turned
on.  The dashed lines in both figures represent the initial wave packet, while
the dotted lines represent the interacting wave packet.

Note that these calculations approximate the incident scattered
wave function at the collision time $(t=0)$.  This why the wave
functions do not look like separated transmitted and reflected waves.
As discussed in section 2, the exact incident wave function at $t=0$
is sufficient to compute the differential cross section.  These
figures also 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.
\begin{figure}
\centering
\includegraphics[scale=.5]{fig_6.pdf}
\caption{Real part of $\langle x \vert \psi_{0i}(-3) \rangle$
and $\langle x \vert \psi_{i}(-3) \rangle$}
\label{figure 6}
\end{figure}

\begin{figure}
\centering
\includegraphics[scale=.5]{fig_7.pdf}
\caption{Imaginary part of $\langle x \vert \psi_{0i}(-3) \rangle$
and $\langle x \vert \psi_{i}(-3) \rangle$}
\label{figure 7}
\end{figure}  

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$.

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 figures 6 and 7 used a time step
$\Delta_t=.1$.  Figures 8 and 9 compare 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}
\centering
\includegraphics[scale=.5]{fig_8.pdf}
\caption{$\Delta t =.1,.05$}
\label{figure 8}
\end{figure}  

\begin{figure}
\centering
\includegraphics[scale=.5]{fig_9.pdf}
\caption{$\Delta t =.1,.05$}
\label{figure 9}
\end{figure}  

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.  Figure 10 and 11 show 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}
\centering
\includegraphics[scale=.5]{fig_10.pdf}
\caption{$t=3.0, 6.0$}
\label{figure 10}
\end{figure}  

\begin{figure}
\centering
\includegraphics[scale=.5]{fig_11.pdf}
\caption{$t=3.0, 6.0$}
\label{figure 11}
\end{figure}  

Due to the oscillating nature of the wave function, the spatial
resolution of the intervals that define the cylinder sets should be
sufficiently small that the wave function is effectively constant on
them.  Figures 12 and 13 show 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}
\centering
\includegraphics[scale=.5]{fig_12.pdf}
\caption{$\delta_x=.01, .005$}
\label{figure 12}
\end{figure}  

\begin{figure}
\centering
\includegraphics[scale=.5]{fig_13.pdf}
\caption{$\delta_x=.01, .015$}
\label{figure 13}
\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.  Figures 14 and 15 compare 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}
\centering
%\includegraphics[scale=.5]{fig16.pdf}
\includegraphics[scale=.5]{fig_14.pdf}
\caption{$\Re, Vol=\pm 25,\pm 50$}
\label{figure 14}
\end{figure}  

\begin{figure}
\centering
%\includegraphics[scale=.5]{fig17.pdf}
\includegraphics[scale=.5]{fig_15.pdf}
\caption{$\Im, Vol=\pm 25,\pm 50$}
\label{figure 15}
\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.  Figures 16 and 17 compare
calculations where the potential is evaluated at the center or left
endpoint of each interval for $\delta x=.01$ Figure 16 shows the real
part of the wave function, where the sample points are at the left
(dash-dot line) and the center (dashed line) of the interval.  Figure
17 shows the imaginary part.  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.  Figures 18 and 19 repeat the
calculations shown in Figures 16 and 17 by increasing the resolution
by a factor of 2.  These figures show a corresponding reduction in the
difference between the two curves.
\begin{figure}
\centering
%\includegraphics[scale=.5]{fig14a.pdf}
\includegraphics[scale=.5]{fig_16.pdf}
\caption{Re left vs center}
\label{figure 16}
\end{figure}  

\begin{figure}
\centering
%\includegraphics[scale=.5]{fig15a.pdf}
\includegraphics[scale=.5]{fig_17.pdf}
\caption{Im left vs center}
\label{figure 17}
\end{figure}
 
\begin{figure}
\centering
%\includegraphics[scale=.5]{fig14.pdf}
\includegraphics[scale=.5]{fig_18.pdf}
\caption{Re left vs center}
\label{figure 18}
\end{figure}  

\begin{figure}
\centering
%\includegraphics[scale=.5]{fig15.pdf}
\includegraphics[scale=.5]{fig_19.pdf}
\caption{Im left vs center}
\label{figure 19}
\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 
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 20.  The solid curve shows
the Gaussian approximation to the 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}
\caption{Born approximation vs smeared Born approximation}
%\label{figure 23}
\label{figure 20}
\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
21.  The imaginary part of $\langle p \vert V \Omega_0 (V=0,t=3) \vert \Psi_{i0}(-3)\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}(-3)\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}
\caption{Exact smeared Born approximation vs V=0 wave function Born approximation} 
\label{figure 21}
%\label{figure 25}
\end{figure}  

The approximate calculations give
$V(x) \langle x \vert \Omega (t) \vert \psi_{i0}(-t)\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 figures 22 and 23.
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
(dash dot curves) compared to the curves which show the corresponding
quantities that replace the final momentum by a Gaussian approximation
to a delta function (solid 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-7
and 10-19.  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}
\centering
%\includegraphics[scale=.5]{fig21.pdf}
\includegraphics[scale=.5]{fig_22.pdf}
\caption{Real(T) - Fourier transform vs wave packet} 
%\label{figure 27}
\label{figure 22}
\end{figure}  

\begin{figure}
\centering
%\includegraphics[scale=.5]{fig21.pdf}
\includegraphics[scale=.5]{fig_23.pdf}
\caption{Real(T) - Fourier transform vs wave packet} 
%\label{figure 28}
\label{figure 23}
\end{figure}  

Figures 24 and 25 compare the smeared transition matrix elements
(using the direct Fourier transform method - dash-dot curves) to sharp
momentum transition matrix elements computed by numerically solving
the Lippmann-Schwinger equation (dashed curves) using the method in
\cite{Rubtsova}.  The Lippmann-Schwinger calculations are represented
by the dashed curve while the path integral calculations are
represented by the dotted curve.  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 curves in fig
20. suggest that this effect is a few percent.  All of the errors can
be reduced at the expense of a larger calculation.
 
\begin{figure}
\centering
%\includegraphics[scale=.5]{2steps_lc_re.pdf}
\includegraphics[scale=.5]{fig_24.pdf}
\caption{Re(T) exact vs path integral}
%\label{figure 20}
\label{figure 24}
\end{figure}  

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

\section{Summary and Conclusions}
 
This work demonstrated the feasibility of performing scattering
calculations using real-time path integrals.  The calculations were
based on a reinterpretation of the path integral as the expectation
value of a potential functional with respect to a complex probability
distribution on a space of paths.

The space of paths was constructed by dividing the total time into a
large number of small time slices.  At each time slice space was
divided up into a large number of small windows, including two 
semi-infinite windows.  The space of continuous 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.

For the convergence of the calculations the equivalence class of paths
must be sufficiently fine that the potential functional has
approximately the same value for all paths in the same equivalence
class.  In addition, the time steps must be sufficiently small for the
Trotter product formula to converge.

The difficulty is the large number of equivalence classes.  Each
equivalence class can be represented as the intersection 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 is the basis of the
illustrated 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.  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 equivalence class 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 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).  Most of the calculations took a few minutes on a desktop
personal computer.  One of the surprising aspects of these
calculations is the stability of the sums over the complex
probabilities.

The one positive is that in all cases of interest the relevant 
volume is finite; it is limited by the range of the interactions
and size of the wave packets. 

The calculations presented were performed in minutes on a desktop
computer. 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.  No attempt was made to compute powers of the product
of the one step probability matrix with the potential matrix.  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.  None of these potential efficiencies
were utilized. 

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 to 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 the 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 quantity 
\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 interesting question is whether this method can be scaled
to more particles or fields?  Considering the information needed to
perform these one-dimensional calculations, it is expected that
few-body calculations will be feasible, however applications to
many-body systems or fields do not appear to be feasible in the
absence of a significant advance in computational strategy.

Some observations may be relevant in this context.  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}
\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 paramters 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}

