# Program Summary

## Program Summary

Title of the Programs: phi3   

Computer: Sun 19   

Operating system: Unix   

Programming language used: FORTRAN 77   

Peripherals used: Laser Printer   

Number of lines in distributed program: 1510   

Keywords: Nonperturbative, bound states, Feynman-Schwinger representation, Monte-Carlo, numerical evaluation.   

Nature of physical problem:   
The program provided here evaluates the mass and distribution probabilities of the fully interacting n-body propagator ($`n3`$) for scalar $`\chi ^2\varphi `$ interaction in 3+1 dimensions. The evaluation takes into account all self energy, vertex dressing, and exchange interaction contributions except those involving matter loops (the quenched approximation).   

Method of solution:   
The Feynman-Schwinger representation approach is used to express the field theoretical Green’s function in terms of a quantum mechanical path integral. The resultant expression is evaluated using a Monte-Carlo simulation.   

Restrictions of the program:   
Only $`n=1,2`$, and $`3`$ body propagators are considered. The extension to $`n4`$ requires straightforward modifications in the program.   

Typical running time:   
About 1 day for the 1-body propagator with the self energy.

## LONG WRITE-UP

## 1 Introduction

In nuclear physics one is often faced by problems that require nonperturbative methods. The best known example is the problem of bound states. Even if the underlying theory may have a small coupling constant (such as in QED), and therefore allows the use of perturbation theory in general, the treatment of bound states are inherently nonperturbative. The n-body bound state is defined by the pole of the interacting n-body propagator. A perturbative approximation of n-body propagator does not produce the bound state pole location. Therefore it is essential that reliable nonperturbative methods that take all orders of interaction into account are developed. For this reason, numerous nonperturbative methods have been developed and successfully used in the literature. Some of the best known examples are lattice gauge theory  (LGT), and relativistic bound state equations .

In a recent paper  we (along with authors J. Tjon, and F. Gross) have discussed yet another method known as the Feynman-Schwinger representation(FSR). The basic idea in the FSR approach  is to integrate out all fields at the expense of introducing quantum mechanical path integrals over the trajectories of particles. Replacing the path integrals over fields with path integrals over trajectories has an enormous computational advantage. The advantage is due to the fact that the path integration over trajectories involves a variation of lines rather than fields in a volume. Therefore the degrees of freedom is considerably fewer. The FSR approach differs from the LGT in that it utilizes a space-time continuum, therefore maintaining the Poincare symmetry. On the other hand it should be pointed out that the FSR approach is not without its drawbacks. In particular, how to extend the FSR approach to include fermions is not known. In the past researchers  have attacked the fermion problem using various approximations. An exact  result involving fermions is an important problem and requires further study.

While being able to calculate a nonperturbative result by itself is interesting, an additional motivation in studying the FSR approach is to determine which subsets of diagrams give the dominant contribution to the n-body propagator. This is particularly important in determining what kind of approximations are reasonable within the context of bound state equations. An example of how the FSR results can be used to compare different nonperturbative approximation schemes was presented in Refs . In those works the emphasis was on the development of the formalism and application to the 1 and 2-body propagators. In Ref , The FSR prediction for the 1-body mass in Scalar QED was compared by the rainbow-Dyson-Schwinger equation prediction. It was found  that while the FSR approach provides a real mass pole for all coupling strengths, the Dyson-Schwinger equation provides a complex mass pole beyond a critical coupling strength. Furthermore it was found that, for Scalar QED in 0+1 dimension, the vertex corrections to the exchange interaction do not contribute to the 2-body binding energy . These examples demonstrate the potential usefulness of the FSR approach. The knowledge about the nonperturbative propagators and vertices provided by the FSR is valuable as an input in testing and improving the modeling of other nonperturbative approaches such as Dyson-Schwinger equations. 

The two and three-body bound state sectors provide possibilities for the application of the FSR formalism. In particular it is important to see how various bound state equations (such as Bethe-Salpeter, Gross (spectator), Blankenbekler-Sugar, Equal-time, and nonrelativistic) compare with the quenched FSR results. Applications of the FSR approach to 2 and 3 body states with comparisons to various bound state equation results in $`\chi ^2\varphi `$ theory is currently under study and will be presented in a separate article. 

