# Monte-Carlo Study of Bound States in a Few-Nucleon System – Method of Continued Fractions –

## I Introduction

The exact treatment of nuclear many-body systems is an important benchmark for investigating the dynamics of multi-nucleon systems. It is a testing ground of potential models and equations of motion, such as those involving multi-nucleon interaction, relativistic dynamics, and possible modification of nucleon properties inside nuclei.

At present, various methods to solve many-body problems have been developed. One of the standard methods is the Faddeev or Faddeev-Yakubovsky (FY) equation. The FY equation has been applied to the investigation of the ground state of <sup>4</sup>He using a realistic nuclear potential. A common problem in attempts to extend the ‘exact’ study of many-body systems is the rapid increase of the number of channels and dimensions as the number of particles increases. One way to overcome this problem is to diagonalize the Hamiltonian in a subspace spanned by appropriate bases chosen using a stochastic method. The other way is the direct calculation of the matrix elements using a Monte-Carlo method. In the Monte-Carlo method, one can perform an effective multi-dimensional integration by taking the importance sampling from the essential part of the phase space.

Extensive works have been carried out by using a variational Monte-Carlo method (VMC) and a Green’s function Monte-Carlo method (GFMC). VMC and GFMC have been applied to the study of the ground state properties of nuclei with mass number $`A7`$ using realistic nuclear potential. The bound state energies and also the electro-magnetic properties of few-nucleon systems have been studied extensively. In VMC, we can obtain the upper bound of the ground state energy for which the matrix elements are evaluated by Monte-Carlo integration. Simultaneously, we obtain an analytic form of the trial function which describes the ground state quite well. In GFMC, we start the calculation with a trial function obtained by VMC, and the imaginary-time evolution $`e^{\tau H}`$ extracts the improved ground state wave function from the trial function. Here, the imaginary-time evolution is carried out by using the quantum Monte-Carlo random walk method. Due to the $`e^{\tau H}`$ projection, we only obtain the properties for the lowest state of the system with a given quantum number. Excited states in the diffusion Monte-Carlo method can be investigated using a diagonalization procedure in the large space shell-model calculation with an auxiliary-field Monte-Carlo approach and in one- and two-dimensional examples of GFMC. For the purpose of conducting systematic studies in few-nucleon systems, it is necessary to investigate excited states as well as the ground state. Therefore we should develop a method for this purpose in the context of few-nucleon physics.

In this paper, we propose a new Monte-Carlo method to study the bound state problems in few-nucleon systems. Here, our purpose is to find an effective method to extract eigenfunctions without referring to the imaginary-time evolution method, but still using the Monte-Carlo random walk. Therefore, the method is not restricted to the ground state problem. We use the method of continued fractions (MCF) as a central formalism in our new calculation method. MCF was originally proposed by J. Horáček and T. Sasakawa  as a fast calculational method to evaluate T-matrix elements and was applied to the three-nucleon system. The wave functions which appear in the MCF iteration are generated so as to be orthogonal to the lower order ones. Then the matrix elements calculated by those wave functions rapidly converge. Moreover, this method is simple and allows us to use any function as a starting trial wave function and to adopt a nuclear potential including a many-body potential.

The direct application of MCF with Monte-Carlo integration, however, causes a few serious problems: First, the original form of the MCF iteration formula  is not suitable to apply a Monte-Carlo random walk, since the modified Green’s functions used in each iteration are written in a separable form. We find with a slight modification that the iteration formula can be written in a form similar to the diffusion equation. There is, however, a crucial difference between them. In our new iteration formula, a higher order MCF wave function is generated from a lower order one. Therefore, the equation cannot be a simple iteration of the same function as in the diffusion equation. Second, in order to guarantee the orthogonalization of each MCF wave function, we must deal the difference of two MCF wave functions. This causes a sort of the “sign problem”. Here, we propose a correct method of counting the weight in the random walk procedure for an equation of a non-iterative type. As a result, we find out that our method is closely related to “the multiple cancellation” method. With our cancellation method, we can take care of the sign problem in our formula.

