- Email: [email protected]
A Monte Carlo method for nuclear many body problems
Nuclear Physics 3 2 (1962) 8 2 - - 9 7 ; ( ~ ) North-Holland Publishing Co., Amsterdam Not to be reproduced by photoprint or microfilm without written permission from the publisher
A M O N T E CARLO M E T H O D FOR NUCLEAR M A N Y B O D Y PROBLEMS E R I C H W. S C H M I D Florida State University, Tallahassee, Florida and
Oak Ridge National Laboratory, Oak Ridge, Tennessee t R e c e i v e d 7 A u g u s t 1961 A M o n t e Carlo m e t h o d is g i v e n to solve n u c l e a r m a n y b o d y p r o b l e m s in t h e f r a m e of t h e R i t z ' v a r i a t i o n a l principle. T h e a n t i s y m m e t r i z a t i o n a n d all i n t e g r a t i o n s a r e p e r f o r m e d b y a r a n d o m process. T h e m e t h o d is f o r m u l a t e d for a t w o c l u s t e r s y s t e m , u s i n g oscillator w a v e f u n c t i o n s a n d h a r d core forces, b u t it c a n be generalized to r a t h e r a r b i t r a r y forces a n d w a v e f u n c t i o n s . As a t e s t e x a m p l e , t h e g r o u n d s t a t e e n e r g y a n d r m s - r a d i u s of t h e ~-particle h a s b e e n calculated. F o r t h e final p a r a m e t e r values, a n a c c u r a c y in t h e e n e r g y of 10 % (2 %) w a s o b t a i n e d a f t e r a c a l c u l a t i o n t i m e of 40 sec (16 rain) w i t h a n I B M 7090 c o m p u t e r . T h e a c c u r a c y in t h e r a d i u s is by a factor 8 b e t t e r . T h e m i n i m u m e n e r g y d e p e n d s r a t h e r s t r o n g l y o n t h e f o r m of t h e t r i a l w a v e f u n c t i o n n e a r t h e h a r d core.
Abstract:
1. Introduction T i l e Ritz' variational principle, which is usually employed to calculate properties of light nuclei, demands that one minimizes the expression
E
=
f 9*HAgdT f 9*A9dz
"
(1)
with respect to parameters in the trial wave function 9. In this equation, H is the Hamiltonian of the system and A denotes the complete antisymmetrization of the wave function. We assume H to be symmetric with respect to an exchange of any two particles (that means A can be exchanged with H), therefore the antisymmetrization of 9* can be dropped. Recent forms for the phenomenological interaction between two nucleons all have a repulsive core of small radius which is also considered to prevent a nucleus from collapsing. The hard core imposes rather complicated conditions on the wave function % since 9 has to be zero within and at the surface of the core. These boundary conditions, together with the antisymmetrization make an analytical evaluation of the expression (1) so complicated and lengthy, that it is extremely difficult to apply it for nuclei more massive than a triton or ~-particle. t \Vork s u p p o r t e d in p a r t b y t h e Office of N a v a l R e s e a r c h . C o m p u t i n g t i m e on t h e I B M 7090 c o m p u t e r a t E g l i n A F B f u r n i s h e d w i t h o u t cost b y A F O S R t h r o u g h c o n t r a c t A F 49(638)-278. 82
NUCLEAR MANY BODY PROBLEMS
83
Therefore, a Monte Carlo method has been developed which allows one to calculate E from eq (1) without any analytical calculation. The integrations as well as the antisymmetrization, are performed b y a random process. The basic idea of the method is the following. The integrands ~ * H A ~ ( : ,¢*AH~v) and ~*A~v are written in the form p ( r l . . , r N ) w ( r l . . , r~), where p is called the probability function and contains essentially ~* A (for details see next paragraph), w is called the weight function and contains essentially H~v or ~, respectively. Random co-ordinates are generated with the probability distribution p ( r ~ . . , rN) and a certain number of terms from the sum represented b y A are selected at random according to a uniform distribution over all N! exchange terms. For these random numbers, together with the exchange orders (which also contain a factor + 1 according to even or odd permutations), the weight functions are calculated. The values obtained in this way are estimates for the integrals S ~v*HA~vdz and S ~v*A~vdr and the averages of a large number of these estimates converge toward the value of the integrals. In principle, any kind of variational function and any kind of interaction force can be used in this method, but the calculation time and also the length of the computer program imposes certain restrictions. Test calculations with the ,t-particle have shown that the method will be useful mainly for the following problems: 1) Three and four body problem. Rather elaborate variational functions and interaction forces can be used to get information on nuclear forces in addition to the information following from the two-nucleon scattering data. 2) Light nuclei up to carbon or oxygen. Some interesting properties of the nuclei (energy levels, radii, 15- and y-transition probabilities, clustering effects) can be investigated, using rather simple central forces of Serber type with a hard core and cluster model wave functions with few, but effective, variational parameters. Emphasis in this work is put on the second problem. The Monte Carlo method will be discussed in sect. 2. Sect. 3 tells something about the computing technique and sect. 4 gives the result of a test calculation with the ~-particle. 2. Method
An integral I = f p ( x l . . . X3N)W ( X l . . . x3N)dxl.., dX3g
(2)
can be evaluated in the following way 1). One generates random co-ordinates x~i, x2 ~..... x~N according to the normalized probability distribution p(xl... XaN). For every set of 3N random co-ordinates the value of the weight function W(Xl. . . X3N) ;~ calculated: Wi . . W. ( X. l i i x3:~). (3)
84
ERICH
W.
SCHMID
The value w~ is an estimate for the integral I, and the average of the estimates converges under a certain condition to the integral: w i
--->I.
(4)
The condition of convergence is, that the variance v, defined as v = f p (X)l.
• . Z3N[W
(X1 .
. . X3N ) - 1 ]
2d x l . . . d x s N
(5)
is finite. The variance gives also the statistical error A I of an average of n estimates (standard deviation)" (6)
AI = V~/n.
From this relation it can be seen already that the error will decrease with the square root of the computing time. The integration method discussed above is called a "Monte Carlo Method". We want to use it to carry out the integrations in a Ritz' variational calculation for light nuclei. In order to give a simple description of the method, we shall restrict ourselves in the following to a two-cluster system *) described by oscillator functions and central exchange forces with a hard core. The trial wave function has the form (7)
W = 919~z H l(r,k)S, i>1¢
where 91 and 92 describe the internal motion of the two clusters and Z describes the relative motion, while S is the spin and isobaric spin function. The functions 91, 95 and Z fulfill the boundary condition at infinity but not in the range of the hard core. The "cut off" factor I I , > k/(r~k) is used to give the right boundary f(r)
r¢
d
F i g . 1. T h e " c u t o f f " f u n c t i o n
1"
/(r).
t I n p r i n c i p l e , e v e r y w a v e f u n c t i o n d e s c r i b i n g a n u c l e u s c a n b e w r i t t e n a s a s u p e r p o s i t i o n of t w o c l u s t e r w a v e f u n c t i o n s z).
NUCLEAR MANY BODY PROBLEMS
85
condition at the surface of the hard core ¢. The function /(r) has the form shown in fig. 1. For values of r smaller than the range of the hard core, ](r) is zero and [(r) is equal to 1, for large values of r. The smallest value of r, for which /(r) is equal to 1 defines the "healing distance" d, that is the separation distance of two nucleons beyond which the hard core no longer influences the wave function b y its boundary condition. The Hamilton operator of the system has the form ?~2 S -- 2M
X
aU particles
v,2+
X
I
e~ ekX
all pairs
(s)
where
V,k ~ V(r,k)Ew+mP*(i, k)+bP~(i, k)--hP'(i, k)l
(da)
is the two body potential and P~(i, k), P*(i, k) and P*(i, k) are the operators which exchange the space, spin and isobaric spin co-ordinates of particle i and k. The numbers e~ are equal to the proton charge if particle i is a proton and zero if it is a neutron. We assume that the wave function ~v does not contain any motion of the total centre of mass; therefore, we can take the centre-ofmass as origin of our co-ordinate system. 3 (N-- 1) co-ordinates (N is the number of particles) are then independent, while 3 co-ordinates are determined b y the condition N
X r, = O.
(9)
The integrands ~v*HA~v and ~p*A~o of eq (1) have to be written now in the form p ( r l . . . r n ) w ( r l . . , rn). We use the same probability function p for both integrands, therefore we do not have to care about the normalization of p. In principle p ( r l . . . rn) could be taken to be different from zero and constant within a rectangular box in all 3N co-ordinates which encloses the integrand at all points where it is essentially different from zero. The function w ( r l . . . rN) would then be the integrand itself (besides a normalization factor). But this would give a large variance, since in a many-dimensional box the main part of the volume consists of "corners" of the box, where the integrand is close to zero. It can be seen from eq. (5) that the variance is small when p represents the main part of the integrand and w is a function which varies as little as possible**. This requirement is most significant in the region where the integrand is falling off to zero, since this region represents the main part of tile integration volume. ? The " c u t off" factor II~>~/(r~) does n o t change the geometrical properties of the wave function; L, S and the p a r i t y r e m a i n good q u a n t u m n u m b e r s . Therefore, the splitting-off of the centre-of-mass m o t i o n is no p r o b l e m for w a v e functions of the f o r m (7), unlike m o s t other m a n y particle w a v e functions (see also refs. 3-~)). it The variance would be zero if p is the integrand and w i s a c o n s t a n t n o r m a l i z a t i o n factor which is in this case the integral itself.
86
ERICH W. SCHMID
The possibilities of choosing p ( r l . . , rN) are restricted by the difficulty of generating random numbers. For a distribution function depending only on one coordinate, the generation of random numbers is easy. A transformation function can be obtained by a one-dimensional numerical integration which transforms random numbers of a uniform distribution into random numbers of the desired distribution (for details see sect. 3). For some distributions like exponential and Gaussian distributions, library programs exist, which generate the random numbers in a more direct way. For a general distribution function depending on m a n y co-ordinates the generation of random numbers is difficult *. We demand therefore, t h a t p is either a product of functions each of which depends only on one co-ordinate or that p is a linear combination a P l + b p ~ + ••• of such products of functions. In the first case every factor for itself is used to generate a random co-ordinate. In the latter case the normalized distribution functions Pl, P, . . . . are used alternatively according to their weights a, b. . . . . As trial wave functions for the internal and relative motion of our clusters we adopt the usual oscillator wave functions, which have the form given below. An unexcited m-particle cluster is described by ~(~) = exp [--~1 a (u 2 + v 2 + w 2 )],
(10)
where r 1 = {(--u--v--w)+R~, r 3 =
r 2 = }(u--v+w)+R~,
½(u+v--w)+R~,
r4 :
(11)
½(--u+v+w)+R~
are the position vectors of the 4 nucleons, R~ is the centre-of-mass vector and a is a variational parameter. A triton trial wave function is described by ~(tr) = exp [--½a 3 2 +½v
)l,
(12)
with rl = ½(u+v)+Rtr,
r2 = ½ ( u - - v ) + R t r ,
r3 = - - u + R t r .
(13)
For a deuteron cluster, a Gaussian function is a good description only when the cluster is strongly bound, If it is loosely bound, an exponential function or a superposition of different Gaussian functions is better: ~o(d) = exp(--½au 2)
~(d) = exp(--~a[u[),
or
(14)
with r I =
1 -flU +
Rd,
r2 = - - l u
+ Rd.
(15)
The relative motion X has the form Z =
R"exp(--aR2)Y~,o(
v~, 9 )
* In order to obtain them by transformation from a uniform distribution, transformation w h o s e J a c o b i a u d e t e r m i n a n t is a g i v e n f u n c t i o n .
(16) one has to find a
NUCLEAR
MANY
BODY
87
PROBLEMS
with R 1 -- -
M S -
R,
MI+M2
Rz -- -
--M1 -
R,
MI+M2
(17)
where M 1, M 2 are the masses of cluster 1, 2). Eq. (17) states, that the centreof-mass of the two-cluster system is at the origin. The functions (10), (12), (14) and (16) are all of the type described above, that is, they are products of one-dimensional functions. This is also true for the oscillator wave functions of bigger clusters, since they can be built up of smaller clusters in the above way. The factor ~ 2 ~ g is contained in ~* as well as in ~. The factor in ~2can be split into a factor (91~2X)s which is symmetric with respect to any exchange of particles and a factor which is not symmetric. As probability distribution p we can take then, if we forget for a moment about the statistical antisymmetrization, P = ~l~$Z(~192Z)8.
(18)
The star coming from ~* has been dropped since the above mentioned functions are real. The weight functions for numerator and the denominator of eq. (1) are 1 wl=~*HA %
1 w2-~/*A %
(19)
respectively. The random co-ordinates generated according to eq. (18) are transformed into cartesian co-ordinates by the transformations (11), (la), (15) and (17). The quantities w1 and w2 are calculated numerically in cartesian co-ordinates. For more complicated trial wave functions than the ones described above (as for instance trial functions which include higher oscillator shells), we can use the same distribution function as above, with parameter values which give about the same or a slightly larger rms-radius as would arise from the trial wave function. Eqs. (19) are also valid in this case. The exact form of the probability function p is not very important, since it has only an influence on the variance. In the method described in this paper the main contribution to the variance comes from the Hamiltonian in the region where two or more particles are close together. A somewhat different probability function would therefore increase the variance only little, while the computing time as well as the preparation time will be smaller when a simple probability function is used. The weight functions (19) are still inconvenient because of the antisymmetrization operator A. This operator splits the weight functions into a sum of "exchange terms". Many of these terms are usually zero because of the orthogonality of the spin functions. Some of them may be equal to each other. But there is usually still a large number of different terms. We want, therefore, to treat the antisymmetrization by a random process.
88
E R I C H W. SCHMID
The antisymmetrization operator A is defined by A ---- X ± P ( 1 , 2 . . . . . N).
(20)
The sum extends over all permutations of the N nucleons, the sign being taken according to whether the permutation is even or odd. Instead of operating with the whole sum of operators on the wave function ~ we take only a certain number m of them (together with the right sign), which we select at random. The m permutation operators could be selected from a uniform distribution over all N! operators. But in this case most of the operators would give zero contributions to the estimates of our integrals. In order to find the permutation operators which give non-zero contributions we have to look at the spin and isobaric spin function S. Since the trial wave function m a y be chosen to have any symmetry, S can be chosen in a very simple form. Usually it is taken as one single product of single-nucleon spin and isobaric spin functions, or a linear combination of few of such products. We assume in the following, that it has the form of one product. Then the only permutations which give a non-zero contribution to the normalization, the kinetic energy and the part of the potential energy which does not involve spin and isobaric spin exchanges, are the ones for which only particles in the same spin and isobaric spin state are exchanged. From a uniform distribution of these permutation operators, we select at random m operators. For the calculation of the Bartlett and Heisenberg part of the potential energy, where we have spin and isobaric spin exchanges, we remember the following: In eq. (1) we are allowed to carry out any permutation of particles in ~v (and change sign if the permutation is odd) and still get the same result because of the operator A. So we use for the calculation of the Bartlett term instead of the wave function ~v the function --P(i, k)~v. If the spin exchange operator of the Bartlett term is applied to this, we get --P*(i, k)P(i, k) = --P*(i, k)er(i, k). That means, in case of different spins of particle i and k, we can replace in the calculation the spin exchange operator P~(i, k) by zero if particles i and k have different isobaric spins and by the operator --P~(i, k) if the isobaric spin is the same. If the particles i and k have the same spin the exchange operator is simply replaced by + 1 . The Heisenberg part of the potential energy is treated in a similar way. We can now write instead of eqs. (18) and (19) in a symbolic way p = (p1(P2%(Vi(p2%)s A : 1
Wl.~y)*Sv) ,
~b'A, 1
W2 -~W*~V ,
(21)
(22)
with the understanding, that A in (21) means statistical antisymmetrization and is applied only on %
NUCLEAR MANY BODY PROBLEMS
89
The statistical antisymmetrization will certainly increase the variance v. H o w many permutations one has to take for every estimate has to be found out b y a test calculation. For small nuclei, or nuclei of high symmetry like Be 8 it is probably better to use eqs. (18) and (19). When the parameter values have been found, which give the minimum energy, the Monte Carlo method can be applied in a similar w a y as for the Hamiltonian to calculate the expectation values of other operators, like the mean square radius, fl- and 7-transition probabilities.
3. Computing Technique For Monte Carlo calculations of the kind described in the last section a fast electronic computer with a large core storage is necessary, although the calculations are rather simple. But in order to get fairly accurate results, a large number of estimates has to be calculated and in general the application of the method is limited b y the length of the computing time. Also the questions arising while writing the computer programs are mainly concerned with the speed of the computation. In the following, some of these questions will be discussed. Random numbers of certain probability distributions like uniform distribution, exponential, and Gaussian distribution are obtained b y subroutines which are distributed b y the S H A R E organization. Time can be saved b y calculating a set of random numbers only once and storing it on tape. Random numbers of other one dimensional distribution functions can be obtained from a uniform distribution b y a transformation. The transformation, which transforms a uniform distribution of x in the interval (0, 1) into a normalized distribution p(y) is obtained b y a numerical integration of the differential equation dy 1 d~ = p(y~)"
(23)
The integration has to give a table of y = y(x), which is stored in the core storage of the computer. Since the table has to be narrowly spaced in order to spare interpolation, a simple integration rule can be used, as for instance h - , Yn+l--Yn _- - -p(~)
(24)
where h is the spacing in x and ~ is a value between y~ and Yn+1 which can be taken to be y,~+½ (Y~--Yn--1) in regions where Y~--Yn-1 and Y,~+I--Yn are not too different ¢. ¢ This c o n d i t i o n is for i n s t a n c e n o t fulfilled in t h e case p(y) ~ c y 4 e x p ( - - a y 4) n e a r y = O. T h u s for s m a l l v a l u e s of y i t is b e t t e r to i n t e g r a t e a n a l y t i c a l l y w h i l e p u t t i n g t h e G a u s s i a n f u n c t i o n e q u a l to 1.
90
ERICH W. SCHMID
Random numbers of a distribution function I5 = a p l + b p 2 +
...
are generated using a fork, which causes the program to generate numbers either from Pl or P2 and so on, according to the weight factors a, b. . . . The generation of random permutation operators can be done in the following way, when the function S is a single product of nucleon spin and isobaric spin functions. The N nucleons are divided into 4 groups of particles, namely protons with spin up, protons with spin down, neutrons with spin up and neutrons with spin down. As discussed in the last section only exchanges among the particles of each group need to be considered and we can therefore treat each group separately. We define two N-dimensional vectors a and b. The vector a has as its components a, the numbers of our N particles in an order, which is arbitrary except that the first N 1 components are the numbers of protons with spin up, the next N~ particles are protons with spin down, and so on. The vector b will contain the permuted numbers. We start by choosing at random one out of the first N1 components of b. This can be done by generating a random number from a uniform distribution between 0 and 1, multiplying it by N1, adding + 1 and truncating it to an integer. The result tells whether the first, second or N1 t~ of the components in question has to be taken. This component we put equal to al. From the still undefined N 1 - 1 components of the first group we choose at random another one and put it equal to a 2. When all N 1 components of the first group are defined, the procedure continues with the next group and so on. Finally b will contain a permutation of the components of a. Every permutation which will give a non zero contribution to our integrals has the same probability in the above generation process. Whether a permutation is even or odd can be found out in the usual way by counting the exchanges which are necessary to get a from b. For the calculation of the kinetic energy, the second partial derivatives of the wave function ~v with respect to the nucleon co-ordinates have to be calculated numerically. A second partial derivative is defined as 32~p _ lim yj(.. " x j - - e . . . ) - - 2 ~ 0 ( . . , ~X t 2
~-*0
xs...)+v2(..,
xj+e...)
(25)
8"2
We use this equation with a proper value of e, but instead of calculating ~v ( . . . 2 j 4 - E . . . ) directly, the following procedure is used: We expand everyone ofthe factors ~1, 92, Z,/(r,~) which build up ~v according to (7) into a Taylor series in e, which we cut off after the e2-term. These Taylor series we multiply together in order to get ~v(... x~-¢-e...). Now, for a computer with 8 digit numbers, we can choose e so, that if the ,0 term starts in the first digit, the el, e~, E8 and e4 terms will start approximately in the 3 ra, 5 th, 7th and 9 th digit respec-
91
NUCLEAR MANY BODY PROBLEMS
tively, that means the # term is truncated b y the machine. Inserting the thus obtained ~0(... x ~ 4 - e . . . ) into eq. (25) we see that the e° term and all odd terms cancel, and tire second partial derivative is obtained with a four place accuracy. If all coefficients of the Taylor series are tabulated and kept in the core storage, the calculation of each Taylor series takes only about two multiplications and two additions. Nevertheless, because of the 3N derivatives and the N - - 1 factors ](rik ) which contribute to every derivative, the calculation of the kinetic energy estimates is the longest part of the calculation. In the calculation of the potential energy estimates, the Wigner part presents no difficulties. For the Majorana Bartlett and Heisenberg terms twoparticle exchanges have to be calculated. When oscillator wave functions of the type (10), (12), (14) and (16) are used, the wave function can be split into a factor which is symmetric and a factor which is not symmetric with respect to particle exchanges. Only the latter one is affected b y the exchange operations. In the computer we write the spin function S (if it is a single product) as an Ndimensional vector of numbers si. The value s~ = 1 means particle i is a proton with spin up, s i = 2 proton with spin down, s, ~ 3 neutron with spin up and s~ = 4 neutron with spin down. The sum and the difference of s t and s, is enough to decide, which operator has to be applied on ~pin order to get a nonzero contribution to the potential energy estimate. Table 1 gives these operators for every possible combination of particles. TABLE 1 Potential energy operators for various combinations of particles i and k s~+s~
Operator
s~
sk
s,--s~,
1
1
0
2
V(r~) ( w + b - - h + m P ~ k ) + - -
2
2
0
4
g(r~) (w+b--h+rnI~i~) + r ~
3
3
0
6
V(r,,) ( w + b - - h + m P ~ , )
4
4
0
8
V(r~,) ( w + b - - h + m P ~ , )
1
2
:#0
3
V (r,,) [w -- h + (m-- b) P~,] + - -
1
3
:#:0
4
V(r,~) [ w + b + (m +h)P~]
1
4
~:0
5
V (r~,) (w + toPOi)
2
3
v~:O
5
V(r,,) (w +raP;,)
2
4
:/:0
6
V(r,,) [ w + b + (m+h)P~,]
3
4
:/:0
7
V(r~,) [w--h:~- (m-- b) P~,]
82
e2
$2
Yik
92
:ERICH W. SCHMID
The normalization estimates arise as intermediate results during the potential energy calculation. For the speed of the calculation it is important that all functions which are used are calculated and stored as tables in the core storage before the Monte Carlo calculation starts. The tables have to be narrowly spaced so that no interpolation is necessary. In this way also the large storage capacity of modern computers can be utilized to gain speed. The calculation time can be estimated by counting the multiplications, divisions, additions and subtractions. All other operations do not count in a first approximation. In order to calculate the statistical errors, we average not only the estimates w,, but also their squares w, 2. The variance is then calculated approximately from v =--
T~ 5=1
w~2-- n
4=1
wi ,
(26)
which follows from eq. (5) if one inserts (1/n)~ w~ as an approximation for I. The error is obtained by eq. (6). When adding up a large number of estimates, it might happen, that the smaller ones get lost by truncation, if the sum is already very large. This can be avoided by calculating intermediate averages after every 104 or 10~ steps.
4. E x a m p l e Monte Carlo methods have been used by G. Derrick and J. M. Blatt 6) to calculate the binding energy of the triton and by L. Cohen and J. B. Willis ~) to calculate the binding energy and rms-radius of the a-particle. Cohen and Willis report for the a-particle calculation with hard core forces a possible error of 10 % in the energy after a calculation time of 40 h (with their final values for the variational parameters). As a test of the method discussed in this paper, the binding energy and rmsradius of the a-particle has been calculated. The calculations should also show whether the shape of the cut-off function [(r) has a strong influence on the result of the variational calculation or not. The two-nucleon potential was assumed to be of the form (Sa) with
V(r) = t + ~ t V exp[--fl(r--rc) ]
if if
r <
re,
Y ~
YC,
(27)
and the Serber relation m =
w =
¼(l+x),
h =
b =
¼(l--x),
('27a)
where x is the ratio of the singlet and triplet well depth. The hard core radius
NUCLEAR
MANY
BODY
03
PROBLEMS
r e was assumed to be re = 0.4 fm. The constants V, fl and x were calculated with formulas given by T. Kikuta, M. Morita and M. Yamanda 8) and T. M. Blatt and J. D. Jackson 9) in order to fit the two-nucleon data: n-p triplet scattering length a t =
5.378 fm,
n-p triplet effective range
1.704 fro,
rot =
n-p singlet scattering length a s ----- --23.69
fm,
n-p singlet effective range
fm,
ros-~
2.4
For the n-p singlet effective range Kikuta et al. used the two values 2.4 fm and 2.7 fm because of the large experimental error of this quantity. In this paper the value 2.4 fm has been preferred since in a previous paper 10) the singlet force constants have been fitted from the more accurate p-p scattering data and a calculated n-p effective range of 2.41 fm was obtained with these constants. Nevertheless, the uncertainty of the n-p singlet effective range was considered large enough to justify the use of the same well parameter/3 for both the singlet and triplet interaction ¢. The following values for the potential parameters were obtained: V -~ 498.02 MeV,
/3 = 2.5885 fm -1,
x = 0.771.
(28)
A numerical integration showed that these constants also give the right deuteron binding energy within the experimental error. The trial wave function was assumed to be
(1) (1)(o)(1)(1)(o)(o)(o)
~p = e x p ( - - l a ~, ri~) I I ](ri~) 0 i ~
i~>k
0
O"l
1
T1
0
0"2
0
72
1
0" 3
1
T3
0"4
(29)
1 'r 4 "
The Gaussian function is the same as in eq. (10). The quantities
(10k(10t are the usual Pauli spin-and isobaric spin vectors; a is a variational parameter. For the cut-off function ](r), two different trial functions ]l(r) and ]~(r) have been used. Both functions are zero within the range of the hard core and both are 1 for values of r larger than the healing distance d. Between these two points, the functions have the following form: t I n t h e c a l c u l a t i o n of t h e force c o n s t a n t s t h e singlet i n t r i n s i c r a n g e w a s c h a n g e d a little in o r d e r to give t h e s a m e / ~ as in t h e triple case. Since t h e well d e p t h p a r a m e t e r is close to 1 in t h i s case, t h e s c a t t e r i n g l e n g t h is n o t i n f l u e n c e d b y t h i s procedure.
~4
ERI CH W. SCHMID
1) The function/1 is a numerical solution of the two-nucleon Schr6dmger equation [
~ dr d~2 + Vt(r) 1 [l(r) = 0 M
(30)
(re < r < d).
The first maximum of the solution/l(r) is normalized to 1 and defines the healing distance d. The potential Vt(r ) which "generates" the trial function t1(') w a s Vt(, ) = --V t e x p f - - f l ( r - - r e ) ].
(31)
The parameter fl was the same as in eq. (28) but the well depth V t was used as a variational parameter. The healing distance d and the well deptL V t depend on each other, therefore we can also call d the variational parameter. 2) The function/2(r) had the form of the first quadrant of a sine wave: /2(r) = sin ~ ( r - r e ) 2(d--re)
(re < r < d),
(32)
where d was considered to be a variational parameter. Because of the high symmetry of the wave function (29), the Bartlett-and Heisenberg part of the interaction cancels. In all other parts of the interaction as well as in tile kinetic energy part and the normalization, the antisymmetrization operator can be dropped since only the direct term gives a non zero contribution. This example, therefore, can give no test for the statistical antisymmetrization. For the probability distribution tile square of expression (10) was taken according to eq. (18): p = exp[--a(u2+v2+w~)].
(33)
The random co-ordinates u, v, and w were transformed into cartesian co-ordinates by eq. (11). The weight functions were (with either/1 or/2 instead of/) according to eq. (19) _he wl = exp(+-~a X r~k) 1-I / ( r , k ) ~ • V, 2 i>k
i>k
j
e x p ( - - } a X r~k) H ~if,k) + H i>k
i~k
/(r,k) 2
(w+m) X V(r,k)+ r
i~k
=
II/('
'
i~k
kV.
(35)
The weight function wa for the mean square radius was w~ = w2rl 2.
(36)
The computer program was written in FORTRAN language and had some
NUCLEAR
MANY
BODY
PROBLEMS
95
170, mostly short, statements. The rather large number of statements was necessary to optimize the speed of the calculation ¢. The computing time was around 4 mill for every 10000 estimates with an IBM 7090 computer (the tabulation of functions takes only a few seconds). TABLE
2
E n e r g y a n d r m s - r a d i u s for c u t - o f f f u n c t i o n ]t a n d different p a r a m e t e r v a l u e s E (MeV)
r m s (fro)
--23.0+0.5 --25.9±0.5 --27.8±0.5 - - 29.3 ± 0.4 --27.8±0.6 --24.7±0.7 --21.4=[_0.7 --25.7 ~ 0 . 7 --29.4±0.4 -- 28.2 J::O.4 -- 1 5 . 6 ± 0 . 6
a (fm -~)
1.53 1.41 1.31 1.23 1.17 1.12 1.07 1.23 1.24 1.24 1.27
0.5 0.6 0.7 0.8 0.9 1.0 1.1 0.8 0.8 0.8 0.8
TABLE
Vt(MeV); d(fm) 1000; 1000; 1000; 1000; 1000; 1000; 1000; 1200; 950; 800; 500;
0.887 0.887 0.887 0.887 0.887 0.887 0.887 0.821 0.908 0.992 1.45
3
E n e r g y a n d r m s - r a d i u s for c u t off f u n c t i o n /2 a n d different p a r a m e t e r v a l u e s E (MeV)
r m s (fro)
a ( f m -~)
d (fm)
1.22 1.23 1.24 1.25 1.40 1.31 1.28 1.24 1.12
0.8 0.8 0.8 0.8 0.6 0.7 0.74 0.8 1.0
0.757 0.821 0.887 0.992 0.86 0.86 0.86 0.86 0.86
--22.8±1.0 --25.3±0.9 --25.4~0.7 21.6±0.6 --24.8±0.6 --25.8±0.5 -- 25.9::k0.4 --25.6i0.5 --21.1 ± 0 . 9
Tables 2 and 3 show the dependence of the energy and rms-radius on the parameter values for the two different cut-off functions. The minimum values are Emt n =
- - 2 9 . 4 i 0 . 4 MeV,
(rms) = 1.24 fro,
d = 0.91 fm
(37)
(rms) = 1.28 fm,
d = 0.86 fm
(38)
for the first case with /1, and Emt n =
--25.9±0.4 MeV,
for the cut off function /2. * A factor 1.5 in speed could p r o b a b l y be gained b y w r i t i n g t h e p r o g r a m in m a c h i n e l a n g u a g e , b u t t h i s w o u l d increase t h e p r e p a r a t i o n t i m e b y a n u n b e a r a b l e a m o u n t .
96
ERICH
W.
SCHMID
The function ]1 gives a deeper minimum and a smaller error. The statistical error in the energy near the minimum is 10 ~ (2 ~ ) after a calculation time of 40 sec (16 rain) with [1 or after a calculation time of 100 sec (40 rain) with ]2. The error increases for smaller healing distances. The accuracy in the normalization and the radius is better by a factor 8. The error in the energy arises mainly from the region between hard core and healing distance, since the energies there are very large. The function ]1 makes a better balance between potential and kinetic energy at every space point in this region and gives, therefore, a lower variance. The result of the variational calculation tells that the form of the short ranged correlation function of the nucleons is quite important. The concept of using a strict solution of a two body SchrSdinger equation for small particle distances was discussed by N. Austern and P. Iano 12). But their "effective range approximation" does not seem to give good results in case of an ~particle. One would expect to get some kind of "effective range approximation" in our case if one puts Vt from eq. (31) approximately equal to V from (28) (or even smaller because of the singlet part of the interaction). This would give a healing distance of 1.45 fro. But the variation gives a healing distance of 0.91 fm and Vt = 950 MeV. The reason for this is, as discussed by H. J. Mang and W. Wild 11), the space dependence of the averaged potential of the 4 nucleons and the large energy gap between ground state and first excited state of the ~-particle. The experimental value of the ~-particle binding energy is around 28.3 MeV and a variational calculation with accurate forces should give a lower limit for this. But the forces (8a) together with (27), (27a) and (28) and the cut off function/1 give already a bigger binding energy. The rather big difference in the results (37) and (38) indicates that a still considerable larger binding energy should be expected from a refinement of the trial function. For definite conclusions about the nuclear forces, further calculations with refined forces and wave functions will be necessary. The shape and radius of the ~-particle have been determined by electron scattering experiments, which suggest a Gaussian charge distribution with a rms-radius of 1.61i0.05 fro. Since the proton is known to have a finite size with an rms-radius of 0.72 ~ 0.05 fm for an assumed Gaussian charge distribution, the ~-particle corresponds to a Gaussian shape with an rms-radius of 1.44+0.07 fm. The rms-radius 1.24 fm obtained by the variational calculation is therefore too small. This m a y be due partly to the fact that the short tail of the Gaussian trial function becomes especially important in the case of a hard core interaction, which increases the probability of finding a nucleon near the surface of the nucleus. For an exponential trial function, which has a longer tail, Cohen and Willis obtained a rms-radius of 2.92 fro.
N U C L E A R M A N Y B O D Y PROBLEMS
97
I wish to thank Dr. K. Wildermuth /or many discussions and Mr. R. R. Coveyou for his advice in setting up the Monte Carlo Method. I thank Dr. J. Lanutti for his interest and for correcting the manuscript. I wish to thank also Mr. D. K. Cavin, Mr. H. T. Gaines and Mr. C. A. Johnson for their kind assistance in carrying out the calculations. References 1) H. Kahn, Rand Comp., AECU-3259 2) K. Wildermuth and Th. Kanellopoulus, Nuclear Physics 9 (1958) 449; CERN Report 59-23 (1959) 3) D. Pearlstein, Y. C. Tang and K. Wildermuth, Phys. Rev. (to be published) 4) K. Wildermuth, Festschrift zum 60. Geburtstag yon W. Heisenberg (to be published) 5) R. Jastrow, Phys. Rev. $7 (1952) 209 6) G. H. Derrick and J. M. Blatt, Nuclear Physics 17 (1960) 67 7) L. Cohen and J. B. Willis in Nuclear forces and the few-nucleon problem Vol. I I (Pergamon Press, London, 1960) p. 399; Nuclear Physics 13 (1959) 125 8) T. Kikuta, M. Morita and M. Yamanda, Progr. Theor. Phys. 15 (1956) 222 9) J. M. Blatt and J. D. Jacksoh, Phys. Rev. 76 (1949) 18 10) E. W. Schmid and K. Wildermuth, Nuclear Physics (1961) (to be published) 11) H. Mang and W. Wild, Zeit. I. Phys. 154 (1959) 182 12) N. Austern and P. Iano, Nuclear Physics 18 (1960) 612
A M O N T E CARLO M E T H O D FOR NUCLEAR M A N Y B O D Y PROBLEMS E R I C H W. S C H M I D Florida State University, Tallahassee, Florida and
Oak Ridge National Laboratory, Oak Ridge, Tennessee t R e c e i v e d 7 A u g u s t 1961 A M o n t e Carlo m e t h o d is g i v e n to solve n u c l e a r m a n y b o d y p r o b l e m s in t h e f r a m e of t h e R i t z ' v a r i a t i o n a l principle. T h e a n t i s y m m e t r i z a t i o n a n d all i n t e g r a t i o n s a r e p e r f o r m e d b y a r a n d o m process. T h e m e t h o d is f o r m u l a t e d for a t w o c l u s t e r s y s t e m , u s i n g oscillator w a v e f u n c t i o n s a n d h a r d core forces, b u t it c a n be generalized to r a t h e r a r b i t r a r y forces a n d w a v e f u n c t i o n s . As a t e s t e x a m p l e , t h e g r o u n d s t a t e e n e r g y a n d r m s - r a d i u s of t h e ~-particle h a s b e e n calculated. F o r t h e final p a r a m e t e r values, a n a c c u r a c y in t h e e n e r g y of 10 % (2 %) w a s o b t a i n e d a f t e r a c a l c u l a t i o n t i m e of 40 sec (16 rain) w i t h a n I B M 7090 c o m p u t e r . T h e a c c u r a c y in t h e r a d i u s is by a factor 8 b e t t e r . T h e m i n i m u m e n e r g y d e p e n d s r a t h e r s t r o n g l y o n t h e f o r m of t h e t r i a l w a v e f u n c t i o n n e a r t h e h a r d core.
Abstract:
1. Introduction T i l e Ritz' variational principle, which is usually employed to calculate properties of light nuclei, demands that one minimizes the expression
E
=
f 9*HAgdT f 9*A9dz
"
(1)
with respect to parameters in the trial wave function 9. In this equation, H is the Hamiltonian of the system and A denotes the complete antisymmetrization of the wave function. We assume H to be symmetric with respect to an exchange of any two particles (that means A can be exchanged with H), therefore the antisymmetrization of 9* can be dropped. Recent forms for the phenomenological interaction between two nucleons all have a repulsive core of small radius which is also considered to prevent a nucleus from collapsing. The hard core imposes rather complicated conditions on the wave function % since 9 has to be zero within and at the surface of the core. These boundary conditions, together with the antisymmetrization make an analytical evaluation of the expression (1) so complicated and lengthy, that it is extremely difficult to apply it for nuclei more massive than a triton or ~-particle. t \Vork s u p p o r t e d in p a r t b y t h e Office of N a v a l R e s e a r c h . C o m p u t i n g t i m e on t h e I B M 7090 c o m p u t e r a t E g l i n A F B f u r n i s h e d w i t h o u t cost b y A F O S R t h r o u g h c o n t r a c t A F 49(638)-278. 82
NUCLEAR MANY BODY PROBLEMS
83
Therefore, a Monte Carlo method has been developed which allows one to calculate E from eq (1) without any analytical calculation. The integrations as well as the antisymmetrization, are performed b y a random process. The basic idea of the method is the following. The integrands ~ * H A ~ ( : ,¢*AH~v) and ~*A~v are written in the form p ( r l . . , r N ) w ( r l . . , r~), where p is called the probability function and contains essentially ~* A (for details see next paragraph), w is called the weight function and contains essentially H~v or ~, respectively. Random co-ordinates are generated with the probability distribution p ( r ~ . . , rN) and a certain number of terms from the sum represented b y A are selected at random according to a uniform distribution over all N! exchange terms. For these random numbers, together with the exchange orders (which also contain a factor + 1 according to even or odd permutations), the weight functions are calculated. The values obtained in this way are estimates for the integrals S ~v*HA~vdz and S ~v*A~vdr and the averages of a large number of these estimates converge toward the value of the integrals. In principle, any kind of variational function and any kind of interaction force can be used in this method, but the calculation time and also the length of the computer program imposes certain restrictions. Test calculations with the ,t-particle have shown that the method will be useful mainly for the following problems: 1) Three and four body problem. Rather elaborate variational functions and interaction forces can be used to get information on nuclear forces in addition to the information following from the two-nucleon scattering data. 2) Light nuclei up to carbon or oxygen. Some interesting properties of the nuclei (energy levels, radii, 15- and y-transition probabilities, clustering effects) can be investigated, using rather simple central forces of Serber type with a hard core and cluster model wave functions with few, but effective, variational parameters. Emphasis in this work is put on the second problem. The Monte Carlo method will be discussed in sect. 2. Sect. 3 tells something about the computing technique and sect. 4 gives the result of a test calculation with the ~-particle. 2. Method
An integral I = f p ( x l . . . X3N)W ( X l . . . x3N)dxl.., dX3g
(2)
can be evaluated in the following way 1). One generates random co-ordinates x~i, x2 ~..... x~N according to the normalized probability distribution p(xl... XaN). For every set of 3N random co-ordinates the value of the weight function W(Xl. . . X3N) ;~ calculated: Wi . . W. ( X. l i i x3:~). (3)
84
ERICH
W.
SCHMID
The value w~ is an estimate for the integral I, and the average of the estimates converges under a certain condition to the integral: w i
--->I.
(4)
The condition of convergence is, that the variance v, defined as v = f p (X)l.
• . Z3N[W
(X1 .
. . X3N ) - 1 ]
2d x l . . . d x s N
(5)
is finite. The variance gives also the statistical error A I of an average of n estimates (standard deviation)" (6)
AI = V~/n.
From this relation it can be seen already that the error will decrease with the square root of the computing time. The integration method discussed above is called a "Monte Carlo Method". We want to use it to carry out the integrations in a Ritz' variational calculation for light nuclei. In order to give a simple description of the method, we shall restrict ourselves in the following to a two-cluster system *) described by oscillator functions and central exchange forces with a hard core. The trial wave function has the form (7)
W = 919~z H l(r,k)S, i>1¢
where 91 and 92 describe the internal motion of the two clusters and Z describes the relative motion, while S is the spin and isobaric spin function. The functions 91, 95 and Z fulfill the boundary condition at infinity but not in the range of the hard core. The "cut off" factor I I , > k/(r~k) is used to give the right boundary f(r)
r¢
d
F i g . 1. T h e " c u t o f f " f u n c t i o n
1"
/(r).
t I n p r i n c i p l e , e v e r y w a v e f u n c t i o n d e s c r i b i n g a n u c l e u s c a n b e w r i t t e n a s a s u p e r p o s i t i o n of t w o c l u s t e r w a v e f u n c t i o n s z).
NUCLEAR MANY BODY PROBLEMS
85
condition at the surface of the hard core ¢. The function /(r) has the form shown in fig. 1. For values of r smaller than the range of the hard core, ](r) is zero and [(r) is equal to 1, for large values of r. The smallest value of r, for which /(r) is equal to 1 defines the "healing distance" d, that is the separation distance of two nucleons beyond which the hard core no longer influences the wave function b y its boundary condition. The Hamilton operator of the system has the form ?~2 S -- 2M
X
aU particles
v,2+
X
I
e~ ekX
all pairs
(s)
where
V,k ~ V(r,k)Ew+mP*(i, k)+bP~(i, k)--hP'(i, k)l
(da)
is the two body potential and P~(i, k), P*(i, k) and P*(i, k) are the operators which exchange the space, spin and isobaric spin co-ordinates of particle i and k. The numbers e~ are equal to the proton charge if particle i is a proton and zero if it is a neutron. We assume that the wave function ~v does not contain any motion of the total centre of mass; therefore, we can take the centre-ofmass as origin of our co-ordinate system. 3 (N-- 1) co-ordinates (N is the number of particles) are then independent, while 3 co-ordinates are determined b y the condition N
X r, = O.
(9)
The integrands ~v*HA~v and ~p*A~o of eq (1) have to be written now in the form p ( r l . . . r n ) w ( r l . . , rn). We use the same probability function p for both integrands, therefore we do not have to care about the normalization of p. In principle p ( r l . . . rn) could be taken to be different from zero and constant within a rectangular box in all 3N co-ordinates which encloses the integrand at all points where it is essentially different from zero. The function w ( r l . . . rN) would then be the integrand itself (besides a normalization factor). But this would give a large variance, since in a many-dimensional box the main part of the volume consists of "corners" of the box, where the integrand is close to zero. It can be seen from eq. (5) that the variance is small when p represents the main part of the integrand and w is a function which varies as little as possible**. This requirement is most significant in the region where the integrand is falling off to zero, since this region represents the main part of tile integration volume. ? The " c u t off" factor II~>~/(r~) does n o t change the geometrical properties of the wave function; L, S and the p a r i t y r e m a i n good q u a n t u m n u m b e r s . Therefore, the splitting-off of the centre-of-mass m o t i o n is no p r o b l e m for w a v e functions of the f o r m (7), unlike m o s t other m a n y particle w a v e functions (see also refs. 3-~)). it The variance would be zero if p is the integrand and w i s a c o n s t a n t n o r m a l i z a t i o n factor which is in this case the integral itself.
86
ERICH W. SCHMID
The possibilities of choosing p ( r l . . , rN) are restricted by the difficulty of generating random numbers. For a distribution function depending only on one coordinate, the generation of random numbers is easy. A transformation function can be obtained by a one-dimensional numerical integration which transforms random numbers of a uniform distribution into random numbers of the desired distribution (for details see sect. 3). For some distributions like exponential and Gaussian distributions, library programs exist, which generate the random numbers in a more direct way. For a general distribution function depending on m a n y co-ordinates the generation of random numbers is difficult *. We demand therefore, t h a t p is either a product of functions each of which depends only on one co-ordinate or that p is a linear combination a P l + b p ~ + ••• of such products of functions. In the first case every factor for itself is used to generate a random co-ordinate. In the latter case the normalized distribution functions Pl, P, . . . . are used alternatively according to their weights a, b. . . . . As trial wave functions for the internal and relative motion of our clusters we adopt the usual oscillator wave functions, which have the form given below. An unexcited m-particle cluster is described by ~(~) = exp [--~1 a (u 2 + v 2 + w 2 )],
(10)
where r 1 = {(--u--v--w)+R~, r 3 =
r 2 = }(u--v+w)+R~,
½(u+v--w)+R~,
r4 :
(11)
½(--u+v+w)+R~
are the position vectors of the 4 nucleons, R~ is the centre-of-mass vector and a is a variational parameter. A triton trial wave function is described by ~(tr) = exp [--½a 3 2 +½v
)l,
(12)
with rl = ½(u+v)+Rtr,
r2 = ½ ( u - - v ) + R t r ,
r3 = - - u + R t r .
(13)
For a deuteron cluster, a Gaussian function is a good description only when the cluster is strongly bound, If it is loosely bound, an exponential function or a superposition of different Gaussian functions is better: ~o(d) = exp(--½au 2)
~(d) = exp(--~a[u[),
or
(14)
with r I =
1 -flU +
Rd,
r2 = - - l u
+ Rd.
(15)
The relative motion X has the form Z =
R"exp(--aR2)Y~,o(
v~, 9 )
* In order to obtain them by transformation from a uniform distribution, transformation w h o s e J a c o b i a u d e t e r m i n a n t is a g i v e n f u n c t i o n .
(16) one has to find a
NUCLEAR
MANY
BODY
87
PROBLEMS
with R 1 -- -
M S -
R,
MI+M2
Rz -- -
--M1 -
R,
MI+M2
(17)
where M 1, M 2 are the masses of cluster 1, 2). Eq. (17) states, that the centreof-mass of the two-cluster system is at the origin. The functions (10), (12), (14) and (16) are all of the type described above, that is, they are products of one-dimensional functions. This is also true for the oscillator wave functions of bigger clusters, since they can be built up of smaller clusters in the above way. The factor ~ 2 ~ g is contained in ~* as well as in ~. The factor in ~2can be split into a factor (91~2X)s which is symmetric with respect to any exchange of particles and a factor which is not symmetric. As probability distribution p we can take then, if we forget for a moment about the statistical antisymmetrization, P = ~l~$Z(~192Z)8.
(18)
The star coming from ~* has been dropped since the above mentioned functions are real. The weight functions for numerator and the denominator of eq. (1) are 1 wl=~*HA %
1 w2-~/*A %
(19)
respectively. The random co-ordinates generated according to eq. (18) are transformed into cartesian co-ordinates by the transformations (11), (la), (15) and (17). The quantities w1 and w2 are calculated numerically in cartesian co-ordinates. For more complicated trial wave functions than the ones described above (as for instance trial functions which include higher oscillator shells), we can use the same distribution function as above, with parameter values which give about the same or a slightly larger rms-radius as would arise from the trial wave function. Eqs. (19) are also valid in this case. The exact form of the probability function p is not very important, since it has only an influence on the variance. In the method described in this paper the main contribution to the variance comes from the Hamiltonian in the region where two or more particles are close together. A somewhat different probability function would therefore increase the variance only little, while the computing time as well as the preparation time will be smaller when a simple probability function is used. The weight functions (19) are still inconvenient because of the antisymmetrization operator A. This operator splits the weight functions into a sum of "exchange terms". Many of these terms are usually zero because of the orthogonality of the spin functions. Some of them may be equal to each other. But there is usually still a large number of different terms. We want, therefore, to treat the antisymmetrization by a random process.
88
E R I C H W. SCHMID
The antisymmetrization operator A is defined by A ---- X ± P ( 1 , 2 . . . . . N).
(20)
The sum extends over all permutations of the N nucleons, the sign being taken according to whether the permutation is even or odd. Instead of operating with the whole sum of operators on the wave function ~ we take only a certain number m of them (together with the right sign), which we select at random. The m permutation operators could be selected from a uniform distribution over all N! operators. But in this case most of the operators would give zero contributions to the estimates of our integrals. In order to find the permutation operators which give non-zero contributions we have to look at the spin and isobaric spin function S. Since the trial wave function m a y be chosen to have any symmetry, S can be chosen in a very simple form. Usually it is taken as one single product of single-nucleon spin and isobaric spin functions, or a linear combination of few of such products. We assume in the following, that it has the form of one product. Then the only permutations which give a non-zero contribution to the normalization, the kinetic energy and the part of the potential energy which does not involve spin and isobaric spin exchanges, are the ones for which only particles in the same spin and isobaric spin state are exchanged. From a uniform distribution of these permutation operators, we select at random m operators. For the calculation of the Bartlett and Heisenberg part of the potential energy, where we have spin and isobaric spin exchanges, we remember the following: In eq. (1) we are allowed to carry out any permutation of particles in ~v (and change sign if the permutation is odd) and still get the same result because of the operator A. So we use for the calculation of the Bartlett term instead of the wave function ~v the function --P(i, k)~v. If the spin exchange operator of the Bartlett term is applied to this, we get --P*(i, k)P(i, k) = --P*(i, k)er(i, k). That means, in case of different spins of particle i and k, we can replace in the calculation the spin exchange operator P~(i, k) by zero if particles i and k have different isobaric spins and by the operator --P~(i, k) if the isobaric spin is the same. If the particles i and k have the same spin the exchange operator is simply replaced by + 1 . The Heisenberg part of the potential energy is treated in a similar way. We can now write instead of eqs. (18) and (19) in a symbolic way p = (p1(P2%(Vi(p2%)s A : 1
Wl.~y)*Sv) ,
~b'A, 1
W2 -~W*~V ,
(21)
(22)
with the understanding, that A in (21) means statistical antisymmetrization and is applied only on %
NUCLEAR MANY BODY PROBLEMS
89
The statistical antisymmetrization will certainly increase the variance v. H o w many permutations one has to take for every estimate has to be found out b y a test calculation. For small nuclei, or nuclei of high symmetry like Be 8 it is probably better to use eqs. (18) and (19). When the parameter values have been found, which give the minimum energy, the Monte Carlo method can be applied in a similar w a y as for the Hamiltonian to calculate the expectation values of other operators, like the mean square radius, fl- and 7-transition probabilities.
3. Computing Technique For Monte Carlo calculations of the kind described in the last section a fast electronic computer with a large core storage is necessary, although the calculations are rather simple. But in order to get fairly accurate results, a large number of estimates has to be calculated and in general the application of the method is limited b y the length of the computing time. Also the questions arising while writing the computer programs are mainly concerned with the speed of the computation. In the following, some of these questions will be discussed. Random numbers of certain probability distributions like uniform distribution, exponential, and Gaussian distribution are obtained b y subroutines which are distributed b y the S H A R E organization. Time can be saved b y calculating a set of random numbers only once and storing it on tape. Random numbers of other one dimensional distribution functions can be obtained from a uniform distribution b y a transformation. The transformation, which transforms a uniform distribution of x in the interval (0, 1) into a normalized distribution p(y) is obtained b y a numerical integration of the differential equation dy 1 d~ = p(y~)"
(23)
The integration has to give a table of y = y(x), which is stored in the core storage of the computer. Since the table has to be narrowly spaced in order to spare interpolation, a simple integration rule can be used, as for instance h - , Yn+l--Yn _- - -p(~)
(24)
where h is the spacing in x and ~ is a value between y~ and Yn+1 which can be taken to be y,~+½ (Y~--Yn--1) in regions where Y~--Yn-1 and Y,~+I--Yn are not too different ¢. ¢ This c o n d i t i o n is for i n s t a n c e n o t fulfilled in t h e case p(y) ~ c y 4 e x p ( - - a y 4) n e a r y = O. T h u s for s m a l l v a l u e s of y i t is b e t t e r to i n t e g r a t e a n a l y t i c a l l y w h i l e p u t t i n g t h e G a u s s i a n f u n c t i o n e q u a l to 1.
90
ERICH W. SCHMID
Random numbers of a distribution function I5 = a p l + b p 2 +
...
are generated using a fork, which causes the program to generate numbers either from Pl or P2 and so on, according to the weight factors a, b. . . . The generation of random permutation operators can be done in the following way, when the function S is a single product of nucleon spin and isobaric spin functions. The N nucleons are divided into 4 groups of particles, namely protons with spin up, protons with spin down, neutrons with spin up and neutrons with spin down. As discussed in the last section only exchanges among the particles of each group need to be considered and we can therefore treat each group separately. We define two N-dimensional vectors a and b. The vector a has as its components a, the numbers of our N particles in an order, which is arbitrary except that the first N 1 components are the numbers of protons with spin up, the next N~ particles are protons with spin down, and so on. The vector b will contain the permuted numbers. We start by choosing at random one out of the first N1 components of b. This can be done by generating a random number from a uniform distribution between 0 and 1, multiplying it by N1, adding + 1 and truncating it to an integer. The result tells whether the first, second or N1 t~ of the components in question has to be taken. This component we put equal to al. From the still undefined N 1 - 1 components of the first group we choose at random another one and put it equal to a 2. When all N 1 components of the first group are defined, the procedure continues with the next group and so on. Finally b will contain a permutation of the components of a. Every permutation which will give a non zero contribution to our integrals has the same probability in the above generation process. Whether a permutation is even or odd can be found out in the usual way by counting the exchanges which are necessary to get a from b. For the calculation of the kinetic energy, the second partial derivatives of the wave function ~v with respect to the nucleon co-ordinates have to be calculated numerically. A second partial derivative is defined as 32~p _ lim yj(.. " x j - - e . . . ) - - 2 ~ 0 ( . . , ~X t 2
~-*0
xs...)+v2(..,
xj+e...)
(25)
8"2
We use this equation with a proper value of e, but instead of calculating ~v ( . . . 2 j 4 - E . . . ) directly, the following procedure is used: We expand everyone ofthe factors ~1, 92, Z,/(r,~) which build up ~v according to (7) into a Taylor series in e, which we cut off after the e2-term. These Taylor series we multiply together in order to get ~v(... x~-¢-e...). Now, for a computer with 8 digit numbers, we can choose e so, that if the ,0 term starts in the first digit, the el, e~, E8 and e4 terms will start approximately in the 3 ra, 5 th, 7th and 9 th digit respec-
91
NUCLEAR MANY BODY PROBLEMS
tively, that means the # term is truncated b y the machine. Inserting the thus obtained ~0(... x ~ 4 - e . . . ) into eq. (25) we see that the e° term and all odd terms cancel, and tire second partial derivative is obtained with a four place accuracy. If all coefficients of the Taylor series are tabulated and kept in the core storage, the calculation of each Taylor series takes only about two multiplications and two additions. Nevertheless, because of the 3N derivatives and the N - - 1 factors ](rik ) which contribute to every derivative, the calculation of the kinetic energy estimates is the longest part of the calculation. In the calculation of the potential energy estimates, the Wigner part presents no difficulties. For the Majorana Bartlett and Heisenberg terms twoparticle exchanges have to be calculated. When oscillator wave functions of the type (10), (12), (14) and (16) are used, the wave function can be split into a factor which is symmetric and a factor which is not symmetric with respect to particle exchanges. Only the latter one is affected b y the exchange operations. In the computer we write the spin function S (if it is a single product) as an Ndimensional vector of numbers si. The value s~ = 1 means particle i is a proton with spin up, s i = 2 proton with spin down, s, ~ 3 neutron with spin up and s~ = 4 neutron with spin down. The sum and the difference of s t and s, is enough to decide, which operator has to be applied on ~pin order to get a nonzero contribution to the potential energy estimate. Table 1 gives these operators for every possible combination of particles. TABLE 1 Potential energy operators for various combinations of particles i and k s~+s~
Operator
s~
sk
s,--s~,
1
1
0
2
V(r~) ( w + b - - h + m P ~ k ) + - -
2
2
0
4
g(r~) (w+b--h+rnI~i~) + r ~
3
3
0
6
V(r,,) ( w + b - - h + m P ~ , )
4
4
0
8
V(r~,) ( w + b - - h + m P ~ , )
1
2
:#0
3
V (r,,) [w -- h + (m-- b) P~,] + - -
1
3
:#:0
4
V(r,~) [ w + b + (m +h)P~]
1
4
~:0
5
V (r~,) (w + toPOi)
2
3
v~:O
5
V(r,,) (w +raP;,)
2
4
:/:0
6
V(r,,) [ w + b + (m+h)P~,]
3
4
:/:0
7
V(r~,) [w--h:~- (m-- b) P~,]
82
e2
$2
Yik
92
:ERICH W. SCHMID
The normalization estimates arise as intermediate results during the potential energy calculation. For the speed of the calculation it is important that all functions which are used are calculated and stored as tables in the core storage before the Monte Carlo calculation starts. The tables have to be narrowly spaced so that no interpolation is necessary. In this way also the large storage capacity of modern computers can be utilized to gain speed. The calculation time can be estimated by counting the multiplications, divisions, additions and subtractions. All other operations do not count in a first approximation. In order to calculate the statistical errors, we average not only the estimates w,, but also their squares w, 2. The variance is then calculated approximately from v =--
T~ 5=1
w~2-- n
4=1
wi ,
(26)
which follows from eq. (5) if one inserts (1/n)~ w~ as an approximation for I. The error is obtained by eq. (6). When adding up a large number of estimates, it might happen, that the smaller ones get lost by truncation, if the sum is already very large. This can be avoided by calculating intermediate averages after every 104 or 10~ steps.
4. E x a m p l e Monte Carlo methods have been used by G. Derrick and J. M. Blatt 6) to calculate the binding energy of the triton and by L. Cohen and J. B. Willis ~) to calculate the binding energy and rms-radius of the a-particle. Cohen and Willis report for the a-particle calculation with hard core forces a possible error of 10 % in the energy after a calculation time of 40 h (with their final values for the variational parameters). As a test of the method discussed in this paper, the binding energy and rmsradius of the a-particle has been calculated. The calculations should also show whether the shape of the cut-off function [(r) has a strong influence on the result of the variational calculation or not. The two-nucleon potential was assumed to be of the form (Sa) with
V(r) = t + ~ t V exp[--fl(r--rc) ]
if if
r <
re,
Y ~
YC,
(27)
and the Serber relation m =
w =
¼(l+x),
h =
b =
¼(l--x),
('27a)
where x is the ratio of the singlet and triplet well depth. The hard core radius
NUCLEAR
MANY
BODY
03
PROBLEMS
r e was assumed to be re = 0.4 fm. The constants V, fl and x were calculated with formulas given by T. Kikuta, M. Morita and M. Yamanda 8) and T. M. Blatt and J. D. Jackson 9) in order to fit the two-nucleon data: n-p triplet scattering length a t =
5.378 fm,
n-p triplet effective range
1.704 fro,
rot =
n-p singlet scattering length a s ----- --23.69
fm,
n-p singlet effective range
fm,
ros-~
2.4
For the n-p singlet effective range Kikuta et al. used the two values 2.4 fm and 2.7 fm because of the large experimental error of this quantity. In this paper the value 2.4 fm has been preferred since in a previous paper 10) the singlet force constants have been fitted from the more accurate p-p scattering data and a calculated n-p effective range of 2.41 fm was obtained with these constants. Nevertheless, the uncertainty of the n-p singlet effective range was considered large enough to justify the use of the same well parameter/3 for both the singlet and triplet interaction ¢. The following values for the potential parameters were obtained: V -~ 498.02 MeV,
/3 = 2.5885 fm -1,
x = 0.771.
(28)
A numerical integration showed that these constants also give the right deuteron binding energy within the experimental error. The trial wave function was assumed to be
(1) (1)(o)(1)(1)(o)(o)(o)
~p = e x p ( - - l a ~, ri~) I I ](ri~) 0 i ~
i~>k
0
O"l
1
T1
0
0"2
0
72
1
0" 3
1
T3
0"4
(29)
1 'r 4 "
The Gaussian function is the same as in eq. (10). The quantities
(10k(10t are the usual Pauli spin-and isobaric spin vectors; a is a variational parameter. For the cut-off function ](r), two different trial functions ]l(r) and ]~(r) have been used. Both functions are zero within the range of the hard core and both are 1 for values of r larger than the healing distance d. Between these two points, the functions have the following form: t I n t h e c a l c u l a t i o n of t h e force c o n s t a n t s t h e singlet i n t r i n s i c r a n g e w a s c h a n g e d a little in o r d e r to give t h e s a m e / ~ as in t h e triple case. Since t h e well d e p t h p a r a m e t e r is close to 1 in t h i s case, t h e s c a t t e r i n g l e n g t h is n o t i n f l u e n c e d b y t h i s procedure.
~4
ERI CH W. SCHMID
1) The function/1 is a numerical solution of the two-nucleon Schr6dmger equation [
~ dr d~2 + Vt(r) 1 [l(r) = 0 M
(30)
(re < r < d).
The first maximum of the solution/l(r) is normalized to 1 and defines the healing distance d. The potential Vt(r ) which "generates" the trial function t1(') w a s Vt(, ) = --V t e x p f - - f l ( r - - r e ) ].
(31)
The parameter fl was the same as in eq. (28) but the well depth V t was used as a variational parameter. The healing distance d and the well deptL V t depend on each other, therefore we can also call d the variational parameter. 2) The function/2(r) had the form of the first quadrant of a sine wave: /2(r) = sin ~ ( r - r e ) 2(d--re)
(re < r < d),
(32)
where d was considered to be a variational parameter. Because of the high symmetry of the wave function (29), the Bartlett-and Heisenberg part of the interaction cancels. In all other parts of the interaction as well as in tile kinetic energy part and the normalization, the antisymmetrization operator can be dropped since only the direct term gives a non zero contribution. This example, therefore, can give no test for the statistical antisymmetrization. For the probability distribution tile square of expression (10) was taken according to eq. (18): p = exp[--a(u2+v2+w~)].
(33)
The random co-ordinates u, v, and w were transformed into cartesian co-ordinates by eq. (11). The weight functions were (with either/1 or/2 instead of/) according to eq. (19) _he wl = exp(+-~a X r~k) 1-I / ( r , k ) ~ • V, 2 i>k
i>k
j
e x p ( - - } a X r~k) H ~if,k) + H i>k
i~k
/(r,k) 2
(w+m) X V(r,k)+ r
i~k
=
II/('
'
i~k
kV.
(35)
The weight function wa for the mean square radius was w~ = w2rl 2.
(36)
The computer program was written in FORTRAN language and had some
NUCLEAR
MANY
BODY
PROBLEMS
95
170, mostly short, statements. The rather large number of statements was necessary to optimize the speed of the calculation ¢. The computing time was around 4 mill for every 10000 estimates with an IBM 7090 computer (the tabulation of functions takes only a few seconds). TABLE
2
E n e r g y a n d r m s - r a d i u s for c u t - o f f f u n c t i o n ]t a n d different p a r a m e t e r v a l u e s E (MeV)
r m s (fro)
--23.0+0.5 --25.9±0.5 --27.8±0.5 - - 29.3 ± 0.4 --27.8±0.6 --24.7±0.7 --21.4=[_0.7 --25.7 ~ 0 . 7 --29.4±0.4 -- 28.2 J::O.4 -- 1 5 . 6 ± 0 . 6
a (fm -~)
1.53 1.41 1.31 1.23 1.17 1.12 1.07 1.23 1.24 1.24 1.27
0.5 0.6 0.7 0.8 0.9 1.0 1.1 0.8 0.8 0.8 0.8
TABLE
Vt(MeV); d(fm) 1000; 1000; 1000; 1000; 1000; 1000; 1000; 1200; 950; 800; 500;
0.887 0.887 0.887 0.887 0.887 0.887 0.887 0.821 0.908 0.992 1.45
3
E n e r g y a n d r m s - r a d i u s for c u t off f u n c t i o n /2 a n d different p a r a m e t e r v a l u e s E (MeV)
r m s (fro)
a ( f m -~)
d (fm)
1.22 1.23 1.24 1.25 1.40 1.31 1.28 1.24 1.12
0.8 0.8 0.8 0.8 0.6 0.7 0.74 0.8 1.0
0.757 0.821 0.887 0.992 0.86 0.86 0.86 0.86 0.86
--22.8±1.0 --25.3±0.9 --25.4~0.7 21.6±0.6 --24.8±0.6 --25.8±0.5 -- 25.9::k0.4 --25.6i0.5 --21.1 ± 0 . 9
Tables 2 and 3 show the dependence of the energy and rms-radius on the parameter values for the two different cut-off functions. The minimum values are Emt n =
- - 2 9 . 4 i 0 . 4 MeV,
(rms) = 1.24 fro,
d = 0.91 fm
(37)
(rms) = 1.28 fm,
d = 0.86 fm
(38)
for the first case with /1, and Emt n =
--25.9±0.4 MeV,
for the cut off function /2. * A factor 1.5 in speed could p r o b a b l y be gained b y w r i t i n g t h e p r o g r a m in m a c h i n e l a n g u a g e , b u t t h i s w o u l d increase t h e p r e p a r a t i o n t i m e b y a n u n b e a r a b l e a m o u n t .
96
ERICH
W.
SCHMID
The function ]1 gives a deeper minimum and a smaller error. The statistical error in the energy near the minimum is 10 ~ (2 ~ ) after a calculation time of 40 sec (16 rain) with [1 or after a calculation time of 100 sec (40 rain) with ]2. The error increases for smaller healing distances. The accuracy in the normalization and the radius is better by a factor 8. The error in the energy arises mainly from the region between hard core and healing distance, since the energies there are very large. The function ]1 makes a better balance between potential and kinetic energy at every space point in this region and gives, therefore, a lower variance. The result of the variational calculation tells that the form of the short ranged correlation function of the nucleons is quite important. The concept of using a strict solution of a two body SchrSdinger equation for small particle distances was discussed by N. Austern and P. Iano 12). But their "effective range approximation" does not seem to give good results in case of an ~particle. One would expect to get some kind of "effective range approximation" in our case if one puts Vt from eq. (31) approximately equal to V from (28) (or even smaller because of the singlet part of the interaction). This would give a healing distance of 1.45 fro. But the variation gives a healing distance of 0.91 fm and Vt = 950 MeV. The reason for this is, as discussed by H. J. Mang and W. Wild 11), the space dependence of the averaged potential of the 4 nucleons and the large energy gap between ground state and first excited state of the ~-particle. The experimental value of the ~-particle binding energy is around 28.3 MeV and a variational calculation with accurate forces should give a lower limit for this. But the forces (8a) together with (27), (27a) and (28) and the cut off function/1 give already a bigger binding energy. The rather big difference in the results (37) and (38) indicates that a still considerable larger binding energy should be expected from a refinement of the trial function. For definite conclusions about the nuclear forces, further calculations with refined forces and wave functions will be necessary. The shape and radius of the ~-particle have been determined by electron scattering experiments, which suggest a Gaussian charge distribution with a rms-radius of 1.61i0.05 fro. Since the proton is known to have a finite size with an rms-radius of 0.72 ~ 0.05 fm for an assumed Gaussian charge distribution, the ~-particle corresponds to a Gaussian shape with an rms-radius of 1.44+0.07 fm. The rms-radius 1.24 fm obtained by the variational calculation is therefore too small. This m a y be due partly to the fact that the short tail of the Gaussian trial function becomes especially important in the case of a hard core interaction, which increases the probability of finding a nucleon near the surface of the nucleus. For an exponential trial function, which has a longer tail, Cohen and Willis obtained a rms-radius of 2.92 fro.
N U C L E A R M A N Y B O D Y PROBLEMS
97
I wish to thank Dr. K. Wildermuth /or many discussions and Mr. R. R. Coveyou for his advice in setting up the Monte Carlo Method. I thank Dr. J. Lanutti for his interest and for correcting the manuscript. I wish to thank also Mr. D. K. Cavin, Mr. H. T. Gaines and Mr. C. A. Johnson for their kind assistance in carrying out the calculations. References 1) H. Kahn, Rand Comp., AECU-3259 2) K. Wildermuth and Th. Kanellopoulus, Nuclear Physics 9 (1958) 449; CERN Report 59-23 (1959) 3) D. Pearlstein, Y. C. Tang and K. Wildermuth, Phys. Rev. (to be published) 4) K. Wildermuth, Festschrift zum 60. Geburtstag yon W. Heisenberg (to be published) 5) R. Jastrow, Phys. Rev. $7 (1952) 209 6) G. H. Derrick and J. M. Blatt, Nuclear Physics 17 (1960) 67 7) L. Cohen and J. B. Willis in Nuclear forces and the few-nucleon problem Vol. I I (Pergamon Press, London, 1960) p. 399; Nuclear Physics 13 (1959) 125 8) T. Kikuta, M. Morita and M. Yamanda, Progr. Theor. Phys. 15 (1956) 222 9) J. M. Blatt and J. D. Jacksoh, Phys. Rev. 76 (1949) 18 10) E. W. Schmid and K. Wildermuth, Nuclear Physics (1961) (to be published) 11) H. Mang and W. Wild, Zeit. I. Phys. 154 (1959) 182 12) N. Austern and P. Iano, Nuclear Physics 18 (1960) 612