In this paper we present a complete numerical algorithm for the evaluation of $`n`$-body masses ($`n3`$) and distribution probabilities for scalar $`\chi ^2\varphi `$ interaction in 3+1 dimensions. By providing this algorithm we intend to facilitate the comparison of various nonperturbative methods with the exact quenched results in $`\chi ^2\varphi `$ theory. The organization of the paper is as follows: In the next section we present a brief summary of the FSR formalism. In the third section, using various 1, 2, and 3-body cases as examples, we discuss how results are obtained. And in the fourth section we explain the components of the program.

## 2 The Feynman-Schwinger representation for scalar fields

We consider the theory of charged scalar particles $`\chi `$ of mass $`m`$ interacting through the exchange of a neutral scalar particle $`\varphi `$ of mass $`\mu `$. The Euclidean Lagrangian for this theory is given by

$$_E=\chi ^{}\left[m^2^2+g\varphi \right]\chi +\frac{1}{2}\varphi (\mu ^2^2)\varphi .$$
(1)   
The two body Green’s function for the transition from the initial state $`\mathrm{\Phi }_i=\chi ^{}(x)\chi (\overline{x})`$ to final state $`\mathrm{\Phi }_f=\chi ^{}(y)\chi (\overline{y})`$ is given by

$$G(y,\overline{y}|x,\overline{x})=N𝒟\chi ^{}𝒟\chi 𝒟\varphi \mathrm{\Phi }_f^{}\mathrm{\Phi }_ie^{S_E}.$$
(2)   
The final result for the two-body propagator involves a quantum mechanical path integral that sums up contributions coming from all possible trajectories of particles

$$G=_0^{\mathrm{}}𝑑s_0^{\mathrm{}}𝑑\overline{s}(𝒟z)_{xy}(𝒟\overline{z})_{\overline{x}\overline{y}}e^{S[Z]},$$
(3)   
where $`S[Z]`$ is given by

$$S[Z]iK[z,s]iK[\overline{z},\overline{s}]+V_0[z,s]+2V_{12}[z,\overline{z},s,\overline{s}]+V_0[\overline{z},\overline{s}].$$
(4)   
where

$`K[z,s]`$ $`=`$ $`(m^2+iϵ)s{\displaystyle \frac{1}{4s}}{\displaystyle _0^1}𝑑\tau {\displaystyle \frac{dz_\mu (\tau )}{d\tau }}{\displaystyle \frac{dz^\mu (\tau )}{d\tau }},`$ (5)   
$`V_0[z,s]`$ $`=`$ $`{\displaystyle \frac{g^2}{2}}s^2{\displaystyle _0^1}𝑑\tau {\displaystyle _0^1}𝑑\tau ^{}\mathrm{\Delta }(z(\tau )z(\tau ^{}),\mu ),`$ (6)   
$`V_{12}[z,\overline{z},s,\overline{s}]`$ $`=`$ $`{\displaystyle \frac{g^2}{2}}s\overline{s}{\displaystyle _0^1}𝑑\tau {\displaystyle _0^1}𝑑\overline{\tau }\mathrm{\Delta }(z(\tau )\overline{z}(\overline{\tau }),\mu ).`$ (7)   
Here the $`V_0[z,s]`$ term represents the self energy contribution, while the $`V_{12}[z,\overline{z},s,\overline{s}]`$ term represents the exchange interaction (Fig. 1). The interaction kernel $`\mathrm{\Delta }(x)`$ is defined by

$`\mathrm{\Delta }(x,\mu )`$ $`=`$ $`{\displaystyle \frac{d^4p}{(2\pi )^4}\frac{e^{ipx}}{p^2+\mu ^2}}={\displaystyle \frac{\mu }{4\pi ^2|x|}}K_1(\mu |x|).`$ (8)   
While we present expressions for the 2-body case, generalization to an arbitrary n-body system is trivial. The bound state spectrum can be determined from the spectral decomposition of the two body Green’s function

$$G(T)=\underset{n=0}{\overset{\mathrm{}}{}}c_ne^{m_nT},$$
(9)   
where T is defined as the average time between the initial and final states