The validity of our method of counting weight and the convergence of the MCF are extensively examined for a simple bound state problem in one-dimension including the excited states. As a more realistic example of a multi-dimensional problem, we study $`A=4`$ nuclei with the Volkov potential. Our method is successful in the treatment both examples.

We briefly review the method of continued fractions in §2. Our random walk method in a non-iterative equation and multiple cancellation mechanism are discussed in §3. We present our results in §4, and a summary is given in §5.

## II Method of continued fractions

The method of continued fractions (MCF) has been developed as an efficient method to study scattering problems in atomic and nuclear physics. We briefly review this method, here.

The Schrödinger equation of nuclear bound state can be written as

$$|\varphi =G_0(E)V|\varphi \text{ ,}$$
(1)   
where $`G_0(E)`$ is a free Green’s function,

$$G_0(E)=\frac{1}{EH_0}\text{ ,}$$
(2)   
and $`H_0`$ and $`V`$ are the kinetic energy of the nucleons and the nuclear potential, respectively.

In this paper, we use the modified Green’s function (MCFG) approach, since it is easy to implement in Monte-Carlo evaluations. We introduce the arbitrary functions $`|F_0`$ and $`|f`$ which are not orthogonal to the eigenfunction $`|\varphi `$. A new function $`|F_1`$ is defined as

$$|F_1G_0(E)V|F_0\text{ ,}$$
(3)   
and a modified Green’s function $`G_1(E)`$ as

$$G_1(E)G_0(E)\frac{|F_1f|}{f|V|F_0}\text{ .}$$
(4)   
$`G_1(E)`$ is chosen so as to operate only on the space which is orthogonal to $`V|F_0`$ as

$$G_1(E)V|F_0=0\text{ .}$$
(5)   
We eliminate $`G_0(E)`$ from the Schrödinger equation (1) using (4):

$$|\varphi =\left(G_1(E)+\frac{|F_1f|}{f|V|F_0}\right)V|\varphi \text{ .}$$
(6)   
Then we define $`|\varphi _1`$ as

$$|\varphi _1\frac{1}{1G_1(E)V}|F_1\text{ ,}$$
(7)   
and (1) can be written as

$$|\varphi =|\varphi _1\frac{f|V|\varphi }{f|V|F_0}\text{ .}$$
(8)   
Multiplying (8) by $`f|V`$ from the left, we obtain the following condition to determine the bound state energy $`E`$ and the wave function $`|\varphi `$:

$$\mathrm{\Omega }(E)f|V|F_0f|V|\varphi _1=0\text{ .}$$
(9)   
The bound state condition can be written in terms of a continued fraction as

$$|\varphi _i\frac{1}{1G_i(E)V}|F_i\text{ ,}$$
(10)   
$$|F_{i+1}=G_i(E)V|F_i\text{ ,}$$
(11)   
and

$$G_i(E)=G_{i1}(E)\frac{|F_if|}{f|V|F_{i1}}\text{ ,}$$
(12)   
where $`|\varphi _i`$ is the $`i`$-th order wave functions. The $`i`$-th Green’s function $`G_i(E)`$ satisfies the orthogonality relation

$$G_i(E)V|F_j=0\text{ .}(i=j+1)$$
(13)   
We rewrite (10) using (12), multiply $`f|V`$ from the left, and obtain the following iteration formula for $`|\varphi _i`$:

$$f|V|\varphi _i=\frac{f|V|F_i^2}{f|V|F_if|V|\varphi _{i+1}}\text{ .}$$
(14)   
As a result, the bound state condition (9) can be written in the form of a continued fraction:

$$\mathrm{\Omega }(E)=x_0\frac{x_1^2}{x_1{\displaystyle \frac{x_2^2}{x_2{\displaystyle \frac{x_3^2}{x_3_{\mathrm{}}}}}}}\text{ .}$$
(15)   
Here, the matrix elements $`x_i`$ are defined as

