XI B HAHY OF THE UN IVLRSITY Of ILLI NOIS S\o.fcA- v\o. 2A3- 2.49 ceo. 2, The person charging this material is re- sponsible for its return 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 AAtf I 3 1372 -».-"*■»•» " SEB 1 6 Kecd L161— O-1096 Digitized by the Internet Archive in 2013 http://archive.org/details/threeparticleeff248jord 5J0. W I £0>^ Report No. 2U8 Tyuudi September 25, 1967 THREE PARTICLE EFFECTS IN THE PAIR DISTRIBUTION FUNCTION FOR HE^ GAS Harry F. Jordan AUG ij UNIVERSITY Of ILLINOIS THE LIBRARY OF THE •< » Report No. 2U8 THREE PARTICLE EFFECTS IN THE PAIR k * DISTRIBUTION FUNCTION FOR HE GAS by Harry F. Jordan September 25, 1967 Department of Computer Science University of Illinois Urbana, Illinois 6l801 This work was supported in part by + .he Atomic Energy Commission and was submitted in partial fulfillment of the requirements for the degree of Doctor of Philosophy in Physics, September 19^7- Ill ACKNOWLEDGMENT The author wishes to thank Professor Lloyd D. Fosdick whose professional guidance and encouragement were indispensable in preparing this work. Thanks for support in preparation of the manuscript is due to the Department of Computer Science and in particular to Mrs. Gayanne Carpenter and Mrs. Frieda Anderson who typed the text. The author also wishes to express gratitude to his wife whose personal support and encouragement did much to lighten the burden of preparation, iv TABLE OF CONTENTS Page ACKNOWLEDGEMENT iii LIST OF TABLES v LIST OF FIGURES vi 1. INTRODUCTION 1 2. THEORY OF THE PAIR DISTRIBUTION FUNCTION 3 2.1 Definition of the Pair Distribution Function in Terras of the Density Matrix 3 2.2 The Effect of Spin and Statistics on the Density Matrix ...... 6 2.3 Cluster Expansion for the Pair Distribution Function . ..... 9 3- PATH INTEGRAL FORMULATION FOR THE DENSITY MATRIX ... 18 3-1 Expression of Density Matrix Elements in Terms of Wiener Integrals 18 3.2 Upper Bounds on Off -Diagonal Matrix Elements . . 23 3.3 Approximation of the Diagonal Matrix Element by a Multiple Riemann Integral. .......... 26 k. EVALUATION OF g^r^) 37 4.1 The Multiple Integral Approximation for g 1 (r 12 ). • • 37 k.2. Monte Carlo Evaluation of the Multiple Integrals 49 4.3 Integration over the Position of the Third Particle ..... 62 5- NUMERICAL RESULTS 75 5-1 Results and Accuracy Considerations 75 5.2 Numerical Comparisons with Theory and Experiment 85 5.3 Conclusion 97 LIST OF REFERENCES 99 APPENDIX A 102 APPENDIX B 105 VITA. 129 LIST OF TABLES Table Page 1.1 Numerical Values of A.~ 3 g (r ;20,000), X " 3s e (r 12 ;20 ' 000) > and W 2 (r 12^ f ° r T = 5 ° K 110 1.2 Numerical Values of X J g (r ;20,000), \~ 3 g (r^^O.OOO), and W 2 (r 12 ) for T = 10°K. ... Ill 1.3 Numerical Values of X J g (r ;20,000), X~ 3 g (r 12 ;20,000), and W (r^) for T = 20°K, ... 112 l.k Numerical Values of A.~ 3 g (r ;20,000), X" 3 g e (r 12 ;20,000), and W^r^) for T = 30°K. ... H3 1.5 Numerical Values of \~ 3 g (r ;20,000), X _3 g e (r 12 ;20,000), and W 2 (r ) for T = 40°K. ... llU- 1.6 Numerical Values of X~ 3 g (r lo ;20,000), _ C Ld A." 3 g e (r 12 ;20,000), and W^r.^) for T = 50°K. ... 115 1.7 Numerical Values of A~ 3 g (r ;20,000), X _3 g e (r 12 ;20,000), and ^(r^) for T = 75°K. ... Il6 1.8 Numerical Values of X~ 3 g (r ;20,000), *-~ 3 g (r 1o ;20,00C), and W^(r no ) for T - 100°K. ... 117 e ' 12 ' ' ' ' 2. 12 ' 1.9 Numerical Values of X~ 3 g (r^;20,000), \" 3 g (r 12 ;20,000), and W 2 (r i2 ) for T = 273.l8°K. . 118 2 Numerical Values of the Third Virial Coefficient C(T) and the Cluster Integral b Q 119 VI LIST OF FIGURES Figure Page 1 Center -of -mass and relative coordinates for 3 particles ................. 38 2 f 13 (r 12 ,g)xf 23 (r 12 ,g) as a function of § for iVg = 1.1a and T = 20°K. 66 3 X~^g (r^;20,000) as a function of r^/a and T 76 k Computed third virial coefficient as a function of temperature 77 5 X g, and sampling error for T = 5 K 79 6 X g and sampling error for T = 20 K. 80 7 X g, and sampling error for T = 75 K 8l 8 kg-, and sampling error for T = 273.l8°K. ... 82 9 Comparison of X g., with classical and .t semi-classical values for T = 5 K. . . 87 10 -3 Comparison of X g, with classical and semi-classical values for T = 10 K. 88 -3 11 Comparisons of X g, with classical and semi-classical values for T = 20 K. ....... 89 -3 12 Comparison of X g, with classical and semi-classical values for T = 50 K. 90 13 The pair distribution function for average density p = .095 gm/cc. 92 ik The pair distribution function for average density p = .18^- gm/cc. . 93 15 The pair distribution function as a function of temperature for p = .184 gm/cc. ........ 9^ 16 Comparison of claculated third virial coefficient with classical and experimental values ...... 96 Vll LIST OF FIGURES (Continued) Figure Page IT A." 3 g x and sampling error for T = - 5°K • . 120 18 >." 3 g 1 and sampling error for T = - 10°K. . . . . . 121 19 *." 3 g L and sampling error for T = = 20°K, . . . . . 122 20 ^" 3 g 1 and sampling error for T = - 30°K. . . . . . 123 21 ^." 3 g L and sampling error for T = - 40°K. . . . . . 124 22 *.~\ and sampling error for T = = 5 0°K. . . . . . 125 23 *.~\ and sampling error for T = . 75°K. . . . . . 126 2k X" 3 g 1 and sampling error for T = = 100°K. . . . . • 127 25 ^" 3 g x and sampling error for T = - 2T3.l8°K. • • 128 THREE PARTICLE EFFECTS IN THE PAIR h DISTRIBUTION FUNCTION FOR HE GAS Harry Frederick Jordan, Ph.D. Department of Physics University of Illinois, 19^7 The diagonal elements of the 3 -particle density matrix h for He have been formulated in terms of Wiener functional integrals in order to compute both that part of the pair distribution function linear in the density and the third virial coefficient. The 3 - particle interaction energy was taken to be a sum of Lennard -Jones (12-6) potentials between pairs of particles. Evaluation of the Wiener integrals was carried out by a Monte Carlo sampling technique to yield the pair distribution function and third virial coefficient for temperatures from 273 K to 5 K. Comparison of the linear portion of the density expansion of the pair distribution function gave qualitative agreement with experimental values for the pair distribution k function in liquid He measured by neutron diffraction. The values of the third virial coefficient calculated by this method are consistently smaller than the experimental values and consistently larger than the classical values of this quantity over the given temperature range . 1 . INTRODUCTION In the past helium gas has "been a common substance for use in the comparison with experiment of the results of theoretical models and techniques in quantum statistical mechanics. At the time quantum theory began to be well understood the statistical mechanical theory of gases had developed further than the theory of either solids or liquids. The use of helium as a test substance arose due to the fact that of the monatomic gases, for which the application of statistical mechanics presents the least complexity, helium is the only one having small enough atomic weight to exhibit sizeable quantum effects in its thermodynamic properties at moderate experimental temperatures and pressures. The second virial coefficient B(T) has received a great deal of attention, and calculations of B(T) have been made using several different forms of the intermolecular potential (for a survey of much of this work see Hirshfelder, Curtis, and Bird [l]). The two methods most frequently applied to the calculation of B(T) are the expression of B(T) in terms of the phase shifts of 2-particle collision theory and the expansion of B(T) about the classical limit as a power series in Planck's constant h (the Wigner-Kirkwood expansion). The latter method is useful only for high temperatures since the series converges slowly, if at all, for temperatures less than about ^0 K in He . Few quantum mechanical calculations have been made for higher order virial coefficients. Larsen [2] and Pais and Uhlenbeck [3] have done work on the third virial coefficient. Both calculations have made use of the "binary collision expansion of Lee and Yang [h] , for which the convergence properties are not well understood. The Wigner- Kirkwood expansion is applicable to systems of more than two particles and can thus be used to compute higher order virial coefficients, but the calculations are more complex, and the convergence difficulties are the same as in the 2 -particle case. The phase shift approach, of course, is only valid for quantities which may be expressed in terms of 2-body collision theory. Recently two calculations of the trace of the 2 -particle density matrix, a 2-body property more general than the second virial, have been made in ways which admit of generalization to the N-particle problem for N > 2. One of these [5] involves the numerical calculation of the eigenf unctions and eigenvalues of the Schro'dinger equation for the relative coordinate and explicit evaluation of the sum over states involved in the density matrix. The other [6] makes use of a path integral formulation for the quantum mechanical density matrix and a Monte Carlo evaluation of the multiple integrals which arise in the evaluation of these path integrals. This thesis concerns the application of the latter method to the 3 -particle problem in order to evaluate both the 3 _ particle term in the pair distribution function and the third virial coefficient, 2. THEORY OF THE PAIR DISTRIBUTION FUNCTION 2.1 Definition of the Pair Distribution Function in Terms of the Density Matrix k To clarify the mathematical model to be used for He we assume that the system of interest is a gas of N identical particles in a volume V in the limit as N, V -+ °° with n = n/v held constant. The particles, the coordinates of which are denoted by i, i=l,2,...,N, are assumed to interact by means of an additive, spherically symmetric pair potential V(r. . ) , where r. . = | i-jj , so that the interaction energy for •^ J -L J the N-particle system is U(1,2,...,N) = | Z V(r ). (2.1.1) In the numerical calculations the potential function will be taken to be the Lennard- Jones (12-6) potential V(r) = k-K. l fe r - (? (2,1.2) where k = ik.Ok x 10 ' erg and a = 2.56 X 10 cm. This potential is taken because it is computationally simple and gives almost as good a fit to second virial coefficient data, at least, as any of the other k potentials proposed for He . The constants k and a are taken after de Boer and Michels [7] to give a good fit to high temperature second virial coefficients. The assumption of an additive pair potential is not essential to the theoretical development below. In order to produce results for more realistic potentials one need only replace the 3-partiele potential energy, U (1,2,3), in the final formulae by the appropriate 3-particle interaction. The pair distribution function n (l,2)d 1 d 2 is defined as the probability that one of the particles of the system is found in the 3 volume element d 1 at 1 while some other particle is simultaneously in 3 volume d 2 at 2. In the gaseous system under consideration, in fact in any homogeneous fluid, n (1,2) must depend only on the separation r, p between the particles except in the region near the boundaries of the system which becomes negligible as V -» <*. In a classical gas of non- interacting particles the probability of finding a particle in a given volume element becomes a constant independent of the positions of the other particles, so in this case /. Q x N(N-l) , , -v n o'l>V = ~p > (2.1-3) V lim n 2 U,2) - n 2 . (2.1.4) N,V-oa n = const , Quantum mechanically the situation is complicated by the exchange correlation between identical particles, but since this effect falls off as (see Section 3*2) exp[ -const, x (r A) ], Eqs . (2.1.3) and (2.1.4) still hold in the limit r, p » X, where is the thermal wavelength of a particle of mass m at inverse temperature £ = l/kT. If P N (1,2, . ..,N)d 3 l d 3 2 ... d 3 N is the probability that 3 particle i is in volume d i at i for i from 1 to N, then the pair distribution function is given by n 2 (l,2) = l(H)J...JP I (l ) 2,... ) N)d 3 3 d\ ... d 3 N. (2.1.5) The conf igurational probability P can be simply expressed in terms of the diagonal elements of the density matrix for the quantum mechanical canonical ensemble P (1,2,..., N) = . (2.1.6) ^ \ Here H^ is the Hamiltonian operator of the N-particle system N £ in Jd H N = ~ Z f^ + Vi^---*3) > ( 2 - 1 ^) Z is the partition function Z = -~= / ... / d 3 l d 3 2 . . . d 3 N , NX" 3 J V J (2.1.8) and the density matrix is given in terms of the energy eigenf unctions d) (1,2, . . . ,N) of the system as <1' ,2* ,...,$* \e |l,2,...,N> _Qtp = NI\ 3N Z M(V, 2 ',..., N')(j) ,(1,2,..., N)e J , (2.1.9) J J J where H„(t> . = E.d).. The pair distribution function can then be written n 2 ( i'S) = " riN~ J '*'/ < a J '2,...,H|e N |1,2,...,N> d (N-2)!^Z N J v J x d 3 ^ d 3 i+ ...d 3 N . (2.1.10) 2.2 The Effect of Spin and Statistics on the Density Matrix It should be noted that for a system of identical particles the sum in Eq. (2.1.9) extends only over properly symmetrized eigen- states. It will be useful in the following to express density matrix elements for Bose and Fermi statistics in terms of matrix elements for Boltzmann statistics and thus to deal with unrestricted sums over eigenstates. If we let the subscripts + denote Bose statistics and - denote Fermi statistics, then the symmetrized N-particle states for both cases can be written + ,2,...,N> = Z p — r |P(1,2,...,N)> , (2.2.1) - * Jul where P is one of the NJ permutations of N elements and a p is defined by + , (J-r. = +1 , if P is an even permutation, (2.2.2) P " \ -1 , if P is an odd permutation. The Bose and Fermi matrix elements can thus "be written in terms of the Boltzmann matrix elements as + = Z pI £ pI a p „a pl , Z p a p , (2.2.3) since Z ;i Z , = NiZ , where P = P"P' . When space averages are taken over the Boltzmann matrix elements on the right of Eq. (2.2.3) only 1/N! of phase space should be taken. To allow the integration to extend over all space, we define the Boltzmann matrix elements without the factor of NJ appearing in the corresponding Bose and Fermi quantities, <1',2',...,N' |e |l,2,...,N> = \^E ty*AV,2\...,W)tyAl,2,...,K)e J . (2.2.4) all j J J k Although the above expressions are correct for He , which has spin zero, modifications are needed for particles having nonzero spin. If the particles each have spin S, then an N-particle state, | l,s ,2,s , . . . ,N,s >, should specify a spin coordinate, s., as well as a space coordinate, i, for each particle. Since the Hamiltonian is spin independent, matrix elements between spin states yield Kronecker deltas, -f3H N - S(s{,s 1 )...6(s ] ^,s N ) . (2.2.5) If one then averages over spin in Eq. (2.2.3), the spin averaged diagonal density matrix element becomes 2,...,N|e S |1,2,...,H> = ^EZ ... ZE a 6(s ,Ps )...5( V Ps^) - (2S+1) 12 S N X . (2.2.6) Noting that Z 5(s.,Ps.) = 2S+1, if s. = Ps. , s. i 1 i i = 1, otherwise , one finds that if the permutation P is decomposed into cycles, the " PH N, term in contains a factor of (2S+1) for each cycle of length i. As a specific example the diagonal element of the 3-particle density matrix is + 2 <2,^,l|e J |l,2,p V2S+1) 1 " PH ^ + P <3»i>2|e 3 |1,2,3> . (2.2.7) (2S+1) 2 .3 Cluster Expansion for the Pair Distribution Function De Boer [8] has shown how the quantum mechanical cluster expansion of Kahn and Uhlenbeck [9] can he used to express the pair distribution function for a gas in terms of a power series in the 10 density. It turns out that the term in the j-th power of the density depends only on density matrix elements for systems of j+2 or fewer particles. The following brief summary of definitions and results draws from the treatment in de Boer's review article [8]. Denote the diagonal elements of the density matrix by W (1,2,...,N) = . (2.3-1) Let r stand for the set (1,2, ...,h) of coordinates of a group of h particles and define Ursell functions U« and modified Ursell functions \j\ by the relations u 1 (i) - w 1 (i) , u 2 (i,4) = w 2 (i,j) - w 1 (i)v 1 ( J l) , u 3 < i,j,k) = w 3 (i, A j,k) - w 2 (i, J i)w 1 (k) - w 2 (i,k)w 1 ( J i) - w 2 (4,^)W 1 (i) + 2W ] _ i)W 1 ( J i)W 1 (k) , etc. (2.3-2) 11 u( h) (r (h) )=W2, — ,N) = U [ h)( £ (h))W N-h ( ^l^'^N ) v (h)/ (h) . T , f u 2 fef ^ i )\_h-i^h + i--'£i-i^i + i---% ) N i=h+l + E UI. \r/ ,r.,r.)W. T ,'(r, .,..., r. . ,r. ,,..., r. ,,r. ...... r„) •j-_ h -i 3 ~i ~j N-h-2 v ~h+l' 7 ~i-l'~i+l/ '~j-1'~j+1' , ~N / + ... (2.3.4) where r. is used for the coordinate of the i-th particle. If we define ~i cluster integrals tA by b to (P oO) _ i r ru(h) (r (n) r r j D i -~ ' £IJ \r J H l ~ '~h+l , " ,, ~h+i-i ; X d 3 r, n d 3 r, ... d 3 r u . . , (2.3*5) ~h+l ~h+2 ~h+f-l * N "* " then using Eq. (2.3»*0 with h - 2 the pair distribution function can be expressed as 13 n 2 (l,2) = ^LLf • • • / W N (1,2, ' • • '^^ d ^ ' ' ' d3 S > \W Z- *-uW 2 >) r ... r H " g(3 1;;';; g) ^ ^ - A \ Z N ^ V ^ (N-2)!\ 2 r 4 2) ( £ (g) ^) , 3 , r f»» #-' s) AA a ^ ••• ." J . ^,_^ a^a-'a- * j h -J- ff n ' 2)(l(2 ' M) AA IV *>rf X Z N V x J " " ' J , , n 3(n-M d £ d S • • ■ d 2 + 4=1 * ~~ Z N Defining the activity z by dinZ M-\ 3 z) - - -^ , (2.3-7) dN ' and assuming N » £ one obtains in(\ 3 z/ = JnZ^ - inZ N , Ik or J - X-3J \L . (2 . 3 . 8) Thus, taking the limit as N -» <», the expansion of the pair distribution function in powers of the activity becomes n (1,2) = z Z < «b ( . 2) (l,2)z i . (2.3.9) Using Eq. (2.3«*0 with h=l the density, n, can be found as a power series in z also 00 n = E Sb Q 7~ , (2.3.10) £=1 l where we drop the superscript on b, (l) and omit the argument since we are dealing with a homogeneous fluid in which X>\ is independent of position,. It should be noted that homogeneity also requires that the 2 -particle cluster integrals depend only on r „ From Eqs . (2.3.9) and (2. 3 '10) a series for n (r, _) in powers of n can be constructed it— cL The first two coefficients of the series are given by 15 y (r ) = b^(r ) r 2 K 12' 1 v 12 y 73 (r 12 ) = 2i 2 \r 12 ) - ^b( 2 )(r 12 ) (2=3.12) The function obtained experimentally, usually by x ray or neutron diffraction, is the Fourier transform of the pair correlation function , , n 2 (r i2 } S(r 12 ) = g (2.3.13) n From the comments following Eqs . (2.1. 3) and (2.1.4) it is seen that the pair correlation function, g(r ), is normalized so that r li 5ooS( r i2 ) = 1 > 12 (2.3.1^) provided that the interparticle potential has a finite range. If one defines & ± ( r 12 ) > i=l,2,3,..., by :(r ) = b^(r ) A 12' 1 K 12' 1 + Z n 1 g (r ) i=l (2.3.15) it is seen from Eq. (2.3.12) that the coefficient of the term linear in the density is 5 i (r i2 } - A2), ; b l (r l2 } 2b^(r ) - kb b^hr ) 2 ■ 12 y 2 1 ^ 12' (2.3.16) 16 which can be written in terras of the diagonal density matrix elements, W, as g l (r i2 } " W (l,2)J .W 3 (l,2,2) - W 2 (1,2)W 2 (1,^) W_(l,2)W (2,^) + W_(l,2) d 3 £. (2.3.17) The linear term in the density expansion of g(r ) (or n (r )) is thus expressed in terms of the density matrix elements of systems of three or fewer particles. Similar manipulations with the cluster expansion of W yield the well known virial expansion for the pressure P (for example see Huang [10]) Z^ _ ! + B(Tl C(Tl (2 3 18) where the third virial coefficient, C(T), is given in terms of the cluster integrals by = -N?(2b_ - 4b Q 2 ) f (2.3.19) '0 V "3 and where N is Avogadro's number and V is the molar volume. Using the definition of the cluster integrals and of g (r p ), C can be written as IT N? f^ 12 )^ 12 ^\- kh 2 2 ] ' (2.3^29) Since W (r ) and b have already been evaluated for the Lennard- Jones (12-6) potential [6] it remains to calculate g (r ) in order to obtain both the pair distribution function to first order in the density and the third virial coefficient. 18 3* PATH INTEGRAL FORMULATION FOR THE DENSITY MATRIX 3 .1 Expression of Density Matrix Elements in Terms of Wiener Integrals The expression of the quantum mechanical density matrix in terms of functional integrals has been treated by several authors [11,12,13,1^,15] • ^ ne mathematically rigorous form of the functional integral in question is the Wiener integral or expectation value, E{F[x(t)]), of a functional, F[x' x t)], over Wiener measure on the set of continuous functions, x(t), over a finite interval in the independent variable To The following treatment will make use of the Wiener integral formulation [11,12,13,1^], which exhibits normalization constants more explicitly, rather than the more intuitive notation of Feynman [15]> which clarifies the transformation properties of the path integral . The quantity directly related to the density matrix is the conditional Wiener integral, which involves Wiener measure on the subset of continuous paths which have fixed values at both endpoints of the T interval. Let the 3N-dimensional vector (the number of dimensions is arbitrary, but 3N is specified for later application to the system of N particles in 3 dimensions), r(r) = (r, !,t), . . . ,r (t) ), be a continuous vector function of the variable T on the interval < t < (3. Let the interval [0,(3] be divided into k equal intervals of length At = p/k, and take the division points to be = T_ < T, < . • <• < T, . < T. = (3. 1 k-1 k Let r(T;k) be a continuous vector function of t consisting of linear segments over the intervals (t„ . ,T. ), i = l,2,.^.,k, with r. = r(x.;k), and satisfying the endpoint conditions r =0, r = R. 19 If the functional F[t(t)] is defined on all continuous functions, r(i:), satisfying the above endpoint conditions, then the conditional Wiener integral can be defined by the limiting process 00 00 E{F[ £ (T)]|r(p) = R} = lira / ... / F[r(T;k) ]d^, ■00 -00 (3.1.1) where the measure d|~u is given "by K. 3N 2 2 -3Nk du R = (2it) '" exp(— )(2itAt) exp k (r. - r. , Y 1=1 2AT n d 3N r. . (3.1.2) i=l ~ X Note that the measure du n is normalized so that k du k = 1 , (3-1.3) and thus E{l|r(p) = R} = 1 (3.1**0 The relationship between the conditional Wiener integral and the density matrix is given by a theorem due to Kac [13, lM" It states that if the scalar field U(r) has a finite number of discontinuities and is bounded below, then 20 ■BE -3N Z.Y. (R)¥.(R')e J = (2*p) exp (B'-fi)' 2P X E iexp B U(r(x) + R)dx r(B) = R'-Rf , (3.1-5) where Y.(R) are the normalized eigenfunctions of J KsiVs> + Ej-U(R) JSTjCr) = (3.1.6) If one transforms the variable, R, in order to make Eq. (3° 1.6) into the time independent Schrodinger equation for a system with the Hamiltonian ii 2 H =-E. 7T- ST. + U(R), 1 2m ~i y ~' ' (3.1.7) and makes the further transformation t -» t/|3, £( t ) "~* £(t)/nP, which normalizes the T interval to [0,1] while leaving the Wiener integral invariant, then the general density matrix element can be written = exp -^(R'-R) 2 X EX exp • e i u fe j(T)+5 r r(l) = "P (R'-5) (3.1.8) 21 This relationship holds only for the Boltzmann matrix elements since the sum in Kac ' theorem is over all eigenstates, regardless of symmetry* As stated in the references [13,1^-], Kac 1 theorem also requires that U(r) be such that the eigenvalue problem, Eq. (3-1.6), yields a discrete spectrum. If the system is enclosed on a box of volume V, then U(r) should have a component due to the walls of the box. This component certainly causes the spectrum to be discrete and gives a natural normalization for the eigenf unctions . However, if one restricts the temperature to values larger than those at which long range order phenomena, such as Bose condensation, become important, then the potential of the box contributes nothing to the problem in the limit V -* oo. One may neglect this component of U(r) provided one considers the wave functions of the continuous spectrum to have delta function normalization and treats the sum over states Z. as representing a sum over the discrete states plus an integral over the continuous spectrum, A slight reformulation of Eq. (3. 1-8) will be useful in the following work. Wiener measure on r(f) gives rise to the interpolation formula [l6] r(T')(T"-T) + r(T")(T-T') r , , w ,, v-|l ,( T) = -,_, + | L (T ( ^»T) T > j 2 , (3.1.9) where t' < t < t" and where £ is a 3N-dimensional random variable whose coordinate random variables are independent and normally distributed with mean and variance 1= From this formula it can be seen that if r'(T) obeys the endpoint conditions, r'(0) = r'(l) = then both r(x) and r'd;) + : — (R'-R) are distributed as the random variable ~ A ~ ~ 22 ^-f 1 (5'-£) +1[t(i-t)] 2 Equation (3. 1.8) can then be written , " PH Ni = exp -4(5' -g) 2 X E-^exp -PJU Wiit r(x) + R + t(R'-R) dT r(l)=Oh , (3.1.10) where the prime on r'(x) has been dropped. A useful intuitive picture of the Wiener integral is obtained by thinking of an imaginary motion of a system of N particles over a unit "time' 1 interval, < t < 1. At time T = the coordinates of the particles are given by R and at 1 = 1 by R' . The quantity (R'-R) is thus the sum of the squares of the distances between the initial and final positions of each particle, and the Wiener integral in Eq. (3.1.10) is the average, over all possible paths of the motion, of the Boltzmann _Q [T factor, e , with the potential energy U replaced by the average of the potential over the path. The paths are weighted by Wiener measure which is related to the average over the path of the total kinetic energy of the system. It should be noted that nS :(t) 2* 23 is a measure of the deviation of the motion from the shortest straight line path from R to R 1 , and as the classical limit is obtained by letting \ -> 0, the deviation vanishes in this limit. The off-diagonal matrix elements vanish in the classical limit since lim exp \-*0 -4(5' -g) 2 ee 6(R'-R) , (3-l.ll) where 5 is the Dirac delta function, and the diagonal matrix element approaches the Boltzmann factor -PH^ -PU(R) lim = e ~ . (3-1.12) \-0 3-2 Upper Bounds on Off -Diagonal Matrix Elements Both diagonal and off- diagonal matrix elements for Boltzmann statistics are involved in the calculation of g (r ) for a Bose or Fermi system. For example, in the evaluation of the 3 _ particle term, -f3H 2,2|e 3 |1,2,3> + , two types of off -diagonal Boltzmann matrix elements arise. One involves the exchange of coordinates of two of the three particles and is of the form 2k ■PH <^2,^h 3 |g,l^> = exp -4 X (2"1) 2 + (1"2) 2 x E^exp -0 / U, -^r r(x) + (2,1,5) + T[(i*2,^)-(2,1,2)] dT 3 US r(l)=0 (3-2.1) The other type involves a cyclic exchange of the coordinates of all three particles and is typified by the matrix element -PH 2,3|e 3 |2,2,1> = exp -±A (2-l) 2 + (^-g) 2 + (1-3) X EJexp •pJU 3 f-^r(T) + (2,3,1) + T[ (1,2, ^-(2,3,1 )] jdT r(l)=o| . (3-2.2) The form of the factors multiplying the Wiener integrals in these expressions yields a significant estimate of the size and temperature dependence of these off-diagonal elements for intermolecular potentials for which the particles are characterized by strong repulsive cores. Suppose the 3 _ body potential satisfies U 3^~l'~2'~3^ = °°' if r ij < ° fQr any pair i ^ ssl » 2 >^' (3.2.3) The Wiener integrals in Eqs . (3-2.1) and ( 3-2.2) must vanish if the initial distance between any pair of particles is less than a. Hence 25 wherever the matrix element is nonzero the factor multiplying the Wiener P P integral is smaller than exp(-2ira /\ ) for the 2-particle exchange term p P and smaller than exp(-3 :rt cj /\ ) for the 3 - particle exchange term. k For He , approximating the repulsive core diameter by the -ft parameter a = 2.56 X 10 cm of the Lennard-Jones (12-6) potential, one finds that the larger of the two multiplying factors, that for the -•6T o 2-particle exchange, is approximately e » Thus at 10 K one expects the off -diagonal elements to be smaller by a factor of at least .0025 than the diagonal element, for which the multiplying factor is unity. -f3H 2 The argument for the 2-particle exchange term in 2> is the same as the argument for the similar 3 - particle term in Eq. (3»2.l) and yields the same result. Due to restrictions of computing time the numerical scheme to be developed below for the calculation of the diagonal matrix element is limited in accuracy at low temperatures . In all cases, the above bound on the exchange terms shows them to be negligible compared to the numerical uncertainty of the diagonal term. Therefore, in subsequent sections the exchange terms will be discarded and the diagonal Boltzmann matrix element will be used as if it were the same as the diagonal element for a Bose or Fermi system,. The above argument was first used in connection with the 2-particle exchange term in the second virial coefficient [17]. The suppression of this 2-particle exchange term with increasing temperature is actually somewhat stronger than is indicated above. This is verified 26 both by numerical calculations [5,6] and by upper and lower bounds obtained analytically by Lieb [18], who shows that at high temperatures the term is asymptotically proportional to -2 Qualitatively the increase in the constant multiplying -X, in the exponent comes about because the Wiener integral in Eq. (3-2.1), as well as the multiplying factor, decreases exponentially with increasing temperature. This is due to the fact that the set of paths for which the "motion'' of the system from time t = to time T = 1 brings the two exchanged particles no closer than a from each other for any T is assigned exponentially smaller Wiener measure as T -» <». 3-3 Approximation of the Diagonal Matrix Element by a Multiple Piemann Integral The calculation of the diagonal matrix elements for Boltzmann statistics is now reduced to the evaluation of a Wiener integral. For the diagonal element of the 3 - particle density matrix Eq. (3*1-8) becomes 27 1 = E -s exp -Pj U 3 h=r(T) + (l,2,3)jdi r(l)=0 f . (3-3.1) Analytical evaluation of the above Wiener integral is possible, however, only in the cases U 3 (r x ,r 2 ,r ) = , and 3 U 3 (r 1 ,r 2 ,r 3 ) = .^^r.) , where f.(r.) is quadratic in each coordinate of r. . Thus for a realistic 1 ~i' ~i interparticle potential it is necessary to rely on a numerical evaluation of the Wiener integral. Numerical schemes for evaluating Wiener integrals have been proposed and applied by several authors [6,19,20,21,22,23]* The simplest of these methods involves taking a finite but large value of k in the limit definition for the conditional Wiener integral, Eqs . (3.I.I) and (3-1.2), and using the result of a numerical evaluation of the k-fold multiple integral as an approximation to the Wiener integral. Since k is large and the weight function du can be conveniently K. expressed in terms of Gaussian density functions the evaluation of the k-fold integral is usually done by the Monte Carlo method. 28 This simple scheme is adequate for one or two particles under the influence of a computationally simple potential and has been applied successfully to such problems [6,20,22], As the evaluation of the integrand of the multiple integral becomes more complex, however, the computing time requirements become unreasonably large, and for this reason methods have been proposed for reducing the multiplicity k of the integral to be performed [19,21,23] . The present calculation makes use of a slight modification of the procedure developed by Fosdick [21]. The Wiener integral in Eq. (3«3«l) can be written as a Riemann integral of a product of Wiener integrals involving a smaller value of 3 (higher temperature)- If the interval < t < 1 is divided into k equal subintervals, (t. -,,t.), i=l,2,...,k, then the integrand in Eq, (3.3=1) can be written as a product of functionals depending on r(x) over only one subinterval F[r(i)] = exp f u 3 (^ £ (t) ♦ (i,s, 3 )j ** ] n fJrfT)] , i=l (3-3.2) where fiteC-O] = exp T. 1 -PJ UJ_\_r(T) + (1,2, 3) dT T. , V2Jt l-l (3.3.3) 29 If each of the k subdivisions is now further subdivided into n equal parts, Eq. (3-1-1) may be written in the form 00 00 s {F[r(T)]|r(l) = 0} = lim / ... / du nk F[r(T;nk) ] , n->oo ^ ^ •00 -00 00 00 00 00 du IT lim i=l n-»oo' /•••J\,n f i [ ^ nk > ] ' 00 - CO -co -00 (3.3.10 where d(i is defined as in Eq. (3-1.2) in terms of the coarse K. subdivision only and du, „ = (24) 2 i,n exp (r. -r. _)' 2(p/nk) 3Nn X exp ■ n (r. .-r. . , Y -y '~i;j ~i;.]-i y . j=1 2(p/nk) (3.3-5) where r. . = r(f . . :nk). The limit on the right hand side of Eq. (3-3-M is just the definition of a conditional Wiener integral 00 00 lim / . . n-» 00^00 du^f^r^nk)] = E^exp -p U, T i-1 \ r(T) + (1,2,3) Ut r(T. .) = r. ,,r(T.) = r.X . (3-3-6) ~ 1-1 ~ X -1'~N X ~1J 30 The time interval in Eq< (3«3-6) m ay be renormalized and the endpoint conditions altered by the transformations t 1 = k (x - x.^) , ' (t) =vk[r(i) -r. , - i(r. -r. n )l ~ -l-l x ~i ~i-l /J Equation (3«3-^) then becomes 00 00 E(F[r(T)]|r(l) = 0} = dn, (3- 3- 7) -00 - X) X H Ei exp i=l .& r(x) U-, -=r ~T=~ + r - 1 +T(r.-r. J + (i*2,3)W T r(l)=0^ , (3.3-8) where the primes on r(x) and t. have been dropped. Since \//k =- 2rth 2 (p/k)/m 2 , the Wiener integrals on the right of Eq. ( .3 • 3 • 8 ) involve a temperature k times that in the Wiener integral on the left. We proceed, then, by finding a high temperature approximation for the Wiener integral and using Eq. (3«3'8) "to extend its range of validity to lower temperatures. 31 The high temperature approximation to be used is the well known Wigner-Kirkwood expansion in powers of |3 and \ [1,2^]. The most useful form of this expansion for the following work will be the expansion for the quantity G = inEjexp -p / U 3 I — £(t) + r + T(s-r)jdi ■(l)=0 (3.3-9) Introduce a parameter t] multiplying the path, r(l) , in the Wiener integral, and define the following functions of r\ f(ti) = exp P / U n M=L r(T) + r + T(s-r)]dT J 3 V^ ~ ~ ~ ~ J ■ F(tj) = E(f(ii)|r(l) = 0} , and G(t]) - ZhF(r\) > (3.3.10) A formal series expansion for G may then be found by expanding G(t|) as a power series in tj about t) = and taking the limit r) -» 1 to obtain G = lim G(t)) . 1^1 (3-3.11) 32 Up to order t) one finds G(n) =inF(0) + ^T^ F(0) drj q=0 2 1 dF(n) _F(0) dr] " 2 1 d 2 F(n) n J F 0) , 2 T]=0 - 1 ' dr) T) = (3-3.12) The integrand of the Wiener integral F(0) is independent of the path so it is seen from the normalization condition, Eq. (3.1.U), that F(0) - f(0) = exp ■P / U (r + T(s-r))di (3.3.13) The first term in the expansion for G is then G - -0 / U (r+T(s-r))dT , (3.3.1M which is just the exponent of the classical Boltzmann factor with the potential replaced by its average over the straight line path from r to So If this approximation is used on the right hand side of Eq. (3- 3»8), the result is precisely that of using Eq. (3.1.1) with a finite value of k as mentioned above. 33 The coefficient of r\ in the expansion for G is 1 dF(n) F(0) d-q n =o F(0) E^-pf(O) -y — r(i ) • W (r+T(s-r))dT >/2rt ~ ~ 3 ~ r(l)=o| , (3-3.15) where V is the 9-dimensional gradient operator. If one reverses the order of the Wiener integration and the integration over 1, the Wiener integral depends on r(x) only for a fixed 1. This Wiener integral may be evaluated by using the interpolation formula, Eq. (3 .1.9)* f° r r('T') with t' = 0, t" = 1 and taking the expectation value over the normal distribution of the random variable £ E{r(T)|r(l) = 0) = <|[t(1-t)] > = , ~ d V (3.3.16) therefore dF(r , ) 4l = TJ=0 (3.3.17) The coefficient of r\ thus vanishes and the coefficient of t) 1 2. becomes 3^ 1 _2 dn n =o r(l)=0 2 9 9 o- + T P J ST * - Z , r i (T)r J (T) 5I5J U 3 ( ^ T( ^ r)) dT i=lj=l r(l)=Oh , (3-3.18) au. where r.(i) is the i-th coordinate function of r(x) and vr^ denotes partial differentiation of U with respect to its i-th argument. The first term on the right of Eq. (3.3.18) is of order \ «T while the second term is of order (3\ «T . Thus to lowest order in l/T 1 d^F(n ) T]=0 9 9 .2 ~ " / ^ ^ ^|U U (r+T(s-r))E{r.(T)r.(T)|r(l)=0) 1=1 .1=1 J 3 J 2* dT . (3.3.19) Using the interpolation formula again one finds Efr, (T)r (f)|r(l) = 0} = , (ifcj) t(1-t) , (i=j) (3.3.20) Therefore to order T 35 G - r 2 r -P / U (r+T(r-s))dx-^- / V U (m(s-r)h(l-T)di (3-3.21) 3 If Eq. (3.3-21) is used for the Wiener integrals on the right of Eq. (3-3-8) one obtains the approximation 00 00 K{F[r(T)]|r(l)=Oj = J . . . J dn k F[r(-r;k) ]C k , (3-3-22) -00 -00 where C = n exp k 1=1 ^i-l^^i^i-l) + (1,2,3)1 T(l-T)dT (3.3.23) If, for each i=l,2,...,k, the exponential in Eq. (3-3-23) is replaced by the first two terms in its power series expansion, it is seen that -2 to order k C k = n k 1=1 1 B\~ kitti F IL \ Hn& r. , +f(r . -r . . ) ~i-l v ~i ~i-l' + (1,2,3) T(l-r)dT (3-3- 2k) Equation (3.3.22) with Eq, (3-3.24) for C is precisely the formula obtained by Fosdick [21] who shows that if U is sufficiently -2 regular, the formula is correct to order k . In the present 3 application, where U (1,2,3) = £ V(r..) with V(r) given by Eq, (2.1.2) the regularity conditions are not satisfied due to the singularity of 36 V(r) at r = 0, It has been found from calculations of the 2-particle density matrix elements that the use of Eq. (3-3.23) -^ or C, S^ ves better K. results for small k than are obtained using Eq. (3- 3*24) due to the fact that at the singularities of U , Eq. (3»3«23) is bounded in absolute value while Eq* (3« 3 -2k) is not. For this reason the exponential form for C, will be used in the following calculations, k The Wigner-Kirkwood expansion can be carried out in the above manner by repeated use of the interpolation formula to obtain approxi- mations which are correct to order k for any m under the appropriate regularity conditions on U. However, the increasing complexity of the C, obtained in this way causes an increase in computing time which must be balanced against the decrease obtained by using a smaller value of k in Eq. (3.3*22). The special case k = 1, which is equivalent to expanding the logarithm of the diagonal matrix element directly, is simplified by the fact that the 1 integrals may be done immediately. The result of the expansion to order T in this case is in Eiexp -^U 3 fj^r[r(T)+(l,2,3)AdlT r(l)=0J 2 2. 2 -pu 3 (i,2,a) -§^^3(1,2,3)+ |£-[vu 3 (1,2,2)]' -^vV (1,2,^) ' (3.3.25) 96O7T ^ 37 k. EVALUATION OF g^r ) l+.l The Multiple Integral Approximation for g (r ) The complexity of the 3 _ P a rticle functions entering into the calculation of g-,(r 1p ) can be reduced somewhat by consideration of the symmetries of the problem. Since there is no fixed external field, the density matrix must depend only on the relative positions of the particles, and hence the first simplification should be to separate out the center of mass motion. Define the usual center of mass and relative coordinates for a system of three particles, Figure 1, by £ CM = fvi-2-1) > E^-g-A . S = 3-|{g+l) . (k.l.l) The Jacobian of this transformation is unity, and furthermore £j^-£(!*.*«£*iii <*•"> so that the coordinates r_„ r,_, and Q behave as independent particles ~CM ~12 ~ of masses 3^, m/2, and 2m/ 3 respectively. Taking U to be the sum of pair interactions as in Eq. (2,1.1) 1 1 and noting that r, _ = 3-1 = Q + o r io an( 3 r_„ = 3~2 = Q-— r,_ it is seen ~13 ^ ~ ~ d -L^- ~2J ~ ~ ~ 2 ~12 38 FIGURE 1. CENTER-OF-MASS AND RELATIVE COORDINATES FOR 3 PARTICLES. 39 that W (1,2,3) can ^ e separated into the product of a center of mass factor which is independent of U and a factor containing only the relative coordinates, r, and Q, W_(l,2,£) = , " PK PM. , " PH rpl, = < W e l£cM >< £l2^h |r 12 ,§>, (4.1.3) where K - *% CM 2 (3m) is the kinetic energy of the center or* mass motion and H rel " -2^27 " iWT) + U 3 ( ^12'S) &**) is the Hamiltonian operator for the relative coordinates. It is to be understood that the normalization factor of X for each particle coordinate appearing in the definition, Eq_. (2.2.4), of the Boltzmann matrix elements is to be replaced by the thermal wave length for a particle of mass corresponding to that associated with the appropriate transformed coordinate (4.1.5) ko The normalization of the density matrix elements yields = 1 ~CM (4.1.6) so that, using Eq. (3-1=8), the diagonal density matrix element can he expressed in terms of a Wiener integral over paths in the 6-dimensional relative coordinate space of r, _ and g: W 3 (i'2,2) EJexp -Pj U 3 (r(T),q(T);r 12 ,3)dTj r(l)=q(l)=o| (^.1-7) where U (r(T), a (T);r ,g) = V — r( T ) + r 3 ^T -12, 12 + V -^- (nT? ^r)+r(r))+^+ Z 2- ) + V -~ (^3 ^T)-r(x))+^ -12 2^7 2 anc *- r po between the particles it can be seen that if r is taken to lie along the axis of polar coordinates, (p n , 6 n , cp ), for Q, then G(r ,Q) is independent of the direction of the polar axis and of cp . Thus one can write Sl (r 12 ) = SnjTsin e^f^ C(r 12 ,p Q ,e Q ) . (4.1. 14) The integral over the position of the third particle which appears in g (r ) can thus be reduced to a twofold integral. Equations (4.1.8) and (4.1.11) show that F can be written as a product F 123 (r^ = W^TTU 4F 123 (r(T)^(T))-F 12 (r(T))F 13 p q(r) + \ r( T ) " F 12 (r(T))F 23 ^ £ (T)-|r(T)lF 12 (r(T))jr(l)= 3 (l)=o} (4.1.18) Due to an oversight the quantity g (r ) rather than g(r ) was originally calculated. The difference between these two quantities arises because, although the distributions of the arguments of F, and F Q _ in Eq. (4.1.18) are the same as the distributions of the arguments of these functions in Eq. (4.1.12), the use of a linear combination of jj(t) and r(x) as the argument of F in Eq. (4.1.18) introduces a spurious correlation between F 13 fl(') - | E(t) and the functional F (r(r)) multiplying it and similarly for the product F, F o The two factors in these products should be 45 uncorrelated since they arise from the product of two independent Wiener integrals. The function g (r ) can be written as the sum of g (r ) and a term which corrects for the unwanted correlation: where and gl (r 12 ) = g c (r 12 ) + g e (r 12 ) , *e (r 12> =/ d3 S G e tr 12 ,§) , (4.1.19) G e (r i 2 ^) = w^y E { F i 2 ^))k 3 fFa(^ +|£(t) - f 13 ( s (t)) + F 2 3PT^) - |r(T)J - F 23 (q(T)) r(l)= 3 (l)=o| (4.1.20) The numerical calculation will thus involve separate evaluations of the terms g Q ( r 12 ) and § e ( r i 2 ^ Using Eq. (4.1.15) the quantity G (r ,Q) can be written in c v 12' the form ° c < r i2'S> = iOT^) E { F i2 (£(l)) [ r i3lT s(^) + 1 iM - 1 f 23 Pf^ t) -|^ t )I- 1 r(l)=q(l)=0 (4.1.21) The results of Section 3-3 nia y now be used in order to express G (r,p,Q) and G (r ,Q) in terms of multiple Riemann integrals over a 6k-dimensional space. Let r(T;k) be the broken straight line path defined in Section 3-1 with breakpoints at r. /i=0,l, . . ..,k, and similarly for q(x;k). Define the functions k6 D 12 (r(T;k)) = F 12 (r(T;k)) ) - Z -&£ I V"[^;[r, ,«(*,-&■•., )] + r 1Q |T(l-T)dT L 1=1 2nk V?" 1 - ■i ~i-l /J ~12 D 13 ( a (T;k)) = F 13 (q(T;k)) 1 2 r, 2*k £12 exP^ &L.J V" ^[ ai . 1+ T( a .- ai . 1 )]^ + ^ T(l-T)dT D 23 (£(T;k)) = F 23 (q(x;k)) X exp k M 2 L^/ v W a - + ^- 3 - )]+s - -12 T(l-T)dT J where V"(r) = V 2 V(r) (V.1.22) (4.1.23) The approximations for G and G are then c e 00 00 G c (r 12 ,g) Vi'S) D 12 (r(T;k)) D 13 \~2 ^ T;k) + \ ~ (T;k) - 1 -00 -00 X D 23 S a ( T;k ) * ( T ;k))-1 rk qk (1+.1.24) and 47 00 00 G e (r 12'S) ~ OT^y /•••/D 12 (r(x ; k)) -oo -oo d 13 I ^ q(x;k) + - r( T |k) D 13 (q(T;k)) + D 23 l — q(T^k) - ^ r(xjk) I D 23 (q(T;k))| d ^ rk d ^ k > (4.1.25) '3 where 3 _3k dn gk = (2«) 2 (2*/k) ' exp f k (s ,-s J 2 "]^! (4.1.26) with s = r or q, From the discussion in Section 3-3 of the order of accuracy of Eq. (3»3»22) as an approximation for the Wiener integral it is seen that the above approximations for G and G cannot be expected to be -2 correct to more than order k . The order of accuracy may be less due to the singularity of V(r). In view of this fact the integrals over T in Eq. (4.1.22) may be replaced by their Simpson's rule approximations without decreasing the order of accuracy of these formulae. If we let r and q be the 3k-dimensional vectors (r, , r ,,°>,r, ) and (qn ) < i^)-'') ( li ) respectively and recall that r„ = q_ = r, = q, - 0, ~1 7 22. ~k ~0 ~0 ~k ~k then k8 00 00 G c (r 12'S) w , ,s/T (k) 1 . (k)' D 13 2 S + 2* -00 -00 '] »J#»""-h""U d ^rk%k ' (U. 1.2. T ) and 00 00 ».('«.«) ■ obr/ • • • / ^ (k) >W? * (k) + I ; (k) - D i 3 )-%3 > du . d(j. . , rk qk (i^.1.28) where D 12 (£ (k) ) = < { * iBfe ** £ 12 I 1^ v 2 ; -12 1 A, I \ _ _ ? t / A, /~i ~i-l\ 23J Hil. 21 \ -12' l^^i.i + 3tTl^ v l^( X 2i + Si-l N „ £12 )+Q+- i=l u W« *2 X 2 Jx ,V*I-1 + VI — Q.+Q+ -~| + £— V' ! — (- wit 0+Q+ -12 2- 2 (^.1.29) where the + sign in the last equation corresponds to D and the - sign to D . h 9 4.2 Monte Carlo Evaluation of the Multiple Integrals It has "been indicated above that the 6(k-l)-fold integrals appearing in Eqs . (4,1 ,27) and (4 „ 1,28) can be evaluated by a Monte Carlo procedure o The quantities d(i and du have the properties of probability densities on 3 (k-l) -dimensional spaces since the integral of either of these quantities over any 3 (k-l) -dimensional Borel set is non-negative and the ' integral over the whole space "is unity (see Eq<- (3«1»2)). The simplest Monte Carlo scheme for eval- uating the multiple integrals is to choose a sequence of independent (k) (k) 6k-dimensional vectors (r , q ), r = q = 0, i = 1, 2, . . . , M ~i ~l ~k,i ~k,i with probability density du du , evaluate the integrands of G and G for each vector, and average over the set of M vectors- If one lets I (r ■ ', q ) and I (r , q ) represent the integrands of G and G respectively, then by the law of large numbers U» i Z I(r (k) q (k) } -/•■■■/ W k) » 3 (k) ) (4,2,1) dn . djj. = G } rk qk c -00 -00 with probability one and similarly for G , 50 (k) Vectors r with distribution du . can "be generated from ~ rk 3(k-l) independent, normally distributed random variables by means of the interpolation formula, Eq. (3«1.9) • The endpoint conditions require that r = r. = 0, and the other 3(k-l) coordinates can be ^ ~o ~k v ' generated by the equation :i-i (1 - l> + k - (i - -) k v k' 1 i-1 1 -I 2 (4.2.2) for 1 a 1, 2, . .., k - 1, where i. is a 3 -dimensional random variable, the coordinate variables of which are normally distributed with mean (k) and variance 1. The vectors q are generated by the same procedure. ~J This straightforward Monte Carlo scheme has been used [6] to compute the 2-particle density matrix element, Wp(l,2). The 3-particle functions G and G are well adapted, to a . n more sophisticated Monte Carlo procedure which should have a smaller variance than the straightforward scheme, thereby reducing the number of samples necessary for a given accuracy. This procedure is the importance sampling technique first introduced by Metropolis and collaborators [25] for evaluating averages over the Boltzmann distribution. The following discussion is, to some extent, an P-L extension to a continuous sample space of the description of this technique given by Hammersly and Handscomb [26]. Let G and G both be typified by 00 00 where ^'^'-v^r'^'s^ ^- h) and let P(r< k ») = V^ (^-5) W,(l,2) This factor, vhich appears in the integrands of both G and G , is the basis for the importance sampling procedure. It is seen from the definition, Eqs ~ (4,1,29), of D,„(r ') that this quantity is non- negative* Furthermore, to the order of the approximation used in Eqs, (4<,1,24) and (4,1,25), one can write 00 00 w 2 (i,2) =/..-/ d^Ce 00 ) *n rk > e^.2.6) 52 and hence 00 00 f,.,fl>(r M ) dn rk - 1 . (4.2.7) -00 -00 (k) Thus P(r ) du has the properties of a probability density on the ~ rK 3 (k-l) -dimensional space of r . If vectors r and q^ are generated with probability density P(r ) dp: du and L(r , q ) ^ ric q^ '**' ^ is averaged over the vectors so generated, the average will approach G with probability one as the number of samples becomes infinite. Qualitatively, the advantage of this scheme over straightforward sampling is that some of the variation in the integrand l(£ / q ) has been absorbed into the probability density. The variance reduction (k) is indicated intuitively by the fact that vectors r for which (k) (k) (k) D p (r ) and hence l(r , q ) are small will be generated with low probability by the importance sampling procedure. (k) It remains to show how vectors r can be generated with (k) probability P(r ) du . This is accomplished by producing a Markov ~ rK (k) chain, {r , t = 1, 2, 3> • - - } ; with a stationary transition •obability function p(r^ , |r| ), satisfying 53 J ../p(r'W)d„ 00 00 •f-If- JW (k) ii- p( ^ (k>) d ^ k d3(k_1) £' - ^- 2 - 8 > -00 -00 where d 3(k_l) r' = d 3 r'd 3 r' ... d 3 r" , (^.2-9) ~1 ~2 ~k-l and where J is any set in the 3(k-l) -dimensional space of r 1 which is measurable with respect to P(r* ) du , . The law of large numbers ~ r k (k) for Markov chains [27] states that if the Markov chain, {r \ )> is such that Eq. (U.2.8) is satisfied and if there is only one ergodic set then M CO 00 £a£^s (k) >-/---JV k) . M~ M ", -«'"' ~^" % ' f-f-**' S^'' P ^ (k) ) d|J rk' (U ' 2 - 10) for almost all realizations of the Markov chain. Thus if L(r^ , q ) (k) (k) is averaged over r produced by the Markov chain and q ' generated by the straightforward Monte Carlo scheme, the result will be an approximation for G . We must now demonstrate a procedure which will generate a Markov chain which satisfies the above conditions. It will be convenient for the rest of this section to deal with only the 54 x -components, x., of each 3 "dimensional vector, r.. All quantities (k) defined above in terms of the 3k-dimensional vector, r v = (r , r , . . . , r ), r = r =0, are assumed to have analogous ~1 ~2 ~k ~0 ~k definitions in terms of the k-dimensional vector x = (x, , x„, ... , x ), x„ = x = Oo The generalization to three space components will be obvious, and the notation will be much simplified if we assume all vectors to be k-dimensional and drop the use of k as a superscript or subscript to indicate dimensionality for the rest of this section. Let w(x) be defined by ,k-l ,k-l du = w(x)d x, d x = dx, dx_ . .. dx. , , x ~ ~ ~ 12 k-1 (4,2.11) so that w(x) = (2ir) k |2 >2it exp S Z (x.-x. .)' 2 1 1-1 1=1 (4.2.12) A transition probability function, p(x' |x), is then needed to satisfy Eq. (4.2.8), which now becomes /;•■/ P(x') w(x')d k " 1 x' oo oo .f ... ff* . . . f j>(x' \x) P(x)w(x) d k_1 x ■_ d k_1 x' . (4.2.13) -00 -00 55 First define the function p*(x'|x) =l±\ k-1 2 exp -k k-1 / x! Z [x! - -1 1=1 V 1 -1 + x i + l (4.2.11+) This function satisfies the criteria for a transition probability- function since p*(x' |x) > (4.2.15) and 00 00 f'"f P*(2'l;X)a k ~V = !' -00 -00 (4.2.16) The function p(x' |x)' i s now defined by P(x') p*(x«|x) p^y- , if P(x') < P(x) P(x'|x) - P(x) , 1 - P(x") pTxT ,k-l „ d x , (4.2.17) •here Gx"|P(x") < P(x)} is the set of all x" such that P(x") < P(x) and 56 5 is the Dirac delta function. The function p(x' |x) is a valid transition function since it is non-negative -and satisfies 00 00 I ... r P (x , |x)d k " i x' -00 -00 {x 1 |P(x') > P(x)) p*(x' |x) + B(x',x) [x M |P(x") < P( J^(-^)* k -v]^ r p(x,) ki {x'|P(x') < P(x)) =/ . • . /W te>* k ~V (x'|P(x') > P(x)} f ... jr p *(x"ix)d k " i x*' , (x"|P(x") < P(x)) = 1. (4.2.18) 57 In order to show that p(x' |x) satisfies Eq. (4.2.13) it will T he useful to define the "transpose", x , of a vector, x = (x , x , ...,x ), ~ ~ • l d k x Q = x k = 0, by x ± = x k _ i , i = 0, L, . .., k. From Eqs . (4.1.29), (4.2.5), and (4.2.12) it can be seen that w (x) - w (* ) , (4.2.19) and P(x) = P(x X ) . (4.2.20) k-1 T T k-1 T Since P(x)w(x)d x = P(x )w(x )d x , and since the integration with respect to x in Eq. (4.2.13) is carried out over all space, it is clear that Eq. (4.2.13) is equivalent to " J ' P(x , )w(x')d k " 1 x' ■/•;•/ P(x')w(x*) + P(x' T )w(x ,T ) , . ~ ~ ~ ~ K.--L , 5 d x . > oo oo ■/■;•//-/ -00 -00 p(x' |x) + p(x' |x ) P(x)w(x) d k_1 x d k_1 x« . (4.2.21) 58 As a lemma for the proof of Eq. (4.2.21) it is shown in Appendix A that p*(x' |x)w(x) = P*(x T |x ,T )w(x ,T ) . (4.2.22) To prove that p(x'|x) satisfies Eq. (4.2.21) it is sufficient to show that 00 00 v(x') = -00 -00 p(x' |x) + p(x' |x ) P(x)w(x)d k_1 x , = 2P(x')v(x') , (4.2.23) where v(x') is defined by the above equation. Using the definition, Eq. (4.2.17), of p(x' |x) and keeping in mind that P(x) and w(x) are invariant under transposition of the coordinates of x it follows that v(x') =f . . . J U|P(x) > P(x')} p*(x'|x) + p*(x' T |x T ) P(x')w(x)d k ' 1 x (x|P(x) < P(x')) p*(x' |x) + P*(x" T |x T ) + 6(x',x) {x"|P(x") < P(x)} p*(x"|x) + p*(x" T |x T ) 1 - pjy I d k ""V \ P(x)w(x)d k_1 x , (4.2.24) 59 where the integral multiplying the delta function can "be extended over the set {x"|P(x") ^_P(x)} since the integrand vanishes for P(x") = P(x) Performing the integral over the delta function v(x') becomes v(x') - P(x') {x|P(x) > P(x')} p*(x'|x) + p*(x ,T |x T ) w(xjd x {x|P(x) < P(x') p*(x*|x) + p*(x ,T |x T ) P(x)w(x)d k " ;L x / {x"|P(x M ) < P(x')} L p*(,x x ) + p*(x X ) ■/ ••• J {x"|P(x") < P(x')} P*(x"|x') + P*(x" T |x ,T ) P(x')v(x , )d k " 1 x' P(x")w(x')d k ' 1 x n . (4.2.25) Applying Eq. (4.2.22) to the first and second terms in the last equation one finds that the second and fourth terms cancel and there remains v(x') = P(x')w(x') (x|P(x) > P(x' )) p*(x T |x ,T ) + p*(xjx') d x + P(x*)w(x') J ... J {x"|P(x") < P(x'j} p*(x"|x') + P*(x'* T |x iT ) ,k-l „ d x (4.2.26) 6o Since the sets over which the two integrations are carried out are complements, the integrals may be combined into an integral over the whole space. Bearing in mind that the variable of integration, x, may T just as well be replaced by x one obtains the result v(x') = P(x')w(x') 00 00 00 00 j ... yV(x ix' )d k_i x + j . . . y x p^(x T ix ,T )d k " i x T = 2P(x')w(x') , (4.2.27) where the final step makes use of the normalization condition on P*(x|x'), Eq. (4.2.16), It can be seen from the definitions of p(x' |x) and p*(x' |x) that there can be only one ergodic class since the probability of the one step transition from any vector, x, to a vector in the k-1 volume d x' about x' is strictly greater than zero for all vectors, x J , satisfying P(x') > 0. Equation (4.2.10) is thus satisfied for the p(x* [x) defined above, and a procedure which generates a Markov chain with this transition function is a valid importance sampling scheme for the present problem. The transition from a vector, x, of the chain to the next vector, x', can be made using k-1 normally distributed random variables and one uniformly distributed random variable. Let 6l it ii x o = x k o , : i = V^ Sj X . . , + X. , 1-1 1+1 i = l,2,...,k-l , (4.2.28) where £., i=l,2, . . .,k-l, are normally distributed random variables with mean and variance 1. If P(x" ) < P(x), generate a random variable, £, which is uniformly distributed on the interval (0,1) and set P(x" ) (4.2.29) If P(x") > P(x) set ~ ' " * ^ p(x) ' x , otherwise. (4.2.30) It can be seen that the probability of x" given x is P*(x" |x) while the probability of x 1 given x is p(x ! |x). 62 4.3 Integration over the Position of the Third Particle Besides the 6(k-l)-fold integration arising from the approximation to the Wiener integral, it is also necessary to perform the twofold integration over the position of the third particle. This is the integration over (Pqj^q) i n Eq. (4.1.14), Since a large number of Monte Carlo samples are necessary for reasonable accuracy in the Wiener integral approximation, it is impractical to compute G (r p ,Q) and G (r, p ,Q) at each point of a 2 -dimensional grid in Q and apply a standard numerical integration technique. If we perform this last twofold integration by a Monte Carlo technique, however, the sampling may be combined with that for the 6(k-l)-fold integrals. The procedure then is to generate M samples {(r; , qj , £.), t = 1,2,...,M] and to form M z t=l ( ri2 ;M) = | Z L(r[ k) , q[ k) , Q t ) r" 1 ^), (l*. 3 .l) (k) where L is defined by Eq. (4.2.4), r^ is generated from the Markov (k) chain, q;. is produced by the straightforward sampling scheme for Wiener paths, and £ = (p_,0 A ,O) is picked according to a 2-dimensional distribution r(g)d 2 S - r p (p Q )r (e Q )dp Q de Q , (M#2) 63 or r(S)a 2 £ = r x (Q x )r y (Q y )d ~£( T ) ^ n § > then F,p(r(T)) is unchanged due to the spherical symmetry of V(r) while F 23 ~2 S (t) " Z ~2~y Fl 3 ( ~2 S (t) + Z ~2~ and F 23 (q(T))-F i3 (c/T)) 69 Thus ^(V = g' 2 ^) (It. 3.12) and g (r ip ) may be written as (r 10 ) - ■ 12' W 2 (l,2) d 3 Q X E^F 12 (r( T )) (t) F l3l~2^ T ) +: - F 13 (q(T)) r(l) = q(l) = O; (^•3.1 In this last form a more natural origin for the integration *12 over Q is'Q=- — ==-, so let the variable of integration he changed to Q'" = Q + -12 (4.3.1*0 Then g e ( r !2) = J d V G e (r 12'S ! ) * (4.3.15) where TO G e (r 12 ,£') - 2E 1i F 22^ r ^ W 2 (l,2) F 13 [ Q', -|c p(t) Fi 3 (S'^(t)) and where r(l) = q(l) . (A , (^•3.16) F 13 (s;a(^)) - ex p P / VI7- q(T) + Q' .o J V* ~ (t.3.17) The prime on Q' will be dropped for the remainder of the discussion of g (r ). In spherical coordinates the cylindrical symmetry in cp n is not altered by the change in origin provided the polar axis remains in the Q direction, x The integrand of g (r _) is the difference of two terms depending on g. This difference will be small if both terms are small or nearly equal so we consider the behavior of F 13 (Q,0) = e ■PV(Q) (4.3.18) with respect to Q. This function is the limit of both of the terms in the square brackets in Eq. (4.3-16) for r(f) = q(i) = 0. For p n < a, 71 F, _(Q,0) goes to zero rapidly. The function has a maximum at p - 1.12a ' 6 \ ^ and approaches one asymptotically as I 1 + — z" 1 for p » a. The *V integrand in Eq. (4.3-15) should thus be zero for p n « o and approach zero again for p n » a as both terms in the difference approach unity. A simple probability distribution with this same qualitative behavior is the x distribution with n degrees of freedom for n > k. For the 2 £ integration in g (r ip ) then a x distribution with 12 degrees of freedom was used for p n and £~ = cos 9 n was chosen from a uniform distribution, -1 < £ n < 1. Thus we define r e (S)d 2 Q = ^2 exp I" -H dp Q d| Q' ( ^ 3 * 19) The number of degrees- of freedom and the value of y were chosen P empirically during preliminary computations to minimize the variance of g (r p ;M) for small M. In initial calculations of g (r 1p ;M) it was found that g (r p ) was indeed small compared to g (r ) but that the variance of g (r p ;M) was quite large. It was thus necessary to apply a further variance reducing technique to the Monte Carlo sampling in the calculation of g . The technique used was the method of antithetic variates described by Hammersly and Handscomb [26], and although it applies to the Monte Carlo sampling as a whole, it is best described in connection with the sampling in the Q integration. The principle of 72 the antithetic variates method is to generate several correlated samples simultaneously in such a way that their sum has a considerably smaller variance than the samples taken independently. In the present case it is recognized that g (r „) is small so that one tries to, correlate the samples in such a way that their sum is near zero. From Eq. (4.3-16) it can be seen that the cylindrical symmetry of G (r 1P >Q) about the Q axis requires that G e (r 12'-S) = G e (r 12' p Q'*Q + f'V' {k ' 3 - 20) Thus the range of 6 Q can be restricted to < 9 Q < .n y (0 < £ Q < l) and the integrand taken as the sum of G (r, p ,Q) and G (r „,-<^). Another pairing of samples can be made due to the fact that the probability of a Wiener path q(x) is the same as the probability of (k) -q(-T). Since successive path approximations q are chosen (k) (k) independently the samples for q and -q may be combined. From Eqs, (4.1.28), (4,1.29), and (4.2,4) it is seen that the integrand of the multiple integral approximation for G (r.p,£) as given in Eq. (4.3.16) is L e(£ w, s « g) = 2 ^ 13 ( S . f a (k) + - 2 -)- \^> s (k) ] > 73 where D 13 (^S (k) ) = exp *i 1=1 { r«^) +kY [r« ./^ + 5i-i h-Q ^ltH + = Y l/. x 3i + Si-n + Q (4.3.22) (k) (k) Then if r is generated from the Markov chain, q) is generated by the straightforward scheme for Wiener paths, and Q, is chosen from the distribution r (Q) defined by Eq. (4.3.19) but with | n restricted to the range < £_ < 1, then with probability one g e (r 12 ) -llB-JZ A(r[ k) , q[ k \ Q t ) r" 1 ^) , (4. 3 .2 3 ) M-*» t=l where the antithetic variates sample ./ (k) (k) ^ N A (£ > <1 ; S) is defined by V£ (k) . S (k) - S) - W°. a'"'' -S) A(r<*\ *« S ) = | (U. 3.2"0 7^ This sum should be more nearly zero than its component terms since the portion of L which is antisymmetric in Q will cancel between terms 1 and 2 and between terms 3 and k on the right side of Eq. (k.3.2k) while similar cancellation will take place for the component (k) of L which is antisymmetric in the variable q between terms 1 and 3 and between terms 2 and h. 75 5. NUMERICAL RESULTS 5 .1 Results and Accuracy Considerations Using the methods developed in the preceding sections g (*%„; 20,000) and g (r ', 20,000) were calculated for temperatures T = 5°K, 10°K, 20°K, 30°K, 40°K, 50°K, 75°K, 100°K, and 273.l8°K and for values of r in the range 0.6 < r.p/a < 4.0. The third virial coefficient C(T) was computed for all "but the lowest temperature, 5 K, using Eq. (2.3°20)° A straightforward numerical integration was performed over r, p with g, (r^p) approximated by g (r,p; 20,000) + g (r, p ; 20,000) and with the values of Wp(r,p) and b p given in [6], The computations were performed on the ILLIAC II computer located at the University of Illinois. The details of the computation and a comprehensive table of results are presented in Appendix Be A summary of the results for g-^r^; 20,000) = g c (r 12 ; 20,000) + g (r 12 ; 20,000) over the range of temperature is shown in Figure 3> and the computed third virial coefficient C(T) is plotted versus temperature in Figure h< i There are two sources of computational error which cause g (r ip ; 20,000) and C(T) to differ from the correct theoretical values of g-,(r,p) and the third virial coefficient for the interaction potential of Eqs . (2.1.1) and (2,1,2). These are the error introduced by the k-fold integral approximation to the Wiener integrals and the (k) Monte Carlo sampling error. The order k of the approximations r (k) and q to the Wiener paths r(t) and q(f) was chosen for each temperature T so that Tk > 260 K, Thus the accuracy of the multiple 76 4.0- 3.0- 2.0- 1.0- ■1.0- ■2.0-- -3.0- -4.0- 12 T = 273.18°K ,-3 FIGURE 3. X g^hz ; 20,000) AS A FUNCTION OF r^/C AND T, 77 100-- r I cm 6 u \ mole2 4 600 ■■ 500-- 400 -• 300 -• + 200 f + + + + 50 100 150 200 250 300 FIGURE 4. COMPUTED THIRD VIRIAL COEFFICIENT AS A FUNCTION OF TEMPERATURE. * T(°K) 78 integral approximations at any temperature is at least equal to the accuracy of the two term Wigner-Kirkwood approximation (first two terms in Eq. (3 • 3<25)) at 260 K. As an estimate of this accuracy the third term in the Wigner-Kirkwood approximation to the second virial coefficient of He at 26o°K amounts to only 0.7$ of the total. It will be seen that an error of this magnitude is negligible compared to the Monte Carlo sampling error. As an estimate of the Monte Carlo sampling error the sample standard deviations of the sample means, g (r, p ; 20,000) and g (r, p ; 20,000), were computed. The numerical results for g, (r-. p ; 20,000) and the associated error estimates are shown plotted against r, p /a for T = 5°K, 20°K, 75°K, and 273 > 0.5 1.0 1 ,1.5 T 2.0 I I 4— » 1—1 1 — +— t — + — I 1 1— ( 1 1 1 i- 2.5 3.0 £l2 (X k -3 FIGURE 5. X g 1 AND SAMPLING ERROR FOR T=5°K. 80 \" 3 4 3.0- 2.5- 2.0- 1.5- 1.0- 0.5- ■0.5- ■1.0 •- i". I I I T T I 11 ] 1 IS 9.0 tor" 1 - ir -t— i — i i i — i t i i t i t t t — i > i i t i i i i i T f 0.5 1.0 1 1.5 2.0 £2.5 3.0 FIGURE 6. X" 3 g AND SAMPLING ERROR FOR T=20°K. '12 81 A X 9i 2.0 1.5 ■■ 1.0- 0.5- -0.5 - -1.0- -1.5 ■ -2.0- -2.5-- 0.5 1.0 1.5 2.0 2.5 3.0 -t— i — i — i | i i — i i I i i i — i — (— t — till — i — I— i — i — | — i—i — i — i — H r i2 I II I II i i h l FIGURE 7. X" 3 g AND SAMPLING ERROR FOR T*75°K. \-\ 4.0 ■ • 3.0- 2.0 ■■ 1.0 ■■ 82 0.5 1.0 1.5 2.0 2.5 3.0 -• — i — i i i — i — i — i — i — i > i — i — i — i i i — •— i — i — i — t— i — i — i— i — i — i — i — l i > lis. a I * I I -1.0- 2.0- I -3.0 ■ -4.0 ■■ I FIGURE 8. X" 3 g AND SAMPLING ERROR FOR T=273.18'K. 83 integration will increase with decreasing temperature. A factor of (3 also multiplies the exponent of the functional integrand of the Wiener integration so that the variance of the Monte Carlo procedure for computing the multiple integral approximation to the Wiener integral will also increase as T decreases as a result of the increase in (3. The variance of the multiple integral approximation is further influenced by the fact that the Wiener random functions r( t) and q(l) in Eqs . (4.1,8) and (4ol.ll) are multiplied by the thermal wavelength \. Since \ increases with decreasing temperature the variance again increases as T decreases. A final contribution to the same behavior of increasing variance with decreasing T results from the fact noted above that k was increased with decreasing T in order to maintain the accuracy of the multiple integral approximation to the Wiener integrals. The multiplicity, 6(k-l), of the integral thus increases with decreasing T resulting in an increasing variance for the Monte Carlo procedure. The standard deviation of the multiple integral results also (k) depends upon the choice of the initial vector r of the Markov chain {r, ' )« Due to the one step memory inherent in the chain the (k) choice of an initial vector r ' which is improbable with respect to the probability density P(r^ k ')du v (Eqs. (h.1.26) and (U.2.5))can (k) ' cause a large variation in the integrand as the vectors r move (k) (k) from r toward the region of maximum P(r' )du , . An attempt was ~o ° ~ ' rk (k) made to reduce the error caused by a bad choice of r for each pair J ~o * (k) of values of r, Q and T. At high temperatures P(r )du , assigns the 81+ (k) highest probability to vectors in the neighborhood of r = provided that P(r ) is slowly varying over the range of r, p + r . Thus for (k) be r = 0. At lower temperatures the initial vector for the largest ~o ° the highest temperature and largest r, p the initial vector was taken t< l1 vi (k) value of r, ~ was taken to be the last vector r.„ generated by the Markov 12 ~M D J chain for the same value of r at the next higher temperature. For successively smaller values of r, p the initial vector was taken to be the last vector generated at the same temperature and the next higher value of r . It can be seen from the tables in Appendix B that g (r, p ; 20,000) is small with respect to g (r, p ; 20,000) as was to be expected from the fact that g (r, p ) vanishes to order T in the Wigner-Kirkwood expansion, \ g (r, p ; 20,000) is, on the average, an order of magnitude smaller than X g (r, p ; 20,000). The standard deviations of these quantities, however, are of approximately the same size so that the effect of adding g to g in the approximation for g, is primarily an increase in the expected error in g, . The quantities g (r, p ; M) and g (r p ; M) were also computed with M < 20,000 at several pairs of values for r, p and T, The standard deviations in these quantities were observed to decrease as M 2 as would be expected of a well -formulated Monte Carlo procedure. 85 5.2 Numerical Comparisons with Theory and Experiment It can be seen from Figure 3 that the position of the minimum in g ,(r, p ) remains fairly constant at r - l„8a over the entire temperature range . The position of the minimum in g (t, ) can "be interpreted in terms of a crude argument about the spatial arrangement of three particles . If one takes the most probable arrangement of particles as a linear arrangement with each particle falling at the minimum of the potentials due to its two nearest neighbors, then for the Lennard-Jones (12^6) potential there should be maximum probability of finding a particle at distances of 1.12a and 2.2ka from a reference particle. Correspondingly there should be a minimum probability of finding a particle midway between these two positions or at 1.69a. If W p (r,p) is slowly varying near this minimum of the pair distribution function, then the corresponding minimum in g (r, p ) should fall at approximately the same value of r, p . The calculated position of the minimum, 1.8a, agrees well enough with this result to indicate that the minimum can be interpreted physically as the gap between the first and second shells of atoms about the reference atom. Classical and two term Wigner-Kirkwood approximations for g (r ) have also been claculated using respectively the first -pa, and first two terms of Eq. (3°3«25) for <1,2,3 |e J | 1,2, 3>. 86 The integration over the position of the third particle was carried out by a straightforward iterated Simpson's rule approximation. Comparisons between the classical, Wigner-Kirkwood, and quantum calculations are shown in Figures 9, 10, 11, and 12. The classical values of g, (r p ) have previously been calculated by Henderson [28], and the results presented agree with his calculations wherever there is overlap. In observing these comparisons it should be noted that when g (r ) is used to calculate a physically measurable quantity it is multiplied by W p (r, p ) which goes to zero rapidly with decreasing r p for r. p /a < !• Thus values of g (r. p ) for r, p < 0.8a are not significant. The two term Wigner-Kirkwood approximation is fairly accurate down to about T = 10 K. The position of the minimum in g,(r. p ) corresponding to the intershell gap remains approximately constant for the Wigner-Kirkwood data while a significant shift toward smaller values of r, p with decreasing temperature occurs in the classical position of this minimum. At 5 K the Wiener integral method gives results which differ markedly from the two term Wigner-Kirkwood approximation. Results were obtained for only a few values of r, p at T = 5 K due to the amount of computation time involved in the Monte Carlo procedure for large values of the order, k, of the Wiener integral approximation (k = 52 at T = 5 K; . Considering the a posteriori information on the accuracy of the two term Wigner-Kirkwood approximation it would seem that reasonably 87 4 X- 3 9l 1.5-- 1.0 ■■ 0.5- 0.5- WIGNER KIRKW00D 2-TERM CLASSICAL- f 4. ) 1 1 1 1 1 I I 1 1 »— I •— I 1 1 1 1 1 1 1— < + 1 1 1 1 1 1 1 1 +■ ri2 0.5 1.0 1.5 + 2.0 2.5 3.0 + * + V QUANTUM CALCULATIONS \-3_ FIGURE 9. COMPARISON OF X'^WITH CLASSICAL AND SEMI-CLASSICAL VALUES FOR T=5°K. 88 ■0.5- -1.0 - WIGNER KIRKW00D 2 -TERM QUANTUM CALCULATIONS FIGURE 10. COMPARISON OF X^WITH CLASSICAL AND SEMI- CLASSICAL VALUES FOR T=10°K. 89 X q 1 3.0-- 2.5 ■ 2.0 - ■ 1.5 ■- 1.0- 0.5- WIGNER KIRKW00D 2-TERM -0.5- -1.0- 12 QUANTUM CALCULATIONS FIGURE 11. COMPARISON OF Xg 1 WITH CLASSICAL AND SEMI-CLASSICAL VALUES FOR T=20°K. 90 \" 3 2.5- 2.0- 1.5- 1.0- 0.5- -0.5- -1.0- -1.5- ■2.0- WIGNER KIRKW00D 2-TERM £12 CLASSICAL QUANTUM CALCULATIONS FIGURE 12. COMPARISON OF X" 3 g 1 WITH CLASSICAL AND SEMI -CLASSICAL VALUES FOR T=50°K. 91 accurate results could be obtained with Tk > 30 K rather than Tk > 260 K so that it would be possible to repeat the 5 K computation in one tenth the time. There seems to be no experimental data on the pair k distribution function in He gas at low temperatures. The pair distribution function for liquid helium, however, has been measured by neutron diffraction (Hens haw [29]) for several values of the density. Figures 13 and lU show W 2 (r 12 ) [1 + n gl (r 12 )] - d ^ = — &- (5-2.1) n o for T = 5 K and number densities n corresponding to mass densities p = 0.095 gm/cc and p = 0.18*4- gm/cc respectively compared with Henshaw's data for these same densities. Since only 3 -particle effects have been taken into account in the computation, agreement cannot be expected beyond the second peak in the distribution function, To indicate the variation of p(r,„)/p with temperature, smooth curves have been fitted to the numerical results for W 2 ( ri2 ) [1 + n g^r^)] - -j±- (5-2.2) for temperatures 10 K, 20 K, and i+0°K at a density of p Q = 0.184 gm/cc. The resulting curves are shown in Figure 15 • 92 p(hz) 2.0- 1.5-- 1.0-- 0.5- THIS WORK , T=5° K , /D = .095 gm/cc HENSHAW [29] , T = 5.04°K , /> = .095 gm/cc 12 FIGURE 13. THE PAIR DISTRIBUTION FUNCTION FOR AVERAGE DENSITY /o =. 095 gm/cc. 93 pihz) t * 2.0- I THIS WORK , T = 5°K , p = . 184 gm/cc HENSHAW [29] , T =4.2° K , p = . 184 gm/cc 1.5 1.0-- 0.5- I 1 1 1 1 i 1 1 4* i I 1 1 1 — 1 1 1 — 1 — 1 — t I 1 — 1 — 1 1 1 — 1 — 1 — 1 — 1 — I — ^ ri2 1.0 2.0 3.0 FIGURE 14. THE PAIR DISTRIBUTION FUNCTION FOR AVERAGE DENSITY p =.184 gm/cc. 9 k p(hz) ♦ P T= 10° K , /> =.184 gm/cc T = 20°K ,p Q =. 184 gm/cc T = 40°K , p =.184 gm/cc Ti2 FIGURE 15. THE PAIR DISTRIBUTION FUNCTION AS A FUNCTION op TEMPERATURE FOR p =.184 gm/cc. 95 Two factors influence the comparison between Henshaw's neutron diffraction data for liquid He and the numerical results of this work. The first is the choice of the Lennard-Jones (12-6) potential, Eq. (2.1.2), as the correct intermolecular potential for helium. The second is the validity of the cluster expansion for temperatures and densities at which helium is a liquid* The influence due to each of these factors cannot be separated without recomputing W p (r ) and g (r, p ) for a different intermolecular potential. However, barring compensating errors in these two approximations, Figures 13 and ~\.K indicate that both approximations are fairly good for the pair k o distribution of liquid He at 5 K. Using Eq. (2.3-20) the third virial coefficient C(T) was evaluated for all temperatures at which g-,(r,p) was calculated except for 5 K. C(5 K) was not calculated partly due to the large sampling errors in g (r ) at 5 K but primarily due to the fact that values of g (r ) for r > 2 =8 a give a large contribution to the integral for C(5 K) and no numerical results were obtained for this range of r 1C) « Figure l6 shows the numerical results for C(T) compared to the values adopted by Keesom [30] for the experimental third virial coefficient and the values calculated classically from the Lennard-Jones (12-6) potential (Hirschf elder [1]). Over almost all of the temperature range the computed values fall consistently below the experimental results. Since Eq. (2.3>20) does not depend on the neglect of higher order terms in the cluster 9 6 r f_cmL\ 4 u Vmole 2 / 600- 500- 400 ■■ 300 ■ 200 - 100 ■ EXPERIMENTAL DATA , KEESOM [30] THIS WORK CLASSICAL RESULTS , HIRSCHFELDER [l] 50 +— *T(°K) 100 150 200 250 300 FIGURE 16. COMPARISON OF CALCULATED THIRD VIRIAL COEFFICIENT WITH CLASSICAL AND EXPERIMENTAL VALUES. 97 expansion the cause of this difference is evidently the choice of the pairwise additive Lennard- Jones (12-6) potential as the correct 3-particle potential for helium. Some recent calculations by Sherwood, De Rocco, and Mason [31] indicate that nonadditive 3 - body forces have a sizeable effect on the third virial coefficient. The ways in which g (r ) and C(T) differ from the experimental values can be reconciled qualitatively. Figures 13 and Ik indicate that the calculated g, (r, p ) is fairly consistently larger than the experimental g (r p ). On performing the integration over r in Eq« (2.3»20) this would lead to a C(T) which would be smaller than its experimental values as is verified in Figure l6„ 5«3 Conclusion The comparisons of Section 5=2 point to two conclusions k regarding the physics of the He system. The first is that quantum mechanical calculations based on the cluster expansion and assuming a pairwise additive Lennard -Jones (12-6) potential can produce a reasonable qualitative fit to the pair distribution function in k liquid He . Secondly, the consistent deviation of the calculated third virial coefficient from the experimental values over the full temperature range indicates that the choice of the latter potential k for the 3-body interaction in He is not good enough to yield quantitative agreement with experiment, 98 From the mathematical viewpoint it has "been shown that the Wiener integral may "be a useful computational tool. The upper bounds of Section 3-2 indicate that the expression of a quantity in terms of Wiener integrals can yield qualitative and semiquantitative information on the size and behavior of that quantity. Finally, the numerical evaluation of a Wiener integral expression has been shown to be computationally feasible in a situation where the amount of computation involved in integrating the original partial differential equation giving rise to the Wiener integral would be prohibitive. With the numerical techniques developed it would be interesting to investigate the effect on the third virial coefficient of a change in the intermolecular potential. A quantum mechanical calculation of the third virial coefficient might offer a way of comparing the potentials which are indistinguishable with respect to the second virial coefficient. Another possible extension of this work would be an investigation of the behavior of the exchange terms in the 3-particle density matrix. Of course, both direct and exchange terms can be computed for systems of more than three particles using the same formulation, but indications are that the necessary computation increases greatly with each added particle. 99 LIST OF REFERENCES [1] J. 0. Hirschfelder, C. F, Curtis, and R t Bo Bird, Molecular Theory of Gases and Liquids , Chap, 6, John Wiley and Sons; 19 5 k. [2] S. Y. Larsen, "Quantum Mechanical Calculation of the Third Virial Coefficient of HeV' The Physical Review , Vole 130, No. k, pp. li+26-lUUO; May 1963. [3] A. Pais, and G, E, Uhlenbeck, "On the Quantum Theory of the Third Virial Coefficient," The Physical Reviev , Vol. Il6, No. 2, pp. 250-269; October 1959- [k] C. N. Yang, and T. D. Lee, ,! Many-Body Problem in Quantum Statistical Mechanics. I. General Formulation," The Physical Reviev , Vol. 113, No. 5, pp. II65-H77; March 1959- [5] S. Y. Larsen, K. Witte, and J t E. Kilpatrick, , "On the Quantum- Mechanical Pair-Correlation Function of He Gas at Low Temperatures," Journal of Chemical Physics , Vol. kk, No. 1, pp. 213 -220; January 1966. [6] Lc D. Fosdick, and Ho F. Jordan, ''Path Integral Calculation of the Two Particle Slater Sum for He , " Th e Physical Review , Vol, 1^3, No. 1, pp. 58-66; March 1966 . [7] J- de Boer, and A, Michels, "Quantum-Mechanical Calculation of the Second Virial-coeff icient of Helium at Low Temperatures . ;; Phy sic a, Vol 6, No, 5 , PP° ^09-^20; March 1939 . [8] J. de Boer, Molecular Distribution and Equation of State of Gases," Reports on Pr ogress m Physic s, Vol, 12, PP. 305-37^, 19^9- [9] B. Kahn, and G, E. Uhlenbeck, "Theory of Condensation, ' Physica , Vol, 5, pp. 399-^16, May 1938, [10] K. Huang, Statistical Mechanics , Chap. lk, Sects. lA.l and lk.2, John Wiley and Sons; 1963 . [11] S. G, Brush, "Functional Integrals and Statistical Physics," Reviews of Moder n Physics, Vol . 33, No- 1, pp. 79~92; January 1961. 100 [12] I. M. Gel 'f and, and A. M. Yaglom, "integration in Functional Space and its Applications in Quantum Physics," Journal of Mathematical Physics , Vol. 1, No. 1, pp. 48-69; January -February i960. [13] M. Kac, Lectures in Applied Mathematics, Vol. I, Proceedings of the Summer Seminar, Boulder, Colorado, 1957 > Chap. IV, Interscience Publishers, Inc.; 1958. [lk] J. Neyman, Editor, Proceedings of the Second Berkeley Symposium on Mathematical Statistics and Probability, Berkeley , University of California Press; 1951 • [15] R. P. Feynman, and A. R. Hibbs, Quantum Mechanics and Path Integrals , Chap. 10, McGraw-Hill, Inc.; 1965. [l6] P. Levy, "Le Mouvement Brownien, " Memorial des Sciences Mathematiques , Fascicule 126, Gauthier-Villars, Paris; 1954. [17] S. Y. Larsen, J. E. Kilpatrick, E, H. Lieb, and H. F. Jordan, "Suppression at High Temperatures of Effects Due to Statistics in the Second Virial Coefficient of a Real Gas," The Physical Review , Vol. 1^0, No. 1A, pp. 129-130; October 1965. [18] E. H. Lieb, 'Calculation of Exchange Second Virial Coefficient of a Hard-Sphere Gas by Path Integrals," Journal of Mathematical Physics , Vol. 8, No. 1, pp. ^3-52; January 1967 ^ [19] R. H. Cameron, A Simpson's Rule for the Evaluation of Wiener's Integrals in Function Space," D uke Mathematical Journal , Vol. 18, pp. 111-130; 1951. [20] L. Do Fosdick, "Numerical Estimation of the Partition Function in Quantum Statistics, 1 '' Journal of Mathematical Physics , Vol. 3, No, 2, pp. 1251-1264; November-December 1962." [21] L. D- Fosdick, "Approximation of a Class of Wiener Integrals," Mathematics of Computation , Vol. 19, pp.. 225-233; 19^5 • [22] I. M. Gel'fand, and N. N. Chentsov, "The Numerical Calculation of Path Integrals," J. Exptl. Theoret. Phys . (U.S.S.R.) , Vol. 31, pp. 1106-1107; December 1956. iCl [23] A. Go Konheim, and Wo L. Miranker, "Numerical Evaluation of Wiener Integrals,'' I,B,M, R esearc h Report , R0-1540; January 1966. [24] T, Kihara, Y, Midzuno, and T, Shizume, "Virial Coefficients and Intermolecular Potential of Helium/ 1 Journal of the Physical Society of Jap an, Vol. 10, No, 4; April 1955. [25] N. Metropolis, A, W. Fosenblu+h, M. N, Rosenbluth, A, H, Teller, and E. Teller, 'Equation of State Calculations "by Fast Computing Machines,' The J ournal o f Ch emic a l Physics , Vol, 21, No, 6, ppo 1087-1092; June 1953- [26] J, M. Hammers ly, and D, C. Hands comb, Monte Carlo Methods , John Wiley and Sons ; 195^° [27] J. L. Doob, Stochas tic Processes , Chap, 5, John Wiley and Sons; 1953- [28] D. Henderson, "Exact and Approximate Values for the Radial Distribution Function at Low Densities for the 6° 12 Potential," Molecul ar Physics, Vol, 10, No, 1, p, 73; 1965 = [29] D. G, Henshaw, "Pressure Effect. In the Atomic Distribution in Liquid Helium by Neutron Diffraction, " The Physical R eview , Vol, 119, No, 1, pp. 14-21; July i960, [30] W. H. Keesom, Helium, Chap, 2, Elsevier; I9U2. [31] A, E, Sherwood, A, G, DeFocco, and E. A., Mason, "Nonadditivity of Intermolecular Forces; Effects or the Third Virial Coefficient,' Journal of Chemical Physics, Vol. h-k, No. 8, pp, 2984-2994 J April I966. 102 APPENDIX A LEMMA FOR THE PROOF OF Eq. (4.2.21) Lemma: p (x'|x)w(x) = p (x |x' )v(x' ) (4.2.20) Proof: From Eqs . (4.2.11) and (4.2.13) it follows that p (x x jw^x ) k-1 k 2 exp k-1 -k Z 1=1 T ,T x. . + x. . i-l i+l i=l J (A.l) T T By the definition of x , x. = x ., i=0,l, ...,k, so Eq. (A.l) becomes ~ 1 K. - 1 p (x x )w(x' ) i k 2 exp [- k 1 (vi x. . , + x n . ,\p k-i+1 k-i-1*^ i=l - 1 (A.2) 103 Setting £ = k-i, one finds that p (x |x' ;w(x' ; k\ 2 k-1 1 2 (2«r 27' 2 k-1 k-1 -k L Xo + k L x pX „ -, £=1 l £=1 l l k-1 k Z £=2 x „x r£-± k k £ x ili + i + x i-i } ? i=i h ^ i=o * + k Jli X £ X £ +1 " 2 f =1 ^ (A. 3) and making use of the fact that x = x, = x 1 = x' = 0, this "becomes & o k o k p (x x' Jw^x ) k-1 1 k ^ 2 ' to \ 2 l*\ 2 (2k) [^) exp k-1 k-1 -k Z (xjr + k Z x'. L i=l j=l J j-l k-1 + k Z x'.x . , .1=1 J J+1 k k z' (x ^i + x i-i V i=l k-1 2 k -k Z x „ + Zx.x. i=l J=l J j-l (AA) where j = i+1. Equating the dummy subscripts and combining the sums, one obtains 10^ P (x |x' )w(x' ) k-1 1 / \ k ^• 2 (-) 2 ^) 2 exp[- k |f x .-!!M4ik X2 2 ^ (x i " X i-1 ) J ' p (x 1 |x)w(x) , (A. 5) which completes the proof, 105 ..APPENDIX B DETAILED COMPUTATIONAL PROCEDURE AND NUMERICAL RESULTS FOR g c (r 12 ;M), g e (r 12 ;M) and C(T) The step by step procedure for computing the approximation g^(r 12 ;M) to g^(r^) is as follows: 1. Generate M independent Wiener path approximations (k) consisting of 3k-dimensional vectors q)_ ' 3 t = 1,2,.,.,M, drawn from a probability distribution du by Eq. (^.2.2), q = q = 0, 3i 'li 1 - 1 - i-i + ii k W\ i - 1 , i = L,2,...,k-1 , (B.l) where the |. are 3 _ dimensional random variables, the coordinate variables of which are independent and normally distributed with mean and variance 1. (k) 2. Generate an M step Markov chain [r)_ t=l,2, . . . ,M) having the stationary probability D 12 (£ (k) ) W Q (1,2) d ^rk by first generating r; using the procedure described by Eq, (^.2.28) io6 ind setting r^*' = r^ if D^r^') > D^fet-l^ while if D 12^t k) ) , v,k-t + °12 (£ t ) r> = iv with probability 7TY~ D (r { } ) and D (r (k) ) (k) (k) ... KM1 „ , °12^t } rl = r : with probability 1 - , \ D (r l k) ) 12 v ~t-l y 3- Generate M independent random vectors Q = (Q ,Q , 0) , ~t x y t 2 t = 1,2, ... ,M, according to the probability T (Q)d Q. h. Form the average M (^JM) .. i £ L c (r« , q< k » , Q t ) r; 1 ^) . (B.2) The equations pertinent to the above procedure are t 1 (k) (k) m L ^3 (k) 1 (k)^ . "I L c ( -r > r > s) = L D i3\^2 5 + 2 £ y J x [ D 2 3 (^s (k) - 2 1 £ (k) )- 1 ] > C^'3-5) 107 V£ (k) > exp (- & L [ "Fi-.^J-wfil^r^ria w« *• *•" „, 2* V » NTT x / Si + Si-i -12 ]}■ (4.1.29) V(r) = kK [(?) 12 I, K = 14.04 x 10" erg, a = 2.56 x 10" cm , (2.1.2) and T (Q)d 2 S = — i— exp ^2*7 x \ 1 r 27 ' / JzH- exp Tt/ 2/ i dQ x dQ y (4.3-9) The procedure for computing the approximation g (r, p ;M) to g (r, ) is the same as that for g (r, ;M) except that steps 3 and 4 e ±-cL c Ld are modified. 108 3. Generate M independent random vectors Q = (p A ,0 A ,O) , ~t 4 > Q) + L (£ * 4 * -Q) r / (k) (k) n s _ , (k) (k) _s " + L e (r v ', -q v ', Q) + L e (r v , -q v ', -Q) , (U.3.210 l.(sW, S 0J^, ■ 2 [, 13 ( S , £ 3 (k) ♦ I ^ (_k ) - "uCS. 2 (k) ) ] > D 13 (Q,q (k!) ) = exp * £[ T te*« + *) ;+ lfV \fSi s^Tjt + v fe^H-OK^H] (^.3.21) W;.) (1^.3.22) and T (Q)d 2 Q = — % exp e ~ ~ 6ly Ct) % d| Q " ? -S 6 ' M) Std. Dev. W 2 (r L2 ) 1.0 0.11161 0.05+33 0.01^93 O.O5076 O.Ul620 1.3 +0.00768 0.05091 0.02307 0.03620 1.39898 1.5 -0.09163 0.0^704 +0.0260+ 0.05+86 1.5252 + 1.7 -0.155^2 0.0^17+ -0.00652 O.O38U2 1.+1138 1.8 -0.16863 O.O3699 -0.02636 0.03781 1.33505 1.9 -0.10887 , 0.03502 , -0.03+88 0.05512 1.26+86 2.1 -0.10791 '■ 0.03155 +0.02816 0.0+026 1.15^33 2.3 +0.06+28 0.02298 -0.01*913 0.039^0 I.O8719 2.5 0.15227 O.OI892 -0.02889 0.0^606 1.0^950 2.8 0.20509 0.01092 +0.0U91+ 0.0++30 1.02273 Ill Table 1.2. NUMERICA1 VALUES OF k'^g (r ;20,000), \" 3 g e (r 19 ;20,000), AND W 2 (r 12 ) FOR T = 10°K, r J 12 a 3 g c (r !2^ M) Std, Dev. 75 g e (r i2' M) K Std. Dev. W ? (r 12 ) 0.6 0.30204 0.05718 -0,10363 O.O5672 0.00000 0.7 0.43065 0.05918 . 00001 0.8 0.22653 0.05944 -O.H826 0,04663 0,00845 0.9 0.20^32 0.05659 0.11992 1.0 0.21626 0.05623 +0.02633 0,03522 0.43502 1.1 +0.00022 0.05535 0,84777 1.2 -0,09561 0.05339 0,11547 0,03143 1,17350 1.3 -0,13564 0.05041 0,00588 0,03157 1,33885 1.4 -0.26533 0.04638 0.02020 0.03433 1,38106 1.5 -0.31497 0.04414 0,03717 0.03018 1=34965 1.6 -O.43696 0.03915 0.00436 0.03249 1,28194 1-7 -0.^3063 0.03585 +0.03937 0,06796 1,21369 1.8 -C, 48012 0,03213 -0.12736 O.O3236 1,15856 1.9 -0.45266 O.02916 1.11587 2.0 -O.41747 O.02496 +0.0106l 0,03334 1 . 08438 2.1 -0,25859 0.02162 1,06163 2.2 -O.I878O 0.01845 +0.04063 0,03334 1.04551 2.3 -0.09946 0.01625 -0,01458 0,03264 1,03408 2.4 +0.00139 0.01320 +0.02194 0,03164 1.02588 2.5 0.07724 0,01100 +0,05166 0.03082 .1 , 01991 2.6 0.10874 0.00953 -0,03222 0.03318 1.01550 2.7 0.l4l66 0.00755 +0.00224 0,03226 1.01220 2.8 0.13137 0.00595 +0.00217 0.032.17 1 . 00970 2.9 O.14632 0,00471 1,00778 3.0 0,11885 O.OO367 -0.00508 0.033^8 I.OO629 3-2 O.08919 0.00214 1.00420 3.* 0,0667k 0.001^9 1,00289 3.6 0.04890 c, 00099 1.00203 3-8 0.03304 0,00072 ' 1,00145 4.0 0,02155 0.000^7 1 . ' 0106 112 Table 1.3. NUMERICAL VALUES OF X 3 g (r ;20,000), \" 3 g e (r 12 ;20,000), AND W (r^) FOR T = 20°K. r !2 ^Jg c (r l2 ;M) Std. Dev. ^3 g e (r 12 ;M) Std. Dev. W^r^) 0.6 1 . 7702 0.08450 -O.62796 0.07301 0.00000 o.t 1.5075 0.08448 -0.21953 0.05l6l 0.00000 0.8 1.4115 0.08353 -O.HO3O 0.04157 0.00827 0.9 1.1673 . 08096 -0.08044 0.03704 0.14649 1.0 O.78071 O.07988 +O.O2763 0.03097 0.52629 l.i 0.53W O.07558 0.02948 0.03080 0.94670 1.2 +O.21673 O.07107 0.08262 0.02909 1.19770 1.3 -0.013^6 0.06457 O.O3672 0.03007 1.26763 1.4 -0. 31^30 0.05738 0.01662 0.03105 1.24238 1.5 -0.48981 0.05180 +0.03585 0.02992 1.18817 1.6 -0.76050 0.04506 -0.02124 O.03089 1.13647 1.7 -0.85461 O.03851 +0.05727 0.02978 1 . 09690 1.8 -0.90454 O.03174 -0.02649 0.03009 1.06884 1.9 -0.84328 O.02588 -O.O2634 0.03156 1.04937 2.0 -0.722.40 0.02248 -0.03736 0.03103 1.03591 2.1 -0.58432 0.01759 -O.06747 0.03031 1.02652 2.2 -0.38289 0.01473 -0.01376 0.03141 1.01986 2.3 -0,23783 0.01188 -0.01845 O.O2985 1.01508 2.1+ -0.08979 0.00993 -0.01048 0.03039 1.01159 2.5 -0.01504 0.00764 +0.01846 0.03114 1.00901 2.6 +0.03631 0.00593 -0.04298 0.03033 I.OO708 2.7 0.06169 0.00461 -0.01156 0.02973 1.00561 2.8 0.06026 O.OO369 +0.03174 0.02972 1.00449 2.9 0.05736 O.OO298 +0.02549 0.03214 1.00362 3.0 0.05746 0.00232 -0.00994 0.03077 1.00294 3-2 0.04262 0.00149 I.00199 3-h 0.03200 0.00101 1.00137 3.6 0.02172 0.00069 1 . 00097 3.8 0.01541 0.00051 1.00070 U-.O 0.01030 0.00038 1.00051 113 Table 1.4. NUMERICAL VALUES OF A~ 3 g (r ;20 ; 000), V " 3g e (r 12 ;20,00 °- ^ AND W 2^ r l2^ P0E T = 3 °° K ' 12 \ ^ g c (.r i2 ;M) Stdc Dev, j A g J v : .r 1? ;M) Std, Dev, w p' r i?.) -4- 0.6 0.7 0,8 0.9 1.0 1.1 1.2 1 = 3 1.4 1.5 1.6 1.7 1.8 1.9 2.0 2.1 2.2 2.3 2.5 2,6 2.7 2.8 2.9 3.0 3.2 3.4 3.6 3.8 4.0 3 . 7089 3.1127 2.9934 2.6:- 2,0920 1.6787 1.0346 +0.50982 -O.I29O8 -0.62747 -1,0382 -1.1272 -1.2507 -I.I69I -0,96546 -0,71601 -0.46636 -0,28688 -0,15784 -0,07918 -c, 03204 +0,01384 0,01602 0,02666 O.C1990 0.01861 0.01 O.OC919 0,00681 0.00479 0,12325 O.II960 0, 11860 0.11824 0.11073 0,10490 0,09486 O.O869O 0.07376 0.06234 0.05072 0.04l8l 0,03376 0.02662 0,02114 0,01693 0,01305 0,01032 0,00818 0,00650 0.00505 O.OO396 0.00310 0.00248 0.00198 0.00131 O.OOO89 0.0006^ 0.C0048 Oo 00037 -0.55485 -0.33827 -0,06397 -O.08559 -O.OL259 +o.oo4o3 -0.01300 +o.o6i48 0.01243 0,01058 +0,03607 -0,01810 -0,03943 -0.00269 +0,06968 0.03424 +0.00379 -0.00871 -0.00372 -0,00138 +0,0284.1 -0,03363 +0,034146 -0,01367 -0,03350 0.07420 0.05606 0.04174 0.03514 0.03C71 0.03031 0.03004 0.03104 0.03083 0,03122 0.03135 0,0295? 0.0306l Oo03l43 0,03079 0.03105 O.O3.II6 0.03214 0.03316 0.03094 O.O3255 O.C3069 0.03112 0.0312.9 0.0309 3 0.00000 0,00001 0,01251 0.17920 0.58907 0,98567 1.18254 1,21059 1.17199 1.12496 1.08763 l,06l4l 1,04.354 1.03135 1 . 02292 1 - 01701 1.01280 1 . OO976 1,0075? 1.00586 1 . 00462 I.OO367 I.OO294 1,002 38 1 , 00193 1.00131 1.00090 1.00064 1.000^6 1.00034 114 Table 1.5. NUMERICAL VALUES OF \" 3 g (r 12 ;20,000) , \" 3 g e (r 12 ;20,000), AND Wg(r^) FOR T = 40°K. 12 o Jj6 c (r 12 ;M) Std. Dev. 3 g e (r 12 ;M) K Std. Dev. W 2 (r 12 ) 0.6 6.5168 O.I6779 -I.34129 0.09346 0.00000 o.T 5.8136 0.16632 -0.30960 0.05813 0.00002 0.8 5.2420 0.16122 -0.12666 0.04290 0.01505 0.9 4.6279 0.15673 -O.O6758 0.03598 0.21818 1.0 3.7103 0.14460 +0.05650 0.03208 0.66868 1.1 2.8058 0.13353 0.05465 0.03152 1.03258 1.2 2.0697 0.12226 0.03208 o. 03006 1.16868 1.3 1 . 1047 0.10556 +0.02770 0.03180 1.16784 1.4 +0.10555 O.O8817 -O.O2226 0.03205 1.12805 1.5 -O.51087 0.07306 -0.00293 0.03048 1.09045 1.6 -1.0742 0.05946 -0.05644 O.03306 1.06301 1.7 -1.3525 0.04512 +0.03024 O.03016 1.04421 1.8 -1.4343 0.03477 0.01522 0.03201 1.03146 1.9 -1.3756 0.02575 0.01866 O.03218 1.02273 2.0 -1.1223 0.02022 O.OO65I 0.03H3 1.01667 2.1 -0.81214 0.01576 0.03297 0,03173 1.01241 2.2 -0.56541 0.01234 0.02112 O.03182 1.00936 2,3 -0.3^673 . 00976 +0.00832 O.03198 1.00715 2.4 -O.20666 0.00754 -0.03604 O.03083 1.00552 2.5 -0.12062 0.00592 +o.ooi46 O.03188 1,00431 2.6 -O.06254 0.00459 -0,04644 O.03171 1.00340 2.7 -O.03496 0.00363 +0.04956 O.03278 1.00270 2.8 -O.OH60 0.00289 0.00766 0.03141 1.00217 2.9 -0.00110 0.00231 +0.02028 0.03253 1.00176 3.0 -0.00212 0.00189 -0.06352 O.03143 1.00143 3-2 •1-0.00275 0.00127 1.00097 3^ 0.00410 0.00087 1.00067 3.6 0.00224 0.00065 1 . 00048 3.8 0.00191 0.00048 1.00034 4.0 0.00125 O.OOO38 1.00025 115 Table 1.6. NUMERICAL VALUES OF \" 3 g (r 10 ;20,000), c^12 j X " 3g e (r 12 ;20 '° 00) > AND W 2 (r l2 ) F0R T = 5 °° K ' 12 1 / M\ ^ g c (r 12 ;M) Std. Dev. \ 3 g e (r 12' M) Std. Dev. w 0.6 0.7 0.8 0.9 l.o l.i 1.2 1.3 1.4 1.5 1.6 1.7 1.8 1.9 2.0 2.1 2.2 2.3 2.1+ 2.5 2.6 2.7 2.8 2.9 3.0 3-2 3.4 3.6 3.8 4.0 9.566^ 8.7736 7.3618 6.6846 5.5319 4.4098 2.8499 I.5876 +O.33265 -O.48695 -1.1023 -1.6007 -1.6056 -1.5481 -I.2836 -O.91036 -0.60499 -0.40800 -0.24114 -0.l4l69 -0.09400 -0.05507 -O.O3636 -0.02176 -0.01912 -0.00742 -0.00353 -O.OOI89 -0. 00182 -0.00144 0.21727 0.21338 O.20678 0.19733 0.18309 0.16844 O.14929 0.12432 0.10413 0.08494 O.06671 0.04780 0.03575 0.02594 0.01960 0.01526 0.01184 0.00940 0.- 00717 O.OO561 0.00442 0.00350 0.00276 0.00221 0.00182 0,00124 O.OOO88 O.OOO65 0.00050 0.00039 -O.8389O -0.41851 -0.22029 -0.02030 +0.02428 0.01908 0.02397 0.00415 +0,04323 -O.OO670 +0.01203 -0.03959 +0.03151 0.00077 +0.02970 -O.OH85 -0.00923 -0.05462 +0,02053 +0.00280 -O.OO663 -0.01^96 +0,03584 +0,03949 -0, 06284 0.08907 O.06036 0.04423 0.03489 0.03138 O.03063 O.03038 0.03142 O.03178 0,03169 O.03228 0,03172 O.03119 0.03222 O.03230 O.03153 0,03129 0,03129 O.03092 O.03235 0.03239 0,03231 O.03251 O.03107 0.03284 ! 1 0.00000 0.00001 0.01763 0.24798 0.71891 1.05350 1.15232 1.13778 1.10136 1 . 07088 1.04933 1.03467 1.02472 1 , 01789 1.01315 00980 1,00740 1,00566 1,00438 1,00342 1,00270 1,00215 1.00172 1.00140 i,oon4 1,00077 1.00053 1 , 00038 1.0002.7 1.00020 116 Table 1.7. NUMERICAL VALUES OF \" 3 g (r ;20,000), X'^S^r^; 20, OOO), ANDW 2 (r 12 ) FOR T = 75°K. r 12 a ^ g c (r 12^ M) Std. Dev. ^3 6 e (r ]2 ;M) Std. Dev. V r 12' 0.6 0.00000 0.7 17.099 0.35252 -0.57189 0.06499 0.00005 0.8 14.592 0.33201 •- 0.08237 0.04288 0.03028 0.9 13.849 0.32565 -O.O3IO2 0.03544 0.34036 1.0 9.9851 0.28127 +0.01446 0.03462 0.82732 1.1 8.1720 O.25669 +0.07193 0.03240 1.07206 1.2 5 . 6260 0.22191 -0.02478 0.03356 1.11255 1.3 3.4497 0.18889 -O.05651 0.033++ 1.09178 1.4 +1.3873 0.14718 +0.02364 0.03500 1.06573 1.5 -O.29434 O.IIO83 0.05378 0.03440 1.04574 1.6 -I.2900 0.08048 +0.04425 0.03376 1.03187 1.7 -I.9094 0.05446 -0.01024 0.03336 1.Q2246 1.8 -2.0339 0.03728 -0.01130 0.03456 I.OI606 1.9 -1.8976 0.02542 +0.04435 0.03305 1.01166 2.0 -1.5165 O.OI976 -0.00107 0,03400 I.OO858 2.1 -1.0853 0.01530 +0.04414 0.03454 1 . 00641 2.2 -0,7H91 O.OII58 -0. 00974 0.03395 1.00485 2.3 -0.48689 O.OO897 +0.00430 0.03431 1.00371 2.4 -O.33031 0.00688 0.00171 0.03486 1.00288 2.5 -O.22195 0.00540 0.02422 0.03323 1.00225 2.6 -0.15578 0.00422 +O.OI697 0.03309 1.00178 2.7 -0.11421 0.00340 -0.00746 0.03451 1 . 00142 2.8 -0.07068 0.00268 -0.00444 0.03284 1.00114 2.9 -0.05547 0.00218 -0.04366 0,03403 1 . 00092 3-0 -0.04714 0.00182 -0,05734 0.03413 1.00075 3-2 -O.02699 0.00124 1.00051 3.4 -0.01723 O.OOO89 1.00035 3.6 -0.01146 O.OOO69 1.00025 3.8 -0.00782 0.00051 1.00018 4.0 -0.00605 0.00041 1.00013 117 Table 1.8. NUMERICAL VALUES OF X _:5 g ( r;L2 ;20,000), \" 3 g (r^;20,000), AND W^r^) FOR T = 100°K. 12 T g c (r 12' M) Std. Dev. K h g e (r 12' M) Std. Dev. W^Pjg) 0.6 0.7 0.8 0.9 1.0 l.i 1.2 1.3 1.4 1.5 1.6 1.7 1.8 1.9 2.0 2.1 2.2 2.3 2.4 2.5 2.6 2.7 2.8 2.9 3.0 3.2 3-4 3.6 3.8 4.0 28.866 26.301 23.154 20.788 16.211 12.226 8.7524 4.9980 2.5440 +0.15175 -1.2950 -2.2048 -2 . 4492 -2.1632 -1.6750 -1.1995 -0.80334 -0.57866 -0.40427 -0.27668 -0.19290 -0.14304 -O.IO312 -0.08112 -0.06210 -0.04i63 -0,02636 -O.O.1764 -0,012.12 -0.00962 0.52539 0.50866 0.47882 0.45434 0.39965 0.35103 0.30203 O.2369O 0.19197 0.14091 O.O9839 0.05886 O.O3686 0.02597 0.02022 0.01549 0.01158 O.OO898 O.OO698 0.00539 0.00426 0.00341 0.00272 0.00223 0.00182 0.00128 0.00092 O.OOO69 0.00054 0.00043 -1.42417 -0.78737 -0.18906 -0.01739 +0.04359 +0.00389 -0.02168 +0.05336 +0.01317 -O.03295 -0.02260 -O.07566 +0.0621.4 -0.01889 -0.01509 -0.06320 -0.06583 +0.02067 -0.00948 +0.01531 -0.00674 -0.02822 -O.OO883 -0.03814 +0.03840 0.11788 O.06773 0.04i4l 0.03523 O.03273 0.03494 0.03340 0.03437 0.03412 O.03369 O.03360 O.03561 0.03393 0.03417 0.03374 O.03469 0.03448 0.03428 O.03366 O.03375 0.03428 O.03388 0.03550 O.03436 03423 0.00000 0.00005 0.04034 0.40579 0.87784 1.06646 I.08760 1 . 06906 1.04907 1.03^11 1.02378 1.01678 1.0.12' 01 1 . 00872 1 , 00643 1,00480 1 . OO364 1.00278 1.00216 1 , OOI69 1.00133 1.00106 I.OOO85 1 . OOO69 1,00056 1,00038 1.00026 .1 , 00019 1.00014 1,00010 118 Table 1.9. NUMERICAL VALUES OF \" 3 g (r ;20,000), X " 3g e( r 12 ;20 ' 00 °)> MD W 2^ r 12^ F0R T = 2 T3-l8°K. 12 ^ g c (r 12^ M) Std. Dev. 3 6 e (r 12 ;M) Std. Dev. W 2 (r 12 ) 0.6 0.7 0.8 0.9 1.0 1.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8 1.9 2.8 2,9 3o0 3.2 3-4 3.6 3-8 4.0 125.44 115.17 96.324 81.748 61.584 45.317 31.259 19.126 9.6566 +2 . 3783 -2 . 0686 -3.8415 -3.7872 -3.1047 -2.2944 -1.6418 -1.1846 -0.84689 -O.60708 -0,45287 -0,35206 -O.26398 -0.20450 -O.15736 -0,12176 -O.08107 -0.05550 -0. 03776 -0.02727 -0, 02110 2.0011 1.9319 1.7268 1.5693 1.3248 I.0947 0.89040 O.66217 O.47061 O.28913 . 14422 0.07019 0.04234 0.03182 0.02359 0.01739 Oo01324 0.01005 0.00770 O.OO605 0.00483 O.OO389 0.00315 0.00256 0.00213 0.00153 0.00113 O.OOO85 0.00068 O.OOO56 -I.81087 -O.58497 -0.04505 -0.04151 -0.03488 -O.04590 +0.01854 +O.06585 -0.02573 +0.07952 +0.02492 -0.04085 +0.00200 -0.02403 +0.00917 +0.03585 -0.00987 +0.11174 -0.00914 +0.01644 O.05927 0.00352 +O.00066 -0.02901 +0.01181 0.11496 O.06149 0.0i+28l O.03984 O.04038 0.04040 0.04181 0.04014 0.04132 0.04200 0.03909 . 04049 0.04228 0.04070 0.04094 0.04239 0.04263 0.04l4l 0.03990 0.04090 0.04145 0.04258 0.04304 0.04089 0.04266 0.00000 0.00119 O.18533 O.72582 O.98181 1.03326 1.03307 1.02481 1.01745 1.0L213 1.00848 1.00600 1.00430 I.00313 I.00231 1.00173 1.00131 1.00101 1.00078 1.00061 1.00048 1 . OOO38 1.00031 1.00025 1.00020 1 . oooi4 1.00010 1.00007 1.00005 1.00004 119 Table 2. NUMERICAL VALUES OF THE THIRD VIRIAL COEFFICIENT C(T) AND THE CLUSTER INTEGRAL b„ . T(°K) c ( cm ^ b 2 (cm 3 ) \mole J 10 515 1,-68 20 286 -k55 30 226 -1.50 ko 195 -^.15 50 178 -7^27 75 162 -16,62 100 150 -27*58 273.18 107 -13O087 120 X q 1 1.5- 1.0- 0.5- H 1 1 1 I I 1 1 1 | I H 0.5 1.0 1 .1.5 T 2.0 i 1 I ♦—i — i > t + — i — + — i — i — i — i — *— i — i — t- 2.5 3.0 £i2 a ■0.5- FIGURE 17. X'^ AND SAMPLING ERROR FOR T=5°K. 121 \" 3 4 X fll 3.0- 2.5- 2.0 - - 1.5-- 1.0- 0.5- 0.5 ■ -1.0-- 1 T j t ,, , , ,, , xMj i -4 — »— ( 1 T I l-H 1 1— I 1 1 1 1 1 — I I I 1 1 1 1— ' 12 -♦— i — i i i i i * t i t | i i i — i i i i — i — »— i — i — x i i i — t— i — i i » -— — 1.0 " L I 1.5 2.0 t J 2.5 3.0 I FIGURE 18. X' 3 g AND SAMPLING ERROR FOR T=10°K. 122 A 3.0- 2.5- 2.0- 1.5- 1.0- 0.5 ■■ ■0.5- -1.0- I I I -» — i t i — i — i — t i i — i i i T > i i i i i — i t i i i T > I — I 1— I — I — I— i — I 1 — I — {if 0.5 1.0 l 1.5 2.0 £2.5 3.0 I I t r 12 FIGURE 19. X" 3 g 1 AND SAMPLING ERROR FOR T=20°K. 9l 3.0- 2.5- 123 2.0- I 1.5-- 1.0-- 0.5- 0.5 1.0 1.5 2.0 2.5 T 3.0 - 4.IT . ^ '12 H 1 1 1— I 1 I | t I I + 1 H — I — I — t— I — I- II*iI - • ^ 0.5- -1.0-- -1.5 + i i I I I V r3 FIGURE 20. X"°g 1 AND SAMPLING ERROR FOR T=30°K. 124 v-3. 2.5- 2.0- 1.5- 1.0 ■■ 0.5 ■• I I L t i i — i — I — i — i t t | i — i — i — J. | i i i t | i i i i | — i x j . . i . i T I ? ■ ■ mLL ■0.5-- ■1.0- 1.5- -2.0-- 0.5 1.0 1.5 2.0 I I 3.0 I I V FIGURE 21. X" 3 g 1 AND SAMPLING ERROR FOR T=40°K. A g, 2.5 2.0- 1.5-- 1.0- 125 0.5- 1 0.5 1.0 1.5 2.0 2.5 H 1 1 1— I 1 1 1 1 1 1 1— 1 1 1 1 ►— I 1 1 1 1 1 1 1 1 1 +— f — | ^ ±± 3.0 12 i* i ■0.5- 1.0 - -1.5- II 2.0 FIGURE 22. X" 3 g 1 AND SAMPLING ERROR FOR T=50°K. 126 A *\ 2.0 - 1.5- 1.0- 0.5- -0.5 -1.0- •1.5- 2.0- -2.5 ■■ 0.5 1.0 1.5 2.0 2.5 3.0 -t — i — i — 1—\ — i — i — i — i— H — i i i — i — (— t — i — i — i — t— t — i — i — l-H — •— • — •— i — I— I- riz 1*11 I >,' FIGURE 23. \- 3 g! AND SAMPLING ERROR FOR T = 75°K 127 X" 3 g 1 6.0- 5.0 •■ 4.0 - ■ 3.0- 2.0- 1.0- -1.0- 2.0- 3.0- I 4- r l2 H — I— I — I— | — I— I — I- — I — h— I — I — •— • — £ — I I I — I — I— i — I — I — I — I — I — I — I — I — f — I I » -i=. 0.5 1.0 1.5 2.0 2 j 5 I I 1 I 3.0 I I I I I U l .-3 FIGURE 24. X"^ AND SAMPLING ERROR FOR T=100°K 128 a X " 3 9l 5.0- 4.0-- 3.0- 2.0- 1.0 ■ ■ •1.0 ■■ 2.0 - -3.0- -4.0 •• H 1 1 1 h— ) 1 1 1 1 I I 1— I 1 I I t I 1 1 1 1 1 k— < 1 1 »— t 1 H 0.5 1.0 1.5 2.0 «1 I * *30 [12 a I L -3 FIGURE 25. \~*q ± AND SAMPLING ERROR FOR T= 273.18 °K 129 VITA Harry Frederick Jordan was born in Takoma Park, Maryland on March 6, 19^-0. He received the B.A. degree in physics and mathematics Cum Laude from Rice University in 1961 and the M.S. degree in physics from the University of Illinois in I963. During the period June 1962 through August 1966, he was a research assistant in the Department of Computer Science at the University of Illinois with the exception of the summer of l$)6k during which he served as a research assistant at the Los Alamos Scientific Laboratory of the University of California. He is currently an assistant professor of electrical engineering at the University of Colorado. Mr. Jordan is a member of Sigma Xi, the American Physical Society, and the Association for Computing Machinery. AUG I 6 1968