$$T\frac{1}{2}(y_4+\overline{y}_4x_4\overline{x}_4).$$
(10)   
In the limit of large $`T`$, the ground state mass is given by

$$m_0=\underset{T\mathrm{}}{lim}\frac{d}{dT}ln[G(T)]=\frac{𝒟ZS^{}[Z]e^{S[Z]}}{𝒟Ze^{S[Z]}},$$
(11)   
While this result in principle is correct, the convergence to the asymptotic mass is slow due to the continuum contribution. The spectrum of the particle involves the mass pole and a cut beyond this pole representing the continuum contribution. Assuming that the constituents of the bound state are restricted to be at equal times in the initial and final states, the Green’s function can be written as a function of time $`T`$, total displacement $`𝐑`$, and relative coordinate $`𝐫`$, G(T,R,r).<sup>1</sup><sup>1</sup>1Dependence on the initial relative final coordinate $`𝐫_0`$ is implicit.

In order to eliminate the contribution of continuum states, introduce the Fourier transform,

$$\stackrel{~}{G}(T,𝐏,𝐩)=d^3𝐑d^3𝐫e^{i𝐏𝐑+i𝐩𝐫}G(T,𝐑,𝐫)$$
(12)   
where $`𝐏`$ is the CM momentum, and $`𝐩`$ is the relative momentum between particles. Setting both $`𝐏=𝐩=0`$ one has

$$\stackrel{~}{G}(T,0,0)=d^3𝐑d^3𝐫G(T,𝐑,𝐫),$$
(13)   
which eliminates the contribution of the continuum and projects out the s-wave state. While an integration over $`𝐫`$ is not necessary for the elimination of the continuum contribution, it is useful in eliminating the contribution of states with nonzero orbital angular momentum.

While the result for the Green’s function Eq. 3 is exact in the quenched approximation, due to its oscillatory behavior it is not appropriate for Monte-Carlo simulation. In Ref.  it was shown that one can perform a Wick rotation in variable $`s`$ to avoid these oscillations. In the limit $`g^20`$ the dominant contribution to the integral in Eq. (3) can be shown, by using the saddle point method, to come from

$$s=is_0=i\frac{T}{2m}.$$
(14)   
Since the large $`s`$ values do not contribute to the integral even without the interaction term, it is a good approximation to suppress the $`g^2`$ term at large s values. While this suppression is done it is important that the integrand is not modified in the region of dominant contribution $`sis_0`$. This can be achieved by scaling the $`s`$ variable, in the interaction term only, by

$$s\frac{s}{R(s,s_0)},$$
(15)   
where

$$R(s,s_0)1(sis_0)^2/\mathrm{\Gamma }^2.$$
(16)   
In the free case, $`(g^2=0)`$, the width $`W`$ of the region of dominant $`s`$ contribution goes as

$$W=\sqrt{\frac{T}{2m^3}}.$$
(17)   
Therefore, in the free case the dominant contribution to the $`s`$ integral comes from $`i(s_0W)<s<i(s_0+W)`$. This claim is supported by the Monte-Carlo simulation results. In Fig. 2 we present the results for $`s`$-distributions for two different coupling strengths for time $`T=40`$, and $`m=1`$ GeV. According to the estimate given above Eq. (17), for $`g=0`$, the dominant contribution to the $`s`$-distribution comes from the region $`15.53<s<24.47`$ which is in agreement with the result presented in Fig. 2. In order to ensure that the scaling given in Eq. (15) does not make a significant change in the region of dominant contribution, $`\mathrm{\Gamma }`$ should be chosen such that

$$\mathrm{\Gamma }W$$
(18)   
As one increases the coupling strength, the value of $`s_0`$ deviates from its free value. Therefore, in general, $`s_0`$ has to be defined self consistently by monitoring the peak of the $`s`$ distribution. In Figure 2 we display the s-distribution for two different coupling strengths. It is seen that as the coupling strength is increased the peak of the s-distribution moves towards higher s values. In general the peak of the distribution can be parameterized as

$$s_0=C\frac{T}{2m}.$$
(19)   
The dependence of $`C`$ on coupling strength $`g^2`$ is determined self consistently. In Fig. 3 $`C`$, which gives the location of the stationary point through Eq. (19), is plotted as a function of the coupling strength $`g^2`$. According to Fig. 3, it is not possible to find a self consistent stationary point beyond the critical coupling strength of $`g^2=31`$ GeV<sup>2</sup>. A similar critical behavior was also observed in Refs.  within the context of a variational approach.

