Model Study of Three-Body Forces in the Three-Body Bound State
H. Liu
Ch. Elster
The Faddeev equation for the three-body bound state with two- and three-body forces is solved directly as three-dimensional integral equation. The numerical feasibility and stability of the algorithm, which does not employ partial wave decomposition is demonstrated. The three-body binding energy and the full wave function are calculated with Malfliet-Tjon-type two-body potentials and scalar two-meson exchange three-body forces. For two and three body forces of ranges and strengths typical of nuclear forces the single particle momentum distribution and the two-body correlation function are similar to the ones found for realistic nuclear forces.
3nfbound.tex
\documentstyle[preprint,prc,eqsecnum,aps,epsf]{revtex}

\tightenlines   %for single spacing in preprint

\begin{document}
\draft

\title{Model Study of Three-Body Forces in the Three-Body Bound
State}

\author{H. Liu, Ch. Elster}
\address{
Institute of Nuclear and Particle Physics,  and
Department of Physics, \\ Ohio University, Athens, OH 45701}

\author{ W. Gl\"ockle}
\address{
 Institute for Theoretical Physics II, Ruhr-University Bochum,
D-44780 Bochum, Germany.}

\vspace{10mm}


\date{\today}

\maketitle


\begin{abstract}
The Faddeev equation for the three-body bound state 
with two- and three-body forces is solved directly as three-dimensional integral equation.
The numerical feasibility and stability of the 
algorithm, which does not employ partial wave decomposition is
demonstrated. The
three-body binding energy and the full wave function are calculated with
Malfliet-Tjon-type two-body potentials and scalar two-meson exchange 
three-body forces. 
For two and three body forces of ranges and strengths typical of nuclear forces the single 
particle momentum distribution and the two-body correlation function are similar  to the ones found
for realistic nuclear forces.

\end{abstract}

\vspace{10mm}

\pacs{PACS number(s): 21.45+v}

\pagebreak

%****************************************************************************

%\narrowtext

%****************************************************************************

\section{Introduction}


Three-body Faddeev equations for the bound state using local forces directly
without finite rank expansion have been solved since the pioneering work by 
Malfliet and Tjon \cite{malfliet69-1} and Osborn \cite{osborn}.
After this, a huge body of work by different groups followed, and 
calculations using the Faddeev equations were
either performed in momentum space, see e.g. 
\cite{pickle,stadler,stadler2,nbench}, 
in configuration space \cite {lavern,payne,schell0}, or in a hybrid fashion using
both spaces \cite{wu}.
However, the technique was always based on an angular momentum decomposition,
which including spin and isospin degrees of freedom lead to quite a large set of states. 
For instance allowing NN forces to act in all states of total angular momenta $j$ up to
$j$=6, which
is necessary to control the $^3$H binding energy within 2 keV require 102 angular momentum
and isospin combinations (in the literature often called channels). 
In view of this very
large number of interfering terms and having in mind that computer resources are  
still increasing steadily, it 
appears natural to give up such an expansion and  work directly with momentum
vectors or position vectors in configuration  space as it is common in
Greens Function Monte Carlo methods \cite{gfmc}. 
We started an approach based on momentum vectors as variables for the 
three-boson bound state \cite{elster98-1},  which then could be extended to the
three-boson continuum in a straightforward fashion \cite{scatt1}.
In configuration space Faddeev equations have been solved applying vector variables 
for pure Coulomb bound-state problems, namely the $e^-e^-e^+$ and the
$pp\mu^-$ systems \cite{kvin}. 
Including spin  and isospin degrees of freedom together with realistic nucleon-nucleon 
(NN) forces,
a scheme based on momentum vectors has been applied to 
NN scattering leading to integral
equations in two variables \cite{imam1}, which can easily be solved on present day computers. 


Angular momentum decompositions even for three bosons require  quite a tedious algebra
and care in a computational implementation \cite{gloeckle-book}
due to often nontrivial cancellations between
partial wave terms.
This is especially true if one includes 
three-body forces \cite{coon,hueber97-1}. 
In contrast, the formulation in terms of momentum vectors of the
one Faddeev equation for identical particles is very straightforward and transparent. 
In this paper we extend the approach from Ref.~\cite{elster98-1}
and include three-body forces for an investigation of the three-body bound state.
This study is meant only to demonstrate the ease in solving the Faddeev equation including three-body
forces, and working directly with momentum vector variables. Since we do not yet include spin and
isospin degrees of freedom we do not address physical questions directly related to the physical three
nucleon system, but the numerical examples given can shed some light on the wave function properties
of three nucleons.   

The paper is organized as follows. In Section II the Faddeev equation including two- and
three-body forces is formulated in terms of momentum vectors, and its solution,
especially the
intermediate integrations are layed out in detail. Section III displays our choice of
forces and provides the necessary numerical insight for achieving an accurate solution.
Wave function properties are displayed in Section IV.
Finally we summarize in Section V and provide an outlook.




\section{Three-Body Bound State Equation with Three-Body Force}

The bound state of three identical particles which interact via pairwise
forces $V^i=V_{jk}$ ($i,j,k=1,2,3$ and cyclic permutations thereof) and
a genuine three-body force $V_{123}$ is given by the Schr\"odinger
equation which reads in integral form
\begin{equation}
|\Psi\rangle =
G_{0}(\sum^{3}_{i=1}V^{i}+V_{123})|\Psi\rangle.
\label{eq:2.1}
\end{equation}
Here the free propagator is given by
 $G_{0}=(E-H_{0})^{-1}$, where $H_0$ stands for the free Hamiltonian
and $E$ for the binding energy of the three-body system. Any
three-body force $V_{123}$ 
can be decomposed into three different pieces as
\begin{eqnarray}
V_{123}=\sum^{3}_{i=1}V^{(i)}_{4},
\label{eq:2.2}
\end{eqnarray}
such that $V^{(i)}_{4}$ is symmetric under the exchange of particles
$j$ and $k$ ($j \neq i \neq k$). The decomposition suggested in
Eq.~(\ref{eq:2.2}) is natural for e.g. realistic $\pi\pi$ 3N forces,
which are considered at present in all available 3N forces.
Introducing  Faddeev components
$|\Psi\rangle=\sum\nolimits^{3}_{i=1}|\psi_{i}\rangle$ with
\begin{equation}
|\psi_{i}\rangle = G_{0}(V_{i}+V^{(i)}_{4})|\Psi\rangle 
\label{eq:2.3}
\end{equation}
leads to three coupled integral equations
\begin{eqnarray}
|\psi_{i}\rangle&=&G_{0}t_{i}\sum_{j\neq i}
|\psi_{j}\rangle+(1+G_{0}t_{i})G_{0}V^{(i)}_{4}\sum_{j}
|\psi_{j}\rangle \nonumber \\
&=&G_{0}\{t_{i}\sum_{j\neq i}
|\psi_{j}\rangle+(1+t_{i}G_{0})V^{(i)}_{4}\sum_{j}
|\psi_{j}\rangle\}.
\label{eq:2.4}
\end{eqnarray}
The operator $t_{i}$ describes the two-body $t$-matrix in the
subsystem $jk$. If we consider identical particles (here bosons,
since we are omitting spin), the three-nucleon wave function
$|\Psi\rangle$ has to be totally symmetric. As a consequence, 
the Faddeev components
$|\psi_{i}\rangle$ are
identical in their functional form, only the particles are
permuted. Thus it is sufficient to consider only one component, e.g. 
\begin{equation}
|\psi_1\rangle=G_{0}t_1P|\psi_1\rangle +(1+G_{0}t)V^{(1)}_{4}(1+P)
|\psi_1\rangle
\label{eq:2.5}
\end{equation}
In the following the index (1) will be dropped.
The permutation operator $P$ is given as
$P=P_{12}P_{23}+P_{13}P_{23}$
and the total wave function reads 
\begin{eqnarray}
|\Psi\rangle=(1+P)|\psi\rangle.
\label{eq:2.6}
\end{eqnarray}
In order to solve  Eq.~(\ref{eq:2.5}) standard Jacobi
momenta are used,
\begin{eqnarray}
{\mathbf{p}}_{i}&=&
\frac{1}{2}\left({\mathbf{k}}_{j}-{\mathbf{k}}_{k} \right) \nonumber \\
{\mathbf{q}}_{i}&=&
\frac{2}{3}\left({\mathbf{k}}_{i}-
\frac{1}{2}\left({\mathbf{k}}_{j}+{\mathbf{k}}_{k}
\right) \right),
\label{eq:2.7}
\end{eqnarray}
where $ijk=123$ and cyclic permutations thereof. For later clarification
we  label the
coordinates with $ijk=123$ as system of `type ($1$)', the ones with 
$ijk=231$ as `type ($2$)' and the ones with $ijk=312$ as `type ($3$)'. 
With the Jacobi momenta from Eq.~(\ref{eq:2.7}) and omitting 
the arbitrarily chosen index 1, Eq.~(\ref{eq:2.5}) reads
\begin{eqnarray}
\langle {\mathbf{pq}}|\psi\rangle=\frac{1}{E-\frac{p^{2}}{m}-\frac{3q^{2}}{4m}}
\langle {\mathbf{pq}}|tP + V_{4}(1+P)+tG_{0}V_{4}(1+P)|\psi\rangle.
\label{eq:2.8}
\end{eqnarray}
Introducing the symmetrized two-nucleon $t$-matrix
\begin{eqnarray}
t_{s}({\mathbf{p}},{\mathbf{q}};E)=
t({\mathbf{p}},{\mathbf{q}};E)+t(-{\mathbf{p}},{\mathbf{q}};E)
\label{eq:2.9}
\end{eqnarray}
and explicitly working out the permutation operator $P$ in the 
first term of Eq.~(\ref{eq:2.8}) leads to
\begin{eqnarray}
\langle {\mathbf{pq}}|\psi\rangle&=&\frac{1}{E-\frac{p^{2}}{m}-\frac{3q^{2}}{4m}}
\Biggl [ 
\int d^{3}q' t_{s}\left( {\mathbf{p}},\frac{1}{2}{\mathbf{q}}+{\mathbf{q}}';E-\frac{3}{4m}q^{2} \right)
\langle {\mathbf{q}}+\frac{1}{2}{\mathbf{q}}',{\bf q'}|\psi\rangle  \nonumber \\
&+&
\langle {\mathbf{p}}{\mathbf{q}}|V_{4}(1+P)|\psi\rangle
+\frac{1}{2}\int d^{3} \tilde{p}
\frac{t_{s}\left({\mathbf{p}},\tilde{{\mathbf{p}}};E-\frac{3}{4m}q^{2} \right)}
{E-\frac{\tilde{p}^{2}}{m}-\frac{3}{4m}q^{2}}
\langle {\tilde{\mathbf{p}}}{\mathbf{q}}|V_{4}(1+P)|\psi\rangle
 \Biggr ] \nonumber \\
