DENSITYFUNCTION FOR SELF-AVOIDIXC RANDOM WALKSON misconstrued, it is not intended that it should detract in any way from the undoubted usefulness of the “Monte Carlo” method as originally applied by Peri. This technique would seem to us to be one of the ways in which future experimental determination of the surface hydroxylation of oxide surfaces may be correlated with realistic theoretical models.
A
3953
LATTICE
Acknowledgments. We thank Miss J. Murray for the TiC14 experimental data. We thank S. R. C. for a personal grant to A. J. T. and Unilever Chemicals Development Centre for a personal grant to C. G. A. We also thank Unilever Chemicals Development Centre and S. R. C. for other generous financial support of part of this work.
The Density Function for End-to-End Lengths of Self-Avoiding Random Walks on a Lattice by Frederick T. Wall1” and Stuart G . Whittingtonlb Department of Chemistry, University of California, S a n Diego, L a Jolla, California 92057
Accurate coordinate density functions are calculated for the end-to-end separations of short self-avoiding random walks on the tetrahedral and face-centered cubic lattices. The results are extrapolated to describe long walks by assuming an approach to spherical symmetry with increase in number of steps. The form of the density functions so obtained can be represented by an equation used by Domb to describe his exact enumeration data on the cubic lattice. It is suggested that the essential form of the density function is independent of the lattice for long walks.
Introduction The number and configurational properties of selfavoiding random walks is of importance in the theory of polymer molecules,2 percolation problems3 and the susceptibility expansion of the Ising model.4 There is overwhelming evidence from exact enumeration resultss and Monte Carlo calculationse that the mean square end-to-end length, (rn2), of an n-step self-avoiding walk on a lattice has the functional form (r,2>
=
AnY
(1)
where 1.18 < y < 1.22 for walks in three dimensions. There is strong evidence from exact enumeration studies that y is exactly 6/6, a conclusion also derived theoretically by Edwards,’ using a self-consistent field approach for a continuum model. Rei& has recently criticized Edwards’ approach to the problem and, using a variation principle, has obtained a value of y = “3 in a continuum. It has since been shown, however, that R e i d solution of the resulting self-consistent field equation was not sufficiently a c ~ u r a t e . ~ Until recently, relatively little work had been carried out on the probability density function of end-to-end separations of self-avoiding walks. Domb’O and coworkers have fitted exact enumeration data on the cubic lattice to a distribution of the form
~ ( x =) A exp { - (x/u)’]
(2)
I n this expression, F ( x ) is the density function of the x coordinate of the end point of the walk, with v = 2.5 for three-dimensional systems. The arguments of Fisher’l are in agreement with this result, and Reiss* also predicts a non-Gaussian distribution, although Edwards” approach suggests that the distribution over r is a displaced Gaussian. It is the purpose of this paper to calculate F ( x ) exactly for short walks on the tetrahedral and facecentered cubic lattices and to use a procedure for longer walks that assumes an approach of the radial distribu(1) (a) American Chemica Society, 1155 Sixteenth St., N. W., Washington, D. C. 20036. (b) Unilever Research Laboratory, The Frythe, Welwyn, Hertfordshire, England. (2) F. T . Wall and L. A. Hiller, Ann. Rev. Phys. Chem., 5,267 (1954). (3) 8. R.Broadbent and J. M. Hammersley, Proc. Camb. Phil. Soc., 53,629(1957). (4) M.E.Fisher and M. F. Sykes, Phys. Rev., 114,45 (1959). (5) C.Domb, J . Chem. Phys., 38,2957 (1963). (6) F.T. Wall and J. J. Erpenbeck, ibid., 30,634 (1959). (7) S. F.Edwards, Proc. Phys. Soc., 85,613 (1965). (8) H.Reiss, J . Chem. Phys., 47,186 (1967). (9) I. McC. Torrens, ibid., 48,1488 (1968). (10) C. Domb, J. Gillis, and G. Wilmers, Proc. Phys. Soc., 85, 625 (1965). (11) M.E.Fisher, J . Chem. Phys., 44,616 (1966). Volume 75, Number 11 November 1969
3954
FREDERICK T. WALLAND STUARTG. WHITTINGTON
tion to spherical symmetry. This procedure has been carried out to about 90 steps, with results consistent with the form of F ( x ) suggested by the exact enumeration studies of Domb, et ~ 1 . ~ 0 The Generating Function of Short Walks. Suppose that N , (x,y,z) and N,"(x,y,z) are, respectively, the numbers of self-avoiding and unrestricted walks which, starting at the origin, have reached the point (x,y,z) after n steps. Let us define two generating functions G n ( ~ , P , rand ) G n " ( 4 , Y )as
Nn(x,YJz)ffxPV-8
G,(a,P,r) =
(3)
X,YIZ
and X,Y,Z
It is clear that e,, the total number of walks of n steps, is equal to Gn(l,l,l) and, as shown in Appendix 2, that
where (xn2)is the mean square x coordinate of the end point of the walk. Hence, for a lattice with cubic symmetry
Similar equations are applicable to the unrestricted walk. For an unrestricted walk on the tetrahedral lattice, the successive vectors representing steps must be chosen alternately from two groups, A and B, each with x, y, and x coordinates as follou1s12 X
Y
i i i
=
1
Gi
=
giGo = gi
Gz
=
gzG1 - 4Go
= (47
+ CY/PY + P / w + r/c@)
+ Pr/a + ar/P + 4 / r )
(7)
(8)
It is now clear that the generating function for the unrestricted tetrahedral walk assumes the form The Journal of Physical Chemistry
-4 gi2gz - 7gi
= gig2
=
gzGzp--I - 3G2p-2
(11)
=
4 ~ , -~ 3 ~ , -for ~ n
3
3
(12)
from which a solution is readily obtained, namely
I n this factor, the exponents of CY, P, and y, occurring in each term, correspond to the x, y, and z components of the corresponding tetrahedral vectors. For group B, a similar factor, g z , is defined g2 = (l/.Pr
giPg2p-I
G3 = giG2 - P G (10) The term -4Go in Gz removes the four walks that return to the origin at the second step, there being four configurations involving second steps that are the exact reverse of the first steps. The term -3G1 appearing in Ga is the generating function of a "tadpole" consisting of a stem of one step and a loop of two. The coefficient 3 reflects the fact that there can be precisely three such "tadpole" configurations associated with each stem; they are, of course, disallowed in self-avoiding walks. If the walks were only restricted by forbidding loops of two, the general recurrence relations for the generating function with p > 1 would be
c,
Each vector in group B is the negative of a vector in group A and each group constitutes a tetrahedral set. For the purpose of establishing a generating function like eq 3 or 4, we first define a factor g1 to symbolize the vectors in group A g1
= gig2
It follows immediately that
i i i 1 1 7
i i i
Gzo
Gzp+1 = glG;p - 3Gzp-1
1 1 1
i i i
= gi
Go
Gzp
i- i i
-
Gi"
1
Gzp" = gipgzp (9) We shall now assume that the generating functions for self-avoiding walks can be derived from the generating functions for unrestricted walks by introducing appropriate correction terms for loops of 2 , 4 , 6,. . . , 2 p , , , . steps. (On a more general lattice, loops of an odd number of steps might also be involved.) Thus, on the tetrahedral lattice, the first three generating functions become
x Y z
X
1 1 1
=
Gzp-l" =
B
A
Goo
c,
=
4
x
3-1
(13)
as required. It is important to notice that the factor 4 appearing in G2 and the factor 3 appearing in G3 and in succeeding G, reflects the difference in the numbers of ring closures attributable to loops of tq-o a t the origin and at points other than the origin. l 3 Presumably higher loops can be corrected for in a similar way, and with this end in mind, we shall employ a method which in our judgment is highly plausible (12) See, e.g., F. T . Wall, S. Windwer, and P. J . Gam, Methods Computat(onaZ Phys., 1, 220 (1963). (13) F. T. Wall, L. A. Hiller, and D. J. Wheeler, J . Chem. Phys., 22, 1036 (1956).
3955
DENSITYFUNCTION FOR SELF-AVOIDING RANDOM WALKSON A LATTICE although not rigorous. Pursuing this line of reasoning, we write that
- IC-1c QzkGzp-2a P-1
G2p = g2G2p-l
~ 2 p G o
(14)
and 4x3512
P
Gap+i = giG2p
-k-1
QzaGzp-zk+i
(15)
where uZp is the number of polygons of 2 p steps (or the number of returns to the origin), and 4 2 k corrects for ring closures to points other than the origin. (Actually, q2k equals the number of tadpole configurations of total length n and loop size 2L divided by Cn-2k.) A further element of complexity enters when star graphsI4 can occur on the lattice. These graphs represent walks containing more than one loop but no cut points. The treatment described above would correct for such configurations more than once, so it becomes necessary to make appropriate compensation. The net result is that q2k is not just a count of tadpoles, but is also weakly dependent upon n. Leaving aside the problem of the treatment of the star graphs, we turn first to the more immediate calculation of the uZp. By successive application of the recurrence relations, G2, can be expressed as a linear combination of the generating functions of unrestricted walks in the form P
Gzp =
i-0
b2aG2i0
- uzpGo0
(16)
in which the term uzpGO0simply represents the number of additional returns to the origin not accounted for in the expression beneath the summation sign. Hence uZp must be so chosen as to make the value of the coefficient of (cy3y) o in the total expression for G2, equal to zero. I n Appendix I1 it is shown that the value of the coefficient of (aby) O in the expression for G2pois given by
where the summation is taken over all nonnegative values of i, j , IC, and I satisfying the condition i j k 1 = p . The first few values of ~ ( 2 p are ) given in Table I.
+ +
+
3 x 2 ~ 6
Figure 1. Schematic representations of loop closures (a) and tadpole configurations (b) with approximate relative weighting factors.
Knowing the values of bzc and ~ ( 2 4we , can calculate uZpusing
c bzm(2i) P
u2p
=
(18)
i=o
Applying this approach to Gq, G5 and GOwe obtain
G4 = g2G3
- 3Gz
G5
- 3Ga
BIG4
u4 = 0
Go = g2G5
- 3G4 - 24Go
GT = giG6
- 3G5
u6 = 24
(19) That u4 = 0 simply means that there are no ways of closing a loop of four steps in the tetrahedral lattice, whereas u6 = 24 means that there are precisely 24 ways of returning to the origin in six steps. For seven steps it is necessary to determine a value of 46. I n forming a loop back to the origin, it is obvious that there are four ways to leave the origin and three ways to return; on the other hand, to form a loop with a stem, it is clear that there are three ways to leave a stem and, at most, two ways to return, without first forming another loop. (See Figure 1.) The ratio of these net counts is (4 X 3 ) / ( 3 X 2 ) or 2 ; hence, in the absence of star graphs, 42p/u2p = l/2. However, it is necessary to take into account the first 0 graphs,I4 0(5,1,1) of which there are u6/4 in each of the positions associated with G1. Hence
-
US/^ - Uo/4)G1
or GT = giGa
- 3G5 - 6Gi
(20)
At eight steps we must take account of polygons of eight steps and the “generalized tadpole”
4
Table I
I P
d2P)
0 1 2 3 4
1 4 28 256 2716 31504
5
which is a 0(5,1,1) with a unit-length stem. At ten steps a further complexity arises because of the star graph 8(6,2,2). Kevertheless, by proceeding as described above in deriving G7, we find that (14) M. F. Sykes, J. Math. Phys., 2,52 (1961).
Volume 73, Number 11 November 1969
FREDERICK T. WALLAND STUART G. WHITTINGTON
3956
- 3Gs - 8G2 - 48Go = giGs - 3G7 - 8G3 - 12Q = 92Gs - 3Gs - 8G4 - lOG2 - 480Go
Gs = g2G7 Gs Glo
F2d-2) (21)
The values of uZpand c, obtained from these generating functions agree with the counts of Sykes and Essam.16 This approach might be continued somewhat further, but it rapidly becomes exceedingly laborious without any assurance of a significant improvement in the ultimate extrapolation. Accordingly, we shall now turn our attention to a procedure leading to an extrapolation for large n. The Generating Function Projected o n One Dimension. To facilitate carrying our method further, we now find it expedient to project our walks onto one axis. It follows immediately from eq 3 that
G,(a,l,l) =
c N n ( z , y , z ) a Z= c F&)a"
Z,Y,Z
G f i o ( a , l , l )= ( 2 a
a statement which should be asymptotically correct for large values of p . Since F,(z) = F,(-x), it is also clear that for large p
FZP(0) = F2PM
I n principle, this should enable us to estimate values of uzP,as follows. From eq 26 we see that U2P =
-Fzp(O)
+ 2 p 2 P 4 - 1 ) + F2p-1(1)] .o-1
and
+
(23)
By subtracting eq 31b from 31a, and recognizing that F , ( l ) = F , ( - l ) , wefind U2p =
{F2p(2) - F2p(o) ]
+ 2{F2p-i(l) - F2p-i(3) ] +
P-1
C ~lzlc{F2p-21~(2)- F2,-2lc(o) ]
k.=l
c Nn0(X,Y,4 Zl,r
(24)
Evidently F," ( x ) satisfies the difference equation
F,"(z) = 2{Fn-i0(z
- 1) + Fn-io(z + l ) ] (25)
For a self-avoiding walk, we can write =
2{F2p-1(z
- 1 ) + F2p-l(X
+ 1)) -
c q2lc~2p-21c(z)- UzpFo(2)
k-1
(26)
and
- 1 ) + F2P(Z
+ 1)) -
Expressed in finite difference form for walks of 2p steps on the tetrahedral lattice, we see that The Journal of Physical Chemistry
(33)
lU2k
where { is assumed to be independent of IC. I n the absence of star graphs, other than ( n , l , l ) o , = uZk/3, so that we expect l < l / 3 . Since spherical symmetry is not realized for short walks, we shall provide for a gradual approach to spherical symmetry b y assuming that
F2P(O) = F2P(2)(1 where u2pand q2k are precisely the same as those appearing in eq 14 and 15, with qZkweakly dependent on the length of the walk. The Approach to Spherical Symmetry. It is well known that the distribution function of an unrestricted random walk approaches spherical symmetry for large n, and Domb, et aZ.,10 have presented numerical evidence that this is also true for a self-avoiding walk. I n Appendix I11it is shown that, for walks that are statistically spherically symmetric
(32)
If eq 30 were valid for all values of 2p, then eq 32 would simplify considerably. Actually the term in the summation of eq 32 cannot be neglected for values of k approaching p . Accordingly, we proceeded with numerical calculations as follows. For long walks, we assumed that13 q2x =
P-1
F2P+l(Z) = 2{F2,(x
(30)
(22)
so that F,"(z) =
(29)
2
where F ( x ) is the total number of walks allowed in three dimensions with a component of the end-to-end distance equal to z. From eq 7 and 8 we see that
F2,(x)
- 2F2,(0) 3. F2p(2) = 0
+V2p)
(34)
for values of p 3 5 , with X chosen to make the equation numerically correct for p = 5 . For smaller values of p , exact enumeration data have been used throughout. I n Table 11, values of ( T , ~ ) calculated by the method described above are compared with Monte Carlo data. Three different values of ( were used for the calculation. The agreement with Monte Carlo data is quite good when l is equal to l / 6 . The dependence of (rn2)on n is shown in Figure 2 and the limiting gradient indicates that y = 1.20 i 0.02. The dependence of c, on n suggests the connective constant, k, is given by ek = 2.876, close to the values derived from exact numeration'5 (2.878) and Monte CarIo studies5 (2.884). (15) M. F. Sykes and J. W . Essam,Physica, 29,378 (1963).
3957
DENSITY FUNCTION FOR SELF-AVOIDING RANDOM WALKSON A LATTICE Table 11: Table of (rn2)on Tetrahedral Latticea It
Monte Carlo
20 30 40
46.51 75.84 106.9 139.4 247.2
50
80 a
1
1
= 0.15
46.38 75.29 106.4 139.5 249.0
46.24 74.75 105.2 137.5 243.5
Calculations made for A
=
= l/6
r
= 0.18
46.49 75.39 106.6 139.9 250.1
0.1286.
256 X
64-
Figure 3. F,(z) us. x for n = 64 on tetrahedral lattice, (a) unrestricted walks; (b) self-avoiding walks calculated for r = 1/6.
43 16-
4-
I
I
4
16
64
n
Figure 2. Log-log graphs of on tetrahedral lattice; = l/6.
r
11s. n
for self-avoiding walks
The density function F,(x) is plotted in Figure 3 for n = 64 and the density function for the unrestricted
walk is shown in the same figure for comparison, both functions being normalized so that the total number of walks is equal in the two cases. Values of F,(s) calculated from eq 2 are shown in the same figure. The number of walks reaching a distance T from the origin after n steps P,(r) is given by (see Appendix 111)
(35) and the dependence of PB((r)on r is shown in Figure 4. Walks on the Face-Centered Cubic Lattice. The genP,(r) = 4nr2fn(r)
rB’,’(y)
erating function of an unrestricted walk on the facecentered cubic lattice is given by
+ l/a)(P + l/P> + (a + l / d ( Y + l / Y )
G n o ( a , P , ~= ) [(a
+ (P + 1/P>(r+ 1/7)1‘
(36) and the one-dimensional generating function of a selfavoiding walk can be written as G,(a,l,l)
=
dJ,-i(a,1,1)
-
n-1
k=2
+ +
where g = 4(a 1 l / a ) and the ak are again weak functions of n. Exact values of c, and u, are available16 up to n = 12 and ( T , ~ ) is known exactly4to n = 7.
Figure 4. P6((r)us. r for self-avoiding walks on tetrahedral lattice; = ’/a.
On this lattice we have made explicit use of the exact data to estimate ak and its dependence on n for n 5 12. We first allowed aR to assume two values at* and ak depending on whether the k loop appears for the second or subsequent time. By choosing the ak to force agreement with exact numeration values of c,, u,, and (r-2) we obtain the values of aRshown in Table 111. Since aR* N ak it seems reasonable to use only the initial values of ua and the limiting values of aR. We have adopted this procedure up to n = 12 using exact ~~
~
Table I11 k
Uk
ab*
2 3 4
12 48 264 1680
36 165.995 904.329
6
11
ak
11
36 166.147 906.496
(16) J. L. Martin, M. F. Sykes, and F. T.Hioe, J . Chem. Phys., 46, 3478 (1967).
Volume 73, Number 11 November 1969
3958
FREDERICK T. WALLAND STUARTG. WHITTINGTON 256r
X
I
4
n
16
Figure 6. F,(z) us. z for self-avoiding walks on face-centered cubic lattice; n = 64.
64
Figure 5. Log-log graph of , us. n for self-avoiding walks on face-centered cubic lattice.
values of c, and un and have fitted the values of as and ut to an equation of the form
ai,
=
+
udo!~ a d k
+4k2)
(38)
After 12 steps we proceed by forcing the approach to spherical symmetry by imposing the condition Fn(O)
=
F,(1)(1
+ h/n + Wnz)
(39)
where X1 and Xz are chosen to make the equation fit best at n = 10,11, and 12. When a0 = 0.2141, 011 = 0.6667, a2 = 7.0, XI = 0.076, and Xz = 0.396, the dependence of (rn2)on n is shown in Figure 5 and the coordinate and radial distribution functions are shown in Figures 6 and 7. With aproper choice of the parameter u, F ( x ) calculated from Domb's eq 2 can be essentially superimposed on the curve of Figure 6.
Discussion The agreement between values of (rn2>calculated by the present method with those obtained by Monte Carlo methods and by exact enumeration coupled with extrapolation techniques, as well as the agreement between the estimated values of the attrition coefficient, lends considerable support to the validity of the density function obtained. The shape of the density function is only weakly dependent on parameters such as I. The form of the density function is qualitatively in agreement with the forms proposed by Domb, et a1.,loand by Reisss and, with a proper choice of parameters, the agreement with Domb's results can be made quantitative except for the region in the immediate vicinity of the origin. This discrepancy was also noted by Domb.lo Since the form (eq 2) can be fitted to numerical results for the cubic latticelo as well as for the tetrahedral and face-centered cubic lattices studied here, it appears that the limiting form of the density function is probably independent of the particular lattice being studied. The Journal of Physical Chemistry
1
Figure 7. Pea ( r ) us. T for self-avoiding random walks on face-centered cubic lattice.
Appendix I The Relation between the Generating Function and the Mean Square x Coordinate of the End Point of an n-Step Walk. The generating function Gn(o!,P,y) is defined as Gn(a,P,y) =
C Nn(x,Y,z)azPVyZ X,Y,Z
Now
sothat,wheno!=p=y=l
= cn(x2)n
Appendix I1 The Number of Returns to the Origin at the dpth Step of a n Unrestricted Random Walk on the Tetrahedral Lattice. The generating function of the unrestricted walk is given by G.2po(a,P,~) =
(g1SdP
HEATSOF MIXINGAQUEOUS ELECTROLYTES
3959 Appendix 111 The Efect of Spherical Symmetry on the Behavior of the Coordinate Density Function Fn(x) around x = 0. Suppose that fn(r)is the’number of walks which are a t point r after n steps. Define the coordinate density function F,(x) by
Expanding by the multinomial theorem
+ + +
+ + +
where i j k 1 = m n q r = p . I n order that the exponents of a, p, and y all be zero, the following additional conditions must be satisfied
i+j-L-1-m-n+q+r=O
Let p2 = y2 symmetry, then
22
y2
=
+ z2.
Assuming spherical
i-j+k-1-m+n-q+r=O
i - j
-k+
1-m+n+q
- r
=O Differentiating with respect to z gives”
To satisfy these conditions it is necessary that i = m, j = n, L = q and 1 = r so that the coefficient of aopoyo in G2p0is =
c (-)p !
i,j,k,2
Fn’(x) F,”(x)
=
=
-2~fn(x)~
-2rlf,(x)
+ xfn’cx,]
For a self-avoiding walk f n ( 0 ) = 0 so that
i!j!k!Z!
Fn”(0) = 0
where the summation is taken over all nonnegative values of i,j , k , and 1 such that i j k 1 = p .
+ + +
(17) This result has been derived by Domb, et al.10
Heats of Mixing Aqueous Electrolytes. VII.
Calcium Chloride
and Barium Chloride with Some Alkali Metal Chlorides by R. H. Wood and M. Ghamkhar Department of Chemistry, University of Delaware, Newark, Delaware
(Received January 80,1969)
The heats of mixing aqueous solutions of calcium chloride and barium chloride with some alkali metal chlorides have been measured at constant molal ionic strength and 25’. The results are similar to those for magnesium chloride with alkali metal chlorides. At lower concentrations, RThobecomes more positive and the mixing with lithium chloride gives the lowest slope. The mixture with sodium chloride gives the largest value of RTh,,. If the contributions of oppositely charged ions are estimated and subtracted from the heat of mixing, the sign of the remainder correlates in most cases with the properties of the water structure around the two cations. The heats of mixing a t constant molal ionic strength, molar ionic strength and equivalents per kilogram of water are calculated and compared. The heats of mixing at constant molal and molar ionic strengths are not very different. The heats of mixing a t constant equivalents per kilogram are as small as the other heats of mixing and this may be the preferred concentration scale a t high salt concentrations.
to the magnitude Of the heat Of heat Of mixing is mixing two alkali metal chlorides. I n addition, the
(1) R.H.Wood, 3. D. Patton, and M. Ghamkhar, J . Phys. Chem., 73, 346 (1969). [Volume 73,Number 11 November 1969