The insensitivity of the dressed mass to the width $`\mathrm{\Gamma }`$ has been investigated  and found that a choice of $`\mathrm{\Gamma }^2=2W^2`$ was satisfactorily large. Results presented here employ the same value of $`\mathrm{\Gamma }`$.

Prescription given by Eqs. (15), and (16) enables one to perform a Wick rotation in variable $`s`$ to obtain a finite and non-oscillatory expression for the fully interacting two-body propagator:

$`G`$ $`=`$ $`{\displaystyle _0^{\mathrm{}}}𝑑s{\displaystyle _0^{\mathrm{}}}𝑑\overline{s}{\displaystyle (𝒟z)_{xy}(𝒟\overline{z})_{\overline{x}\overline{y}}}`$   
$`\times \mathrm{exp}\left[K[z,s]K[\overline{z},\overline{s}]+V_0[z,s_r]+V_0[\overline{z},\overline{s}_r]+2V_{12}[z,\overline{z},s_r,\overline{s}_r]\right],`$   
where

$$s_r\frac{s}{R(s,s_0)}.$$
(21)   
The path integral is discretized using

$$(𝒟)_{xy}(N/4\pi s)^{2N}\mathrm{\Pi }_{i=1}^{N1}d^4z_i,$$
(22)   
where the s-dependence is critical in obtaining the correct normalization. Note that the integration over the final coordinates is not included in this expression.

The discretized versions of kinetic and interaction terms are given by

$`K[z,s]`$ $``$ $`(m^2+iϵ)s{\displaystyle \frac{N}{4s}}{\displaystyle \underset{i=1}{\overset{N}{}}}(z_iz_{i1})^2,`$ (23)   
$`V_0[z,s]`$ $``$ $`{\displaystyle \frac{g^2s^2}{2N^2}}{\displaystyle \underset{i,j=1}{\overset{N}{}}}\mathrm{\Delta }({\displaystyle \frac{1}{2}}(z_i+z_{i1}z_jz_{j1}),\mu ),`$ (24)   
$`V_{12}[z,\overline{z},s,\overline{s}]`$ $``$ $`{\displaystyle \frac{g^2s\overline{s}}{2N^2}}{\displaystyle \underset{i,j=1}{\overset{N}{}}}\mathrm{\Delta }({\displaystyle \frac{1}{2}}(z_i+z_{i1}\overline{z}_j\overline{z}_{j1}),\mu ).`$ (25)   
We next address the regularization of the ultraviolet (short distance) singularities. The ultraviolet singularity in the kernel $`\mathrm{\Delta }(x,\mu )`$ Eq. (8) can be regularized using a Pauli-Villars regularization prescription. In order to do this one replaces the kernel

$$\mathrm{\Delta }(x,\mu )\mathrm{\Delta }(x,\mu )\mathrm{\Delta }(x,\alpha \mu ),$$
(26)   
where $`\alpha `$ is in principle a large constant. The ultraviolet singularity in the interaction is of the type

$$𝑑zz\mathrm{\Delta }(z,\mu ).$$
(27)   
At short distances the kernel $`\mathrm{\Delta }(z,\mu )`$ goes as $`1/z^2`$. Therefore, the self energy calculation involves a logarithmic type singularity. The Pauli-Villars regularization takes care of this singularity. The Pauli-Villars regularization is particularly convenient for Monte-Carlo simulations since it only involves a modification of the kernel. In order to achieve an efficient convergence in numerical simulations we use a rather small cut-off parameter $`\alpha =3`$. This choice leads to a less singular kernel. The value of $`\alpha `$ can be increased arbitrarily at the cost additional computational time. It should be noted that while we employ a Pauli-Villars regularization throughout this paper, in general, the bound state problem without self energies does not have any UV singularities. UV singularities are only associated with the self energies of the particles. After this brief summary of the formalism, in the next section we present some of the results and discuss the details of the algorithm.

## 3 Running the code: applications