$$x_if|V|F_i\text{ .}$$
(16)   
In the MCF calculation, we only need to generate $`|F_i`$ and $`G_i(E)`$ iteratively, following the iteration procedure (11) and (12). In this procedure, $`G_i(E)`$ projects out the states which are orthogonal to the lower order wave functions, as indicated in (13). Therefore, $`x_i`$ is expected to become smaller for larger $`i`$, and rapidly converges into zero. The advantage of this continued fraction series based on Ref.  for the Monte-Carlo calculation is that the $`x_i`$ are simply given by the matrix element of the potential between $`|F_i`$ and the analytically known function $`|f`$, in contrast to those given in Ref. , where we need the matrix elements between $`|F_i`$.

The bound state energy is determined by solving the equation

$$\mathrm{\Omega }(E)=0$$
(17)   
numerically. In actual practice, we perform a parameter search for $`E`$. Since we do not adopt any projections onto the ground state, we can also apply MCF to study the excited states without any restriction such as the $`e^{\tau H}`$ projection.

Once the energy of the bound state has been determined, the relevant eigenfunction is given as

$$|\varphi =𝒩\left[|F_1+y_1|F_2+y_1y_2|F_3+y_1y_2y_3|F_4+\mathrm{}\right]\text{ ,}$$
(18)   
where $`y_i`$ is defined as

$$y_i\frac{x_i}{x_if|V|\varphi _{i+1}}\text{ ,}$$
(19)   
and $`𝒩`$ is a normalization constant.

## III Continued fractions and Monte-Carlo method

The original form of the MCF iteration formula (11) is not suitable for Monte-Carlo evaluations. The difficulty is caused by the Green’s functions $`G_i(E)`$. Since the $`G_i(E)`$ obey the iteration formula (12) and are of a separable form,

$$G_i(E)=G_{i1}(E)\frac{|F_if|}{x_{i1}}\text{ ,}$$
(20)   
it does not seem practical to generate the random numbers according to the probability of $`G_i(E)`$. In order to resolve these problems, we use the orthogonality relation (13) and modify the iteration formula (11) as follows:

$$|F_{i+1}=G_0(E)V\left[|F_i|F_{i1}\frac{x_i}{x_{i1}}\right]\text{ .}$$
(21)   
The advantage of this formula is that a higher order wave function $`|F_{i+1}`$ is always generated by the free Green’s function $`G_0(E)`$, which has a simple analytical form. Therefore, we can use the Monte-Carlo random walk technique in the MCF iteration procedure.

However, it should be noted that the iteration formula (21) is not of the same form as the original Schrödinger equation. The Schrödinger equation,

$$|\varphi =G_0(E)V|\varphi \text{ ,}$$
(22)   
can be regarded as a diffusion equation. $`|\varphi `$ will converge to a ground state by iterating (22) with a Monte-Carlo random walk. In contrast, (21) can be written as

$$|\psi _1=G_0(E)V|\psi _0\text{ ,}$$
(23)   
where $`|\psi _1`$ and $`|\psi _0`$ represent $`|F_{i+1}`$ and $`[|F_i|F_{i1}x_i/x_{i1}]`$, respectively. In our iteration formula (23), the wave function $`|\psi _1`$ has different functional form from $`|\psi _0`$ and is generated from the two terms $`|F_i`$ and $`|F_{i1}`$. Therefore, we cannot naively apply the random walk technique of the diffusion equation. Otherwise, a destructive interference of $`|F_i`$ and $`|F_{i1}`$ in (21) may cause the “sign problem”.

In this section, we extend the random walk technique for (23) and propose a new method of counting weight. The destructive interference of two terms in (21) can be taken care of by our proper method of counting weight.

### A Random walk in MCF

We study an algorithm to apply the random walk method with the free Green’s function $`G_0(R,R^{})`$ to our following iteration formula of MCF:

