LIBRARY OF THE UNIVERSITY OF ILLINOIS AT URBANA-CHAMPAIGN 510.84 ho.445-450 The person charging this material is re- sponsible for its return to the library from which it was withdrawn on or before the Latest Date stamped below. Theft, mutilation, and underlining of books are reasons for disciplinary action and may result in dismissal from the University. UNIVERSITY OF ILLINOIS LIBRARY AT URBANA-CHAMPAIGN FIB 1 9 !3 7 FES 1 9 ', L161 — O-1096 Digitized by the Internet Archive in 2013 http://archive.org/details/pathintegralcalc445jaco ■**p Report No. kk5 PATH INTEGRAL CALCULATION OF THE THREE- PARTICLE DIRECT AND E^Sf CONTRIBUTIONS TO THE PAIR DISTRIBUTION FUNCTION AND OF THE THIRD VIRIAL COEFFICIENT FOR HE^ GAS by Robert Clarence Jacob son June, 1971 This volume Is bound without no. hkG THE LIBRARY C7 THE - b 1973 ssaas *, N -• » in such N a way that — r -* n. In order that the above assumptions be valid, it will L be sufficient that an actual physical gas be confined to a container the dimensions of which are much greater than both the thermal wavelength, , ( 2TTQti 2 \ 1/2 . . _ . . , TT 4 A. = ( — K — I , and the range of the interaction between any two He V "He / atoms. We also assume the gas is in a state of thermal equilibrium and as such can be described by the canonical density operator: -(% e where N ""'i i=l P i + V ~ N) (2 ' l,1) N r = (r,,...,r ) is the 3N- dimensional position operator, p. is the associated 3-dimensional momentum operator for particle i (=1,2,...,N) and N V N (r ) is the interaction potential for the N particles. Eventually we N will assume V„(r ) to be in the form of a superposition of two-body Lennard- Jones potentials, where V(r. .) = 4a ij &) " - ft) (2.1.2) and r . . = r . - r . -16 ° a = 14.04 X 10 erg a = 2.56 A This assumption is not needed until Chapter 4. The values for a and a were chosen so that a comparison could be made with calculations of Jordan [9]. They were obtained by deBoer and Michaels [3], and provided good agreement" with experiment for the second virial coefficient [12]. In all following equations, we will assume units such that 4 ft (Planck's constant/2r0 = irv, (mass of He molecule) = a = 1. In these 1 18 5 units, B = t - ^ = ~z — where k is the Stephan-Boltzmann constant and T is the kT T o 40.97 temperature of the gas in K, y = 4a|3 = — - — and the thermal wavelength, X = /2tt{3 , is 3.41//T. We will generally be interested in the conf igurational matrix elements of the density operator, where If ) S Illf2' ' • • '5n > = |l,2,...,N> is the eigenvector of the position operators r., i = 1,...,N. The vectors are normalized so that ,Ni N N _ 3N, ,N N, / , IN • INv JIN. , IN IN. (r |r ) = 6 (r -r ) Calculations by Boyd, Larsen, and Kilpatrick [12] actually use a/k = 10.22°K in accordance with [3]. The value adopted here and by Jordan yields a/k = 10.17°K, the discrepancy apparently being due to improved measurements of k. The superscript, N, will generally be omitted for convenience expect in cases where it is necessary for clarity. Occasionally, the superscript, op, will be provided to distinguish the operator from its eigenvalue. In thermodynamic equilibrium, the quantum statistical average of any observable, 8, is given by -r-r -1 "^ TeT = Z N L tr(9e 1N ) (2.1.3) - p v . where Z = tr(e ) is the partition function of the system. In the case of identical particles, the trace must be taken over a complete set of 4 properly symmetrized eigenstates. Since He atoms are spin-zero Bosons, the states must be simultaneous eigenstates of the permutation operator, P, with eigenvalue +1. Accordingly we define the symmetrized configurational eigenstates: |l,...,N> h -L. Z P|1,...,N> /NT p where these kets are normalized so that the associated wave functions are normalized to unity. In terms of symmetrized configurational states, the partition function is Z N = nTJ dr 3N {l,...,N|e ^|l,...,N> + (2.1.4) where the N! is necessary so that equivalent states related by a permutation of configurational values will be counted only once. In terms of the un- symmetrized kets we have Z N "hS dr3Nz (1 N|Pe" PHN |l M> ■ ^7 J d3N V£ N > • < 2 -^> X 3N N , J N whe re W (r N ) = X 3N £ (1 n|Pe ^| 1 N> (2.1.6) N- p is known as the Slater sum. The diagonal contribution is the Boltzmann- _Q U Slater sum, lOr ) = X (l,...,N|e |l,...,N>. In the classical limit N B N -P V N ( rN ) Aim W (r ) = Aim W~(r ) = e B - B -» The "generic" conf igurational distribution function for h particles, n, (r ), is defined as the probability density of the occurrence of any h particles in the element dr about the point r = (r. ,r« , . . . ,r, ) . Thus n. (r h ) = (Z n 6 3 (r° P - r, )) (2.1.7) h ~ {L} k=i ~\ A where X represents the sum over all permutations of h particles A J„,. . , J, U k ] l 2 h N! taken out of the total N, there are ... , ^ , such permutations, and each term in the sum is equal. Hence according to equations (2.1.3) and (2.1.5), equation (2.1.7) can be written, n (r h ) = ^r ^ - f d 3 r d 3 r W (r N ) . (2.1.8) h ~ (N-h)!X 3N Z N J h+1 NN ~ In the case h = 2, n, (r ) = n„(r,,r ) is the pair distribution h ~ 2 ~1 ~2 r function. For h = 1 we have the particle density, n(r) . In the limit of large volume, V, and number of particles, N, it is well known that for interactions of finite range, n(r) — n, and n 9 (r.,r 9 ) -• n„ ( | r_-r.. | ) . 2.2 The Modified Cluster Expansion The generic distribution function n(r ) can be expressed as an expansion in powers of the fugacity, z = e p , where u, is the chemical potential. The so-called "modified cluster expansion" due to deBoer [13], an extension of the cluster expansion by Kahn and Uhlenbeck. [14], depends N on a transformation of the Slater sum, W.,(r ), to the U functions defined by: W(r h ) = U(r h ) W(r\r ) = U(r h ,r„) + U(r h )U(r ) W(r h ,r K ,r x ) = U(r\r K ,r x ) + U(r h ,r R )U(r x ) + U(r x ,r K )U(r h ) + U(r h ,r x )U(r K ) + U(r h )U(r K )U(r x ) (2.2.1) where we have dropped the subscripts on the W's for convenience. Note that the U-functions are defined by the number of arguments and not the dimensionality as with W functions. Thus, U(r ,r .) and U(r, ,r„ , . . . ,r, , 1 ) are different functions, e.g. for h = 2: U(r ,r 3 ) = W(r 1 ,r 2 ,r 3 ) - W(r 1 ,r 2 )W(r 3 > U(r 1} r 2 ,r 3 ) = W^.r^^) - W(r 1 ,r 2 )W(r 3 > - w (f 2 , 53 ) W( 5l ) - w (5 1 '5 3 ) w (f 2 ) + 2W( Il )W( 52 )W( 53 ) ' There are thus two types of U-functions implicitly defined above. Those which do not have the r argument are the functions of Kahn and Uhlenbeck. Those which do include the r argument are useful in obtaining the h- particle generic configurational distribution function, n, (r ), since they become small whenever the particles are too distant from one another or from the group of h particles at location r to interact. This property is due to the fact that the W function tends to factor into products of functions involving clusters of interacting particles. The "interaction" distance depends upon the thermal-wavelength, X, as well as the range of the two-body potential, r_ , and thus increases with lower temperatures. Thus the expansion is expected to converge most rapidly for sufficiently high temperatures and for short-ranged potentials. The actual derivation of the cluster expansion can be found in de Boer's article. The result is:' n K (r h ) = \ 3h z h - 1 I ib< h) (r h )z X (2.2.2) h 1=1 l - where 1 x:x h T' { l ) = ,,.3i-3 I U( I •5w-i'"--»Iwi-i ) x d3r h + r-- d3r h+x-i (2 - 2 - 3) is the "modified cluster integral." Equation (2.2.2) assumes N is large; in the limit as V becomes large so that V » r and V » X 3 and N/V - n the cluster integrals do not depend on the size or shape of the vessel containing the gas, and the volume, V, may be regarded as infinite. Under The discrepancy in the factor of X is due to a difference of definition of W(r^) these being dimensionless . 10 such conditions the functions b* (r) are independent of r so that b^Cr) - b^ and n,(r) -\ 2 Xb,z X = n (2.2.4) 1 ~ X 3 X=l X where the interchange of the limit and the summation is justified as long as the density is small enough that the particles are in the gas phase [15]. In the case of h = 2 we have the pair distribution function which depends only on the relative distance between particles, r 19 = | r_ - r 1 | , 00 n (r r)-.!, Z ib[ 2) (r. _) Z i+1 . (2.2.5) 2 „1 ~2 x 6 x=1 SL 12 The fugacity, z, can be eliminated between these two equations yielding CD where n 2 (r 12 ) = 7^ Vr 12 )(\ 3 n) £ Y l = °' Y 2 = b l 2)(r 12 ) ' Y 3 = 2( b 2 2) (^ 12 > " 2b 2 b[ 2) (r 12 )) Y 4 = 3b^ 2) (r 12 ) - 12b 2 b^ 2) (r 12 ) + (4b 2 - 6b 3 )b< 2) (r^) etc. Rewriting the above, we have n 2 (r 12 ) = n 2 (n^ 0) (r 12 ) + n^Cr^n + ...) , (2.2.6) where and 4° >(r ) - " 2 ( 5l ,r 2 ) 11 / I \ o n 2 (r 12^ = J d r 3^ W 3 ( 5l'52'53' ) " w 2 ( r i' r 2 )(W 2 (r 2' r 3 ) + Vf3'fl } " l) ^ • ( 2 -2.7) The first term n~ (r) involves only two-particle interactions whereas the second term no (O involves two- and three-particle interactions. 2.3 Expansion of the Equation of State The equation of state of a gas may be obtained via the grand canonical partition function, 00 S = E z Z (2.3.1) N=0 The pressure of the gas is then given by PV = kT log E . (2.3.2) We can easily expand log S in powers of the activity, z, 00 CO CO N N 1 N 2 log S = log(l + I Z z iN ) = S z\ - i ( 2 z iN Z ) Z N=l N=l N=l + 3 ( S z Z ) -... N=l so 00 i, PV = kT Z q,z i-1 l where A. X q 2 = Z 2 " 1 Z l = ~^6 I ^l^^Vfl'^ ' W l ( Il )W l ( 52^ = ^el d3r i d3r 2 u 2 ( 5i'l2 ) 12 and it can be shown that in general [16] .3X 1 r a 3JL , T < l \ Thus, for finite-ranged potential, .. ^ h * *im rr~ — ~~ r v - « v x J and CO A kT ^ , X Aim P = ^ 2 b.z . (2.3.3) As with the pair distribution function expansion, we eliminate z with the expansion for n (2.2.4) and obtain a power series in the molecular density n. In order to obtain the familiar virial expansion, we express P as a power series in V : |^- = A + B/V + C/V 2 + ... (2.3.4) Substituting from equation (2.2.4) CO . V" 1 = ~ S Xb.z NX 3 1-1 X and comparing the resulting expansion with (2.2.10) we find A = 1, B = - NX 3 b 2 , C = N 2 X 6 (4b 2 - 2b 3 ) . The coefficients A, B, C are the first, second, and third virial coefficients respectively. Higher order virial coefficients may be obtained in the same fashion. It is to be noted that the second virial coefficient involves no more than two-particle interactions and the third virial coefficient in- volves no more than three-particle interactions. We may express them in 13 terms of the distribution function expansion terms: 2 N N 2 B = ^ a Q C = - — (a 1 -a ) (2.3.5) whe re 00 a Q = 4tt J r 2 dr(l-n^ 0) (r)) (2.3.6) CO a L = 4tt J r 2 dr n^ 1} (r) . (2.3.7) Thus, knowledge of the first two terms in the density expansion of the pair distribution function is sufficient to determine the third virial co- efficient. 2.4 Separation of Direct and Exchange Terms We now use (2.1.6) to express n„ (r) and n„ (r) in terms of diagonal and off-diagonal configurational density matrix elements. For N = 2 and N = 3 (2.1.6) becomes W 2 ( fl'f2 ) = X L<12|e Z |l2) + <2l|e Z |l2>] = n^ 0) (r 12 ) (2.4.1) 9 ~^ H 3i ~^ H 3 W 3 (r 15 r 2 ,r 3 ) = \*[<123|e | 123> + <213|e J |l23> -pH -PH + <32l|e | 123) + <132|e | 123> -pH -PH + <23l|e | 123> + <312|e |l23>]. (2.4.2) The off -diagonal elements represent exchange terms. In W_ there is only one (pair) exchange. In W there are three pair exchange and two cyclic exchanges. We can represent them diagrammatically: 14 12 12 12 12 12 (213| <32l| (132| (23l| <312| It will be useful to arrange the integral, n^ (r) , as follows: n^Cr) = n<» - X 6 <12|e" 2 |l2> X (X b (23|e J |23) + \ b <3l|e Z |3l) - 1)} (2.4.4) n^ 1} (r 12 ) = 2 J d 3 r 3 X 9 <23l| ' 3 |l23> (2.4.5) n^ 1} (r 12 ) = J d 3 r 3 {X 9 <213|e 3 |l23) - \ 6 <2l|e 2 |l2) X (X b <23|e Z |23) + \°<3l|e Z |31> -1)} (2.4.6) n p 1)(r 12 ) = 2 I d ^ 3 U 9 (l32|e 3 |123) - \ 6 <32|e 2 |23) X (X b <12|e 2 |12) + X b (2l|e Z |l2>)} (2.4.7) and where the factor 2 in n^, and n!, accounts for terms in which r,«-» r which are equal. The terms n , xC , and n p (= n p + n ) involve 3-particle direct, cyclic exchange, and pair exchange contributions to the pair distribution function, respectively; hence the subscripts. Similarly, the corresponding contribution to the second and third virial coefficients (2.3.6-7) may be broken up into direct and exchange terms 15 D , C , P a, = a, + a, + a^ (2.4.8) where D r .3 (1), v a i = f drn c (r) (2.4.9) 3_ _(D a i = J d r n c (r) (2.4.10) J-rd 3 r(4 l) (r) + 4 l) (r» J P l P 2 (2.4.11) 16 3. FUNCTIONAL INTEGRALS AND THE DENSITY MATRIX 3.1 The Definition of the Wiener Integral The subject of integration over function space with respect to Wiener measure has been treated in both rigorous and non-rigorous ways by many authors [17-22] and can be approached either in the context of functional analysis or probability theory. The present approach which draws from a paper by Kac [20] represents a comprise between the abstract, but illuminating, probabilistic definition and a simpler, but less flexible, analytic definition. The Wiener integral, E[F[x(*)]}> can be regarded as an average or expectation of the functional, F[x(*)]> over the sample space, Q = [continuous functions x(T):Te[0,t] satisfying x(0) = 0], with the probability measure for subsets of functions being assigned in the following way. Let < t, < t < ... t < t. The subsets I = [x(-r):a n < x(tJ < 12 n u 1 1 b,,...,a < x(t ) < b } H Q (referred to as "measurable rectangles" or 1 n n n J 6 "quasi-intervals") are assigned the probability b. b 1 n P(I) = J dx x ... J dx n P(x 1 |0;t 1 )P(x 2 |x 1 ;t 2 -t 1 )... a a 1 n P(x x , ;t -t , ) , n 1 n-1 n n-1 ' (3.1.1) where P(x, x. ;t, -t .) = k 1 j k y 2tt(t,- t .) exp ( v x i ) r 2 is a continuous function of x. According to the definition of the probabilities (3.1.1), the average or Wiener integral of this functional is given by: E{F n [x(.)]} - J dx x ... J dx n P(x 1 |0;i) ... P(x J X n-l ; n )exp[ - n ;= V(X 1^ ' (3 - 1>5) 18 Now consider the limit as n -♦ ». Since V(x(t)) is a continuous function of t (almost everywhere) it follows that n Xim £ Z V^c(t^ = J V(x(t)) n -» co i=0 q dT a.e. Further, since F [x( # )] is bounded and measurable it follows by the Lebesque dominated convergence theorem [23] that Xim E{F n [x(-)]} - E{ Xim F^xC')]} n — » co n ~* co or t co co n-1 E{exp[- J V(x(T))dr]] = Xim J ... J dx . ..dx II n "* " -co -oo i=0 n x p < x i+ ilVn )exp[ - n E v ( x i>] i=0 (3.1.6) where x„ = 0. Of particular interest in quantum statistical mechanics is the conditional Wiener integral which will be defined as the expectation over the subset of Wiener paths satisfying x(t) = R. It can be shown that this conditional expectation can be expressed: t E{exp[- J V(x(T))d T ]|x(t) = R] = n-1 n P(x., 1 |x.;-) _• n i+l ' in 03 » i= i+ii i'n' n Xim J ... J dx 1 ...dx n _ 1 " - exp[- " L V(x )] n -» co_ ro _« n l > ' 1= (3.1.7) where x = R, x n = and the factor P(x 0;t) is necessary so that the n n 1 J measure is normalized to unity. It can be shown that (3.1.7) holds for all V(x) bounded from below and possessing only a finite number of 19 discontinuities. The appropriate conditional measure for measurable rectangles , I = [x(T):a 1 < xCt^ < b 1 ,...,a n _ 1 < xO^) < b^} < t,...t , < t, 1 n-1 is given by n-1 P(D = Vi . n p (- i+ il-i^i> -[ **VJ d Vl X P(x|0;t) (3.1.8) n-1 where x_ = and x = R. The continuous paths are referred to as conditional On r Brownian paths. It can be shown that they satisfy the following inter- polation formula [24]: X(T) • «< T I )(T"-T?_ +T X (T ")(T- T -) ^ (T-T'y-T) ^ 2 ^^ where t" < t < t ' and § is a random variable, normally distributed with mean and variance 1. The above definitions have an obvious generalization to the space of 3N-dimensional functions, x(t). In equation (3.1.7) we simply replace x. by x. etc., and define the conditional probabilities as P(x. Ix. ;t -t .) = oxt/o ~ kLj k J [2TT (T - T .)] 3N/2 exp <„v_V 2(T k -T.) 2 1 (3.1.10) 2 . where (x. -x.) is the square of the distance between x, and x. in 3N- dimensional Euclidean space. I'j 3.2 Relation Between the Density Matrix and the Wiener Integral We will now establish the important connection between the con- figurational density matrix and the conditional Wiener integral. The following heuristic argument elucidates the relationship. Consider the Boltzmann density matrix element: where H ■ . 2 , t. + V? s Ve> + V- r) - 1=1 1 and we are allowing for different masses. Now we can write -fli n -fig e - P H = (e n } n m n e n (32l) 1=1 since H commutes with itself. Then introducing a complete set of eigen- states of r , |x.), between each operator, we obtain n-1 "££ = T dxf ... f dx 3N . "a = ^i+J 6 ^- I (t n + V%> • Since the operators T and V N do not in general commute, we cannot simply express the exponential of the sum as the product of exponentials: e N e P N # However, we can make an expansion in terms of the commutators, assuming |v(r)l < <= and (3 < »: 21 e *p£- ! k=1 .---' N < 3 - 2 - 6 > where x. = (x. ,x ,...x. ), equation (3.2.5) then becomes »** JL ^ i j 1 •- 1 j t. *- I j * ' 22 n-1 3N/2 ? X II {(§-) exp[- f (x -x ) Z ]} j=0 *" ~ J ~ J ft n " 1 i 1 X exp[- * 2 V N (/R y + r + ±. The error term 0(— ) is assumed to approach zero, and comparison with equation (3.1.7) shows: / it -PH. \ ,-3N 3? r 3/2 r "He , , N 2 n . = X n [m^ exp[- ^ (&& ]} k=l 1 X E{exp[-ft J V - e -e v( ~ r) . 3.3 Approximation to the Wiener Integral 3.3.1 The Mixed Integration Formula The Wiener integral in equation (3.2.8) cannot in general be 2 evaluated exactly (there are some exceptions, e.g. for V(r) = r ). There- fore an approximation will be used. We begin by using (3.2.8) to express equation (3.2.2) in terms of Wiener integrals: "X 2 N ' — (r'-r ) 1 X" 3N n [m? /2 e 2P ~ k ~ k }E{exp[-3 \ VG/&§(t) + r + T(r'-r))d T ] k=l o ~ ~ ~ ~ X |T)(1) = 0} — 2 _ r , 3N , . 3n n : : L » , («\y . - J dx i ••• J d vi n n . n _t [i^] e *pl-\ 2g/n i , — X E{exp[- | J V(y| |(t) + x. + T(x i+1 -x i ))d T ]|T)(l) = 0} where x n = r, and x = r' which are 3N-dimensional vectors and dx. is ~u ~ ~n ~ i 3 3 shorthand notation for dx. -...dx.. ... We make the transformation i,l 1,N x. , - x. . /E— + r, + - (r,' - r, ) ~i,k ~i,k J m, ~k n ~k ~k 25 and finally obtain, 1 E{exp[-0 J V(/P|(t) + r + t (r ' -r))d T ] | T]( 1) - 0} n-1 n p( x |x -±) OXT TXT ' 1 ~1+1'~1 n - r d X 3N r dx 3N — — , J ax l '•• J a Vl P(0 0;1) n-1 1 §(t) X n E{exp[-£JV N (/p(~-+z i + T(z i+1 -z.)) i=0 o + r + (^i)(r'-r))dT]|THl) = 0} (3.3.1) where §(t) s (t=- T| ( T ),...,~t- \(t)) / m l ~ L >n ~ N z ■ (— — x ,(t),..., — - x „(t)) i = 0,1,..., n and x n , - x . = k = 1,2,. ..,N . ~0,k ~n,k Equation (3.3.1) is referred to as the "mixed integration" formula. It has the appearance of a truncation of (3.1.7) except that the function of the points X-.X.....X is replaced by a product of Wiener integrals of functionals over paths between the points. The variables z. are referred to as "path break points." The presence of — multiplying the random function §(t) /n suggests that as n becomes large, the statistical variation of the path (over which V is "time" averaged) becomes small. Similarly as (3 -♦ 0, this path approaches a straight line. The idea of expanding the conditional Wiener integral in a functional Taylor's series about /(3/n T^t) = was suggested by Fosdick [26], and a modification of this idea was employed by * r Jordan in his doctoral dissertation |_ 9] . "Jordan expanded the log of the Wiener integral. 26 3.3.2 The "Simpson's Rule" Approximation Another approach of Fosdick and Jordan [27], which is an extension of work done by Cameron [28] to the case of conditional Wiener integrals, is the so-called "Simpson's rule" approximation. We consider, for simplicity, the one-dimensional Wiener integral, E[f[ ( • )] | x(l) = 0} . A function P(u,t) is chosen so that 1 E[P[x(-)]|x(l) = 0] = \ I du P[p(u,-)] (3.3.2) J -1 where P„[x(0] is any third degree polynomial functional of the form: 1 1 1 P 3 [x(-)] = h Q (T) + J x(T)h 1 ( T )dT + J I x(T 1 )x(T 2 )h 2 (T)dT 1 dT 2 1 1 1 + J I I x(T l )x(T 2 )x(T 3 )h 3 (T)dT l dT 2 dT 3 ' (3 - 3 ' 3) It can easily be shown that equation (3.3.2) is equivalent to the conditions: (letting r^ < t 2 < t~) 1 \ J du p(u, T;L ) = E[x(t 1 )|x(1) = 0} = (3.3.4a) -1 1 |j du p(u, Tl )p(u,T 2 ) = E[x(t 1 )x(t 2 )|x(1) = 0] -1 = t 1 (1-t 2 ) ^ < t 2 (3.3.4b) 1 J -l 2 J du p(u,T 1 )p(u,T 2 )p(u,T 3 ) = E[x(t 1 )x(t 2 )x(t 3 )|x(1) = 0] = (3.3.4c) n where the moments, E[ II x(t,)|x(1) = 0] , are readily calculated by direct i=l L application of the probability measure (3.1.8) or by the recursion relation 27 (3.1.9). It is not difficult to verify that the choice: P(u,t) = < 1-T T > u (3.3.5) -T T < U u > 1-T T > U -P(-u,t) u < satisfies conditions (3.3.4). Assuming that F[x(*)] has a functional Taylor's expansion and that terms of order higher than three are small, we might then approximate the Wiener integral by a Rieman integral: 1 E{F[x(-)]|x(l) = 0} S \ J du F[p(u,-)] (3.3.6) -1 where p(u,-r) is given by equation (3.3.5). This method of choosing an appropriate family of parameterized functions so that the functional integral may be replaced by an ordinary integral for polynomial functionals of degree less than four is analogous to the Simpson's rule for ordinary integrals in which the integral can be replaced by a finite sum for poly- nomial functions of degrees less than four. Hence the approximation (3.3.6) is referred to as the "Simpson's rule" for Wiener integrals. The "Simpson's rule" can easily be extended to the case of functionals in more than one variable. We simply apply the following rule: replace each independent path variable T|. (t) of the functional in the Wiener integral by the function p(u. ,t) (3.3.5), and instead of taking the average over paths, take the average with respect to the independent random variables u. , i = l,2,...,n where each u. is distributed uniformly over [-1,1]. Konheim and Miranker [29] and Fosdick and Jordan [27] showed how to generalize the approximation for Wiener integrals so that higher order rules could be obtained. By taking linear combinations of the functions (3.3.5), 28 m (u.t) = m (u,t) = S c p(u ,t) , (3.3.7) i=l and replacing 9(u,t) for the path and integrating as before over u.,...,u , the approximation would be exact for polynomials of degree 2m+l, providing 2 the c. are the roots of the polynomial m m-1 , m-1 , , lN m ] n , _ Q s z-z +TT Z - ... + (-1) — 7 = . (3.3.8) z . m. It can be shown that the application of the above approximation to the functional, exp[- | J V(x(t))<1t] yields an error of the form (8/n) K where K is a rather complicated m m integral expression involving the term 1 exp[- | J V(s9 m (u,T))dT] where s is integrated over the range [0,1]. In the case of the Lennard- Y/(4n) Jones potential, the exponential term is bounded by e , and, because of the strongly repulsive part of the potential, is probably reasonably small providing 9 (u,t) is real. Unfortunately, for m > 1 the roots of (3.3.8) lead to complex values of c. so that the K can be quite large, especially for calculations in which the repulsive part of the potential is important. Actual path integral calculations carried out by this author have verified this problem. For this reason, the present work does not employ these higher order approximations. Applying the "Simpson's rule" to each of the Wiener integrals on the right side of equation (3.3.1) yields: 29 where and and E{exp[- | J V N (/ 3 (i- |( T ) + fi + T(z. +1 -z.)) /n + r + ( I±L -)(r'-r))dT]|Tl(l) = 0] du 3N t 2 /n + t(z. +1 -z.)) + r + (I±i)(r'-r))dT] + 3 (3.3.9) p'(u . ,t) = ( p(u . ,t),. . . ,— p(u ,t)) i = l,...,n, v/m 1 /nijj £ ( "ij ,T) = (P (u iji' T) 'P (u ij2' T) ' p(u ij3' T ^ J = 1 "--» N z. = ( x. .,..., x.. T ) i = l,2,...,n-l - 1 / mi ~ U Vmj, ~ lN z = z = ~o ~n u.., x.. i = l,...,n; 1 = 1, . . . ,N are 3 -tuples It can be shown that the error, S, is 0(— r). There are n such integrals in n the mixed integration formula (3.3.1); hence the total error in calculating , ,N, " 0H N, Nv . . , nf l v (r |e |r ) is of order C (— r) . n 30 4. THE COMPUTATIONAL APPROACH 4.1 Transformation of the Matrix Elements to the Center of Mass System The density matrix elements used in equations (2.4.4-7) can be expressed via equation (3.2.8) in the form of Wiener integrals which can be approximated by the mixed integration - "Simpson's rule" combination. The matrix elements are of dimensionality 6 and 9. Since we are assuming the gas is contained in a vessel large enough so that its size and shape will not enter into the problem, it is possible to transform to a center of mass coordinate system so that we can reduce the dimensionality of the density matrix elements to 3 and 6. We consider now 3-particle conf igurational density matrix elements. Accordingly we make the transformation to the R, r, Q coordinate system defined by: 5 = I ( £l + l2 + 53 ) !r = £l + Zl + £3 I = ll - 5l !r = ¥ll~l\> S = 13 - 2 ( £l + f2 ) ! Q = 3 £3 " ? ' (4 ' 1 - 1) With this choice of position and momentum coordinates, the corresponding operators satisfy the usual commutation relations, and the Jacobian, B(R ,r Q ) r , — ^ — ^ — r , is unity (for each q component). Hence, the choice r 2q' r 3q ; I R, r, Q) = 1 1 2 3) yields properly normalized eigenkets of the operators, N R , r , Q P . Assuming V (r ) = Z V(r. .), the Hamiltonian becomes i + V R r Q 1 A 2 l R "» m r " 2 and m Q " 3 whe re m R = 3 , m = -r and r 31 We next consider the eigenvalues of the permuted kets, P|l 2 3). We obtain them by operating on p| 1 2 3) by R , r , Q P . For example let P|l 2 3) = | 2 3 1); then R° P |2 3 1> -£ (r ° P + r ° P + r° P ) | 2 3 1) ■ 5 (r 2 + r 3 + r x )|2 3 l) = r| 2 3 l) r° P |2 3 1) - (r° P - r° P ) 1 2 3 1> - (^-r^ J2 3 l) = (Q - \ r)|2 3 1) Q° P |2 3 1) =[r° P -\ (r° P + r° P )]|2 3 1) [r l " 2 (r 2 + r 3 } ^ 2 3 l) • C-lj -|r]|2 3 1) . Hence , | 2 3 1> = c|R, Q - \ r, - \ Q - | r) (4.1.3) where c is determined by the normalization condition and is equal to unity since | R, r , Q) = | 1 2 3) . Similar manipulations yield the other corres- ponding permuted kets. |l 3 2> « |R, Q + j r, - \ Q + | r) (4.1.4) | 2 1 3> = |R, - r, Q> . (4.1.5) Since the center of mass coordinate, R, does not appear in the Hamiltonian, the Hamiltonian may be broken up into two commuting parts: H 3 " H R + "s" 1 32 where and 4 2 2 P P h" 1 = ■—- + t 5 - + V(r) + V(Q - -^ r) + V(Q +fr) , 3 2m 2m^ ~ I ~ ~ ^ ~ r Q Thus, (1 2 3|e 3 |l 2 3) = rel = (r,Q|e 3 |r,Q> -pH rel = 3 3/2 X" 3 (r,Q|e ' 3 |r,Q> . (4.1.6) Similarly, rel <2 1 3|e 3 |l 2 3> = 3 J// \' J (-r,Q|e J |r,Q> (4.1.7) -(3H 3 (2 3 l|e | 12 3) -pH rel = 3 3/ V 3 (Q - \ r, -|q - | r|e 3 |r,Q> (4.1.8) (1 3 2|e 3 |l 2 3) -6H rel = 3 3/ V 3 . (4.1.9) Similarly, the two-particle density matrix elements become: rel -j3H_ _ ,„ _ "BH,, (4.1.10) where H" 1 = P 2 + V(r) . (4.1.11) 2 r 33 4.2 Multiple-Integral Approximation of the Density Matrix 4.2.1 Three -Dimensional Matrix Elements We now combine equations (3.2.8), (3.3.1), and (3.3.9) to obtain the approximating integrals for the three-dimensional density matrix element (4.1.10): -pH rel ~( r '- r ) 1 2 3/2 \- 3 = e 4p E{exp[-P J V(/2p T1(t) + r + T(r , -r))dT]|T)(l) = 0} (4.2.1) -(r'-r) 2 3n J d ^xn J ^T exp[ " n t Z =l X f V(/2£ (— p(u t) + x + T (x.-x. )) yn + r + (^i)(r'-r))d T ] + O(^) n whe re d^ xn = [P(0|0;1)]" 1 n LP(x i+1 |x.; ^)]dx 3 . . .dx 3 ^ , (4.2.2) and x = x = ~o ~n P(U ± ,T) S (p(u il ,T),p(u i2 ,T),p(u i3 ,T)) , and we have expanded the product of n 3-dimensional integrals over u as one 3n-dimensional integral over u.,u ,...,u , the product of exponentials becoming now the exponential of a sum. du 3n Since ' du. i — - — = 1, the multiple integral (4.2.1) can be regarded as an average of the exponential, and the Monte Carlo method can be employed 34 for its evaluation. The integral over t can be evaluated by a quadrature formula or by an exact expression. Since p(u.,T) is a discontinuous function of t, the point of discontinuity for the q component being u. , it will first be necessary to break up the integral into a sum of integrals over regions in which p(u,t) is a continuous function of t. To do this, we consider a particular value of u. = (u.., u, „, u.»). Then, defining the time partitions t. by: T n = 0; t. = |u..|, j = 1,...,3; t, = 1; we let k. be distinct integers e {0,1,2,3,4} which satisfy: k~ = 0, k = 4, and S T k < T k < T k < \ < T k H l ' 1 2 3 4 (We need not consider the case of equalities since events in which u's are equal occur with probability zero.) It is to be noted that t. and k. both J J depend on i; we will not explicitly express this dependence, however, for notational simplicity. The integral over t is then partitioned in the following manner: 1 4 k J +1 /55 J V(/2p p(u, T ) + ...)dT = 2 J V( Jf- p(u., T ) + ...)d T . J=0 _ k. u. . Let S. be the sign function of u.., i.e. S. = i 1 ^> , j = 1,2,3. Then, J ° ij' j ju I » J » « by (3.3.5), p(u ,t) = p(S.t.,t) = S.p(t.,t). Consider now the interval (0,t ]. Here T < t., j = 1,2,3; so from (3.3.5), p(t.,t) = -t; hence K l J J p(u. ,t) = (-S t,-S t,-S_t) . Next we consider the interval, (t, ,t, 1: ~ ~i 1 z j k ' k„ J we have p(T ,t) = 1-t p(u ,t) = S (1-t) K l k l k l P(T k ,T) = -T pd^ ,T) = -S k T P(T. ,T) = -T p(u , T ) = -S T R 3 k 3 k 3 35 (i) U 1 l V , 11 WC UCLIUC tllC LU1CC U Lllli: US LUUd i vctiurs (J j = 0,1,2,3 by = U ' Wj 'q = ^n~ x t L V^'V 5 q = 1,2,3 (4.2.3) and so on. Finally, if we define the three dimensional vectors a ~J ° -J J or 1 3 k J +1 r--^ (^ J dT ... = I J V(|uJ lj; + U^\|)dT (4.2.4a) J =0 T k where U, (ij) = a? i) + w. , + r + — (r'-r) -1 ~j ~i=l ~ n ~ ~ U^ X) = - a. (l) + A^. + - (r'-r) (4.2.4b; ~z ~j ~i n ~ ~ and U) . = /2g x. , AUJ. - U) - UU l K -l ~i ~i ~i-l The integral over t is now in a form where it can be evaluated exactly (in the case of the Lennard-Jones potential) or by means of a quadrature method. 4.2.2 Six-Dimensional Matrix Elements In like manner, we can express approximations for the six-dimensional density matrix elements: , - BH 3 el , (Q',r'|e |Q,r> . V, There is an asymmetry between the coordinates Q and r due to the different 2 1 masses -r and •= associated with them; accordingly, the associated paths will be given different labels, q(T) and T](t) respectively, and the corres- ponding "Simpson's rule" parameters will be denoted by u. and v. re- spectively. It can be seen from (3.2.8) that the associated Wiener integral can be expressed in the form: „3/2.6 r 2 tt , .2 , 1 rr . , N 2 n 3 \ exp[r- — r (Q'-Q) + "o ~o (r'-r) ] * x- ~ ~ \ X = 0} = r du. r d^ r ^ r ^1 ~ J J ^xn J ^yn J 3n J 3n R n - X exp[- * E (V(p(u.,T),x.,r,r') i=l + V(p + (u.,y.,T),x^,Q + ,Q 1+ ) + ^ ( £~ ( ~i'~i' T) '*i'2~'3'~ ) J (4.2.9) where V(p(u t),x r,r') = f V(/2j3 (7- p(u., T ) + x. , T+i - 1 + T( ~i^i-1 }) + I + ( ^~^ )( f ''f ))dT (4.2.10) and + 1 *i = 2 ( ^ 3 Zi - *i^ £ ( ~i'~i ,T) "2 (v ^ 3 £ ( yi' T) ± P( u i» T )> • We now express the integral over t as a sum of integrals over continuous functions. Proceeding as before we fix u. ,v. and partition the integrals over t. The integral, V(p(u. ,t) ,x. , r ,r') , is identical to that for the 3-particle density matrix (4.2.1). The remaining two integrals require six partitions: 38 T q "i, q |. ! 3< ^ 6 T Q - T - 1 . Now, let k_,...,k- be defined so that U D k_ = 0, k = 7, and ! s /2B xT , Aw? = uu7 - uT . ~i K ~i ~i ~i ~i-l It is to be noted that equations (4.2.4) have the same form as equations (4.2.12). In both cases, we must calculate integrals of the form: k k 1 - J v( iHi + H2 T i )dT = I v ( x ( T >) dT » V = T . T. J X(t) = /u* + 2U-U 2 t + \] 2 2 t 2 or " k 1 1 i. X( T ) 12 X(T) 6 in the case of a Lennard- Jones potential. Although this integral can be evaluated exactly, the resulting exact expression is a rather involved function of at least three parameters: — y— , -r , and T. . Since, in an actual computation, the values of the parameters involved in V 1 would correspond to a large number of Monte Carlo events, it would be time consuming to evaluate the exact expression each time V is needed; furthermore, a table of integrals would be prohibitively large. By approximating V 1 by a quadrature formula, the time of calculation of V 40 can be kept small without sacrificing accuracy. The two point Gauss- Legendre quadrature has the advantages of being simple and having an error term proportional to (t,-t.) V " (as in Simpson's rule). Using this k J approximation we have, V' = A[V(t") + V(t + )] + 0(A 5 ) where T.-T. + A = 2 - 1 and t" = t. + A(l + .57735027) 4.3 Strategy for the Monte Carlo Approach We now discuss the general procedure to be followed in the evaluation of the multiple integrals used in the calculation of the pair distribution function by the Monte Carlo method. Equations (2.4.4-7), (4.2.1) and (4.2.9) suggest that in general we will be concerned with multile integrals of the form, where the integration "over the third particle," r„ , in equations (2.4.4-7) has been transformed to the relative coordinate, Q = r„ - t (r 1 + r») . The dimensionality of the Q-space integration can be reduced by considering the following coordinate system: Let us take r = r„ - r 1 along the Q -axis as in Figure 2 with the center of mass of particles 1 and 2 at the origin. Then Q is the radius vector. Since the integrand of 7? is actually a function of two- and three-particle density matrix elements of particles 1, 2, and 3, and since the interactions are spherically symmetric, there is 41 azimuthal symmetry with respect to Q, i.e. the density matrix elements shall be independent of 0. Hence in (4.3.1) we are free to set = — Figure 2. Coordinate System for Calculations ;o that Q = (0,0 ,Q ). The integral over Q is then restricted to the V ■Q right half-plane and takes the form: , 3n , 3n du r dv * - 2TT J V% I dQ Z I %n I ^xn J 3T J T31T NS'M*£> (4.3.2) To convert the multiple integral, 7?, into a form useful for Monte Carlo sampling we must first consider a probability distribution for ,Q . We 2 assume for now that a suitable distribution function, F(Q ,Q ), has been y z chosen. The integral then becomes i ■ 1 1 «£*». ra£ v j %n j d , xn i if j ^ 0-oo 2 2 x [tt F(Q,x,y,u,v) ttF, r(Q y ,Q z ) ■] = • (4.3.3) 42 The multiple integral, 71, could at this point be estimated by a Monte Carlo procedure. The variance can still be reduced somewhat, however, by further importance sampling. It will turn out that part of the function, F, can be factored out, normalized, and used as a probability distribution for the selection of certain variables. Let these variables be symbolized by the vector § and the rest by T\. Thus F(§,T]) = U(5)G(|,T)), (U(?) > 0) and 71 = J d|w(|) J dT]aj(Tj)F(|,Tj) (4.3.5) = J d§w(|)U(|) J dTp(T))G(|,Tp (4.3.6) where w(§) and ^(§) are weighting functions. Now define the function U(§) P(|) = -f- (4.3.7) where Z = J d|w(|)U(|) (4.3.8) so that 7? = Z J d|w(|)P(|) J dTjiu(Tj)G(|,7j) . (4.3.9) The variables ^ a * e thus weighted by the probability distribution w(§)P(§). It is intuitively apparent that it is advantageous to use P(§) in the probability distribution for ^ since it biases samples of § in favor of those for which the factor U(^) in F(^ ,Tj) is relatively large. In general, it is difficult if not impossible to generate random variables with distribution P(?)w(^) directly. Instead of taking independent samples, it will prove easier to generate a Markov chain of random variables £(k) whose distribution converges to P(§)w(^) as k -» co . In the case of a discrete sample space, Metropolis et al. [30] showed how to do this in 43 calculating averages with the classical Boltzmann distribution, and the method was formalized by Hammersley [31]. Jordan [9] applied the method to the evaluation of Wiener integrals by extending Hammersley' s argument to the case of a continuous sample space. The present approach stands some- where between that of Hammersley and that of Jordan. Let us assume that in a Monte Carlo scheme, the first choice § (0) = ^ is determined by some arbitrary distribution (e.g. if it is always the same number, the distribution of |(0) would be a 6 function). The next choice £,(1) = 5' would be determined by a Markov transitional probability density p(§'|5) which satisfies the conditions J P(|'||> d 5' = l (4.3.10) and w(|')P(|') = J p(|'|?) p (p w (p d | • (4.3.11) Equation (4.3.10) is the usual normalization condition for conditional probabilities. Equation (4.3.11) prescribes that if § has distribution w(^)P(£)> then the new random vector ^', determined by the Markov transition probability p(§'|§), will have the same distribution. The next random vector %(2) ■ 5" would be determined by the same transitional distribution, p(§"|C') and so on. The fact that the distribution for §" depends on §' and not on ^ is the Markov property of the process. Continuing in this way, we would obtain a sequence of vectors, ^(k), the distributions of which converge to P(?)w(^), equations (4.3.10) and (4.3.11) being necessary and sufficient conditions for this convergence. Another condition, that the states § be ergodic is also required and holds in the cases studied here. 44 Furthermore, if equations (4.3.10) and (4.3.11) hold and the states %(k) and T|(k) form an ergodic set and if 71 = J d§w(§)P(^) ! ' dTp(Tl)G(5 ,11) t M is estimated by — 2 G[^ (k) ,T|(k)] where each Tj (k) is distributed in- M k=1 ~ dependently by W(\\) (^(T]) can also be regarded as a Markov transition probability density for T| since it satisfies conditions similar to (4.3.10-11)), then by the law of large numbers for Markov chains [32] the estimate will converge to the integral, 71 > with probability one as M -• «, In order to find a transition probability satisfying (4.3.10-11) we introduce the "trial transition density" P*(|'||) = w(|')p**(|',|) (4.3.12) where p**(?'>?) is a function satisfying P**(|',p = P**(|»|') (4.3.13) and J w <5 , )P**(§ , f|)d|' = 1 • (4.3.14) We then have the following ■k Theorem: The choice, HV) p(5') p*P(fr Lf wr <1 P*(|'||) + 6 <|'-|) J d§" P(|'|?)= HZ") (4.3.15) ^Pcfr^ P(S") P(S') X P*(|"||)(l - p^-) if p(f)-> l /V This theorem resembles that of Jordan [9]. The choice of p*(§'|§) and the proof, however, follow more closely the approach of Hammersley and Handscomb [31]. 45 satisfies conditions (4.3.10) and (4.3.11). Proof: (1) We show I = J d§'p(§'|§) = 1: Let P<§') r( i ,)H 7(fr ; then from (4.3.15) and (4.3.12) i = J "OpP^CI'^Mpdi' +J w(|')p**(|'|)d5' {|':r(|') < 1} U' :r ?') > l ) + J w(5")p**(|",|)(l-r(|"))d|" U":r(|") < 1} the first and the second half of the third terms cancel leaving I = j w(| , )p**(| , ,|)d|' = 1 . (2) We prove (4.3.11) first by showing p(|'||)P(|)w(|) = P(?||') p (i') w 0p • (4.3.16) Case 1: If P(§') < P(£): PCS 1 ) p(|'U)p(|)w(|) = p ,v *(r»|)w(|*) p^y- p (|)w(|) = P**(|,|')w(|)P(|')w(|') - p(§|? , )P(?')w(?') • 46 Case 2: If P(§') > P(§), §' * §: (if ?' = ? (4.3.16) holds trivially) p(|'||)P(|)w(|) = p**(|',|)w(|')P(|)w(|) = P**(|,|*)w(|)P(|)w(|') PCS) = p*(|»|') pfry p(|')w(|') = pCglpPCpwCi') . Thus (4.3.16) holds in all cases. Hence J p( rip p (|) w (p d i = j p ( iir )p ?) is chosen to be a function of the path variables only. Further, in the case that ^ involves two paths x and y (as with n p (r)) we assume ~ 2 the "trial" paths to be chosen independently with the same distribution, i.e. P**(J',|) = f(x',x)f(y',y) where x and y are 3n-tuples representing the path break points. The "trial transition probability" for the joint distribution then factors: 47 p*(|'|p = f ^ , ^ )w i ( ^ ,)f( y , »y )w i ( y ) --- where w. (x) is the weight function for the path x and is given by the relation w.(x)dx 3(n " 1) s du. . (4.3.17) 1 ~ r nx We determine f(x',x) by requiring for simplicity that the "trial transition probability" be Gaussian. We thus choose the form for f(x',x) n 2 w 1 (x)f(x',x) = C exp[- | Z (ax^ + px^ + yx + &x L > ] i=l 3n = (2tt) 3/2 (^) 2 exp[-| Z (x 1 -x i _ 1 ) 2 ]f(x',x) (4.3.18) i=l where we have used equations (4.2.2) and (3.1.10) in expressing w..(x). The symmetry requirement (4.3.13) implies f (x' ,x) = f(x,x') which from (4.3.18) yields the conditions: 2 . 2 a - 1 = y 2 2 P - 1 - 6 ap + 1 = Y6 a6 - PY • The second of these equations is implied by the first, third, and fourth equations so that only three are independent, and a, (3, and Y can be determined in terms of the arbitrary constant 6. The normalization condition (4.3.14) determines C in terms of this same constant. We thus finally obtain: 48 3n 2, 2 / xir/ i \ /o \ 3 /2 , n(l+6 K / ,, ,.2. n w 1 (x)f(x',x) = (2rr) ( — ^ — «■) exp{-(l+6 ) -^ X S [x!-x. . --~=7 (x -x J] 2 ). f4.3.19) i=l - 1 - 1 " 1 /i+6 ~ ~ If we make the transformation to the 3n-tuple, x = (x, ,x 2 ,...,x ): x' = —=- (X + 6x) ~ /ST 2 ~ ~ then X has the distribution w (x) and x = X- = 0- In summary, the random sequence (^ (k) ,T|(k)) is generated by the following iterative procedure. Assuming we have chosen %, (k-1) ,T^(k-l) and have determined U(§(k-1)): (1) Choose Tl(k) according to the distribution w(T|) . (2) Choose the "trial" vector §' in the following way: Paths : The 3-vectors x-ijXoj'-'jX -i are generated by the weight function, w. (x) , i.e. according to the "interpolation" formula (3.1.9) (x = X =0)- The "trial" path coordinate is then determined by the formula: x! = * (x, + 6x.(k-l)) i = 0,1,. ..,n . (4.3.20) ~ x /l+6 Other variables in £ ' : "Trial" choices for variables which are not path variables are made according to their corresponding weight functions. P(S') U(§') (3) The ratio, r(|') - p(^(k-l)) = u(%(k-l)) ' is calculated - If r(§') > 1 then 5(k) - §*. If r(§') < 1 then 49 r I' with probability r(^') |(k) =/ ^ 5(k-l) with probability l-r(§') . The problem of choosing the initial vector £,(0) still remains. In theory, any choice will yield a sequence converging to a random variable which has distribution, P(^)w(§). The number of steps required for con- vergence depends (in the probabilistic sense) on the distribution of the initial vector, ^(0)- The selection process will be discussed later, however . n The multiple integral, —, is thus estimated by the sum, 1 M — E G(^ (k) ,T](k)) . The normalization function, Z, can be estimated by M k=l ~ Monte Carlo sampling. If 6 = 0, the estimate for Z comes out as a by- product: 1 M Z = ± Z UO-'Ck)) (4.3.21) k=l ~ where Z ' (k) are the "trial" vectors used to determine §(k). Another advantage of taking 6=0 arises when the same random numbers, and, therefore, the same trial vectors, are used for calculations involving only different values of r. Consider a particular event, k, in a Monte Carlo sampling process. Assume a trial vector, 5'> has been selected and is tested as in step (3) above. Whether it "succeeds" or not depends on its value and the value of r. Suppose that the calculations are performed for a set of values of r. There will in general be a group of calculations over a set of r values such that the trial vector succeeds and ^(k) = ?'• There will be other groups involving r values such that the appropriate vectors are given by perhaps ^(k-1), or perhaps §(k-2) and so on. Each 50 group involving a different value for % (k) will in general refer to a sub- set of one or more r values. By carefully keeping track of these groups and performing calculations which are independent of r only once for each group, considerable computational time can be saved. 51 5. EVALUATION OF THE PAIR DISTRIBUTION FUNCTION TERMS 5.1 The Cyclic Exchange Term The simplest three-particle contribution to the pair distribution function (2.2.6) is the cyclic term (2.4.5). As suggested in the previous section, by applying equations (4.1.8) and (4.2.5) to equation (2.4.5) we obtain 3 2 Q* n£ 1} (r) = 2e 4p J d 3 Qe P E{f[T)(t) ,r ,Q _ ] X F[l]"(T),Q",-Q + ]F[T] + (T),Q + ,-r]|T)(l) = q(l) = 0} (5.1.1) + where F, T] (t) were defined in (4.2.6-8). The multiple integral approxi- mation is then employed, and, as in equation (4.3.1), we have yn - 2 ' where * l$r «*- 5 .=, I (v " + v 23 + »S>*1 + "(-S) 2 i=l o n (5.1.2) v" = V(/23(t- P( u ' t > + x -,- i + t(x.-x, ,))■ + r + iz y n «, ~ ~i-i ~i ~i-i ~ ( I± r li )(Q"-r)) (5.1.3a) V 3 J 2 V(/23(t- P "(u,v,t) + x" + T(x"-x7 .)) + Q" zj y n ~ ~ ~ ~i-i ~i ~i-i ~ - 2 ( I± n ^ i ) Q) (5.1.3b) 52 V 2 * h V(/2p(7- p + (u,v, T ) + x| + T (x+-x+ ,)) + Q + (I: T :i)( S + + 5 )} ' (5.1.3c) Consider now, the integration over the third particle. According to the discussion in Section A. 3: Jd 3 Q ... = tt J dCyiQ z ... where Q = (0, Q y , Q z > . Choosing spherical polar coordinates, Q = Q|sin 9|, Q = Q cos 6 and Q = |q| W e get 1 J d Q ... = 2tt J Q dQ j d(cos 6) ... -1 2 2 -0 /B 2 2 The factor Q e K suggests that we choose according to a ^ distribution with three degrees of freedom: -I -1 1 2 T(Q 2 ) = 2tt 2 B 2 (Q 2 ) 2 e" Q /P (5.1.4) J T(Q 2 )dQ 2 = 1 . Equation (5.1.2) then becomes: n (» (t) .£e"*e r2 jr«wJ*<£pl /2 -1 _____ Q O'i where X = /2nB is the thermal wavelength. The factor, exp[- *- Z V dT] n i = 1 « 12 1 L can be taken as a weighting function for importance sampling with a Markov process. We therefore define the probability distribution: 53 n 1 P c (x,u,Q 2 ,cos 6) = [Z c (r)] _1 exp[- | S J v" d T ] (5.1.5) i=1 where the normalization condition requires: 0-1 2 X exp[- | 2 J V^2 d T ] . (5.1.6) 1=1 o The multiple integral finally takes the form: 3 -^-r 2 n C 1>(r) = 7" e 4? Z C (r)8 C (r) (5.1.7) where 3n n l «C (r) 2 J dP C S %n I 3T «*rf- I * I ^Xs^ (5 - 1 - 8) and we have used the shorthand notation: J dP c ... » J r«W J f J * J if P c (x,u,Q 2 ,a)... 0-1 2 dv 3 where a = cos 9. Since | dP [* d|jl P — ^— = 1, g r (r) can be regarded as a weighted average and can be estimated by Monte Carlo sampling as outlined in Section 4.3. Thus M n 1 g c (r) - ^ I exp[-^ 2 J (V^(k) + V^(k))d T ] (5.1.9) k=l 1=1 54 where the argument (k) implies that the variables of integration are now 2 the realization of a stochastic process: x(k) , y(k), u(k), v(k), Q Ck), 2 a(k) where x(k), u(k) , Q (k) , and a(k) are components of a Markov process and y(k) and v(k) have independent steps. 5.2 The Direct Term 5.2.1 Separation of the "Hard-Core" and "Correlation" Terms The three-particle contribution to the pair distribution function which does not depend on the exchange of particles will be referred to as the "direct term," n (r) , and is defined by equation (2.2.4). Choosing our usual center of mass coordinate system and making use of the azimuthal symmetry, the integral becomes CO CO n^ = tt J dQ 2 J dQ z [W 3 -W(W + + W"-1)] (5.2.1) where we use the shorthand notation: W 3 = \ y r -pH -p V W = X 6 (l 2|e 2 |l 2) = W + = X 6 <1 3|e 2 |1 3) = , -PH -PV W' = X b <2 3|e Z |2 3> = 0, Q > 0], S. Then let us partition S by a circle of radius R n about particle 2. The set S is then considered to be made up of two subsets S n and S' where s o = *V V % >0 > Q Z > °' Q y + (Q Z " 2 r)2 * % ] and S 1 = compliment of S in S = S/S n = S^.. Thus we may express the integral over Q as CO CO 2TT J d i! d % ••• = 2tt JT d i d % ••• o o s - 2 *n d i d * z ••• + 2n n dQ y dQ z ••• (5 - 2 - 3) s Q s> Recalling that the interaction potential, V(r), is positive for r < 1 and increasing rapidly as r becomes smaller we expect that the average inter- action, V.., should be quite large for small values of r.-r. < 1. Thus, we expect that there should be some maximum value of L, < R„ < 1, such 56 -0V 23 that e is negligible for Q over the region S for most paths so that -0V 12 -pV 23 -pv 13 -£V 23 (e e e ) and (e ) can be neglected. We will call R n the hard core radius and S the hard core region since it is in this region that the hard core (strongly repulsive) interaction between particles 2 and 3 dominates the density matrix elements. Having chosen such an R, , the direct term becomes: n£ 1} (r) = 2tt W JJ dQ 2 dQ z (l-W + ) y s o + 2tt JJ dQ dQ [W -W(W + +W"-1)] (5.2.4a) S* = n_ (r) + n_,(r) . (5.2.4b) The first integral involves only a two particle term, W , and can be evaluated by quadrature over a table of values for W . The details will be described in Appendix A. The second term can be expressed as one Wiener integral: 2 Jn , -^12, "P V 23 -P V 13 . -P V 23, , "P V 13. I dQ y z 2tt JJ dqj d Q z (e 12 (e 23 e 13 - - (e 13 > + 1)) -pv 23 -pv 23 where in the second term the (e ) cannot be replaced by e since -pv 12 -pv 23 there is a correlation between e and e so that -pv -pv -pv -pv -pv -pv (e ll e 23 > f , and similarly, ^ -pv 12 -pv 13 (e )(e ). It will prove useful, however, to employ this approximation and split the integral in the following way: 57 1 2 n gl (r) - n s , (r) + tig, (r) whe re and i .. 2 " pv i2 " pv 23 _ P v n n s ,(r) = 2tt JJ dQ^ dQ z . (5.2.5b) That the two integrals should converge is shown by the following argument. Let us consider the asymptotic behavior of the integrand of n Q ,(r). According to our interpretation in Section 3.2, V.. is the time averaged interaction between particles i and j whose uncertainty in position is of order \ =/2n(3. As |Q| -* « we expect V and V to approach zero like 1 5 9 (—7) so long as JO J » CT MAy , the maximum deviation of the particles from locations 1, 2, and 3 for a particular path. According to the conditional Brownian measure, these paths which deviate from their initial positions a distance more than X do so with exponentially small probability. Thus, for each X there is a maximum distance, IL. , such that if |Q | > RwAy then |Q j » (TwAy for "most" paths and V ~ 0(— 7-), i = 1,2. Now under r 13 these conditions -pv - P V _ (e - l)(e - 1) = P V 23 V 13 + 13 23 for "most" paths. Thus we expect R Q ] (IS where S = [Q + ,Q": Q + > Q" > 0; Q + + Q" > r > Q + - Q - }. Let us approximate n , (r) by an integral over the truncated region, S" = [Q + ,Q~: r° < q" < r„] n s , where it is assumed that R. is large enough so that the truncation error is negligible. We then partition this region by the sets, S., i = 1,...,7 formed by the interactions of pieces of concentric rings of radii, t + - t R, , R„, R • S. = (Q > Q : R. > Q > R. n 1 OS, i = 1,2,3 which will be 1' 2' 3 i L i — — i-l J called "basic sets." We let s. = s. n s~ s. = s_ n s7 111 4 J 1 s„ = st n si s. = st n s~ 2 2 1 5 3 2 3 " S 2 ° S 2 S 6 " S 3 " S 3 s 7 = s- u s; u s; U S+ where S is the compliment of S with respect to S . It is easy to verify 7 that S. OS. =0 (i^i) and U S. = S". An example is illustrated in l j J i=l i Figure 6. The sets, S., have been chosen so that within each region, (D (?>Q )(D (§,Q ) - 1) does not change appreciably. Thus, in the example of Figure 3 if we choose perhaps R_ = .8, R, = 1.6, R„ = 2.4 and R„ = 3.2, there would be relatively little overall change of F* with respect to Q between these radii. Likewise 1-D (Q ,^) should have little structure with respect to Q within the same intervals. Hence the product (D (c,Q ) - 1) (D (£,Q ) - 1) will not change appreciably over the inter- sections of "basic sets." Figure 6. Example of Region Partitions for S": The Regions S.,...,S_ In summary, the integral in equation (5.2.14) over the infinite region, S ! , is approximated by an integral over the finite region, S" , which is partitioned into subregions, S., and we have Vj.Cr) ~ L n SI1 (r) = 2 W E JJ d(Q + ) 2 d(Q~) 7 Z i=l S. l .+ ..+ X J dP(D (Q ,|)-1)(D"(Q",|)-1) (5.2.15) where + q" % - (0,Q y ,Q z ±2 -> h^ 2 (Q + ) 2 (Q') 2 ] (Q z+ |r) 2 (5.2.16) J 67 Since it is assumed that there is no real structure of the integrand with respect to Q , Q in each region, it will be reasonable to assume uniform ■K 2 -2 distributions of (Q ) and (Q ) over the regions S. for Monte Carlo sampling. Thus, we estimate (5.2.15) by 7 M i 1=1 1 k_i X (D'O^i'lk)" 1 ) (5.2.17) +.2 7 where A(S.), the normalizing constant, = ff d(0 ) d(Q ) , and M = S M. . i J J i=1 i b i Equation (5.2.17) suggests that the same random variables § are used for each region S.. Although this simplification is justified in the sense that the estimate n QM (r,M) -• n Q ,,(r) as M -* » (provided M. -* «, VS. f 0), correlation between the regional estimates over each S. does slow down the ° l overall convergence. It is assumed, however, that the time saved in computation makes up for this fact. The method of generating the variables §, was discussed in Section 4.3. Details such as the generation + ~ k of Q random variables over the various regions and the determination of A(S.) are covered in Appendix B. It remains to determine the sample sizes, M. . We would like to choose them so as to minimize the total variance of the estimation, 1- 7 fi ,,(r,M), subject to the condition, Z M. = M, where M. = if S. =0. b i=l l ii We consider the general problem of evaluating integrals of the form: 68 g = J dQ J dpi f (|,Q) where J* dp, - 1 S" K = E A( V mttJ 1 d ? J" d * f( i'2> J -1 J s J where K S" = U S '; S. OS. =0, i^j; andA(S.)=j dQ . J- 1 s. S J Then let the total Monte Carlo estimate be given by K L = 2 A(S ) (5.2.18) where M. "I" j? "Sk-i^' J J k=l (5.2.19) is the Monte Carlo estimate over the region S . , and IL and Q, are random variables with distributions d\i and uniform over S., respectively. Neglecting correlations between regions S. and S., we obtain for the total variance : K var^) ? S [A(S )]" var( M ) (5.2.20) and ■k.' h A varf£( Sk-8k J) M. k=l J var«f> M ) = — - S var{f(5.,Q. (j) )} + 2 S cov{f(§ Q (j) ),f(? Q ( J } )} kS )]2 * (5.2.21) Thus, equation (5.2.20) becomes K [A(S )] 2 var^) r E ^ var(f) . (5.2.22) j=l J As an approximation to the minimization of the actual variance of g^, then, we choose M. so as to minimize (5.2.22) subject to the condition that K J E M. = M. Regarding the corresponding continuum problem of minimizing j=l J K T = E C./x. (5.2.23) subject to the constraint 1-1 J J K E x. = M , (5.2.24) 70 we find that /C x. = M - 2 /C. i=l L is a necessary condition for the minimization of T subject to (5.2.24). Evidently, for the integer case, if we let A(C.) /var(f.) M. = int 1 M , K (5.2.25) S /var(f.) A(S.) i=l X X then (5.2.22) will be close to its minimum. 5.2.3 Multi-Stage Sampling Scheme We desire to employ equation (5.2.25) for the evaluation of n_,,(r,M) in equation (5.2.17). Unfortunately, the variances (5.2.21) are not known. We therefore perform sampling in stages. In each pass, an estimate of the variance in each region S. is made, and from this inform- ation a new distribution of M. is assigned and the process continues, the estimations of M. approaching (5.2.25). The total variance will not be at a minimum for any particular M but is expected to be close to a minimum especially for large M. This multi-stage sampling procedure will now be described in detail: Consider the L-stage sampling procedure for the estimation of K 1 Z A(S j ) A(S~yJ 1 d Qj d ^ f ?>Q) (5.2.26) J — - 1 - J c J 71 by K L J (*) (*) n( X J) g " S A(S ) — E S f(^ ',Q j=l J M j X= k=l ~ k ~ k (5.2.27) where M. = £ M (X) 4=0 J We calculate g by the following algorithm: 1) Choose the initial sample distribution, K<°>, I M <°>=M <0 > . 1 i-1 J 2) Begin sampling with the initial sample distribution and calculate M (0) j k=1 ~* ~* and 0<» - Z J [f( 5k (0) ,q k <°'J>)] 2 3) Estimate the standard deviation for each region by a. = A(S.) J J M. 1 r (2) 1/2 4) Calculate the new sample distribution by: M (1) .. ** 'i j L K £ CT i-1 : 5) Make the next pass (X = 1,2,...L): M<*> j j k=1 ^ „k 72 M (I) g. (2 >->g? 2 > + i CfCS^.Q^)] 2 k=l k '^k a. A(S ) " l G U) _U'=o j ' 2n G. Z M i'=0 J (X') 1/2 M <*-l) M A L K S a. 1=1 ] 6) Repeat 5) until 4 = L and finally let K G. A(S.) g - s M, a = g K S j-l a. M. 1/2 (5.2.28) where O is the approximate standard deviation of g. O The error of the estimation of a has many origins: 1) Finiteness of M sample (variance of variance estimate). 2) Lack of independence of regional estimates G.. 3) Lack of independence of path variables due to the Markov chain. Hence a is not to be taken too seriously but serves as a rough indicator of the accuracy of the estimate 73 5.2.4 The "Correlation" Term, n , (r) The procedure outlined in Section 5.2.2 for the calculation of 1 2 n , (r) does not yield satisfactory results for n , (r) as defined in equation (5.2.13). The following procedure yields smaller variance in the result, particularly as r becomes large. We may re-express (5.2.5b) in the following form: n , (r) = 2TTW^(r) JJ d(T dQ < (^— - 1) (e ZJ +e iJ )> s , ' Z W*(r) (5.2.29) where we have used the relations \e \e )) = {{e )e >, i = 1,2, W*(r) = . We recall that the region S' = SVS^, is in the first (Q > 0, Q > 0) U y z quadrant. It will prove useful to extend the integral to include the image with respect to the Q -axis. Let us use the following notation: S S [ % > °' Q z > °> Q y + (Q z " I r)2 ^ V <" S ) S H KJy > °'% < °' Q y + (Q z + 2 r)2 ± V S = compliment of S n with respect to the first quadrant in the Q ,Q plane y z r S^ = compliment of S n with respect to the fourth quadrant Then, S 1 = S and 74 2 n s ,(r) = 2TTW JJ ... = nW JJ since the integrand is symmetric with respect to the Q -axis. The integral -3V 23 (5.2.29) may be separated into two equal integrals involving e and e respectively. These two integrals may be combined to form: -pv 12 _ p - 2 n (r) = 2TT W* . b Z — — y Z W*(r) s o us o + The region of integration excludes S . Over the region, S n , the term, -pv 13 e , is considered to be zero for most paths so that the integral will be unaffected by the inclusion of S n . Hence we can integrate over the right half -plane except for the region S n . Furthermore, by choosing the origin for Q to be located at the position of particle 2, the integral assumes the simpler form: 1 2 / x o t T B/ n r ^2,_ r» d(cos 8) n , (r) = 8tt W (r) J Q dQ J — ^-z L R -i X E{(P(Tj(T),r) - l)FCn + ( T ),Q + r)|Tl(l) = q(l) = 0} where F (Tj (t) ,Q + r) was defined by equation (4.2.6), and P(T|(T),r) = B -1 [W (r)] F(T](T),r). Finally, applying the multiple integral approximations, we have 75 ! „ s ,(r) = 8n U *( r )J dQ/Mf^Jd, R o - 1 , 3n 3n , . X J ^rJ ^r (p D <*•»•*> - 1} X D + (x,y,u,v,Q + R) + 9 (-5) (5.2.30) n where D (x,y,u,v,Q ) was defined by equation (5.2.10), and Q = (0,Q /l - cos 6, Q cos 6) 2 Equation (5.2.30) is suggestive. It implies that n„, (r) is the integral over Q of the average difference between D with and without the influence of the probability weight P . The variable, Q, cannot be considered as part of the average, since the integral over the product space which includes Q does not exist. This is seen by changing the order of integration so that Q is integrated before any of the paths. The resulting integral is infinite for an infinite number of paths, and hence the iterated integral diverges. By Fubini's theorem, the integral over the product space does not exist. However, the space of cos 9 may be considered as part of the product space, (cos 9) XxXyXuXv, since it is bounded and the integrand of (5.2.30) is finite. Let us now express (5.2.30) as a Monte Carlo estimate: R- r R. 2 n s ,(r) r 2 n s „(r,M) = 8tt J dQG M (Q,r) (5.2.31a) where ,, G (Q,r) = - I [D (x(k),y(k),u(k),v(k),Q(k) + r) k=1 - D + ( X (k),y(k),u(k),v(k),Q(k) -1- r)]Q 2 , (5.2.31b) 76 and the subscript, S" , indicates that we are approximating S' by a truncation of the upper limit on Q at R_. The distributions of the random variables x(k), y(k), u(k), v(k) were discussed in Section 4.3. The random variable, Q(k), is formed from ,Q /l - a Q(k) = (0,Q /l - a(k)% Qa(k)) where a(k) is uniformly distribution on [-1,1]. xOO has distribution, du , and can be taken to be identical to the "trial" path of Section 4.3. nX The advantage of this procedure is that if we take the value of 5 in equation (4.3.20) equal to zero, the value of xOO will equal the value of x(k) for all "trial" paths which "succeed." In such cases, the summand need not be evaluated since it will be zero. It is apparent that as r becomes large (relative to X) the probability density, P (x,u,r) -» 1 , and more and more paths will succeed. Hence relatively little computational effort is expected for larger values of r. 5.2.5 The "Superposition Approximation" Finally, we discuss an approximation which reduces the calculation of rL (r) to an integral of two-particle density matrix elements only. This approximation was used by Fosdick and Jacobson [33,34] and has proved to be useful for large r and high temperatures. In this approximation, the three-particle term: W 3 = (5.2.32) in equation (5.2.1). The factorization is expected to be valid providing 77 the correlations between the terms are negligible. To see qualitatively how this might come about let us consider the particles 1, 2, and 3 moving in random paths for the various values of Q in the integral. As the temperature becomes large, the paths shrink to the points 1, 2, and 3 and the correlations become small simply because the particles are not moving much. Also as r becomes large, particle 3 tends to interact with either particle 1 o_r particle 2 but interacts weakly with both simultan- eously. That is, due to the exponential nature of the Brownian paths one can consider spheres about each point of a certain radius, R, = 0(X), in which "most" paths occur. Thus the chances that the particles move out of these spheres are small. If the range of the two-particle potential is r_ we may consider concentric spheres about these spheres of radius r n + R, . (J A. This bigger sphere is the region within which two particles can interact. Now if r > 2r n + 4R, the sphere of particle 3 cannot overlap with both that of particle 1 and of particle 2 for most paths. Hence the important -8 (V + V + V ) -8(V -!- V ) ,. fc . P ^ 12 23 13 ; . t . P ^ 12 23 ; contributions to e are either e or -6(V 12 + V 13 ) e and V- _ is close to zero for most paths so that correlations -3v 12 -3V. 3 between e and e , i = 1,2 are small. Replacing (5.2.32) in equation (5.2.1) we have in the "superposition" approximation, CO CO nP-\r) r n (r) = tt f dQ 2 f dQ W(w"'"-1) (W _ -l) . d sp J y J z Let us consider, y " z - 00 g s P (r) s -wfer = n ? dQ y ? dQ z f (Q+)f (Q_) (5 - 2 - 33) -co 78 where + + f(Q') = 1 - W" . As in Section 5 . 2 we can restrict the integration to the first quadrant of 2 the - Q plane, consider a "hard core" region S Q in which f (Q ) = 1, transform to the Q ,Q coordinate system, and truncate the integral beyond + a circular region of radius R, in which f (Q ) ~ 0: g (r) = — JJ Q + dQ + Q"dQ" f (Q + ) + £- JJ Q + dQV dQ~f (Q + )f «}") SP r S (5.2.34) = g. (r) + g q (r) s o s i where S Q = {Q + > Q": Q" < R Q ,Q + + Q" > r > Q + - Q _ ] and S-l = 'Q + > Q": R q r > Q + - Q - } . Nov, a table of f (x) can be calculated by Monte Carlo sampling as suggested in Section 4.2. The first integral, g (r) can be expressed as a one dimensional quadrature, and the second can be evaluated by an iterated quadrature method. The details are discussed in Appendix A. The effectiveness of the "superposition" approximation is demonstrated in the following table of g (r) compared with path integral calculations sp done by Jordan [9] not employing the "superposition approximation," for T = 5 K and 10 K. Other comparisons are found in L33]. Table 1 Comparison of n (r) Obtained by Superposition Approx- imation with Path Integral Calculations of rL (r) 79 T = 5°K T = : 10°K r n sp tf } < r) n sp r, (1) n D (r) 1.0 .22 .05 + .02 .30 .10 + .02 1.3 .20 .01 + .07 .04 -.18 + .07 1.5 -.02 .14 + .07 -.40 -.42 + .05 1.7 -.12 .22 + .06 -.51 -.52 + .04 1.9 -.11 .14 + .04 -.45 -.50 + .03 2.1 -.04 .12 + .03 -.29 -.27 + .02 2.3 .06 .06 + .02 -.10 -.10 + .02 2.5 .14 .16 + .02 + .05 + .08 + .01 2.8 .19 .21 + .01 .130 .132 + .006 BO 5.3 The Pair Exchange Term According to equations (2.4.6-7), the contribution to the pair distribution function due to the exchange of two of the three particles is split into two integrals n (r) and n (r) . The second involves the 1 2 exchange of particles 2 and 3; and in the Wiener integral formulation a (Q~) 2 Q factor of e occurs which can serve as a weighting function for the integration over the third particle as in the cyclic exchange term. The n p (r) integral however involves the exchange of particles 1 and 2 and 1 _ r 2/g hence contains the factor e , and so the convergence of the integral over the third particle (Q) depends upon a cancellation between two- and three-particle terms. Thus n (r) is calculated by a method similar to 1 that for il (r) , and il \r) is calculated by a method similar to that for n (r) . It is clear that il (r) should be more difficult to cal- 1 (1) culate than il, (r) . Moreover for the purpose of calculating the third 2 virial coefficient it is not necessary to have both pair terms. For these reasons, we will not be concerned with the calculation of n^ (r) . 1 We eliminate ru, (r) by considering the third virial coefficient 1 term, a, = d r (n^ (r) + tl, (r)). Combining the integrals over r , 12 3 P using equations (2.4.6-7) , a, can be expressed c£ = J* d 3 r 2 d 3 r 3 {[\ 9 <2 1 3|e 3 |l 2 3> - \ 6 <2 l|e" 'l 2> X (\ b <2 3|e l \2 3) + X b <3 2|e 2 |3 2»] + 2[\ y X (X b r«.* Since the resulting integrals over r and r of each bracket are independent of r.. , r , and r , the variables r.. ,r ,r are dummy variables and the integral of each bracket is invariant under the trans- formation r.. < > r . Hence the two brackets have equal integrals and we may write a l = U d ' r ^^^ + ^ d3r 2 d3r 3 W E (r 12 )(W D (r 31 ) - w E (r 23 ) " 1) (5.3.2) where ,6, - PH 2 W n (r..) h \° —23 —31 —21 -SV -SV -SV becomes \~ 9 (5.3.5) then a^ = 3 J d 3 r(v(r) - a) + 2a(2b-a) (5.3.6) where a = J d 3 rW £ (r) (5.3.7) b = J d 3 r(l - W D (r)) . (5.3.8) In terms of Wiener integrals the three-particle term is ,3 " R (Q } V(r) = J d J Q e P X E { F [5»Q + ;T]]F[Q",-Q";7]"]F[Q + ,r;5 + ]|5(l) = q(l) = 0}. (5.3.9) Since the integral is over all-space, we can integrate with respect to Q = Q - — r just as well. Using equation (4.2.9) and expressing (5.3.9) in terms of spherical coordinates: Q , a - cos 8, and we have (dropping the superscript on Q ) V(r) - 2n J q 2 dQ J da ." » ' Jd , xnJ%nI ^!j^! 0-1 '22 xexp[-| S J (v\l + V^ + V^)dr] +0(^) (5.3.10) i=l o n 83 where V^ h V (/23 (i- p(u, T ) + x , + t(x -x )) /n ~~ ~i i 111 + r + ( T+ f" L )Q) (5.3.11a) V^ = V(/2p (7- P "(u,v,t) + xT + t(x' - xT )) /j /n ~ ~ ~ ~i-i ~i ~i-i + Q - 2 ( :L V i )Q) (5.3.11b) V 13 ~ V(/2P ( 7~ £ +( "'y' T) + *i-l + r( *i " *i-l )} + Q + r - ( I± | li )Q) (5.3.11c) and Q = (0,Q/l-a 2 , Qa) . The Q integration can be normalized by expressing the exponential factor 2 in terms of the x distribution with three degrees of freedom: _ 1 _ 3 1 2 T(Q 2 ) = 2 tt 2 p 2 (Q 2 ) 2 e" Q /P . (5.3.12) Then 2 „ J Q 2 d Q J da e"^ . . . . ^TJ J r(Q 2 )dQ 2 f f» . . . 0-1 2 -1 We next choose as a weighting function for importance sampling, n 1 P^ n) (x,y,u,v,Q 2 ,a) = \~ exp[- | 2 J V^ 2 d T ] (5.3.13) P 1=1 84 where 0-1 ll n 1 X exp[- | S J V^3 dr] . (5.3.14) n i=1 Hence and J T(Q 2 )dQ 2 ... P< n) (x,y,u,Q 2 ,a) ^ J dP^ n) = 1 *<*> ' 372 z p I dp p n> «pC" I .=, J < v u + v n' dT ] 2 1=1 + 0(-~) . (5.3.15) n 2' Now, it is easy to verify so and where 3/2 Z p = ^-3— a '+ O(^) (5.3.16) \ n V(r) = a U(r) + (^) (5.3.17) n a* = 12 tt a J r 2 dr[U(r)-l)] + 2a(2b-a) (5.3.18) n 1 U(r) = J dP< n) exp[- J E J (v[* + vJ*)d T ] (5.3.19) i=l M n 1 ~. ± S exp[- I 2 J (Vj~(k) + V^(k))d T ] (5.3.20) k=l i=l 85 is calculated by Monte Carlo sampling as described in Section 4.3. In this 2 case all random variables x(k), y(k), u(k), v(k), Q (k) , a(k) are obtained by importance sampling via a Markov process based on dP . 5.4 The Path Integral Programs Path integral computations were performed on an IBM System/360 model 75 computer by means of four FORTRAN programs which we shall identify as C, D , D 2 , and P corresponding to equations (5.1.9), (5.2.17), (5.2.31a) and (5.3.20) respectively. These programs were based on the methods developed in the previous sections. In all programs, the Markov chain parameter, 6 (see equation (4.3.20)), was taken to be zero, and the same random numbers were used for all values of r for a given run. Although this tended to correlate the results with respect to r, it presumably saved con- siderable computational time. The initializing vector §(0) for the Markov chain mentioned in Section 4.3 was determined in the following way. First, the arbitrary value of zero was assigned to all path components; next a subroutine used for determining the next vector in the chain was called repeatedly for a number of times, M„, until it could be assumed that: (a) a "trial" vector had "succeeded" for all values of r (U(§(M )) would then be defined in the program); (b) the probability distribution of §(M n ) would be close to P(^(M )) (equation (4.3.7)). The chain was then continued, being used now in the Monte Carlo sampling. The initial process in which the chain is generated without sampling is referred to as the "initial relaxation" of the chain and tends to reduce the variance of the result. Since the states generated by a realization of a finite Markov chain tend to distribute about states corresponding to a particular relative maximum of the probability distribution it is possible, in the case that 86 there is more than one relative maximum, that the state space will not be covered adequately by only one Markov chain. To avoid this situation, the averages over at least ten independent realizations of finite Markov chains were taken, each chain being allowed an initial relaxation before sampling was taken. This method has been suggested by several authors [31,35]. An estimate of the standard deviation of the result based on the independent chains was calculated in programs C and P. In the case of programs D. and D 9 , each independent realization corresponded to a stage of the multi-stage sampling procedure described in Section 5.2.3, and the standard deviation was estimated according to equation (5.2.28). 87 6. RESULTS AND DISCUSSION 6.1 Tests of the Path Integral Programs 6.1.1 Harmonic Oscillator Tests The programs C, D , D , and P were tested by temporarily assuming the two-body interaction to be a harmonic oscillator rather than the Lennard- Jones potential (2.1.2). That is, we let V(r) = | r 2 . (6.2.1) With this assumption, the appropriate three-body density matrix elements can be expressed as products of one-dimensional density matrix elements of the form: = p(X,X',3> (6.2.2) z which are solutions of the Bloch equation, * 1 * 2 "vX ? (fo - ■£- Z-z + -r— X ) P(x,x' ;P) - ** 2m o * x subject to the initial condition: £im d(x,x' 53) - S(X-X') (6.2.3) 3 - and the boundary conditions lira o( x ,x' ;P) = o . X "• °° The solution to the above equation can be determined and is given by °(x,x';3) = [^ sinh(^)]- 1/2 ™ ^ 2 2 X exp{- — [coth(u£)( X + x ' ) - 2XX' cosech(UUg) ]} . (6.2.4) 88 Therefore the pair-correlation functions can be evaluated exactly as a check against the path integral calculations. Some of the results for the potential (6.2.1) are listed in Tables 2 through 8. Program D^ was altered slightly so as to calculate I/ d 3 r 3 W^(l,2,3) = J d 3 Q E{|fV|T|(1) = q(l) = 0} rather than 2 n (t } -f — = J d 3 Q e{| (F + -i)(F'-i)|-n(i) = a (i) = o} W (r) S' to facilitate the analytical calculation. The value of the "hard core radius," R n , was set to zero. The numbers after the + signs in the tables are the computed standard deviation. It will be noted that in Table 7 the two independent calculations for n = 4 agree to within the estimated standard deviation, whereas in Table 3 for n = 8 the numbers agree to with- in two standard deviations. The reason could be due to the fact that the standard deviation estimates in programs D.. and D are based on correlated samples and are therefore inaccurate whereas estimates in programs C and P are based on independent samples . In each table except Table 4, the effect of increasing the number, (n-1) , of "break points" of the paths is shown to be a definite decrease in error as expected. Results for largest n are seen to agree with the exact values to within the estimated standard deviation with the exception of those from program D-. Errors in program D have a third source besides the finiteness of n and the variance of the sample: namely, the quadrature error in the integration over Q. (A trapezoidal rule was used in this integration.) The function, G (Q,r), defined by equation (5.2.31b) is a 89 Table 2 Comparison of Calulations of X~ [1^(1,2)]" f d 3 r w!?(l,2,3) Obtained from Program D. and by Exact Expression Assuming Harmonic Oscillator Potential. (T = 0.6, y = .3/T) r n = 4 exact 0.1 (.136 + .004) x 10" 1 .132 x 10 -1 0.5 (.129 + .004) X 10' 1 .129 X 10 _1 1.0 (.114 + .003) X 10' 1 .118 x 10" 1 2.0 (.85 + .02) X 10" 2 .86 X 10" 2 3.0 (.47 + .02) X 10" 2 .50 X 10" 2 Table 3 Comparison of Calculations of X [W (1,2)] d r„ W (1,2,3) Obtained from Program D.. and by Exact Expression Assuming Harmonic Oscillator Potential. (The two columns for n = 8 represent independent samples.) (T = 0.4, y = .3/T) r n=2 n=4 n=8 n=8 exact 0.1 (.74+.03)Xl0 _3 (.81+.03)X10" 3 .802X10" 3 1.0 (.57+.02)xl0 _3 (.64t.02)XlO" 3 ( .66+.03)xl0" 3 ( . 71+.02)xl0 _3 .727X10" 3 2.0 (.46+.02)XlO -3 (.42+.02)XlO -3 ( .48+.02) XlO" 3 ( .51+.02)xl0" 3 .539xl0" 3 V- Table 4 2 B Comparison of Calculations of g (r) = n , (r)/W (r) Obtained from Program D with Exact Expression Assuming Harmonic Oscillator Potential. (The two columns for n = 2 represent independent samples.) (T = 0.6, y = ,3/T) g D (r)A 3 r n = 2 n = 2 exact (,38+.03)Xl0" 2 (.24±.03)xl(f 2 .436xl0' 2 2.5 (.11-K1)X10~ 2 (.21+.07)Xl0" 2 .285X10" 2 91 CO H I W < .-I CO cfl • > 11 > 4-1 o « U3 -i «rl o c i-i 01 PL, J_> o o W O 4-1 4J >— I CO 3 ^ CO i— i CN 00 o m i—i O CO uo i—i i— i X co CO CN QJ o o O 1—1 o O 1 .—i i—i o X X 1—1 s-^ /— V X r-~ m — i i— i CO 11 U •H CO i-> >- CJ C QJ T3 0) 4J C 4-1 O CO c PU o -a- 2 i-i • o o I 4J o (0 II M 1—1 U-l — i H •i-l v -' T3 O 0> CO c O '-N •i-i • CO cj en •u •r-l QJ J2 C -I o O CU P E ^N J-i CO u cO co ^-^ X ? 4-1 r^ ,— i 00 c i C QJ QJ CO •H 13 i— I E c J3 T3 3 QJ CO C CO (X H CO CO QJ < -o /— \ c ^i ^ -H v— ^ i—l ID QJ 4-1 > C 14-1 •H QJ O 4-1 CO O QJ CO QJ U c & CX o CO QJ •H QJ U 4J pei CO sfr i— i C D O II o •H 1—1 co C CO CO O QJ U U O 14-1 CX H-t o c3 CO c g o *-> 6 CO O 3 •H CO i—i M X o CO W CJ a g "O o o C 5 u CO 4J CNJ CM CM cm CM 4J o o o O o CJ 1—1 r-l 1— 1 i—i 1—1 m X X X X X X o\ 1—1 1—1 vD ro 01 00 CNJ ^o s-^. 11 h* r- CM CM CM CM CM 1 CM 1 o 1—1 O o O O o X i—l 1—1 i—l i—l 1—1 /~\ X X X X X vO /~*\ e-N •-N /"■" \ •— s vO m -N /*"S o m St r*- m 1.5 at both temperatures. There seems to be a consistent difference between the two functions for T = 1.0. This is perhaps not entirely due to the correlation effect J.1 98 A « D (r)/x3 " l - 7 ° K - «. p /X 3 ■J l W 3 7.00 — ♦ 8.00 9.00 — I 10.00 Figure 9. Direct and Exchange Contributions to n^ 1} (r)/\ 3 at 1 . 7°K 102 The three-particle term, n (r), is plotted in Figure 10 for T = 1°K and T = 1.7°K and compared with values at T = 5°K and T = 10°K [9]. It is observed that the function becomes rapidly large with de- o o * creasing temperature, and the relative minimum exhibited at 5 K and 10 K at r ~ 1.8 becomes a slight inflection near the relative maximum at 1.7 K and 1 K. The relative minimum has been interpreted as the valley between the two peaks which correspond to preferred positions of three particles arranged linearly to be within each other's potential well [9], (Particles located at r = 0,1.2,2.4) As the temperature becomes lowered, we expect the distribution of the positions of the particles to become more spread out (X becomes larger) so that the region between the preferred positions becomes an overlap region perhaps accounting for the appearance of the peak at lower temperatures. One would expect the peak in n (r) to correspond, on the other hand, to a triangular arrangement of three particles since it would require the least kinetic energy for the cyclic exchange of the particles. This arrangement could also account for the peak in ru (r) at the low temperatures . The computed values of the third virial coefficient, C(T), obtained via equation (2.3.5) are in qualitative agreement with experimental measure- ments as shown in Figure 11, i.e. they are negative and decrease with de- creasing temperature. The experimental data however are rather uncertain and marginal at best. The computed value at 1.7 K agrees with a smooth extrapolation from Keller's data [37,38], although values extrapolated from Keesom's data [39] seem to disagree with the calculation by more than an order of magnitude. The large relative error in this calculation is 2 due to a large cancellation from 0L - cc n . It has been shown [9] that n 2 (r) exhibits this dip at r ~ 1.8 for the temperature range 5°K < T < 273°K. 103 10 6.0 7.0 8.0 — 1 — 9.0 10.0 (1) 3 Figure 10. n (r)/X as a Function of Temperature 104 8 T (°K) - 80 — I 1.20 1.60 2.00 X2.40 2. 80 X 3.20 3.60 TT -** 4.00 4.4G 4.80 °g o o Path Integral calculation* ■Caller (1955) [37] Keller (1955) [38] Keetom and Walatra (1940) re-evaluated by Keller (1955) [37] KiestemaWer and Keeiom C1946) [39] Keeaom and Kraak H935) t40] Figure 11. Computed and Measured Values of the Third Virial Coefficient, C(T), as a Function of Temperature 105 A theoretical calculation due to Larsen [7] based on the binary collision expansion method of Lee and Yang [5] yielded the result, C (1.7) = -2.2 X 10 cm /mole which deviates by over an order of magnitude from the result of the present path integral calculation, C (1.7) = -8.8 x 10 cm /mole . In addition to using a low temperature expansion, Larsen' s calculations were based on the following assumptions: (1) that the two-body interaction potential is given by a square well with a finite repulsive core, (2) that there exists a three-body bound state, and (3) that the binding energy could be approximated by neglecting the repulsive core. Assumption (1) can account only for a minor difference between C (1.7) and 1j C (1.7). Assumption (2) is quite important, requiring the addition of a tern proportional to e where E is the binding energy; and it is possible B that the error resulting from assumption (3) is considerable. 6.3 Conclusions In summary, the results of the previous section suggest that the exchange contributions to the pair distribution function for He gas are almost negligible at as low a temperature as 1.7 K but are significant at 1 K . The third virial coefficient, however, due to a large cancellation, is sensitive to three particle exchange effects at T = 1.7 K as well as at 1 K. It has also been shown that three-particle contributions to the pair distribution function not involving exchanges can be reasonably approximated in terms of two-particle diagonal conf igurational density matrix elements through the "superposition approximation" for temperatures as low as 1 K. I 06 The success of the superposition approximation suggests the feasibility of a so-called "expanded superposition approximation" of the form: = J d x x ... J d x^ n + E „ where and x. (j ' k) = (X. ., X. .) ~i ~i,j' ~i,k N is a 6- tuple consisting of the components of x. corresponding to particles j and k. By expressing the density matrix elements in terms of Wiener integrals and expanding the exponential functionals in powers of the paths, one can show that the error term, E , of this approximation is 0(l/n 2 ). There seems to be no obstacle to the attainment of further path integral calculations of the third virial coefficient at lower temperatures except an increase in the necessary computing time. It might be possible that, even below 1 K, a certain reduction of effort could be attained through approximation of the direct term by the "superposition 107 approximation." It should be relatively straightforward to extend the path integral approach to the case of non-additive potentials; and in principle, the path integral formalism could be employed for calculations involving four or more particles, although computing time would be a severe limiting factor. With the use of coming parallel processing computers, this limitation can probably be considerably eased. APPENDIX A QUADRATURE FOR g_ (r) , g Q (r) s o b i We have, by equation (5.2.31), 108 g SP (r) r g q (r) + g q (r) where g s (r) = ^U Q + dQ + Q"dQ'f(Q + ) g s (r) = & JJ Q + dQ + Q"dQ"f(Q + )f(Q*) 1 S, and S„ = .+ .+ [Q > Q" : Q" < R Q , Q -!- Q" > r > Q - Q - } = {Q + > Q" : R Q < Q" < R L , Q + + Q _ > r > Q + - Q - } We consider first g (r) . It can be verified that r+R r |2 J" xdx(RQ-(x-r) 2 )f(x) + 2tt[| Ro-Rq(|) + |(f) 3 ] g s (D = < R r+R, 2TT r xdx(RQ-(x-r) 2 )f(x) if 2 < R if 1 * R o (A.l) r-R, where use is made of the fact that f (x) = 1 for x < R_. The integrals in equation (A.l) can be evaluated by a modified trapezoidal rule making 2 2 use of the weighting function x(R_-(x-r) ). Consider then, the integral over the interval [x, ,x ]: I(X V X 2 ) - J xdx(R 2 -(x-r) 2 )f(x) and let us Linearly interpolate f(x) within this interval, i.e., 109 x - x \ f x - x, f (x) r | — f (x.) + 1 x„ - x„ 1 X„ - X, '2 1 '2 1 f(x 2 ) Integrating, we obtain, I(x 1 ,x 2 ) = R " "' ^(x^xp +|0 2 ( Xl ,x 2 ) 20 3 (X 1' X 2 } X (f( Xl ) - f(x 2 )) + R 2 - r 2 10 1_ , 2 2, 2 , 3 3, 1,4 4. - (x 2 ~x 1 ) + - r(x 2 -x 1 ) - 4^ x 2~ x l' ) f(x 2 ) where 2 2 J 1 (x 1 ,x ) = x 1 x 2 + x 2 - 2x L 2 2 3 3 2 (x is x ) = x^ + x x x 2 + x 2 - 3x L * r ^- 3 ^ 22 j 3 ^ 4 3 l' X 2^ = X 1 X 2 X 1 X 2 X 1 X 2 X 2 - 4x, Then assuming a table x., f(x.) for i = 1 , . . . ,N is known, we can express the integral over an arbitrary interval [a,b] as: b I(a,b) = J xdx(R 2 - (x-r) 2 )f(x) V 1 E I(x. ,x ) + I(a,x. ) + I(x. ,b) i=j J l J 2 110 where x. , < a < x. and x. < b < x. ,, and f(a) and f(b) can be deter- mined by linear interpolation. Next we consider g (r) . There are three cases. b l Case I: ^ < R Q R 1 +R Sc; (r) = ^ J Q"dQ"f(Q") J Q + dQ + f(Q + ) - 1 R Q * Case II: R < f < R L . r/2 r-hQ" < r > - f 1 J Q'dQ'f (Q") J" QdQ + f(Q + ) R r-Q" 4tt n 1 - - - V r + + + + f 1 J Q dQ f(Q ) J Q + dQ + f(Q ) r/2 q- g s. Case III: R. < | R l r+Q~ g s (r)-~J Q'dQ'f (Q") J Q + dQ + f(Q + ) 1 R r-Q" The above integrals are of the form r 2 x+r ^J" xdxf(x) J* ydyf(y) . (A. 2) r, a (x) Define accordingly, i|f(x 1} x 2 ) = [x f(x 2 ) + Xj^f (x 1 )](x 2 - x^ x 2 which is the trapezoidal rule approximation to 2 J xf (x)dx. Then (A. 2) X l can be approximated by: Ill J2" 1 J { S [(g(x.)Y(x.) + 8(x 1+1 )Y(x i+1 ))(x i+1 -x t )] i=J + (g(r )Y(r ) + g(x )Y(x ))(r,-x ) ^ ^ Jo J 9 ^ Jo + (g(x. )Y(x. ) + g(r )Y(r ))(x. -r )} J l J l J l where g(x) = xf (x) , X J 1 " 1 - r l r > Q + - Q"; Q + > Q"} . Now the sets, S. , can be expressed in terms of two parameterized sets as follows S l ~ S I (R 0' R 1 ) S 4 " S II ( R » R 1 ; R 2' R 3 ) $2 = Sjj (Rq>R^'jR-^)R2^ ^5 = ^ii ( R ^> R 2' R 2' R 3^ S 3 = SjC^.Rj) S 6 = S l( R 2 ,R 3 ) Sy = S (R ,R ;R ,r + R«) where if R x > R 2 or R 2 < | S I (R 1 ,R 2 ) = [Q < Q" < Q 3 ; Q 2 < Q + < Q 3 1 n S otherwise where and % = h> Q 2 = Q 3 = R 2 113 !•*! R x or R > R or R^R, > r or R3+R, < r S II (R 0' R 1 ;R 2' R 3 ) " + [Q < Q" < Q 1 ; Q 2 < Q < Q 3 ] n s o therwise where r - R 3 , R Q + R 3 < r % = \ R 2 " r ' R 2 " R >r R_, otherwise Qi- \ Q 2 = *3 = r - R^ R x h R 2 < r R , otherwise r + R 1S R_ - R. > r R„, otherwise Th +.2 e normalization factors, a{S.} = J d (Q ) d(Q), can be calculated by determination of A{s (R ,R )} and a{s (R^LjR R ) 114 As an example of how the A's can be determined we consider S (R. ,R ) as shown in Figure 12. The shaded area, S (R- ,R ), can be regarded as the intersection of the rectangular region (S) and the isoceles right triangle in Q ,Q space. All possible "areas" can be built up by a superposition of "elementary areas," A{cr.}, e.g. in the figure, A{S I (R 1 ,R 2 )} = Ata^r-R^} + A{a 2 (R 1 ,R 2 )} where Figure 12. The Regions, S, S and a. (S s (Q + ) 2 , Tl ■ (Q") 2 ) ,2 aCqjCR')} = J d| J dT| = rR' Z (| R'-r) - \ 1 4 r ,r.2 , .1/2.2 2 R 2 F ? R - (r-R ) A{a 2 (R 1 ,R 2 )} ■ J 2 d§ | dTl ^ — « 1 R 2 (R 2 (r . Ri) 2 (r-R x ) R 2 1 v 2 115 The calculation of all possible regions S (R R ) and S (R R ;R R ) is achieved by the determination of a total of nine "elementary areas." They are (including the two just defined), f K /2 > 2 8 3 Ata,(R)} - J dT| J = f rR Cr-lf /2 > 2 2 2 (R 2 -rT R 2 A{a (R R )1 = J" dTl J* d| = (R 2 -r 2 )[(R -r) 2 -R 2 ] 2 1/2 2 R^ (r+Tl 17 ) |[(R 2 -r) 4 - Bj] - | r[(R 2 -r) 3 - rJ] 2 2 R 2 R 2 A{a 5 (R 15 R 2 )} -J dTl J d? = \ (R 2 - R^) 2 R 2 1 r//2 R/2-y A{CT (R)} = 2 j dy J dx(x -y Z ) ://2 + 1 - 1 where Q = — (x + y) Q = — (x - y) /2 /2 = Afj^R)} 2 2 R 2 R 3 A{a 7 (R 15 R 2 ,R 3 )} = J dTl J 2 1/2 2 R^ (r-Tf'V = (R 2 - r 2 )(R 2 - R 2 ) + | r(R^ - R^) - | (R* - rJ) 116 2 2 R 2 R 4 A{a 8 (R 1 ,R 2 ;R 3 ,R 4 )} -J dT\ J d§ = (R* - R^CR* - R*) R l R 3 2 2 R 3 R l A{a 9 (R 19 R 2J R 3 )} = J* d§ J" = A{a 7 (R 2 ,R 3 ,R 1 )} . R* (? 1/2 -r) 2 Flow charts for the calculation of A{s (R ,R )} and A{S (R ,R, ;R_,R~)} in terms of A{a.} are exhibited in Figures 13 and 14 117 c* II cr pU c f> _ cr i .k 6 VI \ m \ < M + / < £/ cr + o + < + O + cr i < + o cr + cr i o cr < c >% o 1—1 (V c o %T c >, o ,0 •r-l 4-1 a) Q y— \ o ■H m CO X a 1—1 o M a u c QJ 0- id P g. .C < o i— 1 X d c 0) s -- / o •r-l o o 4-1 c ■r-l ^w* W o fr 0) a- 3 Oh C/j 00 r-l u CM O P-> a o CT> o vO O O CTv O CO O ■ +1 o o r- O ON i—l + 1 o o o o O O o o o o as to 00 00 00 i—l 00 vO • +1 ' + i—i 1—1 00