The Monte-Carlo simulation starts by choosing an initial configuration for the trajectories of the particles. The choice of the initial trajectory is arbitrary (except at the end points). In Fig. 4 the evolution of the action as a function of number of updates for two different initial conditions is displayed. Starting with a random initial trajectory is analogous to a high temperature system and the termalization (reaching to the ground state) takes a long time (about 1000 updates). However if one starts by a classical free trajectory, the initial configuration is analogous to a frozen system without any fluctuations and the termalization results in an increase of the action. However asymptotically both results should converge as shown in Fig. 4. We usually start with an orderly system and disregard the first 1000 updates. Step sizes of the random walker in configuration space and in s-space should be chosen such that the average acceptance rate of Monte-Carlo updates are about % 50. In each run we typically make 500000 updates of trajectories. In order to reduce the statistical errors runs must be repeated (usually more then 10 times) using different random number seeds. The correlation function $`X(n)`$ of sampled configurations in each run is defined as

$$X(n)\frac{m(i)m(i+n)m^2}{m^2},$$
(28)   
where $`m(i)`$ is the mass measurement at the $`i`$’th update. The correlation function $`X(n)`$ measures how the information about a given configuration is lost as a function of the number of updates $`n`$. In Figure 5 we show the correlation function for the 1-body problem. According to Fig. 5 the number of updates necessary for the correlation between sampled trajectories to vanish is around $`n=1000`$. The number of configurations sampled in our simulations, which is around 500000, is well above the correlation length 1000, thereby insuring that uncorrelated trajectories are sampled during the Monte-Carlo integration.

The time $`T`$ required to reach to the asymptotic limit increases as $`g^2`$ increases. For $`g^2=25`$ the asymptotic value of the mass, given by Eq. (11), is obtained around $`mT=40`$. In particular for the result shown in Fig. 6, as $`g^2`$ ranged from $`g^2=0`$ GeV<sup>2</sup> to $`g^2=31`$ GeV<sup>2</sup> the asymptotic time values used were increased from $`mT=35`$ to $`mT=45`$. As time $`T`$ is increased the number of steps $`N`$ should also be proportionally increased so that the step size of the particle trajectory remains the same. As one increases the coupling strength $`g^2`$, trajectories of the particles deviate from the classical trajectory to a greater degree. This increase in fluctuations requires that one uses a higher $`N`$ value. The number of steps $`N`$ particles take between the initial and final coordinates is typically chosen to be $`35<N<45`$.

In Fig. 6 the $`g^2`$ dependence of the one-body dressed mass is presented. 1-body masses presented here are lower than those in . This is due to the fact that in that work $`C`$ was assumed to be approximately 1. However it was subsequently realized that $`C`$ deviates from 1 significantly (see Fig. 3) as coupling strength is increased, and a self consistent determination of $`C`$ is important to pin down the critical point. According to Fig. 3 there is no self consistent stationary point Eq. (19) beyond the critical coupling strength $`g^2=31`$ GeV<sup>2</sup>, and the 1-body dressed mass becomes unstable for $`g^2>31`$ GeV<sup>2</sup> (Fig. 6).

Next two applications involve the two and three-body bound states of equal mass particles. While the algorithm provided is capable to take into account all self energy and vertex dressing corrections to the bound state, in the bound state applications provided here the self energy contributions of particles were not taken into account. This is controlled by a switch in the input file (see Table 4.1 for input options). In Fig. 7 the $`g^2`$ dependence of the two-body bound state mass is presented. Beyond the critical coupling strength of $`g^2=100`$ GeV<sup>2</sup> the 2-body mass becomes unstable. During the Monte-Carlo simulation the radial distributions of particles are stored in histograms. In general for an n-body bound state there are $`n(n1)/2`$ relative distances. Since the particles are assumed to have equal time coordinates in the final state, essentially the relative distance is equal to the spatial distance. Let $`r_{ij}=|r_ir_j|`$ be the distance between particles $`i`$ and $`j`$. In this case the histogram $`P_{ij}(r_{ij})`$ stores the number of final state configurations sampled in the interval $`(r_{ij}\delta /2,r_{ij}+\delta /2)`$, where $`\delta `$ is the bin size. The range of $`r_{ij}`$ values is specified in the program using