$$\psi _1(R)=G_0(R,R^{})W(R^{})\psi _0(R^{})𝑑R^{}\text{ .}$$
(24)   
Here we re-scale the variables so that the free Green’s function $`G_0(R,R^{})`$ is normalized as

$$\left[_R^2+1\right]G_0(R,R^{})=\delta (RR^{})\text{ ,}$$
(25)   
where $`R`$ denotes an $`M`$-dimensional vector, i.e. $`R=\{z_1,\mathrm{},z_M\}`$ defined in the standard manner given in Ref. . The potential $`W(R)`$ is defined as

$$W(R)V(R)/B\text{ ,}$$
(26)   
where $`B`$ is the binding energy.

For simplicity, we assume that the wave functions $`\psi _1(R)`$ and $`\psi _0(R)`$ and the potential $`W(R)`$ are positive definite. Note that the positive sign of $`W(R)`$ implies that this is an attractive potential. We assume $`\psi _1(R)`$ and $`W(R)\psi _0(R)`$ are represented by the sum of points $`\{R_{k,i}\}`$ distributed with the probability $`\psi _k(R)`$ and the weights $`\{W_{k,i}\}`$, where $`k=0,1`$. The question is how $`\{R_{1,j}\}`$ and $`\{W_{1,j}\}`$ should be generated from given sets of $`\{R_{0,j}\}`$ and $`\{W_{0,j}\}`$. The relation can be obtained by studying the following overlap integral of $`\psi _1(R)`$ with an arbitrary function $`F(R)`$:

$$=F(R)\psi _1(R)𝑑R\text{ .}$$
(27)   
$``$ is given by using $`\{R_{1,j}\}`$ and $`\{W_{1,j}\}`$ as

$$=\underset{j}{\overset{N}{}}W_{1,j}F(R_{1,j})\text{ ,}$$
(28)   
while the same integral is expressed by $`\{R_{0,j}\}`$ and $`\{W_{0,j}\}`$ using (24) as

$$=F(R)\underset{i}{\overset{N}{}}W_{0,i}G_0(R,R_{0,i})dR\text{ ,}$$
(29)   
where $`N`$ is the number of sampling points.

We multiply the integrand of (29) by unity in the form $`G_0(R,X)/G_0(R,X)`$ and obtain

$$=A(R,X)G_0(R,X)𝑑R\text{ ,}$$
(30)   
where $`A(R,X)`$ is defined as

$$A(R,X)F(R)\underset{i}{\overset{N}{}}W_{0,i}\frac{G_0(R,R_{0,i})}{G_0(R,X)}\text{ .}$$
(31)   
Here, $`X`$ can be an arbitrary vector in $`M`$-dimensional space. We impose hyper spherical boundary conditions on $`G_0(R,X)`$ for the bound state problem,

$$G_0(R,X)=0\text{for}|R|,|X|\mathrm{}\text{.}$$
(32)   
$`G_0(R,X)`$ becomes a function of distance $`||`$, where $``$ is defined as $`RX`$. Since $`G_0(R,X)`$ is normalized as

$$G_0(R,X)𝑑R=G_0(||)𝑑=1\text{ ,}$$
(33)   
$`G_0(||)`$ can be regarded as a probability function.

We change the integral variable from $`R`$ to $``$ in (30) and generate a set of random numbers $`\{_j\}`$ with probability $`G_0(||)`$. The algorithm to generate $`\{_j\}`$ in $`M`$ dimensional space is given in Ref. . Equation (30) is rewritten as follows:

$$=\underset{j}{\overset{N}{}}F(X+_j)\frac{_iW_{0,i}G_0(X+_j,R_{0,i})}{NG_0(X+_j,X)}\text{ .}$$
(34)   
Since (34) holds for any values of $`X`$, we choose $`X`$ as $`R_{0,j}`$ for the $`j`$-th random number $`_j`$. Comparing (34) with (28), we can define the new points $`R_{1,j}`$ and the weight $`W_{1,j}`$ as