\label{eq:2.10}
\end{eqnarray}
A three-body force (3BF) with two  scalar meson exchanges and a 
constant meson-nucleon amplitude
can be written in the form of  Eq.~(\ref{eq:2.2}) with
\begin{eqnarray}
V_4 \equiv V_4^{(1)}\propto\frac{F(Q^{2})}{Q^{2}+m^{2}_{s}} \: 
\frac{F(Q'^{2})}{Q'^{2}+m^{2}_{s}}
\label{eq:2.11}
\end{eqnarray}
and a cutoff function
\begin{equation}
F(Q^{2}) = \left(\frac{\Lambda^{2}-m_s^{2}}{\Lambda^{2}+Q^{2}}\right)^{2}.
\label{eq:2.12}
\end{equation}
The momentum transfer ${\mathbf{Q}} \: ({\mathbf{Q}}')$ is given by 
\begin{eqnarray}
{\mathbf{Q}}&=&{\mathbf{k}}_{3}-{\mathbf{k}}'_{3} = 
{\mathbf{p}} -{\mathbf{p}}' -
 \frac{1}{2}({\mathbf{q}}-{\mathbf{q}}')  \nonumber \\
\mathbf{Q}'&=&{\mathbf{k}}'_{2}-{\mathbf{k}}_{2} =  
{\mathbf{p}} -{\mathbf{p}}' +
\frac{1}{2}({\mathbf{q}}-{\mathbf{q}}'),
\label{eq:2.13}
\end{eqnarray}
as indicated in Fig.~\ref{fig1}.
\par
For the evaluation of  Eq.~(\ref{eq:2.10}) matrix elements of the form
$\langle{\mathbf{pq}}|V_{4}(1+P)|\psi\rangle$ need to be calculated.
From Fig.~\ref{fig1} we see that $V_4$ can be considered as a sequence of meson
exchanges in the subsystem (12), called for convenience subsystem 3, and
subsystem (31), called 2. Inserting a complete set of states of the type 3
between $V_4$ and $(1+P)|\psi\rangle$ and another complete set of states of type
2 between the two meson exchanges (see Eq.~(\ref{eq:2.11})), leads after
a straightforward evaluation  to 
\begin{eqnarray}
\langle {\mathbf{pq}}|V_{4}(1+P)|\psi \rangle &=&\int
d^{3} q'
\frac{F((-{\mathbf{p}}-\frac{1}{2}{\mathbf{q}}-{\mathbf{q}}')^{2})}
{(-{\mathbf{p}}-
\frac{1}{2}{\mathbf{q}}-{\mathbf{q}}')^{2}+m^{2}_{s}}
 \nonumber \\
&\times& \int d^{3} p'
\frac{F((-{\mathbf{p}}+\frac{1}{2}{\mathbf{q}}-
\frac{1}{2}{\mathbf{q}}'-{\mathbf{p}}')^{2})}
{(-{\mathbf{p}}+\frac{1}{2}{\mathbf{q}}-\frac{1}{2}{\mathbf{q}}'-{\mathbf{p}}')^{2}+m^{2}_{s}}
\langle {\mathbf{p}}'{\mathbf{q}}'|(1+P)|\psi\rangle  
\label{eq:2.14}
\end{eqnarray}
The propagators in Eq.~(\ref{eq:2.14}) contain linear combinations of three
or four momentum
vectors and thus the two integrations would involve magnitudes of vectors and angles
between them. 
Realizing that both meson-exchange propagators in the 3BF term
only depend on the momentum transfer in a two-body subsystem, one can rewrite
 Eq.~(\ref{eq:2.14}) as
\begin{eqnarray}
\langle{\mathbf{p}}{\mathbf{q}}|V_{4}(1+P)|\psi\rangle &=& \int
d^{3} p'd^{3} q'\
\langle{\mathbf{p}}{\mathbf{q}}|{\mathbf{p}}'{\mathbf{q}}'\rangle_{2}
\nonumber \\
&\times&\int d^{3} p''
\frac{F(({\mathbf{p}}'-{\mathbf{p}}'')^{2})}{({\mathbf{p}}'-{\mathbf{p}}'')^{2}+m^{2}_{s}}
\nonumber \\
&\times&\int d^{3} p'''d^{3} q'''\ _{2}
\langle{\mathbf{p}}''{\mathbf{q}}'|{\mathbf{p}}'''{\mathbf{q}}'''\rangle_{3}
\nonumber \\
&\times&\int
d^{3} p''''\frac{F(({\mathbf{p}}'''-{\mathbf{p}}'''')^{2})}{({\mathbf{p}}'''-{\mathbf{p}}'''')^{2}+m^{2}_{s}}\
_{3} \langle{\mathbf{p}}''''{\mathbf{q}}'''|\Psi\rangle.
\label{eq:2.15}
\end{eqnarray}
Here the subscripts $1,2,3$ of the bra and ket vectors stand for the different types
of coordinate systems described by Eq.~(\ref{eq:2.7}). Though the integrations
of  Eq.~(\ref{eq:2.15}) look more complicated, the
meson-propagators show a very simple form. In fact, the integrations over ${\mathbf{p}}''$
and ${\mathbf{p}}''''$ have an identical functional form. 
Defining
\begin{eqnarray}
F_{3}({\mathbf{p}}''',{\mathbf{q}}''') = \int
d^{3} p''''\frac{F(({\mathbf{p}}'''-{\mathbf{p}}'''')^{2})}{({\mathbf{p}}'''-{\mathbf{p}}'''')^{2}+m^{2}_{s}}\
_{3} \langle{\mathbf{p}}''''{\mathbf{q}}'''|\Psi\rangle, 
\label{eq:2.16}
\end{eqnarray}
the integration of the meson exchange between particles 2 and 1
in Eq.~(\ref{eq:2.16}) is carried out completely
in the coordinate system
of type $3$. Once $F_{3}({\mathbf{p}}''',{\mathbf{q}}''')$ is obtained, it needs to
be expressed in terms of momenta in 
a coordinate system of type 2 in order to carry out the integration
over the remaining meson exchange.  This transformation, labeled
$F_{32}({\mathbf{p}}'',{\mathbf{q}}')$  is explicitly given as 
\begin{eqnarray}
F_{32}({\mathbf{p}}'',{\mathbf{q}}') &=& \int
d^{3} p'''d^{3} q''' \ _{2}\langle
{\mathbf{p}}''{\mathbf{q}}'|{\mathbf{p}}'''{\mathbf{q}}'''\rangle_{3}
\:
F_{3}({\mathbf{p}}''',{\mathbf{q}}''') \nonumber \\
&=&F_{3}(|-\frac{1}{2}{\mathbf{p}}''-\frac{3}{4}{\mathbf{q}}'|,
|{\mathbf{p}}''-\frac{1}{2}{\mathbf{q}}'|,
\frac{(-\frac{1}{2}{\mathbf{p}}''-\frac{3}{4}{\mathbf{q}}')\cdot({\mathbf{p}}''-\frac{1}{2}{\mathbf{q}}')}
{|-\frac{1}{2}{\mathbf{p}}''-\frac{3}{4}{\mathbf{q}}'||{\mathbf{p}}''-\frac{1}{2}{\mathbf{q}}'|}).
\nonumber \\
\label{eq:2.17}
\end{eqnarray}
Here we used that $F_3({\bf p'''},{\bf q'''})$ is a scalar function due to the
total wave function $\Psi({\bf p},{\bf q})$ being a scalar in the ground state.
The integration over the second  meson exchange between particle 1 and 3
in the coordinate system of type $2$ is now given by
\begin{eqnarray}
F_{2}({\mathbf{p}}',{\mathbf{q}}')=\int
d^{3} p''\frac{F(({\mathbf{p}}'-{\mathbf{p}}'')^{2})}{({\mathbf{p}}'-{\mathbf{p}}'')^{2}+m^{2}_{s}}
F_{32}({\mathbf{p}}'',{\mathbf{q}}'). 
\label{eq:2.18}
\end{eqnarray}
The matrix element  
$\langle {\mathbf{p}}{\mathbf{q}}|V_{4} (1+P)|\psi\rangle$ is
finally obtained by integrating 
$F_{2}({\mathbf{p}}',{\mathbf{q}}')$ 
over ${\mathbf{p}}'$ and ${\mathbf{q}}'$, i.e. carrying out the final coordinate transformation
from the system of type $2$ back to the one of type $1$, 
\begin{eqnarray}
\langle{\mathbf{p}}{\mathbf{q}}|V_{4}(1+P)|\psi\rangle &=& \int
d^{3} p' d^{3} q' \langle
{\mathbf{p}}{\mathbf{q}}|{\mathbf{p}}'{\mathbf{q}}'\rangle_{2}
F_{2}({\mathbf{p}}',{\mathbf{q}}') \nonumber \\
&=&F_{2}(|-\frac{1}{2}{\mathbf{p}}+\frac{3}{4}{\mathbf{q}}|,
|-{\mathbf{p}}-\frac{1}{2}{\mathbf{q}}|,
\frac{(-\frac{1}{2}{\mathbf{p}}+\frac{3}{4}{\mathbf{q}})\cdot(-{\mathbf{p}}-\frac{1}{2}{\mathbf{q}})}
{|-\frac{1}{2}\mathbf{p}+\frac{3}{4}\mathbf{q}||-\mathbf{p}-\frac{1}{2}\mathbf{q}|}).
\nonumber \\
\label{eq:2.19}
\end{eqnarray}
Thus, the integration of Eq.~(\ref{eq:2.15}), explicitly given in separate steps
from Eq.~(\ref{eq:2.16}) to Eq.~(\ref{eq:2.19}) contains only integrations over one
vector variable at a time. It should be pointed out that 
Eqs.~(\ref{eq:2.17}) and Eq.~(\ref{eq:2.19}) are 
only three dimensional interpolations. 
Clearly, the in Eq.~(\ref{eq:2.15}) suggested
method is the preferred one for practical calculations.

The Faddeev amplitude $\psi({\mathbf{p}},{\mathbf{q}})$ is given as function of
vector Jacobi momenta and obtained as solution of a three dimensional integral
equation, Eq.~(\ref{eq:2.10}). For the ground state $\psi({\mathbf{p}},{\mathbf{q}})$
is also a scalar and thus only depends on the magnitudes of ${\mathbf{p}}$ and
${\mathbf{q}}$ and the angle between the two vectors. 
In order to solve the Eq.~(\ref{eq:2.10})
directly without introducing partial wave projection, we have to define a
coordinate system. We choose the vector ${\mathbf{q}}$ parallel to the $z$-axis
and express the remaining vectors with respect to ${\mathbf{q}}$.
For the first term in Eq.~(\ref{eq:2.10}) the relevant vectors are ${\mathbf{p}}$
and ${\mathbf{q}}'$. Thus one has aside from the magnitudes the 
following angle relations
\begin{eqnarray}
x&=&\hat{\mathbf{p}}\cdot\hat{\mathbf{q}}=\cos\theta \nonumber \\
x'&=&\hat{\mathbf{q}}'\cdot\hat{\mathbf{q}}=\cos\theta' \nonumber \\
y&=&\hat{\mathbf{p}}\cdot\hat{\mathbf{q}}'=\cos\gamma
\label{eq:2.20}
\end{eqnarray}
where
\begin{eqnarray}
\cos\gamma=\cos\theta\cos\theta'+\sin\theta\sin\theta'\cos(\phi-\phi')
=xx'+\sqrt{1-x^{2}}\sqrt{1-x'^{2}}\cos\phi'
\label{eq:2:21}
\end{eqnarray}
Since the $\phi'$ integration is over the full $2\pi$ interval we have a 
freedom of choice for the azimuthal angle angle and set $\phi=0$. With these
choice of variables  Eq.~(\ref{eq:2.10}) with only the first term, i.e. two-body
forces alone, was solved successfully in 
Ref\cite{elster98-1}. 
For the evaluation of the second term in Eq.~(\ref{eq:2.10}),
\begin{eqnarray}
\langle {\mathbf{pq}}|V_{4}(1+P)|\psi\rangle =\langle {\mathbf{pq}}|V_{4}|\Psi\rangle,
\label{eq:2.22} 
\end{eqnarray}
we start with calculating first $F_{3}({\mathbf{p}}''', {\mathbf{q}}''')$, 
Eq.~(\ref{eq:2.16}), and realize, that for this integration we can choose 
${\mathbf{q}}'''$ parallel to the $z$-axis with corresponding simplifications for one of the azimuthal
angles.
This leads to the explicit expression 
\begin{eqnarray}
\lefteqn {F_{3}(p''',q''',x''') = }   \nonumber \\ 
& & \int^{\infty}_{0}
dp''''p''''^{2}\int^{+1}_{-1}dx''''\int^{2\pi}_{0}d\phi''''
\frac{\left( \frac{\Lambda^{2}-m^{2}_{\alpha}}{\Lambda^{2}+(p'''^{2}+p''''^{2}-2p'''p''''y'''')^{2}} \right)^{2}} 
{(p'''^{2}+p''''^{2}-2p'''p''''y'''')^{2}+m^{2}_{\alpha}}
\Psi(p'''',q''',x''''). \nonumber \\
\label{eq:2.24}
\end{eqnarray}
The evaluation of $F_{32}({\mathbf{p}}'',{\mathbf{q}}')$, Eq.~(\ref{eq:2.17}) is
not an integration but rather a three dimensional interpolation and explicitly given by
\begin{eqnarray}
\lefteqn {F_{32}(p'',q',x'') =} \nonumber \\
& & F_{3}\left(\frac{1}{2}\sqrt{\frac{9}{4}q'^{2}+p''^{2}+3p''q'x''},
                   \sqrt{\frac{1}{4}q'^{2}+p''^{2}- p''q'x''},
        \frac{\frac{3}{8}q'^{2}-\frac{1}{2}p''^{2}-\frac{1}{2}p''q'x''}
        {\left|-\frac{3}{4}{\mathbf{q}}'-\frac{1}{2}{\mathbf{p}}''\right|
         \left|+{\mathbf{p}}''-\frac{1}{2}{\mathbf{q}}'\right|}
       \right). \nonumber \\ 
\label{eq:2.25}
\end{eqnarray}
with
\begin{eqnarray}
\left|-\frac{3}{4}{\mathbf{q}}'-\frac{1}{2}{\mathbf{p}}''\right|&=&
\frac{1}{2}\sqrt{\frac{9}{4}q'^{2}+p''^{2}+3p''q'x''}
\nonumber \\
\left|+{\mathbf{p}}''-\frac{1}{2}{\mathbf{q}}'\right|&=&\sqrt{\frac{1}{4}q'^{2}+p''^{2}-
p''q'x''}. 
\label{eq:2.26}
\end{eqnarray}
For the interpolation we apply the cubic splines introduced in Ref\cite{hueber97-1}.
The integration over the second meson exchange, i.e. the calculation of 
$F_{2}({\mathbf{p}}',{\mathbf{q}}')$ of Eq.~(\ref{eq:2.18}) is functionally the same as 
Eq.~(\ref{eq:2.24}), since we can choose the variable ${\mathbf{q}}'$ parallel to the 
$z$-axis. Thus we have the same expression as Eq.~(\ref{eq:2.24}) with 
$p'$, $q'$, $x'$, $p''$, $x''$ and $\phi''$
instead of $p'''$, $q'''$, $x'''$, $p''''$, $x''''$ and $\phi''''$. Finally, the
matrix element $\langle {\mathbf{pq}}|V_{4}|\Psi\rangle$ is explicitly 
obtained by a second interpolation as
\begin{eqnarray}
\langle{\mathbf{p}}{\mathbf{q}}|V_{4}|\Psi\rangle& \equiv &
V_{4}\Psi(p,q,x) \nonumber \\
&=&F_{2}\left(\frac{1}{2}\sqrt{\frac{9}{4}q^{2}+p^{2}-3pqx},
                   \sqrt{\frac{1}{4}q^{2}+p^{2}+ pqx},
        \frac{-\frac{3}{8}q^{2}+\frac{1}{2}p^{2}-\frac{1}{2}pqx}
        {\left|+\frac{3}{4}{\mathbf{q}}-\frac{1}{2}{\mathbf{p}}\right|
         \left|-{\mathbf{p}}-\frac{1}{2}{\mathbf{q}}\right|}
       \right)\nonumber \\ \label{eq:2.27}
\end{eqnarray}
The last term of Eq.~(\ref{eq:2.10}) requires an additional integration
of the matrix element $\langle {\mathbf{p}}{\mathbf{q}}|V_{4}(1+P)|\psi\rangle$  
and the half shell two body $t$-matrix. Again,
with choosing ${\mathbf{q}}$ parallel to the $z$-axis we only have 
three vectors to consider, $\tilde{\mathbf{p}}$, ${\mathbf{p}}$ 
and ${\mathbf{q}}$, thus the integration is of a similar type as
the one of the first term in Eq.~(\ref{eq:2.10}),
\begin{eqnarray}
&&\frac{1}{2}\int d^{3}\tilde{p}
\frac{t_{s}\left({\mathbf{p}},\tilde{{\mathbf{p}}};E-\frac{3}{4m}q^{2} \right)}
{E-\frac{\tilde{p}^{2}}{m}-\frac{3}{4m}q^{2}}
\langle {\tilde{\mathbf{p}}}{\mathbf{q}}|V_{4}(1+P)|\psi\rangle \nonumber \\
&=&\frac{1}{2}
\int^{\infty}_{0} d\tilde{p}\tilde{p}^{2} \int^{+1}_{1}d\tilde{x}
\int^{2\pi}_{0}d\tilde{\phi}
\ \frac{t_{s}\left(p,\tilde{p},\tilde{y_{p}},;E-\frac{3}{4m}q^{2} \right)}
{E-\frac{\tilde{p}^{2}}{m}-\frac{3}{4m}q^{2}}
V_{4}\Psi(\tilde{p},q,\tilde{x}) 
\label{eq:2.29}
\end{eqnarray}
with
\begin{eqnarray}
\tilde{x}&=&\hat{\tilde{\mathbf{p}}}\cdot\hat{\mathbf{q}}\nonumber \\
\tilde{y}_{p}&=&\hat{\tilde{\mathbf{p}}}\cdot\hat{\mathbf{p}}=
x\tilde{x}+\sqrt{1-x^{2}}\sqrt{1-\tilde{x}^{2}}\cos\tilde{\phi}.
\label{eq:2.30}
\end{eqnarray}
We obtain the energy eigenvalue $E$ and Faddeev component $\psi(p,q,x)$ of 
the three-body system by solving 
Eq.~(\ref{eq:2.10}). 

Finally, we want to give the explicit expression for the full wave
function from Eq.~(\ref{eq:2.6}), which is 
\begin{eqnarray}
\Psi(p,q,x)&=&\psi(p,q,x) \nonumber \\
&+&
\psi\left( 
\frac{1}{2}\sqrt{ \frac{9}{4}q^{2}+p^{2}+3pqx},
           \sqrt{ \frac{1}{4}q^{2}+p^{2}-pqx },  
\frac{\frac{3}{8}q^{2}-\frac{1}{2}p^{2}-\frac{1}{2}pqx}
     {\left|-\frac{3}{4}{\mathbf{q}}-\frac{1}{2}{\mathbf{p}} \right|
        \left|-\frac{1}{2}{\mathbf{q}}+{\mathbf{p}} \right| } 
\right) \nonumber \\
&+&
\psi\left( 
\frac{1}{2}\sqrt{ \frac{9}{4}q^{2}+p^{2}-3pqx},
           \sqrt{ \frac{1}{4}q^{2}+p^{2}+pqx},  
\frac{-\frac{3}{8}q^{2}+\frac{1}{2}p^{2}-\frac{1}{2}pqx}
     {\left|+\frac{3}{4}{\mathbf{q}}-\frac{1}{2}{\mathbf{p}} \right|
        \left|-\frac{1}{2}{\mathbf{q}}-{\mathbf{p}} \right| } 
\right). \nonumber \\
\label{eq:2.31}
\end{eqnarray}
The wave function is normalized according to
\begin{eqnarray}
\int d^{3}p\:d^{3}q\:\Psi^{2}({\mathbf{p}},{\mathbf{q}})=1.
\label{eq:2.33}
\end{eqnarray}
  

%----------------------------------------------------------------------------------------------

\section{Calculation of the three-body Bound State}
Our model calculations are based on
Yukawa interactions. As two-body force (2BF) we 
employ a Malfliet-Tjon \cite{malfliet69-1} type  potential
\begin{eqnarray}
V({\mathbf{q}},{\mathbf{q}}')=
&-&\frac{g^2_A}{(2\pi)^3}\frac{1}{({\mathbf{q}}-{\mathbf{q}}')^{2}+m^{2}_{A}}
\left(\frac{\Lambda^{2}_{A}-m^{2}_{A}}{({\mathbf{q}}-{\mathbf{q}}')^{2}+
\Lambda^{2}_{A}}\right)^{2}
\nonumber \\
&+&\frac{g^2_R}{(2\pi)^3}\frac{1}{({\mathbf{q}}-{\mathbf{q}}')^2+m^2_R}
\left(\frac{\Lambda^2_R-m^2_R}{({\mathbf{q}}-{\mathbf{q}}')^2+\Lambda^2_R}\right)^{2},
\label{eq:3.1}
\end{eqnarray}
which is modified by a cutoff function of dipole type. With this choice, the 2BF
and the 3BF have similar functional forms for the scalar meson exchanges. The force
in Eq.(\ref{eq:3.1}) is a superposition of a short-ranged repulsive and 
long-ranged attractive 
Yukawa interactions. The coupling constants and exchanged meson masses are characterized 
by subscripts R and A respectively. The exchanged masses $m_{A}$ and $m_{R}$ are those 
from the original Malfliet-Tjon model, the coupling constants are chosen so that the
two-body force gives a binding of the three-body system which is slightly smaller than 
the experimental value of the triton binding, which is $8.48\:{\mathrm{MeV}}$. 
The cutoff masses
$\Lambda$ have values typical for one-boson-exchange models. 
The parameters for the 2BF,
named MT2-II, are given in Table~\ref{table1}. They lead to a two-body binding energy of 0.284 MeV, and
the S-wave phase-shift roughly follows the shape of the 
experimental one for the state $^3S_1$, though being less attractive. Since we
neglect spin degrees of freedom, any closer adjustment to real phase shifts would be meaningless
anyhow.

With this interaction we first solve the Lippmann-Schwinger equation for 
the fully-off-shell
two-body $t$-matrix directly as function of vector variables as described in 
Ref.~\cite{elster98-2}.
This $t$-matrix is then symmetrized to obtain $t_{s}(p',p,x,E-\frac{3}{4m}q^{2})$.   
The eigenvalue equation, Eq.~(\ref{eq:2.10}), for the three-body bound state is solved 
iteratively by a Lanczo's type algorithm described in detail in Ref.\cite{stadler}.
Using the 2BF alone, the binding energy of the three-body system is calculated 
as $E=7.6986\:{\mathrm{MeV}}$.
\par 
The simplest 3BF we want to apply in our study has the functional form
\begin{eqnarray}
V_{4}=\frac{1}{(2\pi)^{6}} \frac{a_{\alpha}}{m_{\alpha}} g^{2}_{\alpha}
\frac{F_{\alpha}(Q^{2})}{Q^{2}+m^{2}_{\alpha}}
\frac{F_{\alpha}(Q'^{2})}{Q'^{2}+m^{2}_{\alpha}}, 
\label{eq:3.2}
\end{eqnarray}
where
\begin{eqnarray}
F_{\alpha}(Q^{2})&=&
\left(\frac{\Lambda^{2}_{\alpha}-m^{2}_{\alpha}}{\Lambda^{2}_{\alpha}+Q^{2}}\right)^{2}.
\label{eq:3.3}
\end{eqnarray}
Choosing $a_{\alpha}$ to be a negative constant makes this force purely attractive. 
The parameters of this force, which we name MT3-I in the following, 
are given in Table~\ref{table2}. They are chosen to give a small attractive contribution 
to the three-body binding energy, such that we end up in the neighborhood of the triton binding
energy. For the force parameters in Table~\ref{table2} 
we obtain $E=8.8732$~MeV. 

In order to solve the eigenvalue equation, Eq.~(\ref{eq:2.10}), for the Faddeev component
$\psi(p,q,x)$, we use Gaussian grid points in $p$, $q$, and $x$. 
The momentum and angle grids
in the integrations of the 3BF in Eqs.~(\ref{eq:2.24}) and (\ref{eq:2.29})
have the same sizes as the ones for $p$, $q$, and $x$.
This is very reasonable, since the integrations over the meson exchange contributions of
the 3BF, Eq.~(\ref{eq:2.24}), require  grids similar in range to the one used to
calculate the two-body t-matrix. The p-grid is defined between 0 and $p_{max}=60 fm^{-1}$,
whereas for the q-grid a maximum value $q_{max} = 40 fm^{-1}$ is sufficient. For the
angle ($x$) integration, the preferred number is 42 grid points. Further
details concerning the grid choices are given in Ref.~\cite{elster98-1}. 
Typical grid sizes are $97 \times 97 \times 42$ to obtain an accuracy  
in the binding energy of 5 significant figures. The convergence of the three-body binding
energy $E$ as function of the number of grid points is shown in Table~III, 
where we see a convergence of the energy eigenvalue and the expectation values within
5 digits. 

The three-body wave function is calculated from the Faddeev component
using Eq.~(\ref{eq:2.31}). Since this wave function enters the eigenvalue equation,
Eq.~(\ref{eq:2.10}), we also need to worry about the quality of the calculation
of $\Psi(p,q,x)$. One check of the overall quality of the wave function
 is a comparison of the expectation
value of the total Hamiltonian $\langle H \rangle$ with the calculated value $E$ from
the solution of the Faddeev equation. Explicitly, we evaluate
\begin{eqnarray}
\langle H \rangle
 \equiv \langle\Psi|H|\Psi\rangle=\langle\Psi|H_{0}|\Psi\rangle+\langle\Psi|V_{II}|\Psi\rangle+
\langle\Psi|V_{123}|\Psi\rangle,
\label{eq:3.4}
\end{eqnarray}
where $V_{II}$ represents the 2BF $\sum^{3}_{i=1}V^{i}$ and $V_{123}$
the three-body defined in Eq.~(\ref{eq:2.2}).    
The expectation value of the kinetic energy $\langle H_0 \rangle$ and the two-body potential
energy $\langle V_{II} \rangle$ are given as \cite{elster98-1}
\begin{eqnarray}
\langle H_0 \rangle& \equiv &
\langle\Psi|H_{0}|\Psi\rangle=3\langle\psi|H_{0}|\Psi\rangle
\nonumber \\
&=& 3\cdot
8\pi^{2}\int^{\infty}_{0}p^{2}dp\int^{\infty}_{0}q^{2}dq\left(
\frac{p^{2}}{m}+\frac{3q^{2}}{4m} \right)\int^{+1}_{-1}dx
\psi(p,q,x)\Psi(p,q,x) \nonumber \\
\label{eq:3.5}
\end{eqnarray}
and
\begin{eqnarray}
\langle V_{II} \rangle & \equiv &
\langle\Psi|V_{II}|\Psi\rangle=3\langle\Psi|V^{1}|\Psi\rangle
\nonumber \\
&=&3\cdot
8\pi^{2}\int^{\infty}_{0}p^{2}dp\int^{\infty}_{0}q^{2}dq\int^{+1}_{-1}dx\int^{\infty}_{0}p'^{2}dp'
\int^{+1}_{-1}dx'
\nonumber \\
&\times&\Psi(p,q,x)v_{1}(p,p',x,x')\Psi(p',q,x')
\label{eq:3.6}
\end{eqnarray}
where
\begin{eqnarray}
v_{1}(p,p',x,x')=\int^{2\pi}_{0}d\phi
V^{1}(p,p',xx'+\sqrt{1-x^{2}}\sqrt{1-x'^{2}}\cos\phi).
\label{eq:3.7}
\end{eqnarray}
The expectation value of the three-body potential energy, $\langle V_{123} \rangle$, 
is given by 
\begin{eqnarray}
\langle V_{123} \rangle& \equiv &
\langle\Psi|V_{123}|\Psi\rangle=3\langle\Psi|V_{4}|\Psi\rangle
\nonumber \\
&=&3\cdot
8\pi^{2}\int^{\infty}_{0}p^{2}dp\int^{\infty}_{0}q^{2}dq\int^{+1}_{-1}dx
\Psi(p,q,x)V_{4}\Psi(p,q,x).
\label{eq:3.8}
\end{eqnarray}
Here the integrations need the evaluation of the matrix element 
$\langle{\mathbf{p}}{\mathbf{q}}|V_{4}|\Psi\rangle$ of Eq.~(\ref{eq:2.14}). The expectation
values of the kinetic and potential energies are listed in Table~III as functions of
the size of the $p-q-x$ grid.  Table~III
shows that the expectation values given above 
converge for 5 significant figures when the grid sizes for $p$ and $q$ exceed $77$ points.
However, Table~III
also shows that the expectation value $\langle H \rangle$ does not exactly 
converge to the calculated value $E$ despite increased grid size. The difference between
the two quantities is of the order of 10 keV. This behavior is in 
contrast to calculations based solely on two-body forces. 
Calculations with $V_{123}=0$
are shown in Table~\ref{table4} as function of the size of the  $p-q-x$ grid. 
Here the convergence
of $\langle H \rangle$ to the calculated value $E$ as function of the grid size is 
much better than in the case where an attractive 3BF is included. 

It is well known that three nucleon forces based on multi-meson exchanges can have
attractive as well as repulsive pieces. Thus we also consider a model of this type
given as
\begin{eqnarray}
V_{4}&=&\frac{1}{(2\pi)^{6}}\frac{a_{\alpha}}{m_{\alpha}} g^{2}_{\alpha}
\frac{F_{\alpha}(Q^{2})}{Q^{2}+m^{2}_{\alpha}}
\frac{F_{\alpha}(Q'^{2})}{Q'^{2}+m^{2}_{\alpha}}\nonumber \\
&+&\frac{1}{(2\pi)^{6}}\frac{a_{\alpha\rho}}{\sqrt{m_{\alpha}m_{\rho}}}
g_{\alpha}g_{\rho} \left(
\frac{F_{\alpha}(Q^{2})}{Q^{2}+m^{2}_{\alpha}}
 \frac{F_{\rho}(Q'^{2})}{Q'^{2}+m^{2}_{\rho}}+
 \frac{F_{\rho}(Q'^{2})}{Q'^{2}+m^{2}_{\rho}}
 \frac{F_{\alpha}(Q^{2})}{Q^{2}+m^{2}_{\alpha}} \right).
\label{eq:4.2}
\end{eqnarray}
Here the first term represents an attractive force, characterized by a negative
coupling
$a_{\alpha}$, whereas the second term represent a repulsive force, i.e.
$a_{\alpha\rho}$
is positive. Since the masses of the exchange mesons are different, the form of the
second, repulsive term guarantees that $V_4$ is symmetric under a permutation of
nucleons
2 and 3. The cutoff functions $F_{\alpha}$ and $F_{\rho}$ have the same functional
form as
given in Eq.~(\ref{eq:3.3}). The parameters of this 3BF, named MT3-II in
the
following, are given in Table~\ref{table5}. They are chosen so that the correction due
to this
3BF to the three-body binding energy calculated with the 2BF
MT2-II is small. The binding energy $E$ with this MT3-II 3BF gives $E =
8.6478$~MeV. 

The expectation values of the kinetic and potential energies are listed in
Table~\ref{table6a} as functions of the size of the $p-q-x$ grid. Again, the expectation
values converge within five significant figures when the grid sizes for $p$ and $q$
exceed 77 points. However, now the difference of the expectation value for the total
Hamiltonian $\langle H \rangle$ deviates from the calculated eigenvalue $E$ only by
5~keV, a number being similar to calculations carried out in a partial wave decomposition
and based on realistic forces \cite{NoggaPhd}.

All these numbers are not meant to provide insight into the physics of three interacting nucleons, but
serve only as a demonstration that this technique allows a very accurate and easy handling of typical
nuclear forces consisting of attractive and repulsive (short range) parts. In addition, they will
serve as benchmarks for future studies.


%------------------------------------------------------------------------------

\section{Wave Function Properties}

Despite the fact that there is no spin dependence and questions about effects of realistic forces can
not be posed, we want to display some wave function properties, which are often studied in the context
of realistic forces. It will turn out that qualitatively they also appear in our simple three boson
model. The probability of finding a nucleon with momentum $q$ in the nucleus is given as
\begin{equation}
n(q)= 2\pi q^2 \int^{\infty}_{0} dp \: p^2 \int^{+1}_{-1}dx \Psi^{2}(p,q,x) .
\label{eq:4.1}
\end{equation}
The total wave function $\Psi(p,q,x)$ is given by Eq.~({\ref{eq:2.31}). 
In Fig. 2 we show $n(q)$ with and without three-body forces. There is hardly any change, a fact with
has been noticed before in the context of realistic nuclear forces \cite{NoggaPhd}. 
In addition, the shoulder of
the distribution around 2-4~fm$^{-1}$ is qualitatively similar to the case when using realistic
forces. The momentum distribution $n(q)$ is shown in Fig.~\ref{fig2} for three different cases,
one based on a calculation with two-body forces alone (dotted line), and the other two for which the
two different 3BFs are included. First we notice that the two 3BF's, though different in character, lead
to essentially the same momentum distribution. Compared to the momentum distribution given by
the two-body force alone, the minimum is shifted to a slightly higher momentum.

Another property often investigated is the probability to find two nucleons at a distance $r$. 
To obtain that quantity we generate the total wave function in configuration space as 
\begin{equation}
\Psi({\mathbf{r,R}})=\int d^{3}p\ d^{3}q\
\Psi({\mathbf{p,q}}) \exp(i{\mathbf{p}}\cdot{\mathbf{r}})
\exp(i{\mathbf{q}} \cdot {\mathbf{R}}).
\label{eq:4.3}
\end{equation}
Here $\Psi({\bf p}$, ${\bf q})$ is the total wave function in
momentum space as  given by Eq.~(\ref{eq:2.6}). The variables ${\bf r}$ and
${\bf R}$ are conjugate to the Jacobi momenta ${\bf p}$, ${\bf q}$ and
given as
\begin{eqnarray}
{\mathbf{r}}&=&{\mathbf{x}}_{2}-{\mathbf{x}}_{3}, \nonumber \\
{\mathbf{R}}&=&{\mathbf{x}}_{1}-\frac{1}{2}({\mathbf{x}}_{2}+{\mathbf{x}}_{3}).
\label{eq:4.4}
\end{eqnarray}
where ${\bf x}_1$, ${\bf x}_2$ and ${\bf x}_3$ are the
coordinates of three nucleons in configuration space \cite{gloeckle-book}.

For the explicit calculation of the double Fourier transformation we first consider
the ${\bf q}$-integration
\begin{eqnarray}
\int d^{3}q\exp(i\mathbf{q}\cdot\mathbf{R})\Psi(\mathbf{p,q}).
\end{eqnarray}
We choose the vector ${\bf p}$ parallel to the z-axis and define the angles
$\hat{\mathbf{q}}\cdot\hat{\mathbf{z}}=x_{q}$ and $\hat{\mathbf{R}}\cdot\hat{\mathbf{z}}=x_{R}$.
Since the integration is carried out over all space, we can set $\phi_{R}=0$, and
obtain
\begin{eqnarray}
&&\int d^{3}q\exp(i\mathbf{q}\cdot\mathbf{R})\Psi(\mathbf{p,q})
\nonumber \\
&&=\int^{\infty}_{0}q^{2}dq\int^{+1}_{-1}dx_{q}\int^{2\pi}_{0}d\phi_{q}\exp(iqR\hat{\mathbf{R}}\cdot\hat{\mathbf{q}})\Psi(p,q,x_{q}),
\end{eqnarray}
where
\begin{eqnarray}
\hat{\mathbf{R}}\cdot\hat{\mathbf{q}}=x_{q}x_{R}+\sqrt{1-x^{2}_{q}}\sqrt{1-x^{2}_{R}}\cos\phi_{q}.
\end{eqnarray}
Thus, the integration over $\phi_{q}$ can be carried out separately
\begin{eqnarray}
&&\int^{2\pi}_{0}d\phi_{q}\
\exp(iqR\sqrt{1-x^{2}_{q}}\sqrt{1-x^{2}_{R}}\cos\phi_{q})\nonumber
\\
&&=2\pi J_{0}(qR\sqrt{1-x^{2}_{q}}\sqrt{1-x^{2}_{R}}\cos\phi_{q}).
\end{eqnarray}
Summarizing the above leads to the intermediate result
\begin{eqnarray}
&&\int d^{3}q\exp(i\mathbf{q}\cdot\mathbf{R})\Psi(\mathbf{p,q})
\nonumber \\
&&=2\pi\int^{\infty}_{0}q^{2}dq\int^{+1}_{-1}dx_{q}
J_{0}(qR\sqrt{1-x^{2}_{q}}\sqrt{1-x^{2}_{R}}\cos\phi_{q})\exp(iqRx_{q}x_{R})\Psi(p,q,x_{q})
\nonumber \\
&&\equiv 2\pi \: \Psi_{p}(p,R,x_{R}).
\end{eqnarray}
Next, we consider the integration over ${\bf p}$, where it is convenient to choose 
the vector ${\bf r}$
parallel to the z-axis. Thus, the following angle, $\hat{\mathbf{R}}\cdot\hat{\mathbf{z}} 
\equiv \hat{\mathbf{R}}\cdot\hat{\mathbf{r}} \equiv x_{R}$,
needs to be considered, and the integration over
$x_{p}$ and $\phi_{p}$ can be carried out separately as
\begin{eqnarray}
\int^{+1}_{-1}dx_{p}\int^{2\pi}_{0}d\phi_{p}\exp(iprx_{p})=4\pi\frac{\sin(pr)}{pr}
\end{eqnarray}
Finally, the Fourier transform of $\Psi(p,q,\hat{\mathbf{p}}\cdot\hat{\mathbf{q}})$
can be calculated as
\begin{equation}
\Psi(r,R,x_R)=
\frac{8\pi^{2}}{r}\int^{\infty}_{0}dp\left[\sin(pr)p\Psi_{p}(p,R,x_{R})\right],
\label{eq:4.5}
\end{equation}
where
\begin{equation}
\Psi_{p}(p,R,x_{R})=\int^{\infty}_{0}q^{2}dq\int^{+1}_{-1}dx_{q}
J_{0}\left(qR\sqrt{1-x^{2}_{q}}\sqrt{1-x^{2}_{R}}\right)\cos(qRx_{q}x_{R})\Psi(p,q,x_{q}).
\label{eq:4.6}
\end{equation}
With this the two-body correlation function $c(r)$ is defined as
\begin{equation}
c(r)=2\pi
r^{2}\int^{\infty}_{0}dR \: R^2 \int^{+1}_{-1}dx_{R}\Psi^{2}(r,R,x_{R}).
\label{eq:4.7}
\end{equation}
The correlation function $c(r)$ describes the probability to find two nucleons within a
relative distance $r$. In Fig. 3 the correlation functions are displayed based on a 
calculation with two-body forces (dotted line) and based on calculations with the two different 3BF's. 
Though our model is very simple, the functions $c(r)$ are similar to the ones obtained with 
realistic forces \cite{NoggaPhd}. We see that the maximum of $c(r)$ is shifted slightly to a smaller
value of $r$ once a 3BF is included, which is consistent with the minimum of $n(q)$ being shifted to a
slightly higher momentum. The position of the maximum of $c(r)$ does not depend on the type of 3BF,
however the actual shape of the function does. 

Since our three-body system consists of three identical nucleons acted on by scalar
forces, the three nucleons form a ground state where the most probable positions of the nucleons
have the 
shape of an equilateral triangle. The expectation values of the Jacobi coordinates $r$ and
$R$ can be calculated as
\begin{eqnarray}
\langle r \rangle =
\int^{\infty}_{0}R^{2}dR\int^{\infty}_{0}r^{2}dr\int^{+1}_{-1}dx_{R}
\ r\Psi^{2}(r,R,x_{R}), \nonumber \\
\langle R \rangle =
\int^{\infty}_{0}R^{2}dR\int^{\infty}_{0}r^{2}dr\int^{+1}_{-1}dx_{R} \
R\Psi^{2}(r,R,x_{R}).
\label{eq:4.8}
\end{eqnarray}
Here the values  $\langle r \rangle$ and $\langle R \rangle$ are the 
length and height of the equilateral triangle. The geometrical relation between the length
and height of an equilateral triangle is given by
\begin{equation}
\langle r \rangle:\langle R \rangle=2 : \sqrt{3}.
\label{eq:4.9}
\end{equation}
We also define a deviation $\delta$ by
\begin{equation}
\delta=\frac{\langle r \rangle/\langle R \rangle-2/\sqrt{3}}{2/\sqrt{3}}\times 100.
\label{eq:4.10}
\end{equation}
to account for deviations from an ideal geometric triangle. 
In Table~\ref{tabledel} we show the corresponding values for $\langle r \rangle$, $\langle R
\rangle$, and $\delta$ for the MT2-II 2BF alone and the cases where the two different 3BF
discussed above are added. In all cases listed, the deviation from an ideal equilateral
triangle is 3\% or less, also indicating that our calculations are very accurate. 


We would like to add a little excursion into a playground with forces. Assume there are only
three-body forces. Can one generate a three-body wave function which has about the same binding
energy, single nucleon momentum distribution, and two-body correlation function as given by two-body
forces alone? This is indeed possible. For that aim we have chosen purely attractive two- and
three-body forces, the parameters of which are given in Table~\ref{table11}. The binding energies and
expectation values of kinetic and potential energies for the two cases are displayed in
Table~\ref{table7}. It
turns out that also $n(q)$ and $c(r)$ are close to each other as shown in Figs. 3 and 4. 
We do not know of a physical realization of such a scenario with pure three-body forces, but maybe the
reader may find that little excursion equally entertaining as we do. 


%------------------------------------------------------------------------------

\section{Summary and Outlook}

We derived and calculated  the Faddeev equation for three identical bosons
interacting by two- and three-body forces. The equation is
formulated in momentum space directly in terms of momentum vectors, i.e.
without angular momentum decomposition. It is demonstrated that
this equation can be solved by integrating over magnitudes
of momenta and various angles. In doing so we encounter interpolations
which are carried out by cubic Hermitian splines \cite{hueber97-1}. The Faddeev equation
is solved by iteration using a Lanczo's type method \cite{stadler,saake}. An accuracy
of 5 digits in the energy eigenvalue can easily be achieved. In
comparison to an angular momentum decomposition which is commonly used
\cite{gloeckle-book},
this direct approach has great advantages. It avoids the very
 involved angular momentum algebra  occurring for
 the permutations and especially  for
the  three-body forces \cite{coon,hueber97-1}.

 All two- and three-body forces we employ are of meson exchange type, either purely
 attractive or attractive and repulsive. The mesons responsible
 for the attraction have masses around 300 MeV , the ones for the repulsion
 around 600 MeV. These forces, serving as a reference,
  are chosen with a view towards nuclear physics. Thus the reference two-body
 forces lead to a three-body binding energy somewhat smaller than
 8.48 MeV (the experimental $^3$H binding energy), and the reference
 three-body forces add about 1 MeV additional binding energy.

   For these type of forces we evaluated the single nucleon momentum
 distributions and the two-body correlation functions.
 These quantities turned out to be qualitatively very similar to what is
 achieved with realistic spin dependent forces \cite{NoggaPhd,NoggaPL}. 

   We also evaluated the expectation values $\langle r \rangle$ for a pair distance and 
$\langle R \rangle$ for
 the distance of a third particle to the c.m. of the corresponding
 pair. In case of an equilateral triangle the ratio of these two
 quantities is $2/\sqrt{3}$.  Since we have three identical bosons
 the position of the three particles in the ground state 
should be with highest probability 
at the corners of an equilateral
 triangle. Indeed the corresponding ratio for the expectation values
$\langle r \rangle$ and $\langle R \rangle$ turned out to be $2/\sqrt{3}$ with an
 numerical error of about 1-2\%.


  Including spin (and isospin) degrees of freedom is an additional
 task for the future, which will increase the space of states
  and will lead to
 coupled equations, but only a strictly finite number of equations.
 The form in which they will appear will depend on the way the
  spin degrees of
 freedom will be incorporated. One possibility will be the extension of  the helicity
 formalism chosen for the NN system in  a three-dimensional
  notation \cite{imam1}
 to three nucleons. The other possible  extension of the here presented formulation
is the incorporation of relativity in the instant form of
 dynamics \cite{coester}. The momentum space formulation seems ideal for
 that. First steps have been already undertaken \cite{kamada,kamada2}. Since
 relativity will be of importance at higher energies our treatment
 with momentum vectors and avoiding angular momentum decomposition will
 be awarding in that respect.


%------------------------------------------------------------------------------

\vfill

\acknowledgements This work was performed in part under the
auspices of the U.~S.  Department of Energy under contract
No. DE-FG02-93ER40756 with Ohio University. The computational support of the 
Ohio Supercomputer Center (OSC) for the use of their facilities under Grant No.~PHS206,
the Neumann Institute for Computing (NIC) under project JIKP01 and the National Energy
Research Supercomputer Center (NERSC) is acknowledged.  The authors want to thank
A.~Nogga for helpful and stimulating discussions.



%------------------------------------------------------------------------------
%\pagebreak

\begin{references}

%\bibitem{sofianos} V.A. Roudnev, S.L. Yakovlev, S.A. Sofianos, `Bound State Calculations for
%Three Atoms without Explicit Partial Wave Decomposition', Preprint
%xxx.lanl.gov/physics/0204025.

\bibitem{malfliet69-1} R.A. Malfliet, A.J. Tjon, Nucl. Phys. A, {\bf 127}, 161 (1969).
\bibitem{osborn} T.A. Osborn, `Faddeev�s Equations for local potentials�, SLAC Report No.
17 (1967).

\bibitem{pickle} A.~Picklesimer, R.A.~Rice, and R.~Brandenburg, Phys.
Rev C{\bf 45}, 2045 (1992), {\it ibid}, 547 (1992), Phys. Rev. C{\bf 44},
1359 (1991).
\bibitem{stadler} A.~Stadler, W.~Gl\"ockle, P.U.~Sauer, Phys. Rev. C{\bf 44}, 2319 (1991); 
\bibitem{stadler2} A.~Stadler and P.U.~Sauer, Phys. Rev. C{\bf 46}, 64 (1992).

\bibitem{nbench} A.~Nogga, D.~H\"uber, H.~Kamada, W.~Gl\"ockle,
Phys. Lett. {\bf B 409}, 19 (1997).

\bibitem{lavern} A. Laverne, C. Gignoux, Nucl. Phys. {\bf A203}, 597 (1973).

\bibitem{payne} J.L.~Friar, B.F.~Gibson, and G.L.~Payne, Z. Phys. {\bf A 301},
309 (1981), C.R.~Chen, G.L.~Payne, J.L.~Friar, B.F.~Gibson, Phys. Rev. C{\bf 31}, 2266 (1985).

\bibitem{schell0} N.W.~Schellingerhout, L.P.~Kok, and G.D.~Bosveld,
Phys. Rev. A{\bf 40}, 5568 (1989); N.W.~Schellingerhout, J.J.~Schut,
and L.P.~Kok, Phys. Rev. C{\bf 46}, 1192 (1992).

\bibitem{wu} Y. Wu, S. Ishikawa, T. Sasakawa, Few-Body Systems {\bf 15},
  145 (1993).

\bibitem{gfmc} J.~Carlson, Phys. Rev. C{\bf 36}, 2026 (1987),
  Phys. Rev. C{\bf 38}, 1879 (1988); J.G.~Zabolitzki, K.E.~Schmidt,
  and M.H.~Kalos, Phys. Rev. C{\bf 25}, 1111 (1982); J.~Carlson and
  R.~Schiavilla, Rev. Mod. Phys. {\bf 70}, 743 (1998).

\bibitem{elster98-1} Ch. Elster, W. Schadow, A. Nogga, and W. 
Gl\"ockle, Few-Body System, {\bf 27}, 83 (1998).

\bibitem{scatt1} W. Schadow, Ch.~Elster, and W.~Gl\"ockle, Few-Body Systems {\bf 28}, 15
(2000).

\bibitem{kvin} A.A.~Kvitsinsky and C.-Y. Hu, Few-Body Systems {\bf 12},
7 (1992).

\bibitem{imam1} I. Fachruddin, Ch.~Elster, W. Gl\"ockle,  Phys. Rev. C{\bf 62}, 044002-1
(2000).

\bibitem{gloeckle-book}W. Gl\"ockle,
The Quantum Mechanical Few-Body Problem, Springer-Verlag (1983).

\bibitem{coon} S.A. Coon and W. Gl\"ockle, Phys. Rev. C{\bf 23}, 1790 (1981).

\bibitem{hueber97-1} D. H\"uber, H. Witala, A. Nogga, W. Gl\"ockle and
H. Kamada, Few-Body System, {\bf 22}, 107 (1997).

\bibitem{Fujita} J. Fujita, H. Miyazawa, Prog. Theor. Phys. {\bf 17}, 360 (1957).

\bibitem{elster96-1} Ch. Elster, E.E. Evans, H. Kamada and W.  Gl\"ockle, 
Few-Body System, {\bf 21}, 25 (1996).

\bibitem{elster98-2} Ch. Elster, J.H. Thomas, W. Gl\"ockle, 
Few-Body Systems, {\bf 24}, 55 (1998).

\bibitem{NoggaPhd} A. Nogga, PhD Thesis, Ruhr-University Bochum, 2001. 

\bibitem{NoggaPL} A. Nogga, D. H\"uber, H. Kamada, W. Gl\"ockle, Phys. Lett. {\bf B409}, 19
(1997).

\bibitem{saake} A. Saake, Diploma Thesis, Ruhr-University Bochum, 1992 (unpublished).

\bibitem{coester} W. Gl\"ockle, T.-S.H. Lee, F. Coester, Phys. Rev. C{\bf 33}, 709
(1986).
 
\bibitem{kamada} H. Kamada, Few-Body Systems, Suppl. {\bf 12}, 433 (2000).
\bibitem{kamada2}  H. Kamada, W. Gl\"ockle, J. Golak, Ch. Elster, Phys. Rev. C{\bf 66}, 044010 (2002).

\end{references}

%\pagebreak
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


% parameters of mt2-ii - tab I
 
\begin{table}\caption{\label{table1}
The parameters of the MT2-II two-body force.} 
\begin{tabular}{cccccc}
\mbox{  }
 $g^{2}_{A}/4\pi$ & $m_{A}$[MeV] & $\Lambda_{A}$[MeV]&
 $g^{2}_{R}/4\pi$ & $m_{R}$[MeV] & $\Lambda_{R}$[MeV] \\ \hline
  3.5775 & 330.2104 & 1500.0 & 9.4086 & 612.4801 & 1500.0  \\ 
  \end{tabular}
\end{table}

% parameters of mt3-i  tab.II

\begin{table}\caption{\label{table2}The parameters of the MT3-I attractive 3BF.} 
\begin{tabular}{ccccc}
\mbox{  }
 $g^{2}_{\alpha}/4\pi$ & $m_{\alpha}$[MeV] & $\Lambda_{\alpha}$[MeV] &
  $a_{\alpha}$\\ \hline
5.0 & 305.8593 & 1000.0 & -1.73  \\ 
  \end{tabular}
\end{table}


% convergence table of mt2-ii + mt3-i  -  tab. III

\begin{table}\caption{\label{table3}The calculated eigenvalue $E$ from the
 the solution of the Faddeev
equation and the expectation values of the kinetic energy $\langle H_0
\rangle$, the two-body potential $\langle V_{II} \rangle$,  the
three-body potential energy $\langle V_{123} \rangle$ and the total Hamiltonian $\langle H
\rangle$ as functions of the number of grid points NP, NQ and NX for the $p-q-x$ grid. The
calculations are based on the MT2-II 2BF and the MT3-I 3BF.}
 
\begin{tabular}{cccccccc}
\mbox{  }

   NP & NQ & NX &$\langle H_0 \rangle$ (MeV)&$\langle V_{II} \rangle$ (MeV)
                &$\langle V_{123} \rangle$ (MeV)& $\langle H \rangle$ (MeV)
                & $E$ (MeV)   \\ \hline
   45 & 45 & 42 & 31.8838 & -39.4069 & -1.3392 & -8.8623  & -8.8715 \\ \hline
   61 & 45 & 42 & 31.8846 & -39.4073 & -1.3396 & -8.8623  & -8.8711 \\ \hline
   77 & 45 & 42 & 31.8848 & -39.4074 & -1.3397 & -8.8623  & -8.8709 \\ \hline
   77 & 61 & 42 & 31.8900 & -39.4133 & -1.3406 & -8.8639  & -8.8726 \\ \hline
   77 & 77 & 42 & 31.8915 & -39.4149 & -1.3409 & -8.8643  & -8.8731 \\ \hline
   87 & 87 & 42 & 31.8919 & -39.4152 & -1.3410 & -8.8644  & -8.8732 \\ \hline
   97 & 97 & 42 & 31.8920 & -39.4154 & -1.3410 & -8.8644  & -8.8732 \\ 
  \end{tabular}
  \end{table}

% energies of mt2-ii  -  tab IV

\begin{table}\caption{\label{table4}The calculated eigenvalue $E$ of the Faddeev
equation and the expectation values of the kinetic energy $\langle H_0
\rangle$, the two-body potential $\langle V_{II} \rangle$, and the total Hamiltonian
 $\langle H \rangle$ as functions of the number of grid points NP, NQ and NX for the
$p-q-x$ grid. The
calculations are based on the MT2-II 2BF alone.} 

\begin{tabular}{ccccccc}
\mbox{  }

   NP & NQ & NX &$\langle H_0 \rangle$ (MeV)&$\langle V_{II} \rangle$ (MeV)& 
   $\langle H \rangle$ (MeV)& $E$ (MeV)  \\ \hline
   77 & 77 & 42 & 28.6408 & -36.3390 & -7.6983 & -7.6984 \\ \hline
   87 & 87 & 42 & 28.6408 & -36.3391 & -7.6983 & -7.6984 \\ \hline
   97 & 97 & 42 & 28.6408 & -36.3392 & -7.6983 & -7.6984 \\ 
  \end{tabular}
  \end{table}


% parameters of MT3-II  -  tab. V

\begin{table}\caption{\label{table5} The parameters of the MT3-II 3BF.} 
\begin{tabular}{cccc}
\mbox {}
 $g^{2}_{\alpha}/4\pi$ & $m_{\alpha}$[MeV] & $\Lambda_{\alpha}$[MeV] &
  $a_{\alpha}$  \\ \hline
  5.0 & 305.8593 & 1000.0 & -2.69  \\ \hline
  $g^{2}_{\rho}/4\pi$ & $m_{\rho}$[MeV] & $\Lambda_{\rho}$[MeV] &
  $a_{\alpha\rho}$ \\ \hline
  9.0 & 650.0000 & 1900.0 &  2.40 \\ 
  \end{tabular}
\end{table}

%  numerics 2BF + MT3 3BF -  tab. VI

\begin{table}\caption{\label{table6a} The calculated eigenvalue $E$ from the
 the solution of the Faddeev
equation and the expectation values of the kinetic energy $\langle H_0
\rangle$, the two-body potential $\langle V_{II} \rangle$,  the
three-body potential energy $\langle V_{123} \rangle$ and the total 
Hamiltonian $\langle H
\rangle$ as functions of the number of grid points NP, NQ and NX for the $p-q-x$ grid.
The calculations are based on the MT2-II 2BF and the MT3-II 3BF.}

\begin{tabular}{cccccccc}

\mbox{  }

   NP & NQ & NX &$\langle H_0 \rangle$ (MeV)&$\langle V_{II} \rangle$ (MeV)
                &$\langle V_{123} \rangle$ (MeV)& $\langle H \rangle$ (MeV)
                & $E$ (MeV)   \\ \hline
   45 & 45 & 42 & 31.1745 & -38.7859 & -1.0404 & -8.6518  & -8.6454 \\ \hline
   61 & 45 & 42 & 31.1767 & -38.7878 & -1.0409 & -8.6520  & -8.6456 \\ \hline
   77 & 45 & 42 & 31.1773 & -38.7881 & -1.0413 & -8.6521  & -8.6466 \\ \hline
   77 & 61 & 42 & 31.1825 & -38.7880 & -1.0478 & -8.6533  & -8.6477 \\ \hline
   77 & 77 & 42 & 31.1837 & -38.7886 & -1.0481 & -8.6530  & -8.6480 \\ \hline
   87 & 87 & 42 & 31.1842 & -38.7887 & -1.0481 & -8.6526  & -8.6478 \\ \hline
   97 & 97 & 42 & 31.1847 & -38.7892 & -1.0481 & -8.6526  & -8.6478 \\
  \end{tabular}
  \end{table}

\begin{table}\caption{\label{tabledel} The binding energy $E$, the expectation values $\langle r
\rangle$ and $\langle R \rangle$ calculated with the MT2-II 2BF alone and with the addition of
the two different 3BF's described in the text. The deviation $\delta$ characterizes the
deviation from the shape of an equilateral triangle and is
defined in Eq.~(\ref{eq:4.10}).  }
\begin{tabular}{cccc}
\mbox{}
          & MT2-II  & MT2-II + MT3-I & MT2-II + MT3-II  \\ \hline 
$E$ (MeV) &  -7.6980  & -8.873 & -8.6478 \\
$\langle r \rangle$ (fm)& 2.521  & 2.382 & 2.9401 \\
$\langle R \rangle$ (fm) & 2.221 & 2.096 & 2.5945 \\
$\delta$ (\%)   &      1.7   & 1.6  &   1.9 \\ 
  \end{tabular}
\end{table}

\begin{table}\caption{\label{table11}The parameters of the
purely attractive potentials, the 2BF MT2-I and the 3BF MMT3-I leading to the
same three-body binding energy. }
\begin{tabular}{lcccc}
\mbox{}
MT2-I &   $g^{2}_A/4\pi$ & $m_A$[MeV] & $\Lambda_A $[MeV] &  \\ \hline
   & -0.7210  & 330.2104  & 1500 & \\ \hline \hline 
MMT3-I &   $g^{2}_{\alpha}/4\pi$ & $m_{\alpha}$[MeV] &
$\Lambda_{\alpha}$[MeV] &
  $a_{\alpha}$   \\ \hline
 &  5.0000 & 60.0 & 660  & -1.90015  \\
  \end{tabular}
\end{table}

\begin{table}\caption{\label{table7}The calculated eigenvalue $E$ from the
 the solution of the Faddeev
equation,  the expectation values of the kinetic energy $\langle H_0
\rangle$, and  the
potential energy $\langle V \rangle$ for the two-body force MT2-I and the three-body
force MMT3-I. 
Both forces give similar binding energies. The expectation values of $\langle R \rangle$ and
 $\langle r \rangle$ are also given for both cases.}


\begin{tabular}{cccccc}
\mbox{}
Model &    $\langle H_0 \rangle$ (MeV)&$\langle V \rangle$ (MeV) &
      $\langle R \rangle$ (fm) & $\langle r \rangle$ (fm) &E (MeV)   \\ \hline
MT2-I &  66.967  &  -74.547 & 1.592  & 1.783   & -7.5803  \\  \hline
MMT3-I  &  67.306   &  -74.895 & 1.521 & 1.698    & -7.5504   \\
  \end{tabular}
  \end{table}



\pagebreak
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\noindent
\begin{figure}
\caption{Diagrammatic representation of the three-body force $V_4^{(1)}$. Here particle
(1) is singled out by the meson-nucleon amplitude described by the blob.
 The three-body force is then given according to Eq.~(\ref{eq:2.2}).
 \label{fig1}}
 \end{figure}

 \begin{figure}
 \caption{The momentum distributions $n(q)$ calculated with the MT2-II
 two-body potential (dotted line). The solid line represents the calculation of $n(q)$
 with the  MT2-II two-body potential and the MT3-I three-body potential, the dashed line
 the corresponding calculation with the MT3-II three body potential.  
 \label{fig2}}
 \end{figure}

 \begin{figure}
 \caption{The two-body correlation function $c(r)$ calculated with 
 the MT2-II two-body potential (dotted line).
 The solid line represents the calculation of $c(r)$ with the  MT2-II two-body potential and
 the MT3-I three-body potential, the dashed line
 the corresponding calculation with the MT3-II three body potential.
 \label{fig5}}
 \end{figure}

\begin{figure}
\caption{The momentum distributions $n(q)$ calculated with the attractive  
2BF MT2-I (solid line). 
The dashed curve represents the calculation with  the purely attractive
3BF, where the parameters are chosen such that the binding energy and the 
momentum distribution are similar to the one given by
the MT2-I potential. 
\label{fig14}}
\end{figure}

\begin{figure}
\caption{The two-body correlation function $c(r)$ calculated with
the  attractive 2BF MT2-I (solid line). The dashed curve represents
the calculation with  the purely attractive
3BF as used in Fig.~\ref{fig14}
  \label{fig15}}
\end{figure}



\newpage

\input{psfig}

{\bf Fig.1}
\centerline{\hspace{15mm}\psfig{file=fig1.eps,width=80mm}}

\vspace{30mm}

{\bf Fig.2}
\centerline{\hspace{15mm}\psfig{file=fig2.ps,width=120mm,angle=-90}}


\newpage

{\bf Fig.3}
\centerline{\hspace{15mm}\psfig{file=fig3.ps,width=120mm,angle=-90}}

\vspace{10mm}

{\bf Fig.4}
\centerline{\hspace{15mm}\psfig{file=fig4.ps,width=120mm,angle=-90}}
\vspace{10mm}

{\bf Fig.5}
\centerline{\hspace{15mm}\psfig{file=fig5.ps,width=120mm,angle=-90}}


\end{document}

psfig.tex
% Psfig/TeX 
\def\PsfigVersion{1.9}
% dvips version
%
% All psfig/tex software, documentation, and related files
% in this distribution of psfig/tex are 
% Copyright 1987, 1988, 1991 Trevor J. Darrell
%
% Permission is granted for use and non-profit distribution of psfig/tex 
% providing that this notice is clearly maintained. The right to
% distribute any portion of psfig/tex for profit or as part of any commercial
% product is specifically reserved for the author(s) of that portion.
%
% *** Feel free to make local modifications of psfig as you wish,
% *** but DO NOT post any changed or modified versions of ``psfig''
% *** directly to the net. Send them to me and I'll try to incorporate
% *** them into future versions. If you want to take the psfig code 
% *** and make a new program (subject to the copyright above), distribute it, 
% *** (and maintain it) that's fine, just don't call it psfig.
%
% Bugs and improvements to trevor@media.mit.edu.
%
% Thanks to Greg Hager (GDH) and Ned Batchelder for their contributions
% to the original version of this project.
%
% Modified by J. Daniel Smith on 9 October 1990 to accept the
% %%BoundingBox: comment with or without a space after the colon.  Stole
% file reading code from Tom Rokicki's EPSF.TEX file (see below).
%
% More modifications by J. Daniel Smith on 29 March 1991 to allow the
% the included PostScript figure to be rotated.  The amount of
% rotation is specified by the "angle=" parameter of the \psfig command.
%
% Modified by Robert Russell on June 25, 1991 to allow users to specify
% .ps filenames which don't yet exist, provided they explicitly provide
% boundingbox information via the \psfig command. Note: This will only work
% if the "file=" parameter follows all four "bb???=" parameters in the
% command. This is due to the order in which psfig interprets these params.
%
%  3 Jul 1991	JDS	check if file already read in once
%  4 Sep 1991	JDS	fixed incorrect computation of rotated
%			bounding box
% 25 Sep 1991	GVR	expanded synopsis of \psfig
% 14 Oct 1991	JDS	\fbox code from LaTeX so \psdraft works with TeX
%			changed \typeout to \ps@typeout
% 17 Oct 1991	JDS	added \psscalefirst and \psrotatefirst
%

% From: gvr@cs.brown.edu (George V. Reilly)
%
% \psdraft	draws an outline box, but doesn't include the figure
%		in the DVI file.  Useful for previewing.
%
% \psfull	includes the figure in the DVI file (default).
%
% \psscalefirst width= or height= specifies the size of the figure
% 		before rotation.
% \psrotatefirst (default) width= or height= specifies the size of the
% 		 figure after rotation.  Asymetric figures will
% 		 appear to shrink.
%
% \psfigurepath#1	sets the path to search for the figure
%
% \psfig
% usage: \psfig{file=, figure=, height=, width=,
%			bbllx=, bblly=, bburx=, bbury=,
%			rheight=, rwidth=, clip=, angle=, silent=}
%
%	"file" is the filename.  If no path name is specified and the
%		file is not found in the current directory,
%		it will be looked for in directory \psfigurepath.
%	"figure" is a synonym for "file".
%	By default, the width and height of the figure are taken from
%		the BoundingBox of the figure.
%	If "width" is specified, the figure is scaled so that it has
%		the specified width.  Its height changes proportionately.
%	If "height" is specified, the figure is scaled so that it has
%		the specified height.  Its width changes proportionately.
%	If both "width" and "height" are specified, the figure is scaled
%		anamorphically.
%	"bbllx", "bblly", "bburx", and "bbury" control the PostScript
%		BoundingBox.  If these four values are specified
%               *before* the "file" option, the PSFIG will not try to
%               open the PostScript file.
%	"rheight" and "rwidth" are the reserved height and width
%		of the figure, i.e., how big TeX actually thinks
%		the figure is.  They default to "width" and "height".
%	The "clip" option ensures that no portion of the figure will
%		appear outside its BoundingBox.  "clip=" is a switch and
%		takes no value, but the `=' must be present.
%	The "angle" option specifies the angle of rotation (degrees, ccw).
%	The "silent" option makes \psfig work silently.
%

% check to see if macros already loaded in (maybe some other file says
% "\input psfig") ...
\ifx\undefined\psfig\else\endinput\fi

%
% from a suggestion by eijkhout@csrd.uiuc.edu to allow
% loading as a style file. Changed to avoid problems
% with amstex per suggestion by jbence@math.ucla.edu

\let\LaTeXAtSign=\@
\let\@=\relax
\edef\psfigRestoreAt{\catcode`\@=\number\catcode`@\relax}
%\edef\psfigRestoreAt{\catcode`@=\number\catcode`@\relax}
\catcode`\@=11\relax
\newwrite\@unused
\def\ps@typeout#1{{\let\protect\string\immediate\write\@unused{#1}}}
\ps@typeout{psfig/tex \PsfigVersion}

%% Here's how you define your figure path.  Should be set up with null
%% default and a user useable definition.

\def\figurepath{./}
\def\psfigurepath#1{\edef\figurepath{#1}}

%
% @psdo control structure -- similar to Latex @for.
% I redefined these with different names so that psfig can
% be used with TeX as well as LaTeX, and so that it will not 
% be vunerable to future changes in LaTeX's internal
% control structure,
%
\def\@nnil{\@nil}
\def\@empty{}
\def\@psdonoop#1\@@#2#3{}
\def\@psdo#1:=#2\do#3{\edef\@psdotmp{#2}\ifx\@psdotmp\@empty \else
    \expandafter\@psdoloop#2,\@nil,\@nil\@@#1{#3}\fi}
\def\@psdoloop#1,#2,#3\@@#4#5{\def#4{#1}\ifx #4\@nnil \else
       #5\def#4{#2}\ifx #4\@nnil \else#5\@ipsdoloop #3\@@#4{#5}\fi\fi}
\def\@ipsdoloop#1,#2\@@#3#4{\def#3{#1}\ifx #3\@nnil 
       \let\@nextwhile=\@psdonoop \else
      #4\relax\let\@nextwhile=\@ipsdoloop\fi\@nextwhile#2\@@#3{#4}}
\def\@tpsdo#1:=#2\do#3{\xdef\@psdotmp{#2}\ifx\@psdotmp\@empty \else
    \@tpsdoloop#2\@nil\@nil\@@#1{#3}\fi}
\def\@tpsdoloop#1#2\@@#3#4{\def#3{#1}\ifx #3\@nnil 
       \let\@nextwhile=\@psdonoop \else
      #4\relax\let\@nextwhile=\@tpsdoloop\fi\@nextwhile#2\@@#3{#4}}
% 
% \fbox is defined in latex.tex; so if \fbox is undefined, assume that
% we are not in LaTeX.
% Perhaps this could be done better???
\ifx\undefined\fbox
% \fbox code from modified slightly from LaTeX
\newdimen\fboxrule
\newdimen\fboxsep
\newdimen\ps@tempdima
\newbox\ps@tempboxa
\fboxsep = 3pt
\fboxrule = .4pt
\long\def\fbox#1{\leavevmode\setbox\ps@tempboxa\hbox{#1}\ps@tempdima\fboxrule
    \advance\ps@tempdima \fboxsep \advance\ps@tempdima \dp\ps@tempboxa
   \hbox{\lower \ps@tempdima\hbox
  {\vbox{\hrule height \fboxrule
          \hbox{\vrule width \fboxrule \hskip\fboxsep
          \vbox{\vskip\fboxsep \box\ps@tempboxa\vskip\fboxsep}\hskip 
                 \fboxsep\vrule width \fboxrule}
                 \hrule height \fboxrule}}}}
\fi
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% file reading stuff from epsf.tex
%   EPSF.TEX macro file:
%   Written by Tomas Rokicki of Radical Eye Software, 29 Mar 1989.
%   Revised by Don Knuth, 3 Jan 1990.
%   Revised by Tomas Rokicki to accept bounding boxes with no
%      space after the colon, 18 Jul 1990.
%   Portions modified/removed for use in PSFIG package by
%      J. Daniel Smith, 9 October 1990.
%
\newread\ps@stream
\newif\ifnot@eof       % continue looking for the bounding box?
\newif\if@noisy        % report what you're making?
\newif\if@atend        % %%BoundingBox: has (at end) specification
\newif\if@psfile       % does this look like a PostScript file?
%
% PostScript files should start with `%!'
%
{\catcode`\%=12\global\gdef\epsf@start{%!}}
\def\epsf@PS{PS}
%
\def\epsf@getbb#1{%
%
%   The first thing we need to do is to open the
%   PostScript file, if possible.
%
\openin\ps@stream=#1
\ifeof\ps@stream\ps@typeout{Error, File #1 not found}\else
%
%   Okay, we got it. Now we'll scan lines until we find one that doesn't
%   start with %. We're looking for the bounding box comment.
%
   {\not@eoftrue \chardef\other=12
    \def\do##1{\catcode`##1=\other}\dospecials \catcode`\ =10
    \loop
       \if@psfile
	  \read\ps@stream to \epsf@fileline
       \else{
	  \obeyspaces
          \read\ps@stream to \epsf@tmp\global\let\epsf@fileline\epsf@tmp}
       \fi
       \ifeof\ps@stream\not@eoffalse\else
%
%   Check the first line for `%!'.  Issue a warning message if its not
%   there, since the file might not be a PostScript file.
%
       \if@psfile\else
       \expandafter\epsf@test\epsf@fileline:. \\%
       \fi
%
%   We check to see if the first character is a % sign;
%   if so, we look further and stop only if the line begins with
%   `%%BoundingBox:' and the `(atend)' specification was not found.
%   That is, the only way to stop is when the end of file is reached,
%   or a `%%BoundingBox: llx lly urx ury' line is found.
%
          \expandafter\epsf@aux\epsf@fileline:. \\%
       \fi
   \ifnot@eof\repeat
   }\closein\ps@stream\fi}%
%
% This tests if the file we are reading looks like a PostScript file.
%
\long\def\epsf@test#1#2#3:#4\\{\def\epsf@testit{#1#2}
			\ifx\epsf@testit\epsf@start\else
\ps@typeout{Warning! File does not start with `\epsf@start'.  It may not be a PostScript file.}
			\fi
			\@psfiletrue} % don't test after 1st line
%
%   We still need to define the tricky \epsf@aux macro. This requires
%   a couple of magic constants for comparison purposes.
%
{\catcode`\%=12\global\let\epsf@percent=%\global\def\epsf@bblit{%BoundingBox}}
%
%
%   So we're ready to check for `%BoundingBox:' and to grab the
%   values if they are found.  We continue searching if `(at end)'
%   was found after the `%BoundingBox:'.
%
\long\def\epsf@aux#1#2:#3\\{\ifx#1\epsf@percent
   \def\epsf@testit{#2}\ifx\epsf@testit\epsf@bblit
	\@atendfalse
        \epsf@atend #3 . \\%
	\if@atend	
	   \if@verbose{
		\ps@typeout{psfig: found `(atend)'; continuing search}
	   }\fi
        \else
        \epsf@grab #3 . . . \\%
        \not@eoffalse
        \global\no@bbfalse
        \fi
   \fi\fi}%
%
%   Here we grab the values and stuff them in the appropriate definitions.
%
\def\epsf@grab #1 #2 #3 #4 #5\\{%
   \global\def\epsf@llx{#1}\ifx\epsf@llx\empty
      \epsf@grab #2 #3 #4 #5 .\\\else
   \global\def\epsf@lly{#2}%
   \global\def\epsf@urx{#3}\global\def\epsf@ury{#4}\fi}%
%
% Determine if the stuff following the %%BoundingBox is `(atend)'
% J. Daniel Smith.  Copied from \epsf@grab above.
%
\def\epsf@atendlit{(atend)} 
\def\epsf@atend #1 #2 #3\\{%
   \def\epsf@tmp{#1}\ifx\epsf@tmp\empty
      \epsf@atend #2 #3 .\\\else
   \ifx\epsf@tmp\epsf@atendlit\@atendtrue\fi\fi}


% End of file reading stuff from epsf.tex
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% trigonometry stuff from "trig.tex"
\chardef\psletter = 11 % won't conflict with \begin{letter} now...
\chardef\other = 12

\newif \ifdebug %%% turn me on to see TeX hard at work ...
\newif\ifc@mpute %%% don't need to compute some values
\c@mputetrue % but assume that we do

\let\then = \relax
\def\r@dian{pt }
\let\r@dians = \r@dian
\let\dimensionless@nit = \r@dian
\let\dimensionless@nits = \dimensionless@nit
\def\internal@nit{sp }
\let\internal@nits = \internal@nit
\newif\ifstillc@nverging
\def \Mess@ge #1{\ifdebug \then \message {#1} \fi}

{ %%% Things that need abnormal catcodes %%%
	\catcode `\@ = \psletter
	\gdef \nodimen {\expandafter \n@dimen \the \dimen}
	\gdef \term #1 #2 #3%
	       {\edef \t@ {\the #1}%%% freeze parameter 1 (count, by value)
		\edef \t@@ {\expandafter \n@dimen \the #2\r@dian}%
				   %%% freeze parameter 2 (dimen, by value)
		\t@rm {\t@} {\t@@} {#3}%
	       }
	\gdef \t@rm #1 #2 #3%
	       {{%
		\count 0 = 0
		\dimen 0 = 1 \dimensionless@nit
		\dimen 2 = #2\relax
		\Mess@ge {Calculating term #1 of \nodimen 2}%
		\loop
		\ifnum	\count 0 < #1
		\then	\advance \count 0 by 1
			\Mess@ge {Iteration \the \count 0 \space}%
			\Multiply \dimen 0 by {\dimen 2}%
			\Mess@ge {After multiplication, term = \nodimen 0}%
			\Divide \dimen 0 by {\count 0}%
			\Mess@ge {After division, term = \nodimen 0}%
		\repeat
		\Mess@ge {Final value for term #1 of 
				\nodimen 2 \space is \nodimen 0}%
		\xdef \Term {#3 = \nodimen 0 \r@dians}%
		\aftergroup \Term
	       }}
	\catcode `\p = \other
	\catcode `\t = \other
	\gdef \n@dimen #1pt{#1} %%% throw away the ``pt''
}

\def \Divide #1by #2{\divide #1 by #2} %%% just a synonym

\def \Multiply #1by #2%%% allows division of a dimen by a dimen
       {{%%% should really freeze parameter 2 (dimen, passed by value)
	\count 0 = #1\relax
	\count 2 = #2\relax
	\count 4 = 65536
	\Mess@ge {Before scaling, count 0 = \the \count 0 \space and
			count 2 = \the \count 2}%
	\ifnum	\count 0 > 32767 %%% do our best to avoid overflow
	\then	\divide \count 0 by 4
		\divide \count 4 by 4
	\else	\ifnum	\count 0 < -32767
		\then	\divide \count 0 by 4
			\divide \count 4 by 4
		\else
		\fi
	\fi
	\ifnum	\count 2 > 32767 %%% while retaining reasonable accuracy
	\then	\divide \count 2 by 4
		\divide \count 4 by 4
	\else	\ifnum	\count 2 < -32767
		\then	\divide \count 2 by 4
			\divide \count 4 by 4
		\else
		\fi
	\fi
	\multiply \count 0 by \count 2
	\divide \count 0 by \count 4
	\xdef \product {#1 = \the \count 0 \internal@nits}%
	\aftergroup \product
       }}

\def\r@duce{\ifdim\dimen0 > 90\r@dian \then   % sin(x+90) = sin(180-x)
		\multiply\dimen0 by -1
		\advance\dimen0 by 180\r@dian
		\r@duce
	    \else \ifdim\dimen0 < -90\r@dian \then  % sin(-x) = sin(360+x)
		\advance\dimen0 by 360\r@dian
		\r@duce
		\fi
	    \fi}

\def\Sine#1%
       {{%
	\dimen 0 = #1 \r@dian
	\r@duce
	\ifdim\dimen0 = -90\r@dian \then
	   \dimen4 = -1\r@dian
	   \c@mputefalse
	\fi
	\ifdim\dimen0 = 90\r@dian \then
	   \dimen4 = 1\r@dian
	   \c@mputefalse
	\fi
	\ifdim\dimen0 = 0\r@dian \then
	   \dimen4 = 0\r@dian
	   \c@mputefalse
	\fi
%
	\ifc@mpute \then
        	% convert degrees to radians
		\divide\dimen0 by 180
		\dimen0=3.141592654\dimen0
%
		\dimen 2 = 3.1415926535897963\r@dian %%% a well-known constant
		\divide\dimen 2 by 2 %%% we only deal with -pi/2 : pi/2
		\Mess@ge {Sin: calculating Sin of \nodimen 0}%
		\count 0 = 1 %%% see power-series expansion for sine
		\dimen 2 = 1 \r@dian %%% ditto
		\dimen 4 = 0 \r@dian %%% ditto
		\loop
			\ifnum	\dimen 2 = 0 %%% then we've done
			\then	\stillc@nvergingfalse 
			\else	\stillc@nvergingtrue
			\fi
			\ifstillc@nverging %%% then calculate next term
			\then	\term {\count 0} {\dimen 0} {\dimen 2}%
				\advance \count 0 by 2
				\count 2 = \count 0
				\divide \count 2 by 2
				\ifodd	\count 2 %%% signs alternate
				\then	\advance \dimen 4 by \dimen 2
				\else	\advance \dimen 4 by -\dimen 2
				\fi
		\repeat
	\fi		
			\xdef \sine {\nodimen 4}%
       }}

% Now the Cosine can be calculated easily by calling \Sine
\def\Cosine#1{\ifx\sine\UnDefined\edef\Savesine{\relax}\else
		             \edef\Savesine{\sine}\fi
	{\dimen0=#1\r@dian\advance\dimen0 by 90\r@dian
	 \Sine{\nodimen 0}
	 \xdef\cosine{\sine}
	 \xdef\sine{\Savesine}}}	      
% end of trig stuff
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\def\psdraft{
	\def\@psdraft{0}
	%\ps@typeout{draft level now is \@psdraft \space . }
}
\def\psfull{
	\def\@psdraft{100}
	%\ps@typeout{draft level now is \@psdraft \space . }
}

\psfull

\newif\if@scalefirst
\def\psscalefirst{\@scalefirsttrue}
\def\psrotatefirst{\@scalefirstfalse}
\psrotatefirst

\newif\if@draftbox
\def\psnodraftbox{
	\@draftboxfalse
}
\def\psdraftbox{
	\@draftboxtrue
}
\@draftboxtrue

\newif\if@prologfile
\newif\if@postlogfile
\def\pssilent{
	\@noisyfalse
}
\def\psnoisy{
	\@noisytrue
}
\psnoisy
%%% These are for the option list.
%%% A specification of the form a = b maps to calling \@p@@sa{b}
\newif\if@bbllx
\newif\if@bblly
\newif\if@bburx
\newif\if@bbury
\newif\if@height
\newif\if@width
\newif\if@rheight
\newif\if@rwidth
\newif\if@angle
\newif\if@clip
\newif\if@verbose
\def\@p@@sclip#1{\@cliptrue}


\newif\if@decmpr

%%% GDH 7/26/87 -- changed so that it first looks in the local directory,
%%% then in a specified global directory for the ps file.
%%% RPR 6/25/91 -- changed so that it defaults to user-supplied name if
%%% boundingbox info is specified, assuming graphic will be created by
%%% print time.
%%% TJD 10/19/91 -- added bbfile vs. file distinction, and @decmpr flag

\def\@p@@sfigure#1{\def\@p@sfile{null}\def\@p@sbbfile{null}
	        \openin1=#1.bb
		\ifeof1\closein1
	        	\openin1=\figurepath#1.bb
			\ifeof1\closein1
			        \openin1=#1
				\ifeof1\closein1%
				       \openin1=\figurepath#1
					\ifeof1
					   \ps@typeout{Error, File #1 not found}
						\if@bbllx\if@bblly
				   		\if@bburx\if@bbury
			      				\def\@p@sfile{#1}%
			      				\def\@p@sbbfile{#1}%
							\@decmprfalse
				  	   	\fi\fi\fi\fi
					\else\closein1
				    		\def\@p@sfile{\figurepath#1}%
				    		\def\@p@sbbfile{\figurepath#1}%
						\@decmprfalse
	                       		\fi%
			 	\else\closein1%
					\def\@p@sfile{#1}
					\def\@p@sbbfile{#1}
					\@decmprfalse
			 	\fi
			\else
				\def\@p@sfile{\figurepath#1}
				\def\@p@sbbfile{\figurepath#1.bb}
				\@decmprtrue
			\fi
		\else
			\def\@p@sfile{#1}
			\def\@p@sbbfile{#1.bb}
			\@decmprtrue
		\fi}

\def\@p@@sfile#1{\@p@@sfigure{#1}}

\def\@p@@sbbllx#1{
		%\ps@typeout{bbllx is #1}
		\@bbllxtrue
		\dimen100=#1
		\edef\@p@sbbllx{\number\dimen100}
}
\def\@p@@sbblly#1{
		%\ps@typeout{bblly is #1}
		\@bbllytrue
		\dimen100=#1
		\edef\@p@sbblly{\number\dimen100}
}
\def\@p@@sbburx#1{
		%\ps@typeout{bburx is #1}
		\@bburxtrue
		\dimen100=#1
		\edef\@p@sbburx{\number\dimen100}
}
\def\@p@@sbbury#1{
		%\ps@typeout{bbury is #1}
		\@bburytrue
		\dimen100=#1
		\edef\@p@sbbury{\number\dimen100}
}
\def\@p@@sheight#1{
		\@heighttrue
		\dimen100=#1
   		\edef\@p@sheight{\number\dimen100}
		%\ps@typeout{Height is \@p@sheight}
}
\def\@p@@swidth#1{
		%\ps@typeout{Width is #1}
		\@widthtrue
		\dimen100=#1
		\edef\@p@swidth{\number\dimen100}
}
\def\@p@@srheight#1{
		%\ps@typeout{Reserved height is #1}
		\@rheighttrue
		\dimen100=#1
		\edef\@p@srheight{\number\dimen100}
}
\def\@p@@srwidth#1{
		%\ps@typeout{Reserved width is #1}
		\@rwidthtrue
		\dimen100=#1
		\edef\@p@srwidth{\number\dimen100}
}
\def\@p@@sangle#1{
		%\ps@typeout{Rotation is #1}
		\@angletrue
%		\dimen100=#1
		\edef\@p@sangle{#1} %\number\dimen100}
}
\def\@p@@ssilent#1{ 
		\@verbosefalse
}
\def\@p@@sprolog#1{\@prologfiletrue\def\@prologfileval{#1}}
\def\@p@@spostlog#1{\@postlogfiletrue\def\@postlogfileval{#1}}
\def\@cs@name#1{\csname #1\endcsname}
\def\@setparms#1=#2,{\@cs@name{@p@@s#1}{#2}}
%
% initialize the defaults (size the size of the figure)
%
\def\ps@init@parms{
		\@bbllxfalse \@bbllyfalse
		\@bburxfalse \@bburyfalse
		\@heightfalse \@widthfalse
		\@rheightfalse \@rwidthfalse
		\def\@p@sbbllx{}\def\@p@sbblly{}
		\def\@p@sbburx{}\def\@p@sbbury{}
		\def\@p@sheight{}\def\@p@swidth{}
		\def\@p@srheight{}\def\@p@srwidth{}
		\def\@p@sangle{0}
		\def\@p@sfile{} \def\@p@sbbfile{}
		\def\@p@scost{10}
		\def\@sc{}
		\@prologfilefalse
		\@postlogfilefalse
		\@clipfalse
		\if@noisy
			\@verbosetrue
		\else
			\@verbosefalse
		\fi
}
%
% Go through the options setting things up.
%
\def\parse@ps@parms#1{
	 	\@psdo\@psfiga:=#1\do
		   {\expandafter\@setparms\@psfiga,}}
%
% Compute bb height and width
%
\newif\ifno@bb
\def\bb@missing{
	\if@verbose{
		\ps@typeout{psfig: searching \@p@sbbfile \space  for bounding box}
	}\fi
	\no@bbtrue
	\epsf@getbb{\@p@sbbfile}
        \ifno@bb \else \bb@cull\epsf@llx\epsf@lly\epsf@urx\epsf@ury\fi
}	
\def\bb@cull#1#2#3#4{
	\dimen100=#1 bp\edef\@p@sbbllx{\number\dimen100}
	\dimen100=#2 bp\edef\@p@sbblly{\number\dimen100}
	\dimen100=#3 bp\edef\@p@sbburx{\number\dimen100}
	\dimen100=#4 bp\edef\@p@sbbury{\number\dimen100}
	\no@bbfalse
}
% rotate point (#1,#2) about (0,0).
% The sine and cosine of the angle are already stored in \sine and
% \cosine.  The result is placed in (\p@intvaluex, \p@intvaluey).
\newdimen\p@intvaluex
\newdimen\p@intvaluey
\def\rotate@#1#2{{\dimen0=#1 sp\dimen1=#2 sp
%            	calculate x' = x \cos\theta - y \sin\theta
		  \global\p@intvaluex=\cosine\dimen0
		  \dimen3=\sine\dimen1
		  \global\advance\p@intvaluex by -\dimen3
% 		calculate y' = x \sin\theta + y \cos\theta
		  \global\p@intvaluey=\sine\dimen0
		  \dimen3=\cosine\dimen1
		  \global\advance\p@intvaluey by \dimen3
		  }}
\def\compute@bb{
		\no@bbfalse
		\if@bbllx \else \no@bbtrue \fi
		\if@bblly \else \no@bbtrue \fi
		\if@bburx \else \no@bbtrue \fi
		\if@bbury \else \no@bbtrue \fi
		\ifno@bb \bb@missing \fi
		\ifno@bb \ps@typeout{FATAL ERROR: no bb supplied or found}
			\no-bb-error
		\fi
		%
%\ps@typeout{BB: \@p@sbbllx, \@p@sbblly, \@p@sbburx, \@p@sbbury} 
%
% store height/width of original (unrotated) bounding box
		\count203=\@p@sbburx
		\count204=\@p@sbbury
		\advance\count203 by -\@p@sbbllx
		\advance\count204 by -\@p@sbblly
		\edef\ps@bbw{\number\count203}
		\edef\ps@bbh{\number\count204}
		%\ps@typeout{ psbbh = \ps@bbh, psbbw = \ps@bbw }
		\if@angle 
			\Sine{\@p@sangle}\Cosine{\@p@sangle}
	        	{\dimen100=\maxdimen\xdef\r@p@sbbllx{\number\dimen100}
					    \xdef\r@p@sbblly{\number\dimen100}
			                    \xdef\r@p@sbburx{-\number\dimen100}
					    \xdef\r@p@sbbury{-\number\dimen100}}
%
% Need to rotate all four points and take the X-Y extremes of the new
% points as the new bounding box.
                        \def\minmaxtest{
			   \ifnum\number\p@intvaluex<\r@p@sbbllx
			      \xdef\r@p@sbbllx{\number\p@intvaluex}\fi
			   \ifnum\number\p@intvaluex>\r@p@sbburx
			      \xdef\r@p@sbburx{\number\p@intvaluex}\fi
			   \ifnum\number\p@intvaluey<\r@p@sbblly
			      \xdef\r@p@sbblly{\number\p@intvaluey}\fi
			   \ifnum\number\p@intvaluey>\r@p@sbbury
			      \xdef\r@p@sbbury{\number\p@intvaluey}\fi
			   }
%			lower left
			\rotate@{\@p@sbbllx}{\@p@sbblly}
			\minmaxtest
%			upper left
			\rotate@{\@p@sbbllx}{\@p@sbbury}
			\minmaxtest
%			lower right
			\rotate@{\@p@sbburx}{\@p@sbblly}
			\minmaxtest
%			upper right
			\rotate@{\@p@sbburx}{\@p@sbbury}
			\minmaxtest
			\edef\@p@sbbllx{\r@p@sbbllx}\edef\@p@sbblly{\r@p@sbblly}
			\edef\@p@sbburx{\r@p@sbburx}\edef\@p@sbbury{\r@p@sbbury}
%\ps@typeout{rotated BB: \r@p@sbbllx, \r@p@sbblly, \r@p@sbburx, \r@p@sbbury}
		\fi
		\count203=\@p@sbburx
		\count204=\@p@sbbury
		\advance\count203 by -\@p@sbbllx
		\advance\count204 by -\@p@sbblly
		\edef\@bbw{\number\count203}
		\edef\@bbh{\number\count204}
		%\ps@typeout{ bbh = \@bbh, bbw = \@bbw }
}
%
% \in@hundreds performs #1 * (#2 / #3) correct to the hundreds,
%	then leaves the result in @result
%
\def\in@hundreds#1#2#3{\count240=#2 \count241=#3
		     \count100=\count240	% 100 is first digit #2/#3
		     \divide\count100 by \count241
		     \count101=\count100
		     \multiply\count101 by \count241
		     \advance\count240 by -\count101
		     \multiply\count240 by 10
		     \count101=\count240	%101 is second digit of #2/#3
		     \divide\count101 by \count241
		     \count102=\count101
		     \multiply\count102 by \count241
		     \advance\count240 by -\count102
		     \multiply\count240 by 10
		     \count102=\count240	% 102 is the third digit
		     \divide\count102 by \count241
		     \count200=#1\count205=0
		     \count201=\count200
			\multiply\count201 by \count100
		 	\advance\count205 by \count201
		     \count201=\count200
			\divide\count201 by 10
			\multiply\count201 by \count101
			\advance\count205 by \count201
			%
		     \count201=\count200
			\divide\count201 by 100
			\multiply\count201 by \count102
			\advance\count205 by \count201
			%
		     \edef\@result{\number\count205}
}
\def\compute@wfromh{
		% computing : width = height * (bbw / bbh)
		\in@hundreds{\@p@sheight}{\@bbw}{\@bbh}
		%\ps@typeout{ \@p@sheight * \@bbw / \@bbh, = \@result }
		\edef\@p@swidth{\@result}
		%\ps@typeout{w from h: width is \@p@swidth}
}
\def\compute@hfromw{
		% computing : height = width * (bbh / bbw)
	        \in@hundreds{\@p@swidth}{\@bbh}{\@bbw}
		%\ps@typeout{ \@p@swidth * \@bbh / \@bbw = \@result }
		\edef\@p@sheight{\@result}
		%\ps@typeout{h from w : height is \@p@sheight}
}
\def\compute@handw{
		\if@height 
			\if@width
			\else
				\compute@wfromh
			\fi
		\else 
			\if@width
				\compute@hfromw
			\else
				\edef\@p@sheight{\@bbh}
				\edef\@p@swidth{\@bbw}
			\fi
		\fi
}
\def\compute@resv{
		\if@rheight \else \edef\@p@srheight{\@p@sheight} \fi
		\if@rwidth \else \edef\@p@srwidth{\@p@swidth} \fi
		%\ps@typeout{rheight = \@p@srheight, rwidth = \@p@srwidth}
}
%		
% Compute any missing values
\def\compute@sizes{
	\compute@bb
	\if@scalefirst\if@angle
% at this point the bounding box has been adjsuted correctly for
% rotation.  PSFIG does all of its scaling using \@bbh and \@bbw.  If
% a width= or height= was specified along with \psscalefirst, then the
% width=/height= value needs to be adjusted to match the new (rotated)
% bounding box size (specifed in \@bbw and \@bbh).
%    \ps@bbw       width=
%    -------  =  ---------- 
%    \@bbw       new width=
% so `new width=' = (width= * \@bbw) / \ps@bbw; where \ps@bbw is the
% width of the original (unrotated) bounding box.
	\if@width
	   \in@hundreds{\@p@swidth}{\@bbw}{\ps@bbw}
	   \edef\@p@swidth{\@result}
	\fi
	\if@height
	   \in@hundreds{\@p@sheight}{\@bbh}{\ps@bbh}
	   \edef\@p@sheight{\@result}
	\fi
	\fi\fi
	\compute@handw
	\compute@resv}

%
% \psfig
% usage : \psfig{file=, height=, width=, bbllx=, bblly=, bburx=, bbury=,
%			rheight=, rwidth=, clip=}
%
% "clip=" is a switch and takes no value, but the `=' must be present.
\def\psfig#1{\vbox {
	% do a zero width hard space so that a single
	% \psfig in a centering enviornment will behave nicely
	%{\setbox0=\hbox{\ }\ \hskip-\wd0}
	%
	\ps@init@parms
	\parse@ps@parms{#1}
	\compute@sizes
	%
	\ifnum\@p@scost<\@psdraft{
		%
		\special{ps::[begin] 	\@p@swidth \space \@p@sheight \space
				\@p@sbbllx \space \@p@sbblly \space
				\@p@sbburx \space \@p@sbbury \space
				startTexFig \space }
		\if@angle
			\special {ps:: \@p@sangle \space rotate \space} 
		\fi
		\if@clip{
			\if@verbose{
				\ps@typeout{(clip)}
			}\fi
			\special{ps:: doclip \space }
		}\fi
		\if@prologfile
		    \special{ps: plotfile \@prologfileval \space } \fi
		\if@decmpr{
			\if@verbose{
				\ps@typeout{psfig: including \@p@sfile.Z \space }
			}\fi
			\special{ps: plotfile "`zcat \@p@sfile.Z" \space }
		}\else{
			\if@verbose{
				\ps@typeout{psfig: including \@p@sfile \space }
			}\fi
			\special{ps: plotfile \@p@sfile \space }
		}\fi
		\if@postlogfile
		    \special{ps: plotfile \@postlogfileval \space } \fi
		\special{ps::[end] endTexFig \space }
		% Create the vbox to reserve the space for the figure.
		\vbox to \@p@srheight sp{
		% 1/92 TJD Changed from "true sp" to "sp" for magnification.
			\hbox to \@p@srwidth sp{
				\hss
			}
		\vss
		}
	}\else{
		% draft figure, just reserve the space and print the
		% path name.
		\if@draftbox{		
			% Verbose draft: print file name in box
			\hbox{\frame{\vbox to \@p@srheight sp{
			\vss
			\hbox to \@p@srwidth sp{ \hss \@p@sfile \hss }
			\vss
			}}}
		}\else{
			% Non-verbose draft
			\vbox to \@p@srheight sp{
			\vss
			\hbox to \@p@srwidth sp{\hss}
			\vss
			}
		}\fi	



	}\fi
}}
\psfigRestoreAt
\let\@=\LaTeXAtSign




nucl-th0207062.extracted_bib
\bibitem{sofianos} V.A. Roudnev, S.L. Yakovlev, S.A. Sofianos, `Bound State Calculations for
%Three Atoms without Explicit Partial Wave Decomposition', Preprint
%xxx.lanl.gov/physics/0204025.


\bibitem{malfliet69-1} R.A. Malfliet, A.J. Tjon, Nucl. Phys. A, {\bf 127}, 161 (1969).

\bibitem{osborn} T.A. Osborn, `Faddeev�s Equations for local potentials�, SLAC Report No.
17 (1967).


\bibitem{pickle} A.~Picklesimer, R.A.~Rice, and R.~Brandenburg, Phys.
Rev C{\bf 45}, 2045 (1992), {\it ibid}, 547 (1992), Phys. Rev. C{\bf 44},
1359 (1991).

\bibitem{stadler} A.~Stadler, W.~Gl\"ockle, P.U.~Sauer, Phys. Rev. C{\bf 44}, 2319 (1991); 

\bibitem{stadler2} A.~Stadler and P.U.~Sauer, Phys. Rev. C{\bf 46}, 64 (1992).


\bibitem{nbench} A.~Nogga, D.~H\"uber, H.~Kamada, W.~Gl\"ockle,
Phys. Lett. {\bf B 409}, 19 (1997).


\bibitem{lavern} A. Laverne, C. Gignoux, Nucl. Phys. {\bf A203}, 597 (1973).


\bibitem{payne} J.L.~Friar, B.F.~Gibson, and G.L.~Payne, Z. Phys. {\bf A 301},
309 (1981), C.R.~Chen, G.L.~Payne, J.L.~Friar, B.F.~Gibson, Phys. Rev. C{\bf 31}, 2266 (1985).


\bibitem{schell0} N.W.~Schellingerhout, L.P.~Kok, and G.D.~Bosveld,
Phys. Rev. A{\bf 40}, 5568 (1989); N.W.~Schellingerhout, J.J.~Schut,
and L.P.~Kok, Phys. Rev. C{\bf 46}, 1192 (1992).


\bibitem{wu} Y. Wu, S. Ishikawa, T. Sasakawa, Few-Body Systems {\bf 15},
  145 (1993).


\bibitem{gfmc} J.~Carlson, Phys. Rev. C{\bf 36}, 2026 (1987),
  Phys. Rev. C{\bf 38}, 1879 (1988); J.G.~Zabolitzki, K.E.~Schmidt,
  and M.H.~Kalos, Phys. Rev. C{\bf 25}, 1111 (1982); J.~Carlson and
  R.~Schiavilla, Rev. Mod. Phys. {\bf 70}, 743 (1998).


\bibitem{elster98-1} Ch. Elster, W. Schadow, A. Nogga, and W. 
Gl\"ockle, Few-Body System, {\bf 27}, 83 (1998).


\bibitem{scatt1} W. Schadow, Ch.~Elster, and W.~Gl\"ockle, Few-Body Systems {\bf 28}, 15
(2000).


\bibitem{kvin} A.A.~Kvitsinsky and C.-Y. Hu, Few-Body Systems {\bf 12},
7 (1992).


\bibitem{imam1} I. Fachruddin, Ch.~Elster, W. Gl\"ockle,  Phys. Rev. C{\bf 62}, 044002-1
(2000).


\bibitem{gloeckle-book}W. Gl\"ockle,
The Quantum Mechanical Few-Body Problem, Springer-Verlag (1983).


\bibitem{coon} S.A. Coon and W. Gl\"ockle, Phys. Rev. C{\bf 23}, 1790 (1981).


\bibitem{hueber97-1} D. H\"uber, H. Witala, A. Nogga, W. Gl\"ockle and
H. Kamada, Few-Body System, {\bf 22}, 107 (1997).


\bibitem{Fujita} J. Fujita, H. Miyazawa, Prog. Theor. Phys. {\bf 17}, 360 (1957).


\bibitem{elster96-1} Ch. Elster, E.E. Evans, H. Kamada and W.  Gl\"ockle, 
Few-Body System, {\bf 21}, 25 (1996).


\bibitem{elster98-2} Ch. Elster, J.H. Thomas, W. Gl\"ockle, 
Few-Body Systems, {\bf 24}, 55 (1998).


\bibitem{NoggaPhd} A. Nogga, PhD Thesis, Ruhr-University Bochum, 2001. 


\bibitem{NoggaPL} A. Nogga, D. H\"uber, H. Kamada, W. Gl\"ockle, Phys. Lett. {\bf B409}, 19
(1997).


\bibitem{saake} A. Saake, Diploma Thesis, Ruhr-University Bochum, 1992 (unpublished).


\bibitem{coester} W. Gl\"ockle, T.-S.H. Lee, F. Coester, Phys. Rev. C{\bf 33}, 709
(1986).
 

\bibitem{kamada} H. Kamada, Few-Body Systems, Suppl. {\bf 12}, 433 (2000).

\bibitem{kamada2}  H. Kamada, W. Gl\"ockle, J. Golak, Ch. Elster, Phys. Rev. C{\bf 66}, 044010 (2002).