$$0<r_{ij}<\frac{20}{m_i+m_j}\mathrm{Fermi}.$$
(29)   
The choice of range is arbitrary as long as a reasonably smooth histogram is produced. The number of bins in each histogram is chosen to be 100. Therefore for two particles of mass 1 GeV, the maximum radius to be histogrammed is 10 Fermi; and the size of each bin is given by $`\delta =0.1`$ Fermi. In Fig. 8 we present the projection of the radial probability distribution onto a two-dimensional surface. The two-body probability distribution, presented in Fig. 8, shows the result for $`g^2=25`$ GeV<sup>2</sup>. Notice that the probability distribution vanishes at the origin. This is due to the phase space factor $`4\pi r^2`$. In order to be able to compare the probability distribution with a wave function the phase space factor should be factored out. In general the radial wave function amplitude $`\mathrm{\Psi }(r)`$ is given by:

$$|\mathrm{\Psi }(r)|=\sqrt{\frac{P(r)}{4\pi r^2}}$$
(30)   
The normalization of the probability distribution histogram presented in Fig. 9 is arbitrarily fixed such that the maximum entry is equal to 1.

In generating the surface plot (Fig. 8), it is assumed that one of the particles is fixed at the origin. The amplitude of the surface gives the probability of finding the second particle at a distance $`r`$. While the surface plot for a two-body bound state is not necessary for a visual understanding of the bound state structure, the three-body bound state demands a three dimensional plot. In the three-body case there are three relative coordinates. In order to be able to plot the three-body probability distribution we fix the location of two of the particles along the y axis, and calculate the probability of finding the third particle in an arbitrary distance from two fixed particles. In Fig. 10 we represent the probability distribution of the third particle for a given fixed configuration of the first and second particles. Assume that the fixed particles are particle 1 and particle 2. The probability distribution of the third particle is given by

$$P_3(|𝐫_\mathrm{𝟑}|)P_{13}(|𝐫_\mathrm{𝟑}𝐫_\mathrm{𝟏}|)P_{23}(|𝐫_\mathrm{𝟑}𝐫_\mathrm{𝟐}|)$$
(31)   
Probability distribution $`P_3(|𝐫_\mathrm{𝟑}|)`$ includes the phase space contribution. Therefore $`P_3(|𝐫_\mathrm{𝟑}|)`$ represents the probability of finding the third particle in a ring shown in Fig. 11. For example in the first plot of Fig. 10 two fixed particles are very close to each other such that the third particle sees them as a point particle. Therefore the probability distribution of the first plot of Fig. 10 is very similar to the two-body distribution given in Fig. 8. However as the fixed particles are separated from each other the third particle starts having a nonzero probability of being in between the two fixed particles (second and third plots of Fig. 10). Eventually when the two fixed particles are kept away from each other the third particle has a nonzero probability distribution only at the origin (the last plot shown in Fig. 10). In Fig. 12 the two-body distribution function $`P_{12}`$ in a three-body system is shown. All of these plots were produced with three equal particles of mass $`1`$ GeV, and the coupling strength of $`g^2=64`$ GeV<sup>2</sup>.

## 4 Program details

In this section we discuss the details of the program. We start by providing the tables of input parameters 4.1 and arrays 4.2 and then summarize the components of the program.

### 4.1 Description of Input parameters

SIG = 1 GeV Arbitrary momentum scale  IPAR(I) I=1,2,3, Particles present (0=n, 1=y)  IEXCH(I,I) I=1,2,3, Self Interactions (0=n, 1=y)  IEXCH(I,J) $`(I,J)=(1,2),(1,3),(2,3)`$ Exch. Int. (0=n, 1=y)  QM(I) Mass of particle I  G The coupling strength  XMU The exchange mass $`\mu `$  ALPHA The Pauli-Villars mass ratio: $`m_{pv}=\alpha `$ $`\mu `$  BETA = 1 The width of $`R(s,s_0)1\beta (siC\widehat{s}_0)^2/T^\gamma `$  GAMMA = 1 The T dependence of the width.  C The peak location of the s distribution s=$`C\widehat{s}_0`$.   C is determined self-consistently.  N Number of steps along the trajectory  NSMPL $`\mathrm{\#}`$ of sampled trajectories in MC integration  NVOID No. of uncounted samples for initial termalization  ZSTEP Max. step size of the random walker in coord. space  SSTEP Max. step size of the random walker in $`s`$ space  EPSA The short distance(UV) cut-off for preventing   explicit occurrence of 1/0.0 type of singularity.   This is not a regulator however. The UV   regularization is done using Pauli-Villars subtraction  IDUM The seed of the random number generator  Z(I,0,J) Initial coordinates of particle I=1,2,3  Z(I,N,J) Final coordinates of particle I=1,2,3  IWRTA,IWRTM Action, mass to be stored in file ? (0=n, 1=y)  INTOUT Integrate over final coordinates ?   