$$R_{1,j}R_{0,j}+_j\text{ ,}$$
(35)   
and

$$W_{1,j}\frac{_iW_{0,i}G_0(R_{0,j}+_j,R_{0,i})}{NG_0(R_{0,j}+_j,R_{0,j})}\text{.}$$
(36)   
As a result, (35) and (36) show how the distribution of points $`\{R_{1,j}\}`$ and the weights $`\{W_{1,j}\}`$ of $`\psi _1(R)`$ should be generated from those of $`W(R)\psi _0(R)`$. A new point $`R_{1,j}`$ is generated by a random walk from $`R_{0,j}`$ similar to the diffusion equation. However, a new weight $`W_{1,i}`$ is generated by taking into account all of the old points $`\{R_{0,i}\}`$ and weights $`\{W_{0,i}\}`$ through the free Green’s function. In the naive random walk method, $`W_{1,j}`$ is simply given by $`W_{0,i}`$, which projects into the ground state wave function. For a non-diffusion type integral formula like (24), weights of new points should be generated by using this method of counting weight.

### B Multiple cancellation

In general, $`\psi _0(R)`$ and $`W(R)`$ are not positive definite, and we must include their sign information in our formula. This can be done by adding the sign $`s_{0,j}`$ to the weight $`W_{0,j}`$ as

$$W_{0,j}s_{0,j}W_{0,j}\text{ .}$$
(37)   
Similarly, adding the sign $`s_{1,j}`$ to $`W_{1,j}`$, (36) is replaced by

$$s_{1,j}W_{1,j}=\frac{_is_{0,i}W_{0,i}G_0(R_{0,j}+_j,R_{0,i})}{NG_0(R_{0,j}+_j,R_{0,j})}\text{ .}$$
(38)   
Equation (38) is quite similar to the “multiple cancellation” method proposed by J. B. Anderson et al.. The multiple cancellation method with a random walk has been applied to atomic physics, where it has provided very accurate results. In the conventional random walk technique, points with positive and negative signs walk independently. Then the weight of each point grows, and a destructive interference of these wave functions induces large numerical errors. On the other hand, our formula (38) shows that weight and sign of new points are generated by taking into account all the weights and signs of old points. Therefore, the points with positive and negative weights are canceled out in the course of counting the new weight. Equation (38) plays an essential role in the MCF calculation, and its importance will be demonstrated in the one-dimensional example given in the next section.

## IV Applications to bound state problems

The procedure of our method is summarized as follows:

1. Assume an appropriate form of the trial function $`|F_0`$.
2. Distribute the points by the Metropolis method according to the functional form of $`|F_0`$.
3. Calculate the matrix element $`x_i`$ of the $`i`$-th iteration step.
4. Generate new points of $`|F_{i+1}`$ by a random walk according to the procedure described in the previous section.
5. Iterate steps 3 and 4 until the $`x_i`$ and $`\mathrm{\Omega }(E)`$ converge.

The binding energy $`B(=E)`$ is determined by searching for an energy that satisfies the condition $`\mathrm{\Omega }(E)=0`$.

### A One-dimensional model

As a simplest example, we study a one-dimensional potential model. In this one-dimensional problem, we can study the validity of our method by comparing with a conventional numerical method. We take the mass of the particle as $`m=940`$ MeV and a simple central potential $`V(x)`$ of exponential form with a one pion range as follows:

$$V(x)=V_0\mathrm{exp}(a|x|)$$
(39)   
with

$$V_0=69.3(\text{MeV})\text{ , }a=0.711(\text{fm}^1)\text{ .}$$
(40)   
The strength of the potential is chosen so as to produce three bound states at $``$38.7, $``$10.6 and $``$2.01 MeV. The binding energy of the second excited state is small, and is compared to the deuteron binding energy. Since this model is one-dimensional and does not contain any other degrees of freedom, such as spin or isospin, the wave function of the bound states can be classified by their parity. The ground and second excited states are even parity states, and the first excited state is an odd parity state. The validity of our method for the second excited state is particularly interesting, since it has the same symmetry as the ground state.

Introducing a dimensionless variable $`z`$ defined as

$$z=\sqrt{2mB}x\text{ ,}$$
(41)   
the free Green’s function is given by

$$G_0(z,z^{})=\frac{1}{2}\mathrm{exp}(|zz^{}|)\text{.}$$
(42)   
The trial function is taken to be a Gaussian form as follows:

$`|F_0=|f`$ $`=`$ $`{\displaystyle \frac{1}{\sqrt{b\sqrt{\pi }}}}\mathrm{exp}({\displaystyle \frac{x^2}{2b^2}})\text{ for even parity,}`$ (43)   
$`=`$ $`\sqrt{{\displaystyle \frac{2}{b\sqrt{\pi }}}}{\displaystyle \frac{x}{b}}\mathrm{exp}({\displaystyle \frac{x^2}{2b^2}})\text{ for odd parity.}`$ (44)   
We assume $`b=3`$ (fm) for the ground and first excited states and $`b=10`$ (fm) for the second excited state.

In the Monte-Carlo calculation, the numbers of sampling points were 20,000, 70,000 and 100,000 for the ground, first, and second excited states, respectively. We needed a larger number of points for the excited states than for the ground state in order to obtain reasonable accuracy. The initial distribution of points with the probability $`||F_0|`$ was generated by the Metropolis method. A typical number for the thermalization step was $`100`$ in this problem.

The fast convergence of the MCF iteration is related to the orthogonality relation in (13), which is implemented in our iteration formula in (21). As an example, we consider the second order MCF wave function $`|F_2`$,

$$|F_2=G_0(E)V\left[|F_1|F_0\frac{x_1}{x_0}\right]\text{ .}$$
(45)   
$`|F_2`$ is generated by two terms, $`|F_1`$ and $`|F_0x_1/x_0`$, which will cancel each other. The higher order $`|F_i`$ decreases quickly with this cancellation mechanism. Therefore, accurate treatment of this cancellation is important. If we generate $`|F_2`$ from two terms independently, the resultant distribution of the points for $`|F_2`$ is given by the large overlap of the points with positive and negative signs shown in Fig. 1 at the energy of the second excited state. $`|F_2`$ generated from these points will include large errors, which is just our “sign problem”. Because of this sign problem, it is difficult to obtain convergence for the higher order MCF wave functions in the standard random walk method. On the other hand, using our method of counting weight in (38), the cancellation between two terms is performed analytically, and the distribution of the positive points for $`|F_2`$ is clearly separated from the negative points, as shown in Fig. 2. Our method of counting weight plays an important role in MCF iteration with the Monte-Carlo method.

MCF wave functions at each iteration step are shown in Fig. 3. They decrease rapidly as the iteration steps increase, even for the second excited state. This agrees well with the results obtained using the standard method of numerical integration.

Consequently, the matrix elements $`x_i`$ defined in (16) and $`\mathrm{\Omega }(E)`$ also converge rapidly as a function of the iteration step as shown in Fig. 4 and 5, respectively. A few steps of iteration is enough to obtain reasonable convergence for the ground and first excited states. Slower convergence for the second excited state is due to our simple Gaussian trial function, even though one expects nodes in the wave function of the second excited state.

The energy dependence of $`\mathrm{\Omega }(E)`$ is shown in Fig. 6 for the ground state. To determine the binding energy from the condition $`\mathrm{\Omega }(E)=0`$, we fitted $`\mathrm{\Omega }(E)`$ by using the linear energy-dependence of $`\mathrm{\Omega }(E)=0`$ near the binding energy with the least square search method. Since values of $`\mathrm{\Omega }(E)`$ fluctuate due to sampling errors, we performed calculations several times with different sets of random numbers. Averaging these results, we obtained the energies of three bound states and errors estimated by the standard deviation, as shown in Table I. Our results with the Monte-Carlo method are in a good agreement with those of the ‘exact’ calculation for all three states.