### 4.2 Arrays

QM(3) Particle masses    IPAR(3) Determines which particles are present.  Z(3,0:400,4) The trajectories of particles  ZNEW(3,0:400,4) The updated trajectories of particles  SMAX(3) The maximum value of s value to be histogrammed.   (not an integration cut-off)  RMAX(3,3) The maximum relative distance between particles   to be histogrammed. (not an integration cut-off)  SHIST(3,200) Histogram of $`s`$ values  RHIST(3,3,200) Histogram of relative distances between particles  SUMK(3) Kinetic energies of particles  SUMV(3,3) Potential energies of particles: SUMV(I,I): self   energy of the I’th particle, SUMV(I,J): the exchange   energy between I’th and J’th particles.  SUMKN(3) Updated kinetic energies  SUMVN(3,3) Updated potential energies   

### 4.3 Main Program

The main program starts by calling the INPUT subroutine. This subroutine reads input parameters described in Table 4.1. Next, subroutine XINTGR is called to perform the Monte-Carlo simulation. In the following the role of each subroutine and function is explained.

### 4.4 Subroutine INPUT

In this subroutine 12 input parameters are read. In addition to reading these parameters, histograms for the radial and s distributions, and initial trajectories of the particles are initialized. Initialization of trajectories is done using the classical free trajectories of particles.

### 4.5 Subroutine XINTGR

This subroutine controls basic steps of the Monte-Carlo sampling process. First, the kinetic and the potential sums corresponding to the initial trajectories of particles are calculated by calling the XSUMS subroutine. Next step is the termalization process. In order to reach termalization the configuration of trajectories are updated NVOID times by calling the UPDATE subroutine. First NVOID updates are not used in the actual calculation of the bound state mass or the probability distribution. After the termalization the sampling is done for NSMPL times. The results for the $`s`$ values and the relative distances of particles are histogrammed. Finally the bound state mass is calculated.

### 4.6 Subroutine UPDATE

This subroutine is responsible for updating the current configuration of trajectories. Each coordinate and $`s`$ parameter is updated once. Samplings are done according to distribution

$$e^{S[z]}$$
(32)   
distribution. If the ratio $`r`$,

$$r=e^{(S[z^{}]S[z])},$$
(33)   
is larger than 1 updates are always accepted. When $`r<1`$ the update is accepted with a probability of $`r`$.

### 4.7 Subroutine XSUMS

This subroutine calculates kinetic and potential sums before the configuration updates start. Since calculation of the kinetic and potential sums is a costly operation, after every update we calculate the shift in the sums.

### 4.8 Subroutine ACTION

This subroutine calculates the action using known values of kinetic and potential sums. The action is stored in a file to monitor the termalization. ( See Fig. 4 ).

### 4.9 Subroutine UPDSUM

Since calculation of the kinetic and potential sums is a costly operation, after every update we calculate the shift in the sums. Using this shift sums are updated.

### 4.10 Function XDERIV

This subroutine calculates the derivative operator $`S^{}[Z]`$ Eq. (11) which gives the mass of the bound state.

### 4.11 Function DELTA, DELTAP

These subroutine calculate, respectively, the interaction kernel Eq. (8) and its time derivative.

### 4.12 Function DLARAN

This is a random number generator obtained from LAPACK package at Netlib.  DLARAN returns a random real number from a uniform (0,1) distribution. In the actual implementation of the program we have used a Numerical Recipes random number generator (ran2). However for copyright reasons here we provide this alternative random number generator.