Once we obtained the binding energies, the eigenfunctions were given by the MCF wave functions $`|F_i`$ and matrix elements $`x_i`$. The wave functions for the ground, first, and second excited states are shown in Fig. 7. The wave functions are also in good agreement with the ‘exact’ calculation. Note that we started with a trial function without nodes for the second excited state. We could use more realistic trial functions with nodes, and if we did so we would obtain faster convergence. However, our result shows that it is not necessary to know the correct nodes of the wave function when we construct a trial function in this method. In summary, our new Monte-Carlo method works for the calculation of binding energies and wave functions including the excited states, and we believe that our method will be useful in studying few-nucleon bound state systems.

### B Four-nucleon bound state as an example of a multi-dimensional problem

As a realistic example, we study a four-nucleon bound state in three-dimensional space. The Schrödinger equation is expressed as follows:

$`{\displaystyle \underset{k=1}{\overset{3}{}}}{\displaystyle \frac{1}{2\mu _k}}_{r_k}^2|\varphi +{\displaystyle \underset{i<j}{\overset{4}{}}}V_{ij}|\varphi =B|\varphi \text{ .}`$ (46)   
Here, the $`𝒓`$<sub>k</sub> are the Jacobi coordinates defined as

$`\{\begin{array}{ccc}𝒓_1\hfill & =& 𝒙_2𝒙_1\text{ ,}\hfill \\ 𝒓_2\hfill & =& 𝒙_3(𝒙_1+𝒙_2)/2\text{ ,}\hfill \\ 𝒓_3\hfill & =& 𝒙_4(𝒙_1+𝒙_2+𝒙_3)/3\text{ ,}\hfill \end{array}`$and the $`\mu _k`$ are reduced masses corresponding to the coordinates $`𝒓_k`$:

$`\{\begin{array}{ccc}\mu _1\hfill & =& m/2\text{ ,}\hfill \\ \mu _2\hfill & =& 2m/3\text{ ,}\hfill \\ \mu _3\hfill & =& 3m/4\text{ .}\hfill \end{array}`$

Here, we adopt the Volkov potential, which is a central potential and has a soft core at short range:

$$V_{ij}=V_1\mathrm{exp}(a_1r_{ij}^2)+V_2\mathrm{exp}(a_2r_{ij}^2)\text{ ,}$$
(47)   
with

$`\{\begin{array}{cc}V_1=83.34\text{ (MeV) ,}\hfill & a_1=1.60^2\text{ (fm}\text{-2}\text{) ,}\hfill \\ V_2=144.86\text{ (MeV) ,}\hfill & a_2=0.82^2\text{ (fm}\text{-2}\text{) ,}\text{ .}\hfill \end{array}`$Since the ground state energy of <sup>4</sup>He using this potential has been studied by various methods such as the Faddeev-Yakubovsky (FY) method taking into account only the $`S`$-wave component, the stochastic variational method (SVM) and the hyperspherical harmonics (HH), it is a good example to examine our method in a four-body system.

The free Green’s function in dimensionless coordinates is expressed as

$$G_0(R,R^{})=\frac{1}{2(2\pi )^4}\frac{1}{^4}\left(1+\frac{6}{}+\frac{15}{^2}+\frac{15}{^3}\right)\mathrm{exp}()\text{ ,}$$
(48)   
where $`|RR^{}|`$ is a relative distance in the nine-dimensional space. As the dimensionality increases, the free Green’s function becomes more singular at $`=0`$. It is rather difficult to treat the singularity of the free Green’s function with the conventional integration method. Since the random walk method inverts the free Green’s function analytically, we can obtain numerically stable results. A standard method to generate points with the probability $`G_0()`$ is given in Ref. .

The ground state of <sup>4</sup>He is expected to be spatially symmetric and anti-symmetric in spin-isospin space. Then we take a trial function $`|F_0`$ given by the product of a Gaussian wave function $`|\phi _s`$, a correlation function $`g(r)`$ and the spin-isospin function $`|\mathrm{\Xi }`$, as

$$|F_0=𝒩\underset{i<j}{}g(r_{ij})|\phi _s|\mathrm{\Xi }\text{ ,}\text{}$$
(49)   
with

$$|\phi _s=\underset{i<j}{}\mathrm{exp}(r_{ij}^2/2b_0^2)\text{ ,}$$
(50)   
where $`b_0`$ is a size parameter. We used the form of correlation function $`g(r)`$

$$g(r_{ij})=[1\alpha \mathrm{exp}(r_{ij}^2/b_g)]^2$$
(51)   
with

$`\alpha =0.11\text{ , }b_g=0.5\text{ (fm}\text{2}\text{)}\text{ .}`$The parameters of the correlation function are determined to simulate the wave function of a two-nucleon bound state in this potential. The size parameter $`b_0`$ of the Gaussian wave function in (50) is simply determined by the variational method with one parameter. For the trial function $`|f`$, we assume the same functional form as $`|F_i`$, and take a slightly larger value of the size parameter $`b_0`$,

$`b_0`$ $`=`$ $`3.0\text{ (fm) : for }|F_0`$ (52)   
$`=`$ $`4.5\text{ (fm) : for }|f\text{ .}`$ (53)   
Since the Volkov potential is a central potential with no exchange terms, the symmetry of the spatial and spin-isospin wave functions are preserved separately.

The convergence of $`x_i`$ is displayed in Fig. 8, which shows that three-times iterations are enough in this example. Similarly, $`\mathrm{\Omega }(E)`$ also converges rapidly, as shown in Fig. 9. Here, we take $`N=200,000`$ sampling points.

The bound state energy and its uncertainty are estimated in the same way as in the one-dimensional case. We assume a linear function for fitting the energy dependence of $`\mathrm{\Omega }(E)`$, as shown in Fig. 10. As a result, we obtain $`E=29.88\pm 0.59`$ (MeV). Our results for the binding energy of <sup>4</sup>He are in good agreement with the results of the other methods within uncertainties, as tabulated in Table II. It is shown that our method also works for higher dimensional examples.

## V Summary and discussion

The Monte-Carlo method is useful in the study of the bound state of the many-nucleon systems because of its simplicity and efficiency for multi-dimensional calculations. In this work, we propose a new method to study the bound state of a few-nucleon system without using the diffusive-type projection method. Therefore it is expected to work not only for the ground state but also for excited states.

We used the method of continued fractions (MCF) as the main formalism of the calculation. The advantages of this method are its ability to calculate the excited states and the rapid convergence of the continued fraction series. We modified the original formula of MCF into a suitable form for the Monte-Carlo random walk method. We generated MCF wave functions by the random walk method. Here, it was shown that the cancellation mechanism of lower order MCF wave functions is important, and we should use our method of counting weight in this calculation to avoid the “sign problem”.

Our method was examined in application to a one-dimensional problem in detail. We demonstrated that our new method works even for the calculation of excited states. We also studied the bound state of a four-nucleon system and obtained results that compare well with those of other methods. These facts suggest that our method is a powerful method to study the bound states of few-nucleon systems.

Application of this method to few-nucleon systems with realistic nuclear potentials including exchange terms is one of our future problems. Here, we have to deal with the coupled channel problem of wave functions with various types of exchange symmetry. It is important to retain the correct symmetry of the wave function in each iteration step. It is also necessary to take into account the correlation of nucleons in the trial function as much as possible. The wave function obtained by a fully variational calculation should be used as a more plausible trial function in our procedure. This will also accelerate the convergence of the continued fraction series.

## Acknowledgements

The authors would like to thank Professor H. Ohtsubo for many critical discussions and comments.