### 4.13 Functions BESK0, BESK1, KZEONE

BESK0 and BESK1 are modified bessel functions $`K_0(x)`$ and $`K_1(x)`$ which are needed in calculation of the interaction kernel and its derivative. BESKO and BESK1 call KZEONE subroutine . KZEONE subroutine returns real and imaginary parts of $`e^xK_0(x)`$ and $`e^xK_1(x)`$. The KZEONE subroutine is considerably slower than the Numerical Recipes subroutines for modified bessel functions. For copyright reasons we provide KZEONE rather than the Numerical Recipes subroutines for the modified bessel functions.

## 5 Conclusions

In this paper, using the Feynman-Schwinger representation, we have presented an algorithm that calculates 1,2,3 body masses and distribution probabilities in the quenched approximation for $`\chi ^2\varphi `$ theory in 3+1 dimension. The FSR approach provides an efficient method to calculate nonperturbative propagators. In this work we have presented results of applications to the 1, 2 and 3 body states. A detailed comparison of quenched bound state results of the FSR approach with the bound state equation predictions is under study and will be presented in a separate physics article . It is hoped that, through comparison, this simple and rigorous nonperturbative method will enhance our understanding of various nonperturbative bound state models and approximations.

Acknowledgements

We are grateful to F. Gross, and J. Tjon for discussions. The support of the DOE through grant No. DE-FG02-97ER41032 is gratefully acknowledged. The Thomas Jefferson National Accelerator Facility is gratefully acknowledged for warm hospitality and for providing computer resources.

## 6 Test run INPUT

| 1.0 | SIG (GeV) |
| --- | --- |
| 1 0 0 | IPAR(1), (2), (3) |
| 1 0 0 | IEXCH(1,1)(2,2)(3,3) |
| 0 0 0 | IEXCH(1,2)(1,3)(2,3) |
| 1.0 | QM(1) |
| 1.0 | QM(2) |
| 1.0 | QM(3) |
| 2.0 | G |
| 0.15 | XMU |
| 3.0 | ALPHA |
| 0.5 | BETA |
| 1.0 | GAMMA |
| 1.0 | C |
| 35 | N |
| 500000 | NSMPL |
| 5000 | NVOID |
| 2.1 | ZSTEP |
| 4.5 | SSTEP |
| 1.0D-4 | EPSA |
| 1 | IDUM |
| 0. 0. 0. 0. | Initial coordinates |
| 1. 1. 0. 0. |  |
| -1. 1. 0. 0. |  |
| 0. 0. 0. 35. | Final coordinates |
| 1. 1. 0. 35. |  |
| -1. 1. 0. 35. |  |
| 0 0 | IWRTA, IWRTM |
| 1 | INTOUT |

## 7 OUTPUT

| SIG | 1.0 GeV |
| --- | --- |
| Particles present | 1 0 0 |
| Self interactions | 1 0 0 |
| Exchange ” 1-2 | 0 |
| ” ” 1-3 | 0 |
| ” ” 2-3 | 0 |
| QM(1) | 1. GeV |
| QM(2) | 1. GeV |
| QM(3) | 1. GeV |
| G | 2. GeV |
| XMU | 0.15 GeV |
| ALPHA | 3. |
| BETA | 0.5 |
| GAMMA | 1. |
| C | 1. |
| N | 35 |
| NSMPL | 500000 |
| NVOID | 5000 |
| ZSTEP | 2.1 1/GeV |
| SSTEP | 4.5 1/GeV<sup>2</sup> |
| EPSA | 0.10E-03 GeV<sup>-1</sup> |
| IDUM | 1 |
| Initial coordinates: | 0 0 0 0 GeV<sup>-1</sup> |
|  | 1 1 0 0 GeV<sup>-1</sup> |
|  | -1 1 0 0 GeV<sup>-1</sup> |
| Final coordinates: | 0 0 0 35 GeV<sup>-1</sup> |
|  | 1 1 0 35 GeV<sup>-1</sup> |
|  | -1 1 0 35 GeV<sup>-1</sup> |
| INTOUT | 1 |
| z update % | 50.33 |
| s update % | 53.63 |
| Bound state mass | 0.97$`\pm `$0.001 GeV |