摘要:

Faraday Discuss. Chem. SOC., 1987, 83, 139-144 Asymptotic Structure of Diff usion-limited Aggregation Clusters in Two Dimensions Fereydoon Family”? and H. George E. Hentschel$ Department of Chemistry, Massachusetts Institute of Technology, Cambridge, Massachusetts 02139, U.S.A. Owing to the symmetry of the underlying lattice or anisotropy in the growth mechanism, large-scale diffusion-limited aggregation (DLA) clusters exhibit a crossover from isotropic to anisotropic growth. A scaling description of this crossover behaviour is presented and various exponents, including those describing the divergences of the lengths along and perpendicular to the direction of anisotropy, are defined. Scaling relations among the exponents and the fractal dimension are derived. An idealized model, in which the asymptotic shape of an m-fold symmetric DLA is approximated by a star- shaped object with m symmetric spokes, is introduced.Using the analogy between DLA and an equivalent electrostatic problem, the m-spoke model is solved exactly via conformal mapping and the exponents are calculated exactly. The critical mass, N,, above which DLA clusters exhibit an anisotropic structure is calculated and is found to grow as mp, where p is an exponent that depends on the fractal dimension. 1. Introduction Diff usion-limited aggregation and growth processes are topics of considerable current interest. Much of this interest has been stimulated by the diffusion-limited aggregation (DLA) model introduced by Witten and Sander.‘ DLA provides a realistic yet simple model for a wide variety of non-equilibrium phenomena, including electrodep~sition,~” pattern formation in Hele-Shaw cells,8-10 two-fluid flow in porous dielectric breakdown” and chemical dissolution of porous Much of the interest in DLA was stimulated by the assumption’ that DLA clusters are self-similar f r a ~ t a l s ’ ~ and exhibit scale-invariance and can be characterized by a unique fractal dimension D, which is a lattice-independent universal critical exponent. In fact, a number of attempts have been made to determine D analytically. 16-*’ For this reason considerable attention has been generated by the recent observations20-22 that asymptotically large DLA clusters grown in two dimensions on lattices with rn-fold or with growth anisotropy,21,’2 show similar m-fold anisotropy in their global structure.However, little progress has been made in a quantitative understanding of this unusual phenomenon, which is intimately related to the role of anisotropy in related problems of dendritic s o l i d i f i ~ a t i o n ~ ~ , ~ ~ and in viscous fingering.” In the first part of this paper we present a scaling description of the structure of DLA clusters in two dimensions, including the crossover from an initially isotropic growth to asymptotically anisotropic shape observed in large-scale computer simulations.21,22 In the asymptotic limit, DLA clusters with rn-fold anisotropy are star-shaped clusters with m branches emanating from the origin. We define two lengths, RII and R , , to measure the maximal length and the caliper width of the arms of a large DLA in the directions parallel and perpendicular to the axis of anisotropy, respectively.t Permanent address: Department of Physics, Emory University, Atlanta, Georgia 30322, U.S.A. $ Permanent address: Dowell-Schlumberger, Tulsa, Oklahoma 74134, U.S.A. 139140 Asymptotic Structure of DLA In the scaling description we define various exponents which may be used to describe the crossover from isotropic to anisotropic growth. Using scaling arguments we obtain explicit relations among these exponents and the fractal dimension. In the second part of the paper, the scaling exponents are calculated exactly for an idealized model of asymptotically large DLA clusters, providing confirmation of the scaling picture and exact results for the exponents.In this model the asymptotic shape of large DLA’s with rn-fold anisotropy is approximated by a star-shaped object with rn symmetric spokes of length RII . We exploit the electrostatic correspondence with DLA and use a conformal mapping to solve analytically the rn-spoke model and to calculate exactly the various exponents defined in the first part of the paper. The outline of the paper is as follows. In section 2 we present the scaling description and introduce the rn-dependent exponents needed to describe fully the structure of two-dimensional DLA clusters, including the crossover from initially isotropic cluster to asymptotically anisotropic shape. Section 3 is devoted to analytic studies of the rn-spoke model where the exponents defined in the first part are explicitly calulated and the critical size for observing anisotropic DLA clusters is determined. In section 4 we discuss the results and present the conclusions.2. Scaling of DLA Clusters In the presence of an rn-fold anisotropy, arising either from the symmetry of the underlying lattice or the anisotropy of the growth mechanism, an asymptotically large DLA cluster resembles a star-shaped object with rn spokes where each arm of the star is a ‘cigar’ of length RII and width R, .21722726 Initially, we expect Ril and R, to grow with the same power of the cluster mass N as (Rll/u)-- N ” and (R,/a)- N” ( 1 ) independent of rn. Here, a is the lattice constant. This isotropic growth continues until RII = R, = 6, where ( ( / a ) - rnP (2) and p is an exponent which remains to be calculated. We shall show 1ater.that this scaling form is supported by our model calculations; see eqn (14).For RII >> 6, on the other hand, and Eqn (1)-(3) can be combined into the scaling forms (RlI/“) - q I “ ” Q / 6 ) - N”J(rn-P”N”) ( R , / a ) - N”f,(N”a/() - N”f,(m-/”N”). (4a 1 (4b) and From eqn (3) and (4) we see that in the limit x + co the scaling functions JI(x) and fi(x) must vary as X ~ I I ’ ~ and x ~ ” ~ , respectively. This gives the following relations among the exponents andE Family and H. G. E. Hentschel 141 r n = 4 Fig. 1. The rn-spoke model for the asymptotic shape of DLA clusters witl, rn-fold anisotropy. In this example, rn = 4. In order to relate the exponents vll and Y, to the fractal dimension, we note that the star-shaped clusters are expected to have a DLA structure on scales r<< R, as N + a.This means that we can cover a 'cigar' with (RII/R,) 'blobs' of size R,, then N = (Rill R J R J a ) ''' or26 D = 1/ Y = 1 + ( 1 - ~ l l ) / Y,. (6) Eqn (4a) and (4b) and the scaling relations ( 5 ) and (6) can be used to provide a complete description of the scaling behaviour of DLA clusters with rn-fold anisotropy. 3. The rn-Spoke Model In the absence of a tractable approach for the DLA problem in general, we will introduce and study in this section an idealized model of asymptotically large DLA clusters and obtain the scaling results. The model consists of rn spokes of length RII, representing the rn arms of a star-shaped DLA with m-fold anisotropy, as shown in fig.1. The approach we use to solve the rn-spoke model is based on mapping the DLA problem onto an equivalent electrostatic p r ~ b l e m . ~ , ~ ' * ~ ' The random walker in DLA satisfies Laplace's equation V2C(r, t ) = 0, where C ( r , t ) is the concentration of the random walkers at the vicinity of point r at time t. Since the perimeter sites are perfect traps and the walkers cannot penetrate the cluster, C = 0 on the cluster. There is also an additional boundary condition that C( r -+ a) = Co, because random walkers are isotropically launched from infinity. The probability that a walker is absorbed at the perimeter site r is proportional to the flux VC. Thus, the DLA problem maps onto an electrostatic problem where C is the potential obtained by solving the Laplace equation with the boundary conditions C = 0 on the conductor and C = Co on a circle with large radius.In fact, the patterns found in the experiments on dielectric breakdown13 and electrodeposition' are very similar to the computer-generated DLA clusters. To solve Laplace's equation for the rn-spoke model we have to find some analytic function f ( z ) whose real part satisfies the required boundary conditions. Then C = f( z ) . Consider f( z ) = A sin-' [ ( z / RII)-"] + B. (7) Note that f ( z ) has the correct symmetry property, because it is invariant under the transformation z --+z exp (27~i/rn). Also, as r + a, f ( z ) -+ B and thus if we choose B = Co we have the correct behaviour as r + CO. To fix the boundary condition C = 0 on the spokes we need only consider one spoke on the x-axis.For O < x < RIl, 1 < (x/RII)-" < a. In this region (8) sin-' ( X / R ~ ~ ) - " = 7~/2--i In { ( X / R ~ ~ ) - ~ + J [ ( X / R ~ , ) - * ~ - 13).142 Asymptotic Structure of DLA Thus, if Ref(z) represents a solution to our problem, we require T A / ~ + B = 0. Then where Now, if the lattice cell size is a, then the probabilities of growth at the tips pt and sides p s of the spokes can be estimated as pt=m d C / a r l r = R , l u , while p s = rn(Rll/a)(l/r) a c / d 4 I r + R l l a . As pS+pt = 1 and ps>>pt as (Rll/a) -+ a, we have asymptotically or This is the result that Turkevich also derived the rn-dependence. In comparing eqn (14) with eqn pt - m-'/2 ( a / R,I) '/*. (R11/a) =: rn-1/3N2'3. (14) (13) and Scher2O would predict with 4 = 0, except we have As d ( R I I / a ) / d N = p , we find ( 3 a ) we find vll= 2/3 and pi1 = 1/3.(15) Note that eqn (15) agrees with our scaling assumption (3). The relation between the exponent Y, and the fractal dimension D is found from eqn (6) and (15) V, = [3(D - 1)I-l. This implies that the dependence of Y, on rn is through 0, and the m-dependence of Y, can be determined once the dependence of D on rn is known. We have obtained explicit results for the various exponents, defined in section 2, for the rn-spoke model. The exponents p, pll and p, are new and have not been discussed or measured before. However, several theoretical attempts 16-20 have been made in determining D, and its dependence on rn, and the exponents vII and v, have been measured n ~ m e r i c a l l y .~ ' , ~ ~ ' ~ " ~ ' The first of the attempts to determine D analytically was made by Muthukumar,16 and later by Tokuyama and Kawasaki17 and by Honda et aZ.,ls who used mean-field type arguments to determine D. They find that in d-dimensions D = ( d 2 + l ) / ( d + l ) , independent of rn. Using this expression in eqn (16) for d = 2 gives v, = 1/2 (independent of rn). (17) Based on a branched polymer topology, Hentschel" has used a Flory theory approach to determine D and finds that D = (5d2+ 8 ) / ( 5 d + 6) for a diffusion-limited aggregate. Substituting d = 2 in this expression and using eqn (16) we get v, = 4/9 = 0.444. . . (independent of rn). (18) In contrast to the mean-field type calculations which predict an rn-independent D, Turkevich and Scher2' find D = (271+4)/(7~+4) for an rn-fold symmetric lattice withF.Family and H. G. E. Hentschel 143 angle 4 = 27~/ rn between the arms. Substituting the Turkevich-Scher expression D = 2( rn + 1)/( rn + 2) for an rn-fold symmetric cluster in eqn (16) we find This gives v, = 5 / 9 , 1/2, 7/15 and 4/9 for rn = 3, 4, 5 and 6, respectively. In the case of uniaxial bias,” (rn = 2), D = 2. Substitution of this result in eqn (16) gives v, = 1/3, in agreement with the results of Ball et aL2’ Meakin22 has found v,=O.45 for an anisotropic 3-fold symmetric DLA cluster and Meakin and Family26 have reported v, = 1/2 for DLA on the square lattice. More recent, larger-scale simulations2’ of DLA clusters of up to 4 x lo6 particles on a square lattice give Y, = 0.48.These numerical results appear to be consistent within the error bars with the predictions of eqn (17)-( 19), but are not of sufficient accuracy to distinguish between the above predictions. In particular, whether v, does or does not depend on rn is still an open question and needs to be addressed in the simulations. The problem is that it appears to be difficult27 to obtain accurate estimates of v, for anisotropic DLA clusters owing to slow convergence to asymptotic behaviour. The exponents p, and pi1 can also be related to the fratal dimension D. From eqn ( 5 ) , (15) and (16) we find v, = (rn +2)/3rn. (19) (20) p,=(20-1)/[3(0-1)(3-2D)]. (21) -1 p = 2 D - 3 and finally p, is found from eqn ( 5 ) , (15) and (16) As mentioned above, there have been no previous attempts to determine these exponents and we hope that in the future they are measured and can be compared with the predictions (20) and (21).We have now obtained expressions for all the exponents in the case of rn-fold symmetry. Let us use these exponents to answer the question: when is the anisotropy first observed and how does it depend on rn? We would expect the anisotropy to be observable when the angle subtended by a ‘cigar’, $ =: R I I / R , , becomes smaller than the intrinsic angle of symmetry, $ = 2 7 ~ / r n . Using eqn (3) we would expect to see anisotropy for N > N, where with N,- rnp /3 = ( 4 0 - 3)/(2D - 3). Note that N, grows as a power law and for rn = 6 reaches values of the order of 10’. This implies that at the present level of computational capabilities lattice anisotropy cannot be observed for rn > 4, in agreement with the available simulations where no anisotropy has been found for DLA clusters grown on hexagonal lattices with 6-fold symmetry.22 4.Conclusion We have studied analytically the scaling properties of two-dimensional DLA clusters with symmetric rn-fold anisotropy and answered several questions regarding their asymptotic structure. If the clusters are asymptotically unstable with respect to anisotropy induced by the growth mechanism or the symmetry of the underlying lattice as suggested by ~ i m u l a t i o n s , ~ ~ * ~ ~ ~ ~ ~ ” ~ then two lengths, Ril and R,, measuring lengths along and perpendicular to the direction of anisotropy, are needed to describe their scaling behaviour.26 We have used conformal mapping to calculate explicitly, for the rn-spoke model, the dependence of the exponents zq and v, on the degree of symmetry rn through the fractal dimension D.If the asymptotic shape of DLA clusters is analogous to the rn-spoke model, then vII = 2/3 independent of nz, but v, = [3(D - 1 ) I - l and may144 Asymptotic Structure of DLA depend on rn through D. We also calculated the critical size N,, above which DLA clusters exhibit fully anisotropic structure. We find that N , grows as mp, where p depends on D. This result points to two possibilities. First, N , grows so quickly with rn that for m > 4 present simulations will not be able to see the true asymptotic behaviour in DLA. The second possible conclusion is that as m increases the screening of the interior increase and the needle-like picture presented by the simple m-spoke models breaks down.In this case, the tip angle a = .n/2 - 4/2 never approaches its asymptotic limit and vll does not reach its maximum value of 2/3. For small a [ a = RL/(R,l/2)], from the formula of Turkevich and Scher,*' (24) Thus, if the rn-spoke model represents the true asymptotic shape of DLA, then the scaling description and the exact results obtained here give a complete picture of DLA clusters in two dimensions. On the other hand, if the asymptotic shape of DLA is more intricate, then the rn-spoke model will provide a simple and tractable phenomenological tool for further investigations of the asymptotic structure of large DLA clusters in two dimensions.We would like to thank John Deutch for helpful discussions. This work was partially supported by the Office of Naval Research and the National Science Foundation. Acknowledgment is also made by F.F. to the donors of the Petroleum Research Fund, administered by the ACS, for partial support of this research. References 1 Kinetics of Aggregation and Gelation, ed. F. Family and D. P. Landau (North-Holland, Amsterdam, 1984). 2 On Growth and Form: Fractal and Non-Fractal Patterns in Physics, ed. H. E. Stanley and N. Ostrowsky (Martinus-Nijhoff, Dordrecht, 1985). 3 H. J. Herrmann, Phys. Rep., 1986, 136, 153. 4 F. Family, P. Meakin and T. Vicsek, Rec. Mod. Phys., in press. 5 T. A. Witten Jr and L. M. Sander, Phys. Rev. Lett., 1981, 47, 1400.6 R. M. Brady and R. C. Ball, Nature (London), 1984, 309, 225. 7 M. Matsushita, M. Sano, Y. Hayakawa, H. Honjo and Y. Sawada, Php. Rev. Lett., 1984, 53, 286. 8 J. Nittmann, G. Daccord and H. E. Stanley, Nature (London), 1985, 314, 141. 9 G. Daccord, J. Nittmann and H. E. Stanley, Phys. Rev. Lett., 1984, 56, 336. 10 H. Van Damme, F. Obercht, P. Levitz, L. Gatineau and C. Laroch, Nature (London), 1986, 320, 731. 11 L. Patterson, Phys. Rev. Lett., 1984, 52, 1621. 12 K. J. Maloy, J. Feder and T. Jossang, Phys. Rev. Lett., 1985, 55, 2688. 13 L. Niemeyer, L. Pietronero and N. J. Wiesmann, fhys. Rev. Lett., 1984, 52, 1033. 14 G. Daccord and R. Lenormand, Nature (London), in press. 15 B. B. Mandelbrot, The Fractal Geometry of Nature (Freeman, San Francisco, 1982). 16 M. Muthukumar, Phys. Rev. Lett., 1983, 50, 839. 17 M. Tokuyama and K. Kawasaki, Phys. Lett. A, 1984, 100, 337. 18 K. Honda, H. Toyoki and M. Matsushita, J. Phys. SOC. Jpn, 1986, 55, 707. 19 H. G. E. Hentschel, Phys. Rev. Lett., 1984, 52, 212. 20 L. Turkevich and H. Scher, Phys. Rev. Lett., 1985, 55, 1026. 21 R. C. Ball, R. M. Brady, G. Rossi and B. R. Thompson, fhjls. Rev. Lett., 1985, 55, 1406. 22 P. Meakin, Phys. Rev. A, 1986, 33, 3371. 23 E. Ben-Jacob, N. D. Goldenfeld, J. S. Langer and G. Schon, Phys. Rev. Lett., 1983, 51, 1930. 24 R. Brower, D. Kessler, J. Koplik and H. Levine, Phys. Rev. Lett., 1983, 51, 1111. 25 E. Ben-Jacob, R. Godbey, N. D. Goldenfeld, J. Koplik, H. Levine, T. Mueller and L. M. Sander, Phys. 26 P. Meakin and F. Family, Phys. Rev. A (Rapid Commun.), 1986, 34, 2558. 27 P. Meakin, R. C. Ball, P. Ramanlal and L. M. Sander, preprint, 1986. Rev. Lett., 1985, 55, 1315. Received 8th December. 1986

ISSN:0301-7249    

DOI:10.1039/DC9878300139

出版商:RSC    

年代:1987

数据来源: RSC

摘要:

Faraday Discuss. Chem. SOC., 1987, 83, 145-152 Hydrodynamic Properties of Fractal Flocs Pierre M. Adler Laboratoire d'Ae'rothermique du C. N. R.S. 4ter Route des Gardes, F92190-Meudon7 France The flow field around two models of two-dimensional fractal flocs is numeri- cally derived in the zero-Reynolds-number limit and macroscopic quantities such as the drag force exerted by the flow on the flocs are determined. In the range of parameters studied and within the precision of the numerical data, the corresponding dimensionless quantities are mostly functions of the gyration radius, the role of the fractal dimension, if any, is very weak and the average drag is very close to that relative to circular cylinders. The fundamental problem associated with coagulation and the characterization of the resulting aggregates have been recently totally renewed by the concept of fractals.' Witten and Sander' studied the formation of aggregates by computer simulation and interpreted the results in terms of fractal dimension, being the first to do so.Later Meakin3 and Kolb et aZ.4 devised the cluster and cluster aggregation model which was found to be in good agreement with experimental data,' thereby confirming the fractality of flocs. Some theoretical questions of importance such as the asymptotic behaviour of the flocs and the influence of the underlying lattice are not totally clear yet and are beyond the limits of this work. The major purpose of the present paper is to determine some hydrodynamic properties of fractal flocs. To the best of our knowledge, we are only aware of two papers637 which address this matter.The translational friction coefficient of three-dimensional Witten- Sander and cluster-cluster aggregates was determined using a far-field approximation for the hydrodynamic interaction between the particles. These results will be detailed below and it suffices here to note that the hydrodynamic radius of the first class of aggregates was found to be significantly different from the gyration radius; this conclusion was not expected on a physical basis. It was thus found useful to determine some hydrodynamic properties of fractal flocs by the direct resolution of the Stokes equations of motion. One interesting issue in this prospect is the potential influence of the fractal dimension on the results; in this respect, the Witten-Sander and cluster-cluster flocs are good candidates since their fractal dimensions are quite different.An important restriction on our numerical results is that they were only obtained in two dimensions; this is due to the large amount of computer time which would be necessary in three dimensions. In the next section, the construction of such flocs as well as their major properties are briefly recalled. The Stokes equations which govern the flow field are presented and the numerical scheme which is very classical indeed is sketched; the flow is computed around a spatially periodic pattern of identical flocs located at the centre of a channel. The numerical results are then exposed and discussed in the last section. The statistical procedure is presented and the role of physically irrelevant parameters such as the size of the numerical mesh etc.is detailed. Two macroscopic quantities are analysed, namely the seepage velocity through the channel and the hydrodynamic force exerted on the flocs by the fluid. First, the dimensionless results are found to depend only on the ratio of the gyration radius to the channel halfwidth for both kinds of flocs; this property is due to the self-similarity of the aggregates. Secondly, the influence of 145146 Hydrodynamic Properties of Fractal Flocs 100 N 'O 1 / I I I I 0.5 1 2 5 10 Rgl 'A Fig. 1. The size N of the aggregates as a function of the gyration radius R,: 1, Witten-Sander; 2, hierarchical model. The broken lines represent the asymptotic behaviour of these two models: 1, D f = 1.70; 2, Df= 1.5.the fractal dimension on this function is very weak if any; this surprising conclusion is discussed at length and compared with the results obtained in three dimensions. It is checked that this unexpected feature cannot be attributed to similarities of the external perimeters of the two models which were used in the simulation; they are compared and found quite different. General Two kinds of two-dimensional flocs are constructed and used in this paper, namely the Witten-Sander model and the hierarchical model. These two models have been selected for the following reasons. First, the asymptotic regime is very quickly reached and they have been more widely studied than more general models such as that of ref.(8); for practical purposes, this feature will turn out to be crucially important here, in relation to numerical aspects of the computations. Second, the fractal dimensions are significantly different, since they are equal to 1.7 and 1.5 for the Witten-Sander and the cluster-cluster models, respectively. The Witten-Sander aggregate is grown in the usual way, but for a detail which accelerates the convergence to the asymptotic regime.' The random particle does not stick to the cluster when it is on a site adjacent to occupied sites; it is moved until it is on a site already occupied by the cluster. It is then assumed to stick to the cluster, when it was in the previous position. A floc is characterized by its gyration radius, which may be defined as where N is the total number of particles of the cluster, n and rn denote the number of any two particles and r, is the position vector of a particle.The fractal dimension Df may be introduced by the relation Some data are shown in fig. 1, where it is seen that the asymptotic regime is very N=R,D'. (2) quickly reached. The value of D,- is very close to the accepted value of 1.70.''P. M. Adler 147 Fig. 2. Two random aggregates of N = 64 particles on a square lattice: ( a ) Witten-Sander; ( b ) hierarchical. I I I I / / I I 2h j % i % j flow I 1 I I I I * T' 2h 4 - x Fig. 3. Spatially periodic array of fractal flocs in a two-dimensional channel of width 2h. The hierarchical model with linear trajectories is applied without any modification.10311 This particular model was chosen because it is less time consuming than the complete cluster-cluster model and the hierarchical model with brownian trajectories.The influence on the fractal dimension is very small; note that two clusters stick when they occupy adjacent sites. Results are shown in fig. 2 and the fractal dimension is very close to 1.50.'' Finally, for computational purposes, aggregates are assumed to be composed of square particles of side a, instead of circular ones; this feature should not influence the asymptotic character of the results. Typical aggregates are displayed in fig. 2; the different values of Df are clearly visible. Since we are dealing with a two-dimensional system, the Stokes paradox around cylindrical bodies shows up at zero Reynolds number for translational problems whatever their particular shape. Two techniques may be employed in order to avoid this difficulty; either the inertial terms may be inserted in the equations, or the flow may be analysed in a bounded region.The latter technique was chosen here. The following situation was thus analysed. A spatially periodic pattern of identical flocs was created. The centre of gravity of these flocs is located at the centre of the channel and the distance between two successive flocs is equal to the width 2h of the channel (see fig. 3). A flow was forced in this channel by means of an external pressure gradient. The low Reynolds number flow of a Newtonian fluid is governed by the usual Stokes equations: vp = pv2v ( 3 4 vv=o (36) where u, p and p are the velocity, pressure and viscosity of the fluid, respectively.148 Hydrodynamic Properties of Fractal Flocs In general, o satisfies the conditions u = 0 at the wall u = 0 at the surface of the floc u is spatially periodic. (4) Finally the pressure drop through one cell is imposed, or equivalently the macroscopic The seepage velocity u, may be expressed as pressure gradient is a given constant vector which is here parallel to the x-axis.f + h since it has only one component on the x-axis. advantage is taken of the fact that the stress u is divergence-free. Hence, The viscous force F exerted by the fluid on each floc may be easily calculated when F = ds u ( 4 ) where d 0 is the external surface of the unit cell, and ds the outer surface element. The determination of F necessitates the resolution of the two-dimensional system (3), for which several standard methods exist.'2 In order to cope with the continuity equation, the so-called artificial compressibility method was used with a staggered marker-and-cell (MAC) mesh.In essence, the problem is replaced by an unsteady compressible one which is assumed to converge towards the steady incompressible situation of interest. An alternating-direction-implicit (ADI) scheme was used in order to meet the stability conditions with an acceptable computation time. A decisive advantage of this scheme is that it may be readily generalized to three dimensions. No further detail will be given on the numerics, which are classical. The mesh spacing A is the same in both the x and y directions (parallel and perpendicular to the channel axis, respectively) and it is a fraction l / n of the size of the particles: A = a l n .(7) Results and Discussion The average seepage velocity (ij) and the average drag force ( F ) can only be determined by Monte-Carlo calculations. A number N of aggregates of a given mass N are generated; iji and Fi are determined for each of them; the average values may then be expressed as In view of the symmetry of the aggregates, the average drag is parallel to the x-axis and in the following its component (F,) will be denoted simply by ( F ) . These two quantities can now be rendered dimensionless. When the channel is empty and is submitted to a pressure gradient dpldx, the seepage velocity and the velocity at the centre are, respectively, Hence, it is natural to compare (fix) to Go, and ( F ) to puma,.Let us define u" and F" byP. M. Adler 149 Q5 V* 01 0 V* Q5 0 0.5 1 0 0.5 1 Rgl h Rgl h Fig. 4. The dimensionless seepage velocity v" as a function of the dimensionless gyration radius R,/h. ( a ) Witten-Sander model, ( b ) hierarchical model. Conditions are: 2 h / a = 9 ( +), 14 ( x ), 19 (e), 24 (O), 34 (0). The bars indicate the range of spreading of the individual values of u * . The number N of particles in the aggregates ranges from 4 to 64 whenever possible. As is usual in low Reynolds number flow,13 these quantities only depend in principle of the geometry of the system, which may be characterized by the number N of particles in the aggregate, by the ratio R,/h and by the fractal dimension of the aggregate.One may think that this latter number is not sufficient to describe fully the aggregate and so we leave the list open. Thus, v* = v*(N, R,/h, Df, . . .) F* = F*( N, Rg/ h, Df, . . .) However, these relations may be further simplified when the self-similarity of the flocs is used; two aggregates, which contain different numbers of particles, are identical within a change of scale. Hence N should not be included any more in the above list of geometrical parameters. Thus, v* = v*( Rg/ la, Df, . . .) F* = F*( Rg/ h, Df, . . .) The numerical results are also functions of the two parameters n and N, which are irrelevant in the sense that they must be chosen in such a way that they do not influence the final results on v* and F*. It is obvious that the larger these parameters, the better the accuracy; however, a compromise must be found with the total length of the computations which increases very rapidly with these parameters and becomes prohibi- tive.N was chosen equal to 10, which is far from being sufficiently large as it will be seen below; with such a choice, the error corresponding to the largest possible mesh, i.e. n = 2, was usually found to be negligible when compared to the statistical error. Let us look first at the results relative to the dimensionless seepage velocity v*. They are given in fig. 4. Examples of statistical errors are displayed by vertical bars; the ends of the bars correspond to the extreme values of v* which have been obtained for 10 aggregates of a given size.Of course, fluctuations are larger when the aggregate is large compared to the channel width. Calculations were performed for a series of ratios 2h/ a ;150 Hydrodynamic Properties of Fractal Flocs Fig. 5. The dimensionless drag F" as a function of the dimensionless gyration radius Rg/ h. Data are for: +, Witten-Sander; x , hierarchical. The number N of particles in the aggregate ranges from 4 to 64 whenever possible. Three channel widths were investigated: 2 h / a =9, 19 and 29. The bars indicate the spreading of the individual values of F". it is seen that the scaling relation (1 1) is extremely well verified, since the dimensionless seepage velocity only depends upon the ratio R,/h. This is even more remarkable when it is remembered that the smallest flocs are of size N = 4.For instance, for Witten-Sander aggregates v" = 0.200 for R,/h = 0.23,2h/a = 9, N = 4; this is very close to an apparently different situation: ( 2 h / a = 34, N = 32) for which R,/h = 0.19 and v* = 0.027. Another result was totally unexpected: it is clearly seen in fig. 4 that the curves for the two kinds of flocs are almost identical and may be almost exactly superimposed. In other words, there is almost no influence of the fractal dimension Df and all the necessary geometric information seems to be contained in the gyration radius R,, though the structures of the two flocs are very different, as was visually demonstrated in fig. 2. The data on the hydrodynamic drag F" may be presented in the same way. In fig. 5, F" is drawn as a function of the dimensionless radius of gyration for various channel widths, various values of N and for the two kinds of flocs under study.As can be seen, F" adheres well to the scaling relation (11) and hardly seems to depend on Df, at least when the statistical precision of the results is taken into account. This precision is very poor and a sample of N = 10 aggregates as used here is certainly not sufficient. Thus, it may be found remarkable that the differences between the two sets of data are well below the fluctuations; moreover, there is no systematic trend in these differences. The drag force may be plotted in another way. Analytical calculations [see ref. (13) for a review] show that F" has a logarithmic singularity in relation to the StokesP. M. Adler 151 0 1 2 Fig.6. The inverse of the dimensionless drag as a function of the natural logarithm of the dimensionless gyration radius R,/h. Notations are the same as in fig. 5. The solid line corresponds to the analytical results of H a r r i ~ o n ' ~ for circular cylinders. paradox. Hence, it should be of the form ' = 1 [ l o g A - k o + k , ( ~ ) 2 + - F* 47~ R, - * ] / [ l - k 3 ( + ) 2 + - -.I where ko, k, and k, are constants of order 1, which were determined by H a r r i s ~ n ' ~ for circular cylinders. Note that the radius of a cylinder is equal to &R,. The constants are then equal to ko= 1.3147, k, =3.750, k2= 1. (13) Numerical data are plotted in fig. 6 in this representation. Note that this seems to be the adequate framework for the analysis of the results and that, surprisingly, the drag on circular cylinders is very close to the drag on the two kinds of fractal cylinders which were studied here.These results on the relatively weak (if any) influence of Df seem to contradict those obtained previously in three dimensions.677 This disagreement should be weakened by the fact that the flocs analysed in ref. (6) and (7) were much larger than the ones here. However, it seems that the drag force on self-similar flocs in three dimensions should scale as F = Upk( Df, . . .)Rg (14) where U is the translational velocity and k(Df, . . .) is a constant which could depend on D,- and on other geometrical quantities, but not on R,. This relation is the counterpart of eqn (1 1) in three dimensions. It should also be noticed that eqn (1 1) and (13) enable us to pass smoothly from fractal objects for which Df < 2 to non-fractal ones for which Df= 2, such as circular cylinders.This disagreement between numerical data in 2 and 3 dimensions is still worse when it is noticed that the confined configuration which is used here should enhance the differences, since the fluid is forced to pass through certain narrow areas which would otherwise be ignored. The reduction of the dimensionality also acts in the same sense. A possible way to understand the lack of influence of Df was to compare quantitatively the external configurations of the two models of flocs. This was done in two different ways. First, the number N, of external sites was calculated on the samples which were used for the hydrodynamic calculations.Secondly, in reference to the classical Carman152 I Hydrodynamic Properties of Fractal Flocs 10 - 5, ? 4* - 1 , Fig. 7. External characteristics of fractal flocs as functions of the gyration radius. 1, Witten-Sander model; 2, hierarchical model. ( a ) External sites N , . ( b ) Dimensionless wetted perimeter PJa. These results are relative to the statistical samples analysed in fig. 4-6. equation,” the perimeter P, of the aggregates wetted by the external fluid was deter- mined. These two quantities are represented in fig. 7. They are seen to be significantly different and they cannot be used to explain the lack of influence of Df on the seepage velocity and the drag force. Finally, it would certainly be very desirable and interesting to apply the present numerical scheme to other situations, such as three dimensions or larger and more numerous aggregates in order to improve the statistics; other values of Df should also be investigated, especially close to I . A priori, there is no fundamental problem associated with these extensions; the only one that has severely slowed down our efforts in this field is the cost of computer time. References 1 B. B. Mandelbrot, The Fractal Geometry of Nature (Freeman, San Francisco, 1982). 2 T. A. Witten Jr and L. M. Sander, Phys. Rev. Letf., 1981, 47, 1400. 3 P. Meakin, Phys. Rev. Lett., 1983, 51, 1119. 4 M. Kolb, R. Botet and R. Jullien, Phys. Rev. Lett., 1983, 51, 1123. 5 D. A. Weitz and M. Oliveria, Phys. Rev. Lett., 1984, 52, 1433. 6 Z. Y . Chen, J. M. Deutch and P. Meakin, J. Chem. Phys., 1984, 80, 2982. 7 P. Meakin, 2. Y . Chen and J. M. Deutch, J. Chem. Phys., 1985, 82, 3784. 8 P. Meakin, Phys. Rev., 1983, B28, 6718. 9 R. Jullien, personal communication. 10 R. Jullien, in preparation. 11 D. Sutherland, Nature (London), 1970, 226, 1242. 12 R. Peyret and T. D. Taylor, Computational Methodsfor Fluid Now (Springer Series in Computational Physics, Berlin, 1985). 13 J. Happel and H. Brenner, Low Reynolds Number Hydrodynamics (Prentice-Hall, Englewood Cliffs, 1965). 14 W. J. Harrison, Trans. Cambridge Philos. SOC., 1924, 23, 71. Received 8th December, 1986

ISSN:0301-7249    

DOI:10.1039/DC9878300145

出版商:RSC    

年代:1987

数据来源: RSC

摘要:

Faraday Discuss. Chem. SOC., 1987,83, 153-165 Properties of Fractal Colloid Aggregates H. M. Lindsay, M. Y. LinJ D. A. Weitz,* P. Sheng and Z. Chent Exxon Research and Engineering, Rt 22 E, Annandale, NJ 08801, U.S.A. R. Klein Department of Physics, University of Konstanz, Konstanz, Federal Republic of Germany P. Meakin$ Department of Physics, Boston University, Boston, M A 02140, U.S.A. We study the influence of the fractal structure on the physical properties of colloidal gold aggregates. The optical properties are strongly influenced by both the long- and the short-range correlations of the aggregate structure, as well as the electronic plasma resonance of the gold particles. The absorp- tion of the aggregates is dominated by the electromagnetic interaction between nearest-neighbour particles.The angular dependence of the polar- ized scattering reflects the long-range fractal correlations of the particles in the clusters, while that of the depolarized scattering reflects the short-range correlations of the local field corrections. The fractal structure also affects dynamic light scattering, so that rotational diffusion effects play an important role in the decay of the autocorrelation function. Finally, we show that the very tenuous nature of the fractal structure makes the aggregates quite susceptible to deformation under shear stress. There have been considerable advances recently in our understanding of the process of colloid aggregation.' Central to all this work is the recognition that the highly disordered structure of colloidal aggregates can be characterized in quantitative terms.It has been shown both experimentally' and through computer simulation374 that the structure of colloidal aggregates exhibits scale invariance and can be quantitatively described as a fractal. This development has been accompanied by significant advances in our under- standing of the dynamics of the aggregation process and the causal relationship with the structure.536 Thus a new and rather elegant picture of colloid aggregation is emerging. Armed with the new knowledge of the process of colloid aggregation, an important set of questions can now be investigated. Since we can quantitatively characterize the structure of colloidal aggregates, we can investigate their physical properties in much greater detail.Of particular interest is the determination of the physical properties unique to fractal objects. Knowledge of these is crucial not only as a matter of fundamental scientific interest, but also because of the technological importance of colloidal flocs. In this paper we examine several different aspects of the physical properties of colloidal aggregates that are uniquely affected by the fact that their structure is fractal. These properties include a detailed discussion of the optical properties as well as the results of a preliminary investigation of the mechanical stability of the fractal structure. We study the properties of colloidal gold aggregates. Colloidal gold has historically been a widely studied system and is a prototypic lyophobic, metal colloid.It has the advantage that the aggregation process has been studied in considerable detail and is t Also at the Department of Physics, City College of CUNY, New York, NY 10031, U.S.A. 3 Permanent address: E. I. DuPont Co., Wilmington, DE, U.S.A. 153154 Properties of Fractal Colloid Aggregates relatively well understood.’ Furthermore, the aggregation of gold colloids can be easily controlled, and the fractal dimension of the resultant clusters can be reproducibly varied by changing the conditions of the aggregation kinetics5 This gives us additional flexibility in determining the effects of the fractal structure on the physical properties of the aggregates. Colloidal gold also possesses an additional very interesting property, in common with several other, widely studied metal colloid systems.This is the electronic plasma resonance that gives the colloid its unique colours: red when unaggregated and blue or grey when aggregated.’ In this paper, we also consider the effect of this resonance, in combination with the fractal structure, on the optical properties of the aggregates. The recent advances in our understanding of the kinetic aggregation process of colloidal gold have been discussed elsewhere,’ and we present here only a very brief summary. The aggregation of colloidal gold can be divided into two distinct regimes5 In the first, the aggregation is very rapid, and the rate is limited solely by the time taken for cluster collisions to occur via Brownian diffusion. This is diffusion-limited cluster aggregati~n’?~ (DLCA).In the second regime the rate is limited by the low probability of actually forming a bond upon collision of two clusters. This is reaction-limited cluster aggregation*-’* (RLCA). These regimes are achieved experimentally by adjusting the charge on the colloid surface, and hence the height of the potential barrier caused by the repulsive Coulomb interaction between approaching colloid particles.6 When the surface charge is completely removed, so there is no repulsive barrier, DLCA results. By contrast, if the surface charge is reduced so the barrier is still several k,T, RLCA results. The clusters in both regimes are fractal, but their fractal dimensions, d f , are diff erent.5 For DLCA, dr= 1.8, while for RLCA, df = 2.1. In addition, both the cluster mass distribu- tions, and their time dependences, are different for each regime.” For diffusion-limited aggregation, the cluster mass distribution exhibits a broad peak, and the average cluster mass grows roughly linearly with time.I2 By contrast, for reaction-limited aggregation, the cluster mass distribution is power-law, CM cc A T T , with T == 1.5, up to some cutoff mass.This cutoff mass (which characterizes the size of the clusters in the distribution) increases exponentially with time. Thus we can use this knowledge of the aggregation of colloidal gold to obtain clusters of a known structure, as characterized by their fractal dimension, and distribution. Detailed knowledge of the aggregation process also provides us with another valuable method for studying the physical properties of colloidal aggregates. The physics of the aggregation can be formulated as a set of rules and then simulated on a computer.This allows one to generate clusters that have the same fractal dimension and cluster mass distribution as that produced by experiment. Moreover, from the simulations we know the exact coordinates of each particle in the clusters, and these can be used for numerical calculations to determine a variety of physical properties of the clusters. Here we use the results of an off-lattice diffusion-limited cluster-cluster aggregation model” for numerical computations. In this paper we investigate a variety of physical properties of aggregates. Since optical scattering techniques are one of the simplest and most informative methods for studying colloids and aggregation in general, an understanding of the optical properties is crucial.We therefore investigate the optical properties of fractal gold aggregates in detail. We consider the influence of both the structure of the aggregates and the electronic plasma resonance on both the absorption and scattering from the clusters. We also consider the effects of the fractal structure of the aggregates on dynamic light scattering. The size of the aggregates is typically much larger than the inverse of the scattering vector used. We show that this results in a significant contribution from rotational diffusion to the decay of the measured autocorrelation of the scattered light intensity, and we develop a general formalism to properly account for this.Finally, we consider the mechanical properties of the aggregates. A fractal aggregate possesses uniqueH. M. Lindsay et al. 155 wavelength/nm Fig. 1. The absorption (solid line) and integrated scattering (dashed line) from unaggregated gold colloids, ( a ) and ( b ); and from DLCA colloids, 25 min after the aggregation was initiated, ( c ) and I d ) . correlations, which have significant consequences on its mechanical and structural stability. An understanding of these effects is essential since most colloidal aggregates of technological importance are rarely prepared under the controlled aggregation condi- tions used in the laboratory. We discuss preliminary measurements of the mechanical stability of colloidal aggregates and the effects of shear forces on their structure.The remaining paper is divided into four sections. The next section discusses the optical properties of fractal aggregates and the effect of the electronic plasma resonance on these properties. The following section discusses quasielastic scattering from fractal aggregates and the contribution of rotational diffusion to the measured decay time. Then we discuss the changes in the structure of the aggregates induced by shear forces. We present a brief summary in the final section. Optical Properties Absorption The optical properties of gold colloidal aggregates are strongly affected by two features: the electronic plasma resonance that is unique to metal colloids and the fractal structure that characterizes clusters formed by irreversible kinetic aggregation.In qualitative terms, the electronic plasma resonance most strongly influences the absorption of the aggregates, while the fractal structure most strongly influences the static scattering properties. However, each of these two phenomena also has some influence on the other property; the plasma resonance on the scattering and the correlations of the fractal structure on the absorption. In this section we discuss these effects. In fig. 1 we plot both the absorption and the integrated scattering as a function of wavelength for a solution of unaggregated gold colloids and for colloids aggregated by DLCA. These data were obtained using a spectrophotometer fitted with an integrating sphere. Thus the absorption represents the light actually absorbed (transmitted intensity minus scattered intensity), while the scattering represents the total light scattered in all156 Properties of Fractal Colloid Aggregates directions (although there is only appreciable intensity scattered in the forward direc- tion).The measured scattering is corrected for absorption effects. The consequences of the electronic plasma resonance are clearly visible in the absorption of the unaggre- gated colloids, whch has a large peak at ca. 530 nm, resulting in the distinct wine-red colour of colloidal gold. The scattered light is so weak that it it almost unmeasurable with this apparatus, reflecting the very small size of the particles and showing how the plasma resonance affects primarily the absorption. As the aggregates grow, an additional absorption resonance develops at longer wavelength, and the scattering intensity increases.The increased absorption at longer wavelengths caused by this new resonance results in the distinctive change in colour of the colloid upon aggregation, from red to blue. An understanding of the absorption resonances for the aggregated colloid requires an appreciation of the unique correlations between the particles in the ~ l u s t e r . ' ~ At later times there are few if any monomers remaining, yet there still remains a pronounced and distinct resonance at essentially the same wavelength as for the unaggregated colloid. This implies that both absorption peaks are due to gold particles incorporated within the clusters. While each cluster is comprised of many particles, the long-range structure is fractal, meaning that the average density, or volume fraction of gold in a cluster is quite low.Nevertheless, since the cluster is connected, each particle has on average two near neighbours, resulting in very strong and well defined short-range correlations. It is these latter correlations which predominantly determine the optical resonances. The shifts and splitting in the electronic plasma resonance evident in fig. 1 arise from the electromagnetic interactions of the resonances on neighbouring particles. An instructive example of this is seen in the calculations for the optical absorption of a dimer.15 It exhibits two resonances, the excitation of which depend on the orientation of the polarization of the field compared to the axis joining the dimer pair.One resonance is at nearly the same frequency as the single-particle resonance and occurs when the field is polarized normal to the axis joining the two particles. The second is at a lower frequency and occurs when the polarization is parallel to the dimer axis. Physically, the second resonance at longer wavelengths arises from the fact that the electric field extends outside the metal in the narrow gaps which exist when two spheres touch. Thus the near-field coupling between the spheres is basically capacitive in nature and the shift in the resonance has the same physics as that underlying the dielectric anomaly in the Maxwell-Garnet theory.16 It is important to appreciate that even if the spheres are touching and hence electrically connected, the dominant coupling between them is still capacitive because the resistance of a metal contact is extremely large at optical frequencies, since it is only at most a few ingstroms in dialneter and hence is substantially smaller than the mean free path of the electrons.Thus the RC time constant between two touching spheres is very large at optical frequencies. We estimate17 that o >> 1/RC by well over an order of magnitude, confirming the capacitive nature of the coupling. The presence of all the remaining gold particles in the cluster can also have some effect upon the absorption. The sum of the fields of all the remaining particles acting on any given one will have the effect of changing the mean dielectric constant in an effective medium sense.However, the fractal structure of the aggregate means that the mean density of neighbouring particles decreases with distance. Thus we might expect the modification of the effective average dielectric constant to be quite small. To investigate this, we have performed a self-consistent calculation of the optical properties of the clusters using the computer simulations to obtain the exact coordinates of the ~articies.~' The caicuiatioii is based O ~ I Ev;.a!d's self-consistent field method'* in which the local exciting field at each particle site, including multipole scattering effects to all orders, is obtained exactly. The scattered field outside the cluster is then the sum of the radiation emitted by all the particles in the cluster while the absorption is obtained by summing the energy absorbed by each particle.Implicit in this approach is that theH. M. Lindsay et al. 157 polarizability of each gold particle is not given by the expression for a sphere in vacuum. Rather, it is calculated by taking into account all the near-field effects arising from higher-order multipole interactions between the nearest-neighbour spheres. These calcu- lations confirm that the change in the effective average dielectric constant within the cluster is rather small, only ca. 5 % , owing to the relatively low density of gold particles in the fractal structure. Thus the absorption resonances of the aggregated gold colloids are dominated by the strong short-range correlations characteristic of the clusters, with the long-range fractal structure having relatively little influence.While the behaviour of the resonances of a dimer can serve as a qualitative guide, the actual behaviour of real aggregates is considerably more complicated. The short-range structure consists of a combination of chains of gold particles of varying lengths, oriented in all directions. The resonances for each sphere would be determined by their exact local environment, and would include effects of at least the nearest neighbours, and presumably the next few spheres as well. Indeed, the spatial extent of the resonance would depend on the excitation frequency, and would be larger as the long-wavelength resonance is approached. The actual absorption spectrum is thus the result of the combination of the resonances of all the spheres, with their different local environments and orientations combining to cause the observed width.Scattering The gold aggregates scatter much more light than the unaggregated colloid. This is due primarily to the increased size of the clusters. The unaggregated gold particles are much smaller than the wavelength of light, and are therefore relatively ineffective scatterers relative to a cluster of size comparable to a wavelength. We can also see from fig. 1 that the resonances have some effect on the scattering, since there are two peaks in the total scattering as a function of wavelength. This results from the increased polarizability of the spheres caused by the optical resonances. However, rather than considering the integrated scattering from clusters, it is more informative to consider the angular, or k vector, dependence of the scattering, where k = 47rn sin (8/2)/A, with n the index of refraction of the water and 8 the angle between the direction of the incident light and the direction of the scattered light, and A the wavelength. For fractal objects, the long-range correlations can result in the now familiar kPdf dependence in the scattering intensity, when kR, >> 1, where R, is the radius of gyration of the c l ~ s t e r .' ~ When kR,a 1 the scattering becomes independent of k, as the internal fractal structure is no longer probed. To include both regimes, we can represent2' the scattering intensity from a single cluster as where SM(k) = 1/(1+2k2Ri/3df)df'2. This form has the required behaviour that it varies as M 2 for small clusters or small angles and as M / k d f for larger clusters at larger angles.The angular dependence of the scattered light intensity is thus very useful in that it can both determine the fractal dimension of the aggregates at large kR,, and the characteristic size of the clusters, averaged over the cluster mass distribution, at small kR,. The k-df dependence of the scattered light reflects the Fourier transform of the fractal correlations of the particles within the clusters. This implicitly assumes that the field at each particle is simply proportional to the incident field at that position. However, the field at each particle is in fact the local field, which is the sum of the incident field158 Properties of Fractal Colloid Aggregates and the radiated field from all the other particles at the given particle position.Neverthe- less, the scattering properly reflects the Fourier transform of the particle correlations for the values of k probed in light-scattering experiments, as we now show. At any given gold particle within the cluster, the total local exciting field, EL, can be represented as the sum of the incident field, Eo, and the total magnitude of the radiated fields of all the other particles within the cluster, ER, at the given particle. Of course, ER must be calculated self-consistently to properly include the effects of all the particles. However, ER will always be proportional to the incident field. Thus we can express the local field at each gold particle within the cluster as EL = xEo, where x must be determined self-consistently. Furthermore, we can divide the response into x = xave + xfluc where xave is the average value, and is the same at each particle, while xfluc expresses the fluctuations in x and varies from particle to particle.The scattering intensity will then be proportional to the Fourier transform of the spatial correlations of this field. These correlations can be divided into three terms. The first term reflects the correlations of xaveEO with itself. This term will reflect the fractal correlations as it is proportional to the incident field at each particle. The constant of proportionality reflects the change in the average dielectric constant of the cluster, and will affect the magnitude, but not the k-dependence of the scattering intensity.The second term reflects the correlations of xaveEO with xflucEO, and because xRuc has random spatial fluctuations, this term will average to zero. The final term reflects the correlations of xflucEo with itself. Since Xnuc varies randomly in space, it will not exhibit the correlations of the fractal structure. Indeed, we expect no long-range correlations whatsoever. Instead there will only be short-range correlations whose spatial extent will be on the order of a single-sphere radius. Thus, at k vectors probed in light-scattering measurements, this term will be independent of k. Furthermore its magnitude will in general be small, and therefore this term will not contribute significantly. Thus the local field corrections do not contribute to the angular dependence of the light-scattering intensity.These qualitative ideas have been confirmed in a more quantitative fashion using the self-consistent field calculations with the computer-generated clusters. l 7 We calculate the scattered field exactly, properly accounting for the corrections to the local field. This can then be compared to a calculation of the scattering using only the incident field. To compare to experimentally measured results, we average over orientations and sum over many different clusters. When this is done, there is no observable difference between the exact calculations and the ones excluding the local field corrections. In both cases the fractal correlations are reflected in the kPdf dependence of the scattered intensity when kR, >> 1.This confirms the results of our qualitative arguments above. There is one feature unique to light scattering that does allow the effects of the short-range correlations of the local field to be probed. This is the depolarized scattering, or the intensity of scattered light polarized normal to the incident field. Depolarized scattering arises because of the anisotropy induced in the polarizability of each sphere due to the presence of the nearest neighbours. Since the incident field is polarized, the depolarized components can only consist of contributions from the corrections to the local field due to the other particles in the cluster. The k dependence of the depolarized scattering reflects the Fourier transform of the spatial correlations of the local field corrections.Thus the only contribution will arise from the correlations of xflucEo with itself, where we now are considering the depolarized component of this quantity. Again, since these terms vary randomly in space, there will be correlations only on length scales of the order of a single-sphere radius, and the depolarized scattering should be isotropic in k over the range probed in light-scattering experiments. This is indeed observed as is shown in fig. 2, where both the polarized and depolarized light scattering are shown as a function of k on a logarithmic plot. The data were obtained from clusters prepared by DLCA, and were obtained using the 633 nm HeNe laser line. The polarized scattering exhibits fractal correlations, with the slope giving df = 1.77. By contrast, the depolarizedH.M. Lindsay et al. I C 159 + + + + + + + + + + + + + t + + t ++ ++ + t t +++++++* 0.001 0.005 0.01 0.02 0.05 k/nm-’ Fig. 2. Scattering intensity from diffusion-limited aggregated gold colloids as a function of k. The polarized scattering (upper curve) exhibits the linear behaviour on the logarithmic plot expected for the fractal correlations with d , = 1.77. The depolarized scattering (lower curve) is isotropic. scattering intensity is independent of k, reflecting the complete absence of correlations in the regime of k probed. Since the depolarized scattering arises from the anisotropy of the polarizability of the spheres induced by the nearest neighbours, we would expect it to be strongly influenced by the electromagnetic plasma resonance.This is indeed observed, as can be seen by the variation in the intensity of the depolarized scattering with exciting wavelength. When the exciting wavelength is 488 nm, the depolarized intensity is very weak, on the order of 1 O/O of the polarized intensity at high k. As the excitation wavelength is increased and approaches the second absorption resonance, the depolarized intensity also increases, and is ca. 6% of the high k polarized intensity at 633 nm and ca. 10% at 752 nm. This behaviour is consistent with the behaviour of the coupled electromagnetic plasma resonance. The long-wavelength resonance involves fields polarized parallel to the axis joining the spheres and will therefore induce the largest degree of anisotropy in the polarizability.Thus the magnitude of the depolarized scattering will increase as this resonance is approached. Finally, we note that Schaefer and collaborators have claimed2’ that the fractal dimension of colloidal gold aggregates cannot be determined from light-scattering measurements because they are distorted by ‘multiple scattering effects’. The results of this section show that Schaefer’s argument is incorrect and confirm that it is based on an erroneous interpretation of his data.* Dynamic Scattering Dynamic, or quasielastic, light scattering22 has proved quite useful in the study of the kinetics of the aggregation process. It measures the time dependence of the intensity fluctuations of the scattered light. Normally these fluctuations arise from the diffusional * Schaefer observed the slope of the scattering to change when he diluted his colloids.This was due to the shear-induced restructuring of the colloids upon dilution (see next section of this paper). Schaefer misinterpreted this to indicate ‘multiple scattering’ from the aggregates.160 Properties of Fractal Colloid Aggregates loo*oo i 10000 1000 A c c1 .- v) u 2 100 8.01 0.02 0.05 0.1 0.2 0.5 1 ka Fig. 3. The scattering intensity calculated for computer-generated DLCA clusters as a function of ka where a is the diameter of a single sphere. The solid line is the result for a single cluster of mass 500 particles in a single orientation, and exhibits large fluctuations when kR,> 1. The dashed line is an average over 10 clusters and 100 randomly chosen orientations for each, and exhibits the expected linear, fractal behaviour when kR,> 1.motion of the particles in the scattering volume. However, for aggregated colloids, the size of the clusters can rapidly become much greater than k-I, so that the scattering intensity from each cluster is no longer proportional to M 2 , as it is for very small clusters. Nevertheless, since we know the structure of the clusters, we also know the mass dependence of the average scattering intensity. Indeed, if we use a scattering angle such that kR,> 1 for most of the clusters, the average scattered intensity from each cluster will scale as M. However, ths dependence requires averaging over all orientations of the cluster. In actual fact, for a given orientation of a cluster, the scattering will not scale this way at all, and indeed will not even have the power-law dependence expected for a fractal. This is illustrated in fig.3, where we have calculated the scattering from a computer-generated DLCA cluster comprised of 500 particles. We compute I C, exp ( - k - r,) 12, where the sum extends over all the particles in the cluster. The large oscillations in the scattering as a function of k are clearly visible, reflecting the highly disordered structure of the cluster. However, if we average the intensities from ten different clusters calculated for 100 randomly chosen orientations, we regain the expected power-law dependence, as shown by the dashed curve in fig. 3. These results show the importance of the orientational averaging in determining the static scattering.However, we cannot ignore these variations in dynamic light scattering, as they will manifest themselves directly as fluctuations as the cluster undergoes rotational diffusion. This will add an additional decay mechanism in the autocorrelation of the scattered intensity. Physically, this arises from the fact that the cluster is larger than k-' and so is comprised of many different blobs of size k-', whose scattering can interfere. The exact resultant intensity is then very sensitively dependent on the orienta- tion of the cluster. One immediate consequence of this is that it is possible to see quasielastic scattering from a single cluster, even in a homodyne experiment. To investigate the effects of rotational diffusion, we have developed a theory for the autocorrelation of the intensity fluctuations of inhomogeneous objects.23 This theory is a generalization of the treatment for a long rod.22 To apply our theory, we need toH.M. Lindsay et al. 161 know the position of all the particles in the clusters and thus use the computer-simulated cluster aggregates. We calculate the field-field correlation function for a set of rigid clusters made up of identical spheres. Experimentally we measure a homodyne intensity- intensity correlation function; in the limit of an infinite number of clusters, our measured correlation function will be the square of the field-field correlation function plus a baseline. If we have a small number of clusters in the scattering volume, this condition no longer holds.In this paper we deal only with the Gaussian limit of a large number of clusters, as it is usually achieved in our experiments. For simplicity, the spheres are taken as point scatterers ( i e . the form factor is 1 ) . Then the field-field correlation function is where r y ( t ) is the position of particle i in cluster a at time t, n, is the number of particles in cluster a, and a0 is the polarizability of each particle. We now make some simplifying assumptions. We assume that our sample is sufficientiy dilute so there are no correlations between clusters. We assume that the clusters are not too anisotropic, so we may represent translational diffusion of a cluster with a single diffusion coefficient, 0, = kBT/67rqRh, where q is the viscosity of water and Rh is the hydrodynamic radius of a cluster.Finally, we assume that translational and rotational diffusion are independent. This means that we may decompose the motion of the particles into the motion of the centre of mass of each cluster, which will depend on translational diffusion, and the motion of each particle about the centre of mass of its cluster, which depends on rotational diffusion. We then separate terms in the autocorrela- tion function and write it as S ( k , t ) = a i C e x p (-Dak2t) C (exp{-ik.[by(t)-b;’(O)]}) ( 2 ) where b y ( t ) is the displacement of particle i from the centre of mass of cluster a. The first exponential reflects the decay in the autocorrelation function due to the centre-of- mass motion, while the second contains the rotational component.To calculate the rotational contribution, we expand the positions of the particles in each cluster as a series of multipole moments, and calculate the scattering from each moment. For each cluster we use a single rotational diffusion coefficient, 8, = kBt/87r77R:. The rotational contribution can then be expressed as a sum over these multipole moments; (Y i, j ‘2.2 C ST(k) exp [-Z(Z+ l)O,t] /=0 where (3) where j , is the fth order spherical Bessel function, YIM are the spherical harmonics and Cnp is the angle of displacement of particle i relative to an arbitrary axis of cluster a. Combining the translational and rotational terms, the field-field correlation function is 02 S ( k , t)=a;Cexp(-D,k*t) C S p ( k ) exp[-Z(Z+1)8,t].( 5 ) a / = 0 The rate of decay due to translational diffusion depends on k2, while rotational diffusion is independent of wavevector; however, the relative magnitude of each multipole term will vary considerably with k. Thus the effects of rotational diffusion will alter the k dependence of the decay of the autocorrelation function. Using the computer-generated clusters, we calculate the S;Y( k ) terms, and sum to obtain the correlation function. In fig. 4 we plot S , ( k ) as a function of k averaged over162 Properties of Fractal Colloid Aggregates 0.01 0.02 0.05 0.1 0.2 0.5 t Fig. 4. Contributions of the spherical harmonics from 1 = 0 to 7 to the scattering intensity calculated for 39 computer-generated clusters of mass 900 particles. The 1 = 0 term only contributes to the translational decay of the autocorrelation function, while the 1 > 0 terms also contribute to the rotational decay.The static scattering is the sum of all the terms. 39 900 particle clusters, for 1 = 0 to 7. For kR, s 1 only the So term contributes, and the decay will be that expected for purely translational diffusion. Physically this reflects the fact that at small wavevectors the cluster scatters as a point, so the orientation of the cluster is irrelevant. As k R , a 1 we can distinguish regions within the cluster, so orientation is important. For these k values the higher-order S, terms become important, and so rotational diffusion makes a greater contribution to the total decay of the autocorrelation function. To compare our calculated results with experiment, we consider the k dependence of the first cumulant, r, which is the logarithmic derivative of the autocorrelation function as t --+ 0.We use an ensemble of computer-generated clusters with a distribution similar to that produced experimentally by DLCA, and calculate r as a function of wavevector. To determine the rotational contributions, we use Rh = 0.8Rg,f as determined numerically for the computer-generated c1uste1-s.~~ We compare the calculations with experimental results in fig. 5 , where we have divided r by k2 to remove the simple dependence on k. For a purely translational system, r/ k 2 = Do, the mean translational diffusion constant for the distribution of clusters. Our theoretical values for T / k 2 approach Do at small wavevectors, and increase by a factor of ca.2 at !arge k, showing the significant contribution due to rotational diffusion. We note that some of the k dependence of r is due to the fact that we are measuring different moments of the cluster mass distribution as S , ( k ) changes as k is varied, but most of the k dependence is due to rotational effects. The data points in fig. 5 are experimentally measured first cumulants for DLCA aggregates with roughly the same size as our model clusters. The agreement of our calculation with the data is quite satisfactory. The calculations suggest that most of the variation in the experimentally measured decay with wavevector is due to rotational effects, although other effects may also contribute. Finally, we note that while a major contribution to r still comes from the exp (-Dk2t) term, due to translational diffusion, as kR, increases rotational diffusion plays an increasingly important role in quasielastic light scattering.this paper (see Discussion). i We thank Dr P. N. Pusey for pointing out an error in the ratio Rh/ R, used in the original version ofH. M. Lindsay et al. 163 0.002 0.005 0.01 0.02 0.05 k/nm-' Fig. 5. The k dependence of the first cumulant, normalized by k2, for a distribution of DLCA clusters. The solid line is the calculated behaviour, including effects of rotational diffusion. The points are experimental values. Mechanical Properties All of the aggregates we have discussed so far have been formed under simple aggregation conditions, without the presence of additional forces of any sort.However, very often subsequent processing of flocs will subject the aggregates to a wide variety of stresses. We might expect the clusters to be quite susceptible to deformation, given their very tenuous nature. Indeed, a recent calculation of the mechanical properties of colloidal aggregates found that they could possess a fractal structure only up to some critical size when they become unstable to either thermal fluctuations or gravitational d i ~ t o r t i o n . ~ ~ In this section we discuss the effects of shear on the structure of the aggregates and the consequences on the static scattering from them. To apply a controlled amount of shear to the clusters, they were first aggregated by DLCA until the characteristic cluster radius was ca.5000 A, as measured by dynamic light scattering. Then they were flowed through a narrow tube with various pressure drops across the length of the tube to change the flow velocity. The k dependence of the static scattering from the clusters was then compared with that obtained from the same colloid not subjected to shear. An example of the results obtained is shown in a logarithmic plot in fig. 6. The lower curve is from the unsheared aggregates, and exhibits the expected linear behaviour giving df = 1.84. By contrast, the scattering from aggregates subjected to a maximum shear of ca. lo4 Hz, shown in the upper curve, no longer exhibits a simple linear shape. Instead, there is a distinct kink in the curve at k = 0.006 nm-', which corresponds to a length of ca.1700 A. At larger k the scattering still has the same slope as that from the unsheared clusters. However, at small k the slope of the scattering has markedly increased, although it is not really possible to even describe it as linear on the logarithmic plot given the limited range of k probed. The change in the scattering results from shear-induced restructuring of the tenuous fractal structure of the clusters. A plausible picture of the effects of the shear is that it causes the fractal aggregates to bend or deform, causing loops to be formed, thus strengthening the structure. These loops would tend to be formed on larger length scales, where the forces are greater and the structure weaker. Indeed, X-ray scattering164 Properties of Fractal Colloid Aggregates 1 0 .1 C .- 0.001 I I I O - 0°01 do2 0.005 0.01 0.02 0 . k/nm-' Fig. 6. Scattering from DLCA clusters before (lower curve) and after (upper curve) being subjected to shear. The shear causes restructuring of the fractal shape, changing the slope of the scattering at low k. The inset shows schematically the expected effect of the restructuring on the scattering. from gold aggregates, which probes higher values of k then those observed with light, suggests that the fractal structure at small length scales is much more robust, and the restructuring occurs primarily at larger length^.'^ The expected effects of this restructuring on the scattering from the aggregates is shown schematically in the insert of fig. 6. The unsheared aggregates exhibit a linear behaviour in the fractal regime, with I E k-df, until kR, == 1, at which point the scattering becomes independent of k.Here, R , represents an average radius of gyration of the clusters in the distribution. The intensity of the scattering at low k will be proportional to M 2 averaged over the cluster distribution. Since this does not change upon shearing, the scattering from the sheared sample will have the same intensity at very low k. Furthermore, if the small length scale structure is robust and remains fractal, the scattering at high k will also be unchanged. However, the formation of loops will result in a decrease in R,, so that the isotropic region of the scattering will have to extend to larger k. Then, to join properly with the scattering from the fractal region, the scattering in the region of k which probes length scales that have been substantially restructured will have to have a larger slope, as shown by the dashed line in the inset of fig.6. If this region matches that probed in our light-scattering experiments, this simple argument would account for the behaviour observed. An important consequence of this picture is that the restructured aggregates are not necessarily fractal at all on the length scales over which they are restructured. The scattering would neither be linear on a logarithmic plot nor would it approach any well defined slope. This is indeed what we observe: the apparent slope of the low k scattering varies with the magnitude of the applied shear. Furthermore, it does not always appear t o be linear, although we rarely can probe over a sufficient range in k to determine this unambiguously.We note that these observations are in contrast with those for silica colloids.26 While our results are still preliminary, we can summarize the qualitative features as follows: as the shear increases, the restructuring extends to shorter length scales, as determined by the position of the change in slope in the scattering. Greater shear alsoH. M. Lindsay et al. 165 apparently leads to an increase in the apparent slope of the scattering in the region of k probing the restructured shape. Finally, the aggregates prepared by DLCA are considerably more susceptible to restructuring than those prepared by RLCA. The Iatter can exhibit some restructuring effects, but require substantially larger shear. This is consistent with the idea that the larger fractal dimension of the RLCA clusters, df= 2.1, results in a more robust structure.Conclusion To summarize, we have found that the optical absorption of gold colloidal aggregates is dominated by the electronic plasma resonances in the gold spheres, which change from a single resonance in the unaggregated colloid to a pair of resonances in the aggregates due to the electromagnetic interactions between neighbouring spheres. The polarized scattering from the aggregates is dominated by the fractal structure of the clusters, giving a static structure factor S ( k ) K k - d f . The local field corrections due to the neighbouring particles have no effect on the k dependence of the polarized scattering.By contrast, the depolarized scattering is independent of wavevector, and reflects the short range of the spatial correlations of the corrections to the local field due to the other particles in the cluster. In dynamic light scattering we find that rotational diffusion can make a significant contribution to the decay of the autocorrelation function of the scattered light for kR, 2 1 , but translational diffusion remains dominant for all wavevec- tors. Finally, we find that large clusters are subject to restructuring under mechanical stress. At large length scales, where the clusters are weakest, diffusion-limited clusters lose their original fractal behaviour. By contrast, reaction-limited clusters also restruc- ture, but to a significantly lesser degree, owing to their higher initial fractal dimension. References 1 D. A. Weitz, M. Y. Lin, J. S. Huang, T. A. Witten, S. K. Sinha and J. S. Gethner, in Scaling Phenomena 2 D. A. Weitz and M. Oliveria, Phys. Rev. Lert., 1984, 52, 1433. 3 P. Meakin, Phys. Rev. Lett., 1983, 51, 1119. 4 M. Kolb, R. Botet and R. Jullien, Phys. Rev. Lett., 1983, 51, 1123. 5 D. A. Weitz, J. S. Huang, M. Y. Lin and J. Sung, Phys. Rev. Lett., 1985, 54, 1416. 6 D. A. Weitz, M. Y. Lin and C. J. Sandroff, Surf: Sci., 1985, 158, 147. 7 J. A. Creighton, C. G. Blatchford and M. G. Albrecht, J. Chem. Soc., Furaday Trans. 2, 1979, 75, 790. 8 M. Kolb and R. Jullien, J. Phys. Lett., 1984, 45, L977. 9 W. D. Brown and R. C. Ball, J. Phys. A, 1985, 18, L517. in Disordered Systems, ed. R. Pynn and A. Skjeltorp (Plenum, New York, 1985), p. 171. 10 R. C. Ball, D. A. Weitz, T. A. Witten and F. Leyvraz, Phys. Rev. Lett., 1987, 58, 274. 11 D. A. Weitz and M. Y. Lin, Phys. Rev. Lett., 1986, 57, 2037. 12 D. A. Weitz, J. S. Huang, M. Y. Lin and J. Sung, Phys. Rev. Lett., 1984, 53, 1657. 13 P. Meakin, J. Colloid Interface Sci., 1984, 102, 491. 14 P. Dimon, S. K. Sinha, D. A. Weitz, C. R. Safinya, G. S. Smith, W. A. Varaday and H. M. Lindsay, 15 P. K. Aravind, A. Nitzan and H. Metiu, Surf: Sci., 1981, 110, 189. 16 R. Landauer, in Electrical Transport and Optical Properties of Inhomogeneous Media (American Institute 17 Z. Chen, P. Sheng, D. A. Weitz, H. M. Lindsay, M. Y. Lin and P. Meakin, to be published. 18 M. Lax, Rev. Mod. Phys., 1951, 23, 287. 19 D. W. Schaefer, J. E. Martin, P. Wiltzius and D. S. Cannell, Phys. Rev. Lett., 1984, 52, 2371. 20 G. Deitler, C. Aubert, D. S. Cannell and P. Wiltzius, Phys. Rev. Lett., 1986, 57, 3117. 21 J. P. Wilcoxon, J. E. Martin and D. W. Schaefer, in Fractal Aspects of Materials (Mat. Res. SOC. 22 B. J. Berne and R. Pecora, Dynamic Light Scattering (Wiley, New York, 1976). 23 H. M. Lindsay, R. Klein, M. Y. Lin, D. A. Weitz and P. Meakin, to be published. 24 P. Meakin, Z. Y. Chen and J. M. Deutch, unpublished. 25 Y. Kantor and T. A. Witten, J. Phys. Lett., 1984, 45, L675. 26 C. Aubert and D. S . Cannell, Phys. Rev. Lett., 1986, 56, 738. Phys. Rev. Lett., 1986, 57, 595. of Physics, New York, 1978), p. 2. Extended Abstracts, 1985), p. 33, and unpublished. Received 5 th January, 1987

ISSN:0301-7249    

DOI:10.1039/DC9878300153

出版商:RSC    

年代:1987

数据来源: RSC

摘要:

Faraday Discuss. Chem. SOC., 1987, 83, 167-177 Brownian-dynamics Simulation of the Formation of Colloidal Aggregate and Sediment Structure Geoffrey C. Ansell and Eric Dickinson" Procter Department of Food Science, University of Leeds, Leeds LS2 9JT Brownian-dynamics simulation is used to investigate the effect of particle volume fraction and interparticle interactions on the structure of colloidal materials formed by irreversible aggregation. Short-range structure of aggre- gates and sediments is represented by the pair distribution function g ( r ) , and longer-ranged structure by the effective fractal dimension. For unstable DLVO-type particles with a secondary minimum, the simulated short-range aggregate structure resembles that in a stable concentrated dispersion; for particles without a secondary minimum, g ( r ) is relatively unstructured either in sediments or gels.Hydrodynamic interactions play a crucial role in the formation of dense sediments at high field strengths, but are much less important when Brownian motion is predominant. Studies of small-floe dissociation give insight into the likely consequences of neglecting multi- body hydrodynamic interactions in simulations of this type. Brownian motion is the primary process whereby particles come together to form colloidal aggregates and gels. Brownian motion plays a crucial role in the dissociation of small colloidal flocs, in the desorption of particles from surfaces, in the unfolding and denaturation of proteins, in the kinetics of biological reactions and in many other diffusive phenomena.In sedimentation or creaming, Brownian motion has a secondary role compared with the action of gravity; but it is nevertheless important in determining the final sediment or cream structure. The term colloidal particle is essentially synony- mous with the term Brownian particle. To a large extent, then, colloid science is the science of Brownian motion. To simulate a system containing solute particles and a very much larger number of solvent molecules (plus ions etc.), where one is interested only in the dynamics of the solute particles, it is convenient to use a stochastic-dynamics approach rather than the more conventional molecular-dynamics approach used in connection with simple liquids. Where inertial terms are negligible, the motion is entirely diffusive, and the technique is called Brownian dynamics.In a Brownian-dynamics simulation, each particle diffuses in a force field caused by the presence of neighbouring particles, with spatial correlations in the motion determined by position-dependent hydrodynamic interactions between the particles. In this paper we report on how the technique of Brownian-dynamics simulation can be applied to the problem of aggregate and sediment structure formation. We consider in particular how the colloidal structure is affected by such factors as particle size, dispersed-phase volume fraction, external field strength, hydrodynamic interactions and the nature of the interparticle colloidal forces as expressed through the classical DLVO theory. Emphasis will be placed on an analysis of short-range structural features which are not adequately described by simple lattice-based simulation models.And the possible importance of multi-body hydrodynamic interactions will be assessed by considering their effect on the kinetics of small floc dissociation. In a low-density sediment or gel formed from irreversibly aggregated spherical colloidal particles, three spatial scales of structure may be identified: ( a ) short-range order associated with packing and excluded-volume effects, ( b ) medium-range disorder 167168 Brownian -dynam ics Sirnula tion associated with the fractal characteristics of diff usion-limited aggregation and (c) long- range uniformity as expected for a material that is homogeneous on the macroscopic scale.These three spatial scales can be described' by three regions of the normalized particle-particle distribution function G( r ) for pairs separated by centre-to-centre dist- ance r : In eqn (1) a is the particle radius, g ( r ) is the short-range liquid-like pair distribution function which extends out to r = 5, D is the characteristic fractal dimension and 6 is the characteristic network length of the connected gel or sediment. Experimental light-scattering measurements on coagulating polystyrene latices2 and a recent X-ray study of aggregated gold colloids' have both shown that short-range structure as measured by g ( r ) is sensitive to whether the aggregation is rapid or slow. The presence of an intermediate fractal scaling regime is indicative of a material produced from aggregate precursors; a recently reported example4 is the so-called silica aerogel ( D ==: 2.0, c / a = 10).In gels or sediments formed at high or intermediate particle concentration, the fractal regime could well be missing altogether ( 5 = 6). Brownian Dynamics with Hydrodynamic Interactions The algorithm used to simulate the diffusive motion of interacting colloidal particles was devised by Ermak and McCammon.' The translational displacements of N particles during a single time-step A t are given by where x, is the coordinate vector of particle i, FJ denotes the non-hydrodynamic force associated with particlej, D, is the diffusion tensor, k is Boltzmann's constant, T is the absolute temperature and the symbol O denotes quantities at the beginning of the time-step.The stochastic term R, has first and second moments: It is implicit that A t is much larger than the particle momentum relaxation time (m,D,,/ kT for a particle of mass m,), but small enough that FJ and D, remain sensibly constant over the interval. Extended forms of the basic algorithm can be written to allow for rotational Brownian motion6 and the effect of an external flow field.' Computer simulations of systems containing hundreds of particles are feasible only if the interactions between individual particles are of a relatively simple form. In practice, this means neglecting interactions altogether when their magnitude is small, making use of good approximations when their magnitude is large, and moving away from pairwise additivity only as a last resort.With charge-stabilized dispersions, the pairwise additivity assumption of the colloidal interactions is valid so long as the electrolyte concentration is not too low (i.e. K a >> 1, where K - ' is the Debye length). Pairwise additivity of hydrodynamic interactions is not such a good approximation, however, and its adoption in concentrated systems can lead to physically absurd predictions (e.g. negative diffusion coefficients') and computational difficulties in Brownian-dynamics simulations due to non-positive-definite diffusion matrices.' Hydrodynamic interactions are included in the simulation through the mobility tensor p, = D,/kT. For a pair of spheres of radius a, the mobility tensor can be written in the form' O (R,(At))=O, (R,(At)R,(At)) =2D;At.(3)G. C. Ansell and E. Dickinson 169 where r is the absolute scalar separation, $12 is the unit vector between particles 1 and 2, A,( r ) and B,( r ) are distance-dependent mobility coefficients, and 7 is the viscosity of the medium. Explicit expressions for the coefficients A, and B, (i, j = 1,2) in powers of ( a / r ) are available" to order r-I2. To order r-3, with stick boundary conditions, eqn (4) reduces to the Rotne-Prager tensor:12 r;'= (6.rrr7a)-"6,1+ ( 3 / 4 ) ( a / N - 6,)(1+ t 1 2 t 1 2 ) - ( 1/21 ( a / r13( 1 - 8, (3 31 2312 - 1 I. ( 5 ) Eqn ( 5 ) is a convenient tensor from the computational viewpoint, since it is algebraically simple and also physically well behaved. It always gives a positive-definite 3 N x 3N diffusion matrix for r > 2 a in N-particle systems, and the divergence term in eqn (2) can be omitted because we have V .r f' = 0. To go beyond Rotne-Prager hydrodynamics, one is compelled to confront the problem of multi-body interactions. Mazur and van Saar100s~~ have presented a general scheme for obtaining the mobility tensor for an arbitrary number of spheres to any desired order in ( a / r ) . Dominant contributions to translation from clusters of N spheres are shown13 to be of order r - ( 3 N - 5 ) , which means that three-body interactions first appear at order rb4, and four-body interactions at order r-'. The results including three- and four-body interactions were, in fact, anticipated by KynchI4 more than 25 years ago, but, because of their apparent complexity, Kynch's results have remained unused, although not ~nnoticed.'~ To order r-4, the mobility tensor has the formI3 Higher-order terms become increasingly more complicated, and therefore expensive to include explicitly in simulations. One way round the problem is to use an effective screened-pair mobility tensor:9716 it gives reasonable predictions of diffusive behaviour, but the approach has been criticized" as lacking in rigour.On the more pragmatic level, recent studies of fractal aggregate formation1s would seem to suggest that aggregate structure is insensitive to the exact form of particle and aggregate diffusion coefficients. So, in simulations of colloid coagulation kinetics, one may in practice be able to avoid the complication of multi-body hydrodynamics. And, if one is interested only in aggregate structure, it may not be unreasonable to neglect inter-aggregate hydrodynamic interactions altogether.Effect of Hydrodynamic Approximation on Small-floc Dynamics To assess quantitatively the effect of three- and four-body hydrodynamic interactions on the average dynamical behaviour of interacting colloidal particles, we describe the process of Brownian dissociation, whose simulation we have considered previ~usly'~ only at the level of two-body hydrodynamics. We consider groups of three or four DLVO-type particles flocculated in the secondary minimum. The van der Waals attraction UA(r) is calculated from UA( r ) = -(AH/ 12) { [ 4a2/ ( r2 - 4a2)] + ( 2 a / r ) 2 + 2 In [ 1 - ( 2 a / r)']} (7) where A, is the Hamaker constant. The screened Coulombic repulsion, UR(r), is calculated from ~ , ( r ) = 2 .r r ~ , ~ ~ a & In [I +exp ( - K S ) ] ; ( ~ a >> 1) (8) where E~ is the relative dielectric constant of the medium, E~ is the permittivity of free space, t+b0 is the surface potential, s = r - 2a is the surface-to-surface separation, and K170 Brownian-dynamics Simulation Table 1. Mean dissociation time f of small flocs as simulated for various levels of hydrodynamic approximation [up to terms of order ( a / r ) " ] equilateral linear n trimer trimer tetramer 0.37 0.49 0.54 0.32 0.42 0.39 0.42 0.5 1 0.58 0.35 0.50 0.58 0.39 0.53 0.47 estimated errorn k0.03 k0.04 *0.07 ' Taken as *2 standard deviations. is the inverse Debye length, defined for a 1 : 1 electrolyte by K~ = 2e*CN/~,.r,kT (9 1 where e is the electronic charge, N is Avogadro's number and C is the ionic strength.Potential parameters in eqn (8) and (9) are set as follows: a = 0.5 pm, q90 = 30.5 mV, T = 300 K, E , = 80 and C = 1.5 x lop4 mol dmP3. Dissociation events are studied, as a function of the level of hydrodynamic approxi- mation, for trimers starting from both linear and equilateral configurations, and for tetramers starting from a tetrahedral configuration. In trimer simulations the parameter AH is set at 0.83 x J, producing a pair potential with a high primary maximum (> 60 kT) and a shallow secondary minimum (0.8kT at rmin = 2 . 3 5 ~ ) . Initial pair separ- ations are rI2 = r23 = rmin and rI3 = 2rm,, for linear trimers, and r,? = r23 = r I 3 = rmin for equilateral trimers, and a trimer is deemed to have dissociated when each of the pair separations exceeds 1.5 p m ( 3 a ) .A total of 300 separate dissociation events are recorded for both linear and equilateral trimers. In tetramer simulations the parameter AH is set at 1.66 x J, corresponding to a potential with secondary minimum of 2.4kT at rmin = 2 . 3 ~ . Initial pair separations are rI2 = ~ 2 3 = r34 = rI3 = r,4 = r24 = rmin, and a tetramer is deemed to have dissociated when one particle gets to be at least 1.5 p m from the rest. Averages are taken over 100 separate tetramer dissociation events. The simulation time-step is set at 20 ps, and the medium viscosity is taken as 0.865 mPa s (water at 300 K). Table 1 lists values of the mean dissociation time t obtained by truncating the series expansion at various values of (air)".(Terms of order F2 and F 5 are identically zero even for liquid particles of finite viscosity.") The results show that the Rotne-Prager tensor ( n =3) leads to faster dissociation than the leading three-body term ( n = 4). Within the statistical uncertainties, however, the other approximations ( n = 1,6 and 7) give essentially the same times. The data in table 1 follow a pattern of behaviour that is consistent with a more detailed analysis of trimer trajectories*' in which the trend of inferred relative mobilities ( r P 3 > r-' > r-'> r-' > F4) is rarely broken. Overall, perhaps the most significant general trend is for the two-body approximations to give rise to more rapid dissociation than the multi-body approximations, although we would not wish to propose this as being a universal rule.Trajectory studies of sets of systems subject to exactly the same Brownian fluctuations (same sequences of random numbers used to compute { R, >) have demonstrated" that relative particle configurations are very insensitive to the value of n ; only the spatial scale of relative diffusion changes (at constant time) [or, alternatively, the time-scale changes (e.g. f) over constant spatial scale].G. C. Ansell and E. Dickinson 171 Fig. 1. Pair potentials of mean force A and B. The energy U (reduced by k T ) is plotted against centre-to-centre separation Y (reduced by a ) . These results on small-floc dynamics can be used to justify the neglect of multi-body hydrodynamic interactions in simulations of aggregation of DLVO-type particles, especially where one is primarily interested in the structure and size distribution of aggregates, as opposed to the coagulation rate constant.Pairs of particles spend most of their time at separations for which either hydrodynamic interactions can be ignored altogether ( n = 0) or the Rotne-Prager tensor [eqn ( 5 ) ] is adequate. Owing to the strong van der Waals attraction U,( Y), particles spend little time at close separations ( Y S 2 . 1 ~ ) before coagulating irreversibly, and so the nature of plJ at close separations is of little structural consequence. Moreover, at these close separations, the higher-order approxi- mations can give non-positive-definite diffusion matrices, thereby leading to computa- tional problems.Probably the best way to allow for multi-body interactions in a concentrated dispersion is to use the Rotne-Prager approximation, but to replace the viscosity 7 of the medium in eqn (5) by the effective viscosity of the dispersion.** In effect, this just means a simple scaling of the coagulation time-scale by the appropriate viscosity ratio. Effect of Colloidal Interaction and Field Strength on Sediment Structure A procedure of single-particle deposition has been used to simulate sediment formation in a very dilute system of coagulable DLVO spheres. The two pair potentials considered here, A and B, are plotted in fig. 1. Both are weakly attractive for moderately large separations ( r 2 3 a ) . The parameter values in eqn (7) and (8) are set as follows: a = 0.25 pm, Go = 20 mV, T = 300 K, E~ = 80 and AH = 1.94 x J.Curve A corre- sponds to C = 3 x lop4 mol dm-3 in eqn (9); it has a well defined secondary minimum at r = 2.3a, and a small primary maximum at r==2.la. Curve B corresponds to C = 6 x lop4 mol dmp3, and is attractive at all separations. With a sedimenting force contribut- ing to the systematic motion in eqn (2), Brownian particles were accumulated one at a time (so-called particle-cluster aggregation) until a sediment of 2500 particles had been deposited. A Rotne-Prager form was assumed for the hydrodynamic interaction between the moving particle and its nearest neighbour in the existing deposit. Statistical averages were derived from sets of three or five runs at each set of simulation conditions.Full details of the methodology have been reported e l ~ e w h e r e . ~ ~ Type A particles gave sediments with volume fractions of 0.12, 0.20 and 0.33 at field strengths of 500g 15OOg and 5000g, respectively, where g is the normal acceleration due172 Brown ia n-dynam ics Simulation I I L 2 4 6 8 01 r l a 0 ’ I I t 2 4 6 8 r / a Fig. 2. Plots of normalized pair distribution function g ( r ) as a function of reduced separation r / a . ( a ) Potential A: (-) sediment formed at 1500g; (- - -) 512-particle cluster formed at 4 = 0.2. ( b ) Potential B: (-) sediment formed at 500g; (- - -) 512-particle cluster formed at 4 = 0.1. to gravity. Fig. 2(a) shows g ( r ) for the sediment formed with potential A at 1500g. The high probability of pairs at r=2a corresponds to nearest neighbours in contact; the nearby peak at r Z 2 .3 ~ arises from non-bonded pairs in the secondary minimum; correlations disappear into the statistical noise for rb6a; there is no evidence for a fractal scaling regime. The sediment volume fraction for potential B at 500g is 0.09, almost exactly the same value as that obtained with pure Brownian trajectories (no colloidal or hydrodynamic forces), but significantly lower than that found with potential A (0.12). Fig. 2 ( b ) shows g ( r ) for the sediment formed with type B particles. The structure near r=2.3a is now lost owing to the absence of a secondary minimum. Sediments simulated with simple ballistic trajectories (no Brownian motion and no colloidal or hydrodynamic forces) were of volume fraction 0.14, in good agreement with the pioneering simulations of V ~ l d .~ ~ The simulations show that sediment structure is sensitive to the form of the coagulat- ing potential. Specifically, the presence of a secondary minimum leads to denser sediments. Also, the volume fraction increases with increasing field strength, as new particles tend increasingly to fill up holes in the existing structure, instead of sticking immediately to protrusions at the surface of the sediment. Hydrodynamic interactions have little effect on sediment structure at low field strengths, but this is not the case at high field strengths, where the Brownian component of the motion is swamped by the systematic component.’ At high field strengths the timescale for the moving particle to ‘slip past’ a fixed particle in the sediment becomes lower than the Brownian coagulation time, and so immediate ‘sticking’ does not occur, unless the line of approach is so close as to lie inside the critical capture t r a j e ~ t o r y .~ ~G. C. Ansell and E. Dickinson 173 Table 2. Half-life, t 1 / 2 , and total coagulation time, tc, as a function of volume fraction 4 and colloidal potential (A or B) potential A potential B 0.05 90* 10 1100 f 300 10*2 1200 f 400 0.10 65*7 520 * 60 1.2k0.3 200 * 30 0.15 43+6 120 * 30 0.36 f 0.04 49*9 0.20 40*4 110*30 0.22 * 0.04 6.4 * 0.5 0.25 28*3 90 * 20 0.13 * 0.03 4.2 * 0.4 0.30 17*3 82 f 20 0.10*0.02 1.8 f 0.2 Effect of Colloidal Interaction and Particle Concentration on Aggregate Structure In the final section of this paper, we report simulation results for the Brownian coagula- tion of a non-dilute system of 512 DLVO particles in a cubic box with periodic boundary conditions.Here we ignore inter-aggregate hydrodynamic interactions and aggregate rotation, and we focus attention on the nature of the coagulating potential and the overall particle volume fraction. (Previously, in a small two-dimensional s i r n ~ l a t i o n ~ ~ using a Brownian dynamics constraints algorithm,27 we found that, except for very small clusters, aggregate rotation is of minor importance,) The two types of DLVO particle investigated in the coagulation study are identical to those examined in the sedimentation study (see fig. 1). Based on the arguments discussed already, one would expect the secondary minimum in potential A to cause pairs to associate loosely at separations r == 2.3~1, before thermal motion either induces dissociation or causes irreversible aggregation into the primary minimum ( r + 2a) after jumping over the primary maximum.We shall assume that the scalar translational diffusion coefficient D, of an aggregate of m particles is given by D,= (kT/6rr7a)mY where the exponent y is taken to be -0.54 as estimated by Meakin et ~ 1 . ~ ’ from a cluster-cluster aggregation model incorporating Kirkwood-Riseman theory at the level of Rotne-Prager hydrodynamics. The use of eqn (10) enables us to move to larger systems and longer simulation times without losing the essential physics of the problem. It does, however, make an a priori assumption about the fractal nature of the aggregates to be formed, since D, scales as m-”” for aggregates of fractal dimension D.29 Nevertheless, we do not expect the structural conclusions drawn below to be significantly affected by the choice of y ; it has been shown elsewhere3’ that, so long as D, is indeed a decreasing function of m (as must be the case in reality), it matters little to the final long-range aggregate structure what values the diffusion coefficient takes.Coagulation has been simulated for potentials A and B at each of the volume fractions 4 = 0.05,0.10,0.15,0.20,0.25 and 0.30. The time-step was varied between 2.5 and 20 ps depending on the conditions. Initially, particles are positioned at random, but subject to the condition r > 2.4a for each pair in the system.A simulation is run until all particles are incorporated into a single 5 12-particle cluster. Averages are taken over five separate runs at each volume fraction. Although we are not primarily interested in the coagulation kinetics here, some useful insight into mechanism can be derived by first considering the set of coagulation times listed in table 2. Quantities tabulated are the half-life t l / * (time to reduce number of particles by 50% ) and the total coagulation time t , (time to produce a final 5 12-particle aggregate). Owing to the neglect of inter-aggregate hydrodynamic interactions, and174 Brownian-dynamics Simulation Table 3. Number distribution n ( m ) of m-particle aggregates after one half-life for potentials A and B at 4 = 0.1 and 4 = 0.3 as compared with Smoluchowski theory; mmax denotes largest simulated aggregate; N ( m 3 8) denotes number of particles existing in aggregates with m 3 8 4 =0.1 4 = 0.3 theory m 1 2 3 4 5 6 7 2 8 mmax N ( m 2 8 ) 147 * 5 54*4 22k3 15*2 6*3 4*2 3*2 6*2 15 5 2 i 12 139*4 54*7 28*4 16*4 5*2 3*2 4*2 14 32* 14 8*2 140*5 47*6 23 *3 8*3 12*4 4*2 3*2 7*2 26 70* 18 147 * 5 49*5 25*4 13*2 8*2 4 5 2 4*2 3*2 15 54* 15 128 64 32 16 8 4 2 2 18 uncertainty about the most appropriate choice of starting configuration, the absolute values of t l l Z and t , are to be treated with caution, but their relative values contain useful information.We note that t l , Z is much more strongly dependent on the form of DLVO potential than is t , , and this is especially true at low volume fraction, where t , appears to become independent of the potential.The reason for this is fairly clearly understood:26 at short times the rate is controlled by the potential-energy barrier to coagulation (‘reaction-limited’), whereas at long times the rate-determining factor is the speed at which the (large) aggregates move around (‘diffusion-limited’). Despite the wide variation in values of t , / 2 , table 3 shows that the distribution of aggregate sizes is relatively insensitive to volume fraction and pair potential, although there does seem to be a definite tendency towards more particles in larger aggregates ( m 2 8) with potential A and at the higher volume fraction. This trend becomes stronger as the coagulation proceeds: after two half-lives at 4 = 0.3, the largest aggregate in the system contains more than twice as many particles with potential A as with potential B (52, 44, 52, 50 and 36 for the five separate runs with A, as compared with 12, 16, 20, 18 and 18 for B).In every case, simulated distributions were broader than calculated from Smoluchowski theory (see table 3), which is expected to hold only as 4 -+ 0. Deviations from theory are of a similar order. to those reported3’ for aggregating polystyrene latex particles. For each set of simulation conditions, effective fractal dimensions were calculated from plots of log m against log RG( m ) , where RG is the radius of gyration of an aggregate of m particles ( R G - rn””). Values of D are collected in fig. 3 as a function of volume fraction.The data show a slight tendency for type A particles to give more compact structures than type B particles, although the difference is hardly significant compared with the statistical uncertainty (k0.05). It is difficult to make a reliable extrapolation to 4 = 0, but the limiting value of D at infinite dilution does seem to lie below that given by the reaction-limited clustering of clusters Fig. 4 shows the effect of colloidal potential on g ( r ) for the final 512-particle aggregates generated at 4 = 0.15. The plot for potential A has a strong peak at r=2.3a corresponding to non-bonded pairs in the secondary minimum, and a broader peak at r =: 4a corresponding to the ‘second shell’ of liquid-like structure characteristic of stable concentrated dispersion^.'^ (In aggregates formed at 4 = 0.3, the ‘third shell’ is clearly evident.) The plot for potential B, on the other hand, has little short-range structure, in qualitative agreement with theG. C.Ansell and E. Dickinson 2.4 D 2.0 175 - / 7’ I I I I 0 0.1 0.2 0.3 d Fig. 3. Effective fractal dimensions D for simulated aggregates at various volume fractions #: 0, potential A; e, potential B. W 2 4 6 8 Fig. 4. Influence of colloidal pair potential on short-range structure of aggregates generated at # = 0.15. The normalized pair distribution function g( r ) is plotted against reduced separation r / a for ( a ) potential A and ( b ) potential B.176 Brownian-dynamics Simulation Table 4. Number Np of final percolating clusters in the simulations for various values of the shell parameter u [ Np (maximum) = 51 4 potential Np (u = 2 .5 4 N p ( u = 3 a ) Np ( u = 3 . 5 ~ ) Np ( D = 4a) ~ ~~ 0.05 A B 0.10 A B 0.15 A B 0.20 A B 0.25 A B 0.30 A B 2 0 4 3 4 5 0 0 0 0 4 5 5 5 4 5 5 5 5 5 earlier sedimentation data. Indeed, it is instructive to compare the short-range sediment structures with the structures of the final aggregates simulated here at the same overall volume fraction. The latter g ( r ) plots are the dashed curves in fig. 2 ( a ) (potential A, 4 = 0.20) and fig. 2 ( b ) (potential B, 4 = 0.10). While agreement is reasonably good, there are clearly differences outside the statistical error. This suggests that a sediment and a gel, although of the same overall volume fraction, may have rather different micro-structures, and therefore different mechanical properties.In part, this represents the difference between particle-cluster aggregation (sediment) and cluster-cluster aggre- gation (gel). A system is said to have gelled when aggregates of the order of the system size have been produced. This will occur if the original particle volume fraction exceeds some critical percolation threshold &. Computer simulation is a blunt instrument for estimat- ing 4,: in small systems with periodic boundary conditions, there are large fluctuations in structure, which means that percolating configurations may be observed both above and below the true percolation threshold. As pointed out recently by Seaton and Glandt,34 it is possible for the infinite system of replicas to be spanned, but not the small system itself, and vice versa.Nevertheless, even in this preliminary study, it still seems worthwhile trying roughly to estimate 4p, if only to give a handle on the approximate size of the parameter 6 in eqn (1). Table 4 records the number of percolating 512-particle clusters produced in the various simulation runs. Percolation is registered as occurring when the cluster in the basic cell ‘connects’ to its nearest image clusters in each of the x, y and z directions, a ‘connection’ between two particles being deemed to exist if r < a, where a is a conductive shell parameter.” Numbers in table 4 indicate that the final clusters generated at C$ = 0.05 and 0.10 are below the gelation threshold; and, while N , is clearly sensitive to a, most of those generated at 4 2 0.15 are above the gelation threshold.The values of Np at 4 = 0.15 suggest that the form of the colloidal potential has an influence on the gel point, but it would be imprudent to make any generalizations based on these few preliminary pieces of information. We note, however, that a value of +,=0.15 in the simulation corresponds to a characteristic network length of == 25a. This means that the value of D calculated from the simulations at C$ = 0.20, 0.25 and 0.30 are based partly on data in the sol-gel crossover region, and should not therefore be interpreted as true fractal dimensions in the self-similar sense.G. C. Ansell and E. Dickinson 177 One factor which probably distinguishes the model calculations reported here from many real colloidal systems is the absence from the simulations of any allowance for consolidation or relaxation of sediment or gel structure following initial sticking together of the constituent particles.Such considerations will have to be included in the next generation of models if Brownian dynamics simulations are to fulfil their potential as a predictive tool in colloid and material science. E.D. acknowledges support from S.E.R.C. G.C.A. acknowledges receipt of a Studentship from the Ministry of Agriculture, Fisheries and Food. References 1 E. Dickinson, J. Colloid Interface Sci., 1987, 118, 286. 2 D. Giles and A. Lips, J. Chem. Soc., Faraday Trans. I , 1978, 74, 733. 3 P. Dimon, S. K. Sinha, D. A. Weitz, C. R. Safinya, G. S. Smith, W. A. Varady and H.M. Lindsay, 4 D. W. Schaefer and K. D. Keefer, in Fractals in Physics, ed. L. Pietronero and E. Tosatti (North-Holland, 5 D. L. Ermak and J. A. McCammon, J. Chem. Phys., 1978, 69, 1352. 6 E. Dickinson, S. A. Allison and J. A. McCammon, J. Chem. Soc., Faraday Trans. 2, 1985, 81, 591. 7 G. C. Ansell, E. Dickinson and M. Ludvigsen, J. Chem. Soc., Furaduy Trans. 2, 1985, 81, 1269. 8 A. B. Glendinning and W. B. Russel, J. Colloid Interface Sci., 1982, 89, 124. 9 J. Bacon, E. Dickinson and R. Parker, Faraday Discuss. Chem. Soc., 1983, 76, 165. Phys. Rev. Lett., 1986, 57, 595. Amsterdam, 1986), p. 39. 10 G. K. Batchelor, J. Fluid Mech., 1976, 74, 1. 11 R. Schmitz and B. U. Felderhof, Physica A, 1983, 116, 163. 12 J. Rotne and S. Prager, J. Chem. Phys., 1969, 50, 4831. 13 P. Mazur and W. van Saarloos, Physica A, 1982, 115, 21. 14 G. J. Kynch, J. Fluid Mech., 1959, 5, 193. 15 J. Happel and H. Brenner, Low Re-vnolds Number Hydrodynamics (Noordhoff, Leiden, 2nd edn, 1973), 16 W. van Megen and I. Snook, Faraday Discuss. Chem. Soc., 1983,76, 151. 17 C. W. J. Beenakker, Faraday Discuss. Chem. Soc., 1983, 76, 240. 18 P. Meakin, J. Colloid Interface Sci., 1984, 102, 491. 19 J. Bacon, E. Dickinson, R. Parker, N. Anastasiou and M. Lal, J. Chem. Soc., Faraday Trans. 2, 1983, 20 U. Geigenmuller and P. Mazur, Physica A, 1986, 138, 269. 21 G. C. Ansell, Ph.D. Thesis (University of Leeds, 1986). 22 C. W. J. Beenakker and P. Mazur, Physica A, 1984, 126, 349. 23 G. C. Ansell and E. Dickinson, J. Chem. Phys., 1986, 85, 4079. 24 M. J. Vold, J. Colloid Sci., 1959, 14, 168. 25 E. Dickinson and R. Parker, J. Colloid Interface Sci., 1984, 97, 220. 26 G. C. Ansell and E. Dickinson, Chem. Phys. Lett., 1985, 122, 594. 27 S. A. Allison and J. A. McCammon, Biopolymers, 1984, 23, 167. 28 P. Meakin, Z-Y. Chen and J. M. Deutch, J. Chem. Phys., 1985, 82, 3786. 29 W. Hess, H. L. Frisch and R. Klein, 2. Phys., Teil B, 1986, 64, 65. 30 P. Meakin, Phys. Rev. Lett., 1983, 51, 1119. 31 P. G. Cummins, E. J. Staples, L. G. Thompson, A. L. Smith and L. Pope, J. Colloid Interface Sci., 1983, 92, 189. 32 M. Kolb and R. Jullien, J. Phys. (Paris) Lett., 1984, 45, L997. 33 W. van Megen and I . Snook, Adu. Colloid Interface Sci., 1984, 21, 119. 34 N. A. Seaton and E. D. Glandt, J. Chem. Phys., 1987, 86, 4668. 35 S. A. Safran, I . Webman and G. S. Crest, Phys. Rev. A, 1985, 32, 506. pp. 276 and 376. 79, 91. Received 2nd December, 1986

ISSN:0301-7249    

DOI:10.1039/DC9878300167

出版商:RSC    

年代:1987

数据来源: RSC

摘要:

Faraday Discuss. Chem. Soc., 1987,83, 179-191 Dynamics of Colloidal Particles in the Vicinity of an Interacting Surface Alec T. Clark, Moti Lal" and Gill M. Watson Unilever Research, Port Sunlight Laboratory, Bebington, Wirral L63 3JW A Brownian dynamics method has been pursued to studythe diffusive motion of a colloid particle leading to its deposition at an interacting surface. The present study considers two models for the diffusion boundary layer. In model I, the layer is bounded by a pair of absorbing surfaces, one real and the other virtual. The real surface has both static and hydrodynamic interac- tions with the particle, but the virtual surface has no such interactions, serving merely as an absorbing boundary. In model 11, the interaction characteristics of the real surface are the same as in model I, but the virtual surface is a reflecting boundary.Calculation of the mean times and prob- abilities of the particle absorption has been carried out for a number of particle-surface potential-energy functions. The results are in excellent agreement with those derived from the application of the first-passage-time approach to the two models. The computed data are used to obtain values of the rate constant for the process of particle deposition. For a given potential function, the rate constant is found not to differ greatly for the two models, implying that the results are almost model-independent. The increase in the depth and range of the potential function is found to result in significant changes in the rate constant. The dynamics of colloidal particles in the liquid/ solid interfacial region constitutes an important fundamental process which underlies the kinetics of a variety of phenomena such as particle adhesion, sedimentation and detergency.The basis of early theoretical developments of this subject resided in the essential assumption of the Nernst theory that under moderate hydrodynamic conditions, i. e. in the low-Reynolds-number regime, the solid surface is separated from the bulk fluid by a stationary layer which acts as a medium for the transport of particles to and from the surface, the transport occurring by the diffusive mechanism.' Nernst's assumption has been vindicated by subsequent theoretical work of Levich which showed that although at high fluid velocities the convective component of the transport would be appreciable, it is possible to formulate the problem such that the total transport can be attributed to diffusion alone, with the thickness of the diffusion (stationary) layer given by the hydrodynamic condition prevailing in the fluid.2 Levich established quantitative relationships between the station- ary-layer thickness and mean fluid velocity for various hydrodynamic conditions.The theory of Levich has met with considerable success in the interpretation of kinetic results for phenomena involving molecular transport, e.g. kinetics of dissolution of solids and electrochemical deposition.' For processes involving colloidal transport, however, the applicability of the theory is somewhat limited, since it takes little account of the particle-surface hydrodynamic and finite-range static interactions which are the essential features of colloidal systems, but are either weak or absent in molecular systems.Inclusion of interactions in the model led to analytical solutions with severe restrictions on the domain of their validity.334 More recent attempts have pursued the numerical approach to the solution of the transport equation, and a good measure of success has been achieved in obtaining solutions for dilute systems under the steady-state ~ o n d i t i o n . ~ 179180 Dynamics of Colloidal Particles Fig. 1. Schematic representation of a colloid particle residing in the vicinity of a plane surface. The origin, 0, of the coordinate system is located at the surface, which is coincident with the xy plane.d denotes the vertical distance between the centre of the particle and the surface and s is the separation between the plane surface and the particle surface. The past few years have witnessed the emergence of an entirely different computa- tional approach for studying the dynamics of particles immersed in hydrodynamic media. This approach is analogous to the well known method of molecular dynamics, used in the simulation of the dynamical behaviour of molecular systems, but differs from it in that instead of using the Newtonian equation of motion the method follows the Langevin equation for determining the time evolution of the motion of the particles. The original development of the method was due essentially to Ermak and McCammon.6 Sub- sequently, one of the authors in collaboration with Dickinson and coworkers pursued this method to study the kinetics of dissociation of two- and three-particle flocs under- going Brownian motion.’ Computations have since been extended by the latter workers to systems of sedimenting particles8 and to colloidal aggregates subjected to shear forces.’ In this paper, we present the application of the method to simulate the dynamics of a colloid particle under the influence of an interacting surface.Details of the Model for Particle-Surface Interactions We consider a rigid, spherical particle of colloidal dimensions located in the vicinity of a rigid plane wall of infinite dimensions, shown schematically in fig. 1. The semi-infinite volume above the wall is occupied by a Newtonian liquid of viscosity 7).It is assumed that, by virtue of its infinite mass, the wall remains absolutely stationary in the presence of forces of finite magnitude. So the particle-wall relative motion is constituted entirely of the motion of the particle in the hydrodynamic medium, which is governed both by the time-dependent hydrodynamic interactions and the static forces existing between the particle and the surface. The Particle-Surface Hydrodynamic Interaction A vital manifestitation of the hydrodynamic interactions between two or more bodies resides in the changes that are brought about in the various components of the diffusionA. T. Clark, M. La1 and G. M. Watson 181 tensor." For an isolated particle of spherical shape, the diffusion tensor is given by the well known Stokes-Einstein relation where a is the particle radius; k is the Boltzmann constant and T the temperature.The tensor To is identified with the unit tensor I: To=I= 0 1 0 . (1 :: :i Thus an element D, of D is succinctly expressible as where 8, is the Kronecker delta function: 6, = 1 for i = j , S, = 0 for i # j . If the particle is located in the vicinity of a wall, the tensor appearing on the right-hand side of eqn ( 1 ) is no longer I but is modified, owing to the particle-surface interaction, to 1/A, 0 0 .=( ; l/;yy l/;zJ. (2) The parameters A,, and A, account for the effect of the hydrodynamic interaction on particle diffusion parallel to wall, whereas A,, corresponds to the effect perpendicular to the wall. These parameters may be represented by more obvious symbols: = Ax, = A, A, = hZZ.Evaluation of Ail in terms of the particle radius, a, and the particle-surface distance, d, was accomplished by Faxen over 60 years ago who, using the method of reflections, solved the Laplace equation under appropriate boundary conditions for the velocity field." He was able to show that Ail is expressable as a power series in ( a l d ) : 45 a 1 a 5 A d l = - 9 (z) + 1 ( z) - - (-) - - (-) + . . . . 16 d 8 d 256 d 16 d (3) More rigorous solutions of the problem were sought by O'Neil12 and by Goldman et all3 The latter workers showed that Faxen's approximation is in reasonable agreement with the exact results of O'Neil for d / a > 1.04. For d / a < 1.04, however, Faxen's theory progressively underestimates the hydrodynamic effect as the particle approaches smaller separations.Since in the present computations we are interested primarily in the motion of the particle normal to the surface, insistence on high accuracy in the values of All is not warranted. Theoretical determination of A, as a function of d seems to have been first carried out by L o r e n t ~ , ' ~ but his result is valid only for d >> a. In order to pursue a more general solution, applicable at all separations, Brennerl extended the bipolar coordinate approach to the solution of the creeping motion equation for the system, originally developed by Stimson and Jeffery" for a system of two spheres. He arrived at the following expression for AL: (4) 00 A,=:sinha 2 2 sinh (2n + 1)a + ( 2 n + 1 ) sinh 2a 4 sinh2 ( n +;?a - (2n + 1)' sinh' a - ') where a = cosh-' ( d l a ) .182 Dynamics of Colloidal Particles Fig. 2.Variation of the hydrodynamic parameters A , , ( a ) and A, ( b ) with 8, the particle-surface separation reduced in terms of a. For d / a >> 1 , Lorentz's result that is recovered by retaining only the first term in the series. In the limit of d / a approaching unity, Cox and B~enner,'~ applying the singular perturbation method, derived the following simple expression: 1 S A, = - - 0.2 In S + 0.97 with We have found that e d - a s a=--- - a a n ( 5 ) gives sufficiently accurat values of A, for s/ a < 0.1. Experimental studies of McKay et al. appear to confirm that Brenner's solution is substantially correct.18 In the present computation, we have utilised the solutions of Faxen and Brenner in calculating the variation of All and A, with d, shown graphically in fig.2. Particle-Surface Static Interactions The static interactions refer to the time-independent interactions (time average of the hydrodynamic interactions being zero), which would comprise the van der Waals and the electrical double-layer forces acting between the particle and the surface. These are traditionally expressed in terms of the DLVO potential function, which assumes that the van der Waals and the double-layer components of the interaction energy, U ( d ) , are linearly additive: W d ) = K W W + Ue(d). ( 6 ) A simple expression for the van der Waals component, which is adequate for the present purpose, derives from Hamaker's theory."A.T. Clark, M. La1 and G. M. Watson 183 where AIz3 is the Hamaker constant for particle (1) and surface (3) immersed in a medium (2). Bell, Levine and McCartney (BLMc) have carried out an extensive investigation of the application of the Debye-Huckel equation in the development of expressions for U, for a pair of particles of unequal radii and dissimilar surface potentials and $z.20 In the limit of one of the particles approaching infinite size, the BLMc treatment yields the following expression for U,: where I is the dielectric constant of the medium and E~ the permittivity of the vacuum; K is the Debye-Huckel reciprocal length parameter and H = [ a / ( a + ~ ) l ~ ~ ~ + [ ( a + s ) / a ] * ’ ’ . Eqn (8) is more general than that derived by Hogg et al.” as it is valid in a much wider range of particle-surface distance.The particle-surface force, F( d ), is simply the negative of the gradient of U, which is readily obtained by differentiating eqn (7) and (8) with respect to d. Computational Approach Theoretical considerations underlying the Ermak-McCammon method have been dis- cussed in depth by the original authors6 A more recent analysis by Tough et al.** established the consistency of the method not only with Deutch and Oppenheim’s2’ form of the Langevin equation but also with the Smoluchowski equation. Here we confine our attention to those operational aspects of the method which apply to the problem under consideration. The computational viability of the particle dynamics simulation follows from the result that in a time interval ( A t ) , much longer than the relaxation time ( T ~ ) for Brownian fluctuations, the motion of particles in a hydrodynamic medium can be written as a sum of three terms: D kT Ar = V .DAt+- F A t + R (9) where A r represents the displacement vector for the system during the time interval At. The first term contributing to A r is called the gradient term, which is caused by the interparticle distance-dependent changes in D. The second term is the drift term, arising from the presence of systematic forces ( F ) in the system. The third term, the stochastic term, is the Brownian motion contribution to Ar. Eqn (9) is readily recast to express the components of Ar in terms of those of D and F: where Arl is the ith component of Ar. The index j in the summations covers all the components of the displacement vector for the system; e.g.for a system of two particles j will vary from 1 to 6. Although the particle at a surface is a two-body system, in respect of its dynamical behaviour it is equivalent to a system of an isolated particle with a diffusion tensor containing only the self-terms which are no longer independent of the position of the particle in space. This means that the summations in eqn (10) will each include just three terms.184 Dynamics of Colloidal Particles The Gradient Term It is well known that for a pair of particles the gradient of the diffusion tensor is identically zero for Oseen ( first-order) and Rotne-Prager (third-order) approximations for the mobility Therefore, the gradient term has been justifiably neglected in the computations on particle cluster di~sociation.~ In the present case Thus, although contribution of the gradient term to the components of motion parallel to the surface is zero, the component perpendicular to the surface will include a non-zero contribution arising from the variation of the hydrodynamic parameter A, with respect to the particle-surface distance.In calculating this contribution, aA,/dd was obtained by differentiating eqn (4) and ( 5 ) and computing the data from the resulting expressions. The Drift Term In the present model, the sole origin of the drift term is the force deriving from the DLVO potential governing the particle-surface interaction. The force has one com- ponent, in the z-direction.Thus the components of the drift term parallel to the surface are zero. From the distance derivative of the potential the force was calculated at each time step in the main simulation program. The Stochastic Term In computing the stochastic term, it is assumed in common with previous s t ~ d i e s ~ - ~ that to first-order in At, the presence of the gradient and the drift terms does not affect the Gaussian behaviour of this component of the particle motion. A theoretical vindication for this assumption has emerged in the analysis of Tough et a1.22 Thus, the calculation of this term simply involves the generation of a sequence of normally distributed numbers whose mean is zero and standard deviation is equal to unity. The value of Ri is then given as where Xi is a number randomly selected from the sequence, and oii = (2AtDii)"2 where Dii is the diffusion coefficient in direction i.The Time Step The length of the time step assumed in the simulations is limited by two constraints: ( a ) A t should greatly exceed the relaxation time ( T ~ ) for Brownian fluctuations, which is of the order of m / 6 7 ~ q a (for a particle of 0.2 p m radius, TB = lo-' s) and (6) A? should not be too large to violate seriously the assumption inherent in eqn (9) that the diffusion coefficient and the force remain constant during the step length. A safe lower limit on the value of A? would be ca. 1 0 ~ ~ . The upper limit would, of course, dependA. T. Clark, M. La1 and G. M. Watson 185 real s u r f a c e virtual surface (absorbing ) (absorbing) I I I I I I @ I particle 1 I I 1 I real surface virtual surface (absorbing) ( reflecting) I I I I particle I I I I I I I I @ I ( b ) Fig.3. Schematic representation of the models assumed for the interfacial region. ( a ) Model I, ( b ) model 11. on how steep the force and the diffusion coefficient gradients are. In our previous study on particles of 0.5 p m radius, At varied between and s . ~ Results and Discussion As implied earlier, the objective of this work is to investigate the effectiveness of the Brownian simulation method in studying the process of particle deposition at a surface. To this end, we consider two models. The interfacial region is bounded by a pair of surfaces, one real and the other virtual. The real surface has both hydrodynamic and static interactions with the particle, but the virtual surface has no such interactions, serving merely as an absorbing boundary.This model represents a particle in a solid/liquid interfacial region which may diffuse towards the solid (real) surface resulting in its deposition, but which also has the freedom to move away from the interfacial region and be lost forever in the bulk liquid (ie. absorbed at the virtual boundary). Model I1 In this model, the interaction characteristics of the real surface are the same as in model I, but the virtual surface is a reflecting boundary. This model corresponds to an interfacial region in which the particle is irreversibly trapped and therefore is bound to deposit on the (real) surface in finite time.The two models are shown schematically in fig. 3. Our reason for selecting these rather simple models lies largely in their amenability to analytical solution leading to exact results which can be compared with the simulation data. In this way the reliability of the method in obtaining accurate information on deposition kinetics can be evaluated. This is a necessary step towards the application of the method to complex systems, e.g. a cluster of particles approaching a surface, where analytical solutions are not available. We study the motion of a particle of 0.215 p m radius in the interfacial region of width equal to 4a, the space in the region being occupied by a hydrodynamic medium (an electrolyte solution) of 7 = lop3 Pa s-'. The particle and the interacting surface are assumed to be oppositely charged; thus the particle diffuses in the presence of an attractive potential whose range is governed by the level of electrolyte in the medium.186 Dynamics of Colloidal Particles 0 20 40 k -.5 I 60 80 100 s 0.00 0.05 0.10 0.15 0.20 0.25 Fig. 4. Potentials of particle-surface interaction calculated from eqn (6)-(8) with a = 0.215 pm, AlZ3 = kT, E ~ = 80, I,!J1 = 27 mV, I,!J2 = -7 mV and KU = 510(1), 160(II), 50(III), 30(IV) and 16(V). The computations are performed for five potentials constructed from eqn (7) and (8) with values of the various parameters as follows: A,23= kT; ~ , = 8 0 ; (particle) = 27 mV; +* (surface) = -7mV; Ka = 510 (potential function I), 160 (potential function I I ) , 50 (potential function III), 30 (potential function IV) and 16 (potential function V).These potential functions are presented graphically in fig. 4, where the potential energy, U, is plotted as a function of 6. In obtaining these potentials, it is assumed that at surface-to-surface separations <lop3 p m a steep repulsion is produced owing to overlap of the Stern layer. This effect is introduced in the potentials by imposing a vertical barrier of infinite height at 6 = 5 x giving rise to a minimum in the potential function at this distance. A program was developed in Fortran to compute the components of the displacement vector for the particle from eqn (10). Values of D and dD/dz, discussed in the preceding section, were stored as functions of d in the computer in tabular form, and were called when required in the computation.For values of d not existing in the tables, the two quantities were obtained by interpolation between the nearest lower and upper values of d featuring in the tabulated data. The sequence of normally distributed random numbers, { X } , required in the computation of the stochastic term, was generated using the NAG subroutine G O ~ D D F ~ ’ which uses the method of Brent.26 In all the computations performed, the length of the time step ( A t ) assumed was 10-6s. This is two orders of magnitude greater than the value of T~ for the particle. The calculations were carried out on a VAX 785 computer. Starting from the initial position at the centre of the interfacial region ( d = 2a), the particle was allowed to move about in the region in discrete steps computed according to the above procedure.If in the course of its motion the particle touched or crossed either surface in model I, or the real surface in model 11, it was assumed to be absorbed (the particle was precluded from crossing the reflecting boundary in model I1 by stipulating that the value of d would not exceed 3a). The time corresponding to this event was recorded and whether the particle was absorbed by the real or the virtualA. T. Clark, hf. La1 and G. M. Watson 187 (absorbing) boundary. The particle was brought back to its initial position when it started to diffuse again, the motion continuing until it again reached the state of absorption. The computation was carried out for (7-34) x lo6 time steps, during which time substantial numbers of absorption events had taken place in most runs allowing mean times of absorption to be obtained.Other quantities of interest include the relative probabilities of absorption by the real and the virtual boundary in model I. The computed data are presented in tables 1 (model I) and 2 (model 11) for all the five potentials considered as well as for the case where there exists no static interactions between the particle and the real surface. For the system with no static interactions, the number of absorption events on the real surface was not large enough to obtain a good estimate of errors in the computed mean values. In the presence of a long-range attractive interaction, however, the number of absorptions achieved was sufficient to divide the data-into sub-samples and thus to carry out the analysis for the magnitude of the error involved in the estimate.On the basis of this analysis, we believe that the uncertainty in our results ranges from ca. 5% for potential V to ca. 10% for potential 1. Mean times corresponding to the events brought about essentially by diffusive motion, e.g. a chemical reaction in a liquid medium, collision of colloidal particles in a sol etc., are known in the theory of stochastic processes as the first passage or survival times (FPT). The FPT (7) for an event can be derived from the Fockker-Planck equation either by directly integrating it2’ or by pursuing the method of o p e r a t o r ~ , ~ ~ ’ ~ ~ the two approaches yielding identical results. The FPT approach has been considered in the treatment of a variety of problems including protein folding,30 geminate rec~mbination~~ and dissociation of colloidal aggregates3’ For a Brownian particle confined between a pair of boundaries, r is expressible as: where p ( z , tlz,,) is the probability that the particle starting at a position zo assumes a position z at time t.The integration of p ( z , tlz,) between the limits B, and BZ, which define the positions of the two boundaries, gives the probability that the particle at time t exists in the region between the boundaries (the particle will cease to exist when it is absorbed by a boundary). The integration of this probability with respect to t will represent the mean time that the particle will spend in the region before it is absorbed.In model I, both boundaries are absorbing and the particle is subjected to a potential field arising from its interaction with the real surface. For this model, it may be shown that:32 I: ?=; [ I: A, exp ( U / k T ) dz exp ( - U / k T ) d l - ( [I A, exp ( U / k T ) dz exp ( - U / k T ) dE I:‘ where 5 is the weighted mean value of the first passage time for the absorption of the particle on both boundaries: p r and pv are, respectively, the probabilities of absorption on the real and the virtualTable 1. Mean times and probabilities of particle absorption calculated by the method of Brownian dynamics (BD) and by the first-passage-time approach (FPT)-Model I Tr/ s Pr T"/ s P V 7/ S potential function BD FPT BD FPT BD FPT BD FPT BD FPT no static forces 0.1 1 0.102 0.28 0.212 0.085 0.072 0.72 0.788 0.091 0.079 I 0.060 0.059 0.39 0.395 0.043 0.048 0.61 0.605 0.050 0.052 I1 0.060 0.059 0.40 0.395 0.045 0.048 0.60 0.605 0.051 0.052 111 0.051 0.054 0.46 0.408 0.047 0.047 0.54 0.592 0.049 0.050 IV 0.048 0.047 0.43 0.433 0.042 0.043 0.57 0.567 0.045 0.044 V 0.033 0.035 0.5 1 0.498 0.032 0.037 0.49 0.502 0.033 0.036A.T. Clark, M. La1 and G. M. Watson 189 Table 2. Mean times of particle absorption at the real surface-model I1 potential function BD FPT no static forces 0.60 0.528 I 0.21 0.197 I1 0.21 0.197 I11 0.17 0.182 IV 0.14 0.159 V 0.12 0.115 boundary and rr and are the respective first passage times: P”= 1 -Pr* The first passage time for absorption on the real surface has been obtained as x 1; A, exp ( U / k T ) dz exp ( - U / k T ) dl r x 1: A, exp ( U / k T ) drn I/.1; A, exp ( U / k T ) dz. A similar expression holds for 7,. Solution of eqn (15)-(19) was pursued by evaluating the various integrals involved with potential function U given by eqn (7) and (8) and A, by eqn (4) and (5). This was accomplished to a good degree of accuracy by following numerical procedures involving the NAG subroutines DOI A H F ~ ~ and DOI D A F , ~ ~ which are designed to perform single and double integrations, respectively, employing the method of Patterson.35 The triple integrals appearing in eqn (19) were evaluated by calling sequentially the two subroutines in the program. The calculated values for the various quantities are presented in table 1, which also contains the Brownian dynamics results.The FPT values quoted are accurate to the third decimal place. For model 11, since the virtual surface is a reflecting boundary, particle absorption will occur only at the real surface. In this case Deutch’s solution applies27 and so we have: exp ( U / k T ) dz exp (- U / k T ) dl. rr values calculated according to ecp (20) are shown in table 2 together with our Brownian dynamics results on this model. Within the limits of accuracy of the present computation, we find that for both models the values of mean times and probabilities for particle capture, obtained by the Brownian dynamics method, are in excellent agreement with the first-passage-time190 Dynamics of Colloidal Particles Table 3. Values of the rate constant calculated from Brownian dynamics and first-passage-time data on T, and pr d,/s-' (model I) k,/s-' (model 11) potential function BD FPT BD FPT no static interactions 2.5 2.08 1.7 1.89 I 6.5 6.69 4.8 5.08 I1 6.7 6.69 4.8 5.08 I11 9.0 7.55 5.9 5.49 IV 9.0 9.21 7.1 6.29 V 15.5 14.22 8.3 8.70 calculations.Such a good agreement with the theory greatly enhances our confidence in the capability of the method in treating the more complex models mentioned earlier. However, in order to produce accurate results on such models, it would be necessary to generate sequences of elementary displacements perhaps an order of magnitude longer than the present samples. This is within easy reach of computers of even moderate power and speed. We should like to note that in the present study the inclusion of the gradient term in the computation produced results which were found to be distinctly in better agreement with FPT than those which did not include this term.Considering that the contribution of the gradient term in a single step is small in comparison to the Brownian term (the latter being, typically, two orders of magnitude greater than the former), it would appear at first sight that the neglect of this term would not introduce any appreciable error. The difference observed between the two sets of results (with and without the gradient term) are attributable to the systematic effect which this term confers on the particle motion in small increments, cumulating over long periods of time. Absorption on the real surface in the present models is an idealised representation of the process of particle deposition which is experimentally known to follow first-order rate kinetics.For such a process the rate constant is simply the reciprocal of the characteristic time corresponding to the rate-determining step. Assuming that the rate- determining step in the deposition process is the diffusion of the particle through the interfacial region, we may write where kd is the rate constant for the process. In calculating the rate constant from model I, we must take account of the fact that there exists a finite probability for the particle escaping irreversibly from the interfacial region, which corresponds to the particle absorbing on the virtual boundary. Thus pv represents the probability that the process has a zero rate constant and pr the probability that the rate constant is l / ~ ~ .Thus the probability-weighted rate constant is P r kd = -. 7, For model 11, since the probability of particle escape is zero, the rate constant is simply given by eqn (21). The kd values for the two models, calculated with the various potential functions, are presented in table 3. For both the models, the rate constants obtained from Brownian dynamics results are found to be in good agreement with those derived from the first passage times. It is interesting to note that for a given potential function the rate constants for the two models are similar in magnitude, differing at the most by a factor of ca. 1.9. The increase in the rate constant bought about by an increaseA. T. Clark, M.La1 and G. M. Watson 191 in the depth and range of the potential function, although significant, is not dramatic. For example, for systems with potential function V, which has a depth of ca. 350kT and a range of ca. 0.5a, the rate constant is less than one order of magnitude greater than for the system with no static interactions. Concluding Remarks The present study strongly affirms the effectiveness of the Brownian dynamics method as a means for studying colloidal dynamical processes at surfaces, which underlie such important phenomena as particle deposition, removal and transport in porous media. The approach seems particularly suitable for investigating the effects of particle dimensions and interactions, among other factors, on the colloid dynamics in the interfacial region.The excellent agreement with the first passage time theory, achieved for calculations on simple models, encourages one to undertake the study of more complex problems such as the significance of blocking effects in particle deposition, and the influence of an interacting boundary on the kinetic stability of colloidal aggre- gates. The authors gratefully acknowledge helpful discussions with Drs R. B. Jones (Queen Mary College, London) and E. Dickinson (Leeds University). References 1 W. Nernst, 2. Phys. Chem., 1904, 47, 55; W. Nernst and E. S. Merriam, 2. Phys. Chem., 1905, 53, 235. 2 V. G. Levich, Physicochemical Hydrodynamics (Prentice-Hall, Englewood Cliffs, N. J., 1962). 3 D. C. Prieve and E. Ruckenstein, J. Colloid Interface Sci., 1976, 57, 547.4 K. Chari and R. Rajagopalan, J. Chem. Soc., Faraday Trans. 2, 1985, 81, 1345. 5 Z. Adamczyk and T. G. M. van de Ven, J. Colloid Interface Sci., 1981, 84, 497; 6 D. L. Ermak and J. A. McCammon, J. Chem. Phys., 1978, 69, 1352. 7 J. Bacon, E. Dickinson, R. Parker, N. Anastosiou and M. Lal, J. Chem. Soc., Faraday Trans. 2, 1983, 79, 8 E. Dickinson and R. Parker, J. Colloid Interface Sci., 1984, 97, 220. 9 G. C. Ansell, E. Dickinson and M. Ludvigsen, J. Chem. SOC., Faraday Trans. 2, 1985, 81, 1269. 91. 10 G. K. Batchelor, J. Fluid Mech., 1976, 74, 1. 11 H. Faxen, Arkiu. Mat. Astron. Fys., 1923, 17, no.27. 12 M. E. O’Neill, Mathematika, 1964, 11, 67. 13 A. J. Goldman, R. G. Cox and H. Brenner, Chem. Eng. Sci., 1967, 22, 647. 14 H. A. Lorentz, Abhandl, Theoret. Phys., 1906, 1, 23. 15 H. Brenner, Chem. Eng. Sci., 1961,16,242; J. Happel and H. Brenner, Low Reynolds Number Hydrodynamics (Noordhoff Publishing Co., Leyden, 1973), chap. 7. 16 M. Stimson and G. B. Jeffrey, Proc. R. SOC. London, Ser. A , 1926, 111, 110. 17 R. G. Cox and H. Brenner, Chem. Eng. Sci., 1967, 22, 1753. 18 G. D. M. MacKay, M. Suzuki and S. G. Mason, J. Colloid Sci., 1963, 18, 103. 19 H. C. Hamaker, Physica, 1937, 4, 1058. 20 G. M. Bell, S. Levine and L. N. McCartney, J. Colloid Interface Sci., 1970, 33, 335. 21 R. Hogg, T. W. Healy and D. W. Fuerstenau, Trans. Faraday Soc., 1966, 62, 1638. 22 R. J. A. Tough, P. N. Pusey, H. N. W. Lekkerkerker and C. van den Broeck, Mol. Phys., 1986, 59, 595. 23 J. M. Deutch and I. Oppenheim, J. Chem. Phys., 1971, 54, 3547. 24 R. Parker, Ph. D. Thesis (University of Leeds, 1984). 25 NAGFLIR: 1451/0, Mk6, May, 1977. 26 R. P. Brent, Commun. ACM, 1974, p.704. 27 J. M. Deutch, J. Chem. Phys., 1980, 73, 4700. 28 K. Schulten, Z. Schulten and A. Szabo, J. Chem. Phys., 1981, 74, 4426. 29 N. Agmon, J. Chem. Phys., 1985, 82, 2056. 30 M. Karplus and D. L. Weaver, Riopolymers, 1979, 18, 1421. 31 D. Y. C. Chan and B. Halle, J. Colloid Intecface Sci., 1984, 102, 400. 32 N. S. Goel and N. Richter-Dyn, Stochastic Models in Biology (Academic Press, New York, 1974), chap. 3. 33 NAGFLIB: 1745/0, Mk 8, 1981. 34 NAGFLIB: 830/0, Mk 11, 1983. 35 T. N. L. Patterson, Math. Cornput., 1968, 22, 847. Received 3rd March, 1987

ISSN:0301-7249    

DOI:10.1039/DC9878300179

出版商:RSC    

年代:1987

数据来源: RSC

摘要:

Faraday Discuss. Chem. SOC., 1987, 83, 193-198 Boundary Conditions on Hydrodynamics in Simulations of Dense Suspensions E. R. Smith Mathematics Department, La Trobe University, Bundoora, Victoria 3083, Australia The problem of the hydrodynamic boundary conditions on a simulation sample of a hydrodynamically dense colloid in Brownian dynamics is described and the necessity of a periodic boundary condition shown. The lattice sum for the leading-order pairwise hydrodynamic mobility tensor in periodic boundary conditions is analysed and shown to be associated with a quite unphysical model periodic system in which infinite suspending fluid velocities certainly occur. A more complicated boundary condition is intro- duced, in which a large spherically shaped array of cubic simulation sample copies is surrounded by a stationary spherical container with stick boundary condition.Analysis of this boundary condition shows that the resulting simulation sample is physically consistent with a properly defined positive definite hydrodynamic mobility tensor. Some effects of this definition of the simulation sample on simulation measurements are discussed. The Stochastic Dynamics Algorithm The system to be considered consists of N spherical particles all of radius a moving in a fluid of viscosity 7. There are systematic forces (due, for example, to electrostatic interactions, van der Waals interactions or steric hindrance effects) acting between the particles. The particles are confined to a cube of side L and so have density p = N / L 3 . As the particles move, they set u a velocity field in the suspending fluid.A typical interparticle distance is R = LA-''. As the particles move, they are subject to a viscous drag. The system considered is hydrodynamically dense, so that the drag on one given particle is affected by the suspending fluid velocity field set up on the given particle by the motion of the other particles. The simplest case, which is considered here, calculates the viscous drag from the quasi-static linearized Navier-Stokes equation for the suspending fluid velocity c : qv2v=vp; v * v = o ( 1 ) where p is the pressure in the suspending fluid. The particles are assumed not to rotate and the solutions to eqn ( 1 ) obey the boundary condition that the fluid velocity at a point on the surface of particle j with velocity u, is 0 = u,.The solutions to eqn (1) are developed in an asymptotic expansion in a / R and here the expansion is taken only out to order ( a / R ) 3 . In the standard Brownian dynamics algorithm, time evolution is deemed to occur in time steps of fixed length At. The particles are at q ( n ) at the end of time step n. During the ( n + 1)th time step, the particles move at constant velocity to rj ( n + 1 ) = yj ( n ) + R,S( n ) + Rf ( n ) where R.s( n ) is the displacement during the time step ( n + 1 ) due to the systematic forces acting and R;(n) is a random displacement. These random displacements are chosen from a joint Gaussian distribution with covariance matrix (R,;( n)R',( n ) ) = 6kTp.(j, k ; n ) At. (3 1 193194 Colloid Simulation Boundary Conditions Here, k is Boltzmann’s constant, T is the absolute temperature and p(j, k ; n ) is the mobility tensor for the particles evaluated when they are at q( n ) , 1 N.When the suspension is hydrodynamically dilute ( a / R << 1) this tensor is given by the Stoke’s Law diagonal form 67rr]ap(j, k ; n ) = 8j,kI, where I is the 3 x 3 unit matrix. This method is known to give reliable descriptions of the stochastic evolution of hydrodynamically j dilute suspension^.'-^ If the particles move at constant velocity uj, velocity, and so the viscous drag Di on particle j constant velocities y j . Thus N Dj = - 2 C ( j 9 k ) k = 1 1 “ - j s N, then the suspending fluid must be a linear combination of the and the matrix C(j, k ) must be evaluated from eqn (1) and the boundary conditions.Let the net systematic force (i.e. that due to non-hydrodynamic effects) on particle j be 4. Since the particles move at constant velocity, FJ + Dj = 0. This gives the forces FJ as a linear combination of the velocities uj and that linear relation may be inverted to give N u j = c p6, k , ‘ F k k = l where p = 6 - l is the mobility tensor. The solution to eqn (1) for an isolated particle at r, with velocity uo is 3a a3 4 4 u ( r ) = - U l ( r - r o ) * u o + - Q ( r - r o ) uo where 1 1 O(r)=-(I+??) and Q ( r ) = 7 ( 1 - - 3 F ? ) Irl Irl (7) are the ‘Oseen’ and ‘quadratic’ terms, respectively. Here ?= r/lrl, The solution to eqn (1) for a set of particles may then be developed by the reflection principle as an expansion in powers of ( a / 15, k I) with q, k = rj - rk.If this expansion be carried out to order ( a / (rJ, 1)‘ and the drag calculated, eqn ( 5 ) gives the celebrated Rotne-Prager hydrodynamic mobility t e n ~ o r . ~ Note that much recent work on calculating higher-order expansions for this tensor uses an induced force method which does not require explicit calculation of the suspending fluid Neverthe- less, these methods are based on eqn ( 1 ) and so make the assumption of small Reynold’s number flow or very small particle velocities which allows linearization of the Navier- Stokes equation. Boundary Conditions for Simulations A finite simulation sample is supposed to be a model of a large but microscopic fraction of a macroscopic system. What happens at the edges can be crucial in ensuring that the sample (or a subsection of it) really does behave as part of a macroscopic system. It is well known that if a hard wall is used to keep the particles at fixed density, then there will be strong layering effects into the sample.The size of an isolated sample big enough so that its interior will not feel such surface layering appears prohibitive from the point of view of available computational power. The only viable alternative is to impose some form of periodic geometry on the system.E. R. Smith 195 The simplest form of periodic geometry which can be imposed is the minimum image convention. The system is wrapped on a 'cylinder' in all three directions. The sample may be viewed as surrounded by 26 exact copies of itself making a 3 x 3 cubic array.If a particle moves through a surface of this array it re-enters at the same velocity through the opposite surface. I n t h e minimum image convention, the mobility tensors are calculated not with rJ,k but rJ,A which is the vector to rJ (in the sample) from the closest copy of particle k (the sample version of particle k being one of the 27 possibilities). This construction is known to work well with other simulations with short-range forces, but does fail with electrostatic forces. In stochastic dynamics simula- tions the construction fails catastrophically: lo-'* the Rotne-Prager mobility tensor in minimum image boundary conditions is no longer positive definite.I3 The simulation rapidly generates configurations for which this mobility tensor has negative eigenvalues and so cannot be used as a covariance matrix.The alternative is to use a full periodic- boundary condition. The cubic sample of side L is replicated periodically throughout space. Here the Rotne-Prager tensor is positive definite. The expression for the mobility tensor in PBC (periodic boundary conditions) might naively be written via eqn ( 5 ) noting that a given particlej has hydrodynamic interactions with every other particle in the period array so that where A is the simple cubic lattice of vectors with integer components and the asterisk on the sum on n means to leave out the j = k term when n = 0 and to take care of the convergence difficulties inherent in the sum. It is important to recognize that this lattice sum arises from a particular physical model.The model consists of using the quasi-static linearized Navier-Stokes equation to describe hydrodynamics in an infinite periodic array. This means that the fluid velocities must be small throughout this array, or the assumption of linearizability is physically inconsistent. Lattice Sums for Mobility Tensors Eqn (8) contains lattice sums of two terms, one of the 0 tensors and one of the Q tensors. The sum of the Q tensors is known from problems of electrostatic energies in dipolar lattices and is conditionally ~0nvergent.l~ The sum of the 0 tensors has worse convergence properties. This means that if a convergent answer is to be obtained, the order of addition must be defined and all subsequent analysis carried out using the physical picture demanded by the chosen summation order.In this work the sum on n is carried out, for some summand f ( n ) , as c f(4= c f(n) nF!\ n E A InlSS where the array is chosen as a large spherical array of all copy cells When this is substituted into eqn (8) the limit of the whole expression as (9) n with I n / s S . S -+ co is taken and an effective mobility tensor for the periodic simulation sample may be identified. Consider then the sum of the Oseen terms in eqn (8). It is important to remember that the result must hold for any copy cell within the array. Thus consider the relation for the particles in cell m. This is N196 Colloid Simulation Boundary Conditions the asterisk now meaning that when n = m, the j = k term is to be omitted. The Oseen tensors may now be approximated by their expansions in powers of r j k / L and the sums of the various terms approximated by integrals.This gives where p = m/S and J(p)= 1, O(p-x)d3x=&[(5-2p2)1+pp]. X I S 1 Eqn ( 1 1 ) displays the inherent divergences of the lattice sum, as S -+ 00. The term in S2 is in fact zero by Newton's third law, if no external forces act. The term in S' simplifies to The first term here should also be zero when there is no external torque on the system, but the second term is not reduced to zero by the constraints of zero net force and torque. Thus, in the bulk of the lattice array, where Iml= O ( S ) , the lattice sum for the mobility tensor diverges. This divergence is not removed by the requirement that the net force and torque on the sample be zero.When details of the O(So) terms are calculated,15 even the m = 0 cell has an extra piece attached to it. This piece is not zero when the net force and torque are zero. Beenakkerl' has attempted to calculate precisely the lattice sum described here and misses the O ( S ) and O(So) terms which are discussed here. This is because he has used the Poisson summation formula where it does not apply. The result is that he has missed the O( S ) and O( So) terms and has only cancelled out the O(S2) terms (using Newton's third law). His results leave out an O(So) term and are thus incorrect for the problem he poses. Since the lattice sum of Rotne-Prager mobility tensors is divergent for off-centre cells in the array, it makes sense to consider the physical picture a little more closely and to examine the suspending fluid velocity within the array. It seems simplest to consider the solution of eqn ( 1 ) in the array by the reflection principle, since that concentrates the mind of the suspending fluid velocity at each state of the reflection expansion.At a point r = Sp in the array, at the first level of reflection, the suspending fluid velocity is The O(S2) term may be removed assuming that the net velocity of the particles is zero, but again even the requirement of zero net angular velocity is insufficient to remove the O( S ) term. Thus, in the periodic array, as S + a, the suspending fluid velocity diverges at the first level of reflection. This means that there is no physical basis for using the linearized form of the Navier-Stokes equation.It may be shown" that similar divergen- ces arise at each level of the reflection principle expansion to the solution to eqn ( 1 ) in the periodic array. All is not lost, however, when it is realized that the O(S*) and O ( S ) terms in eqn (13) are in fact interior solutions of eqn ( 1 ) inside the sphere of radius S. The particular representation of these solutions required is described in ref. ( 1 5 ) . If the surface of the sphere is considered to be a fixed spherical surface with a stick boundary condition forE. R. Smith 197 the suspending fluid velocity on it, then these terms will reflect off that surface. The reflected contribution to the suspending fluid velocity exactly cancels the diverging terms in eqn (13) and this cancellation occurs at every level of the reflection principle expansion.Notice especially that this outer sphere stick condition means that the net velocity of everything inside the sphere, suspending fluid and array of particles, is zero. This is the condition required to force finiteness in the lattice sum for the fluid velocity. This boundary condition also removed the divergences from the expression for the mobility tensor in off centre cells and the corresponding O(So) term from the mobility tensor in the central cell. The details of the O(So) terms are excessively complicated and so we refer the reader to ref. (15) for them. Discussion The major point of this work is to identify the diverging terms in the lattice sums for suspending fluid velocity and mobility tensors in a periodic array of particles, and to specify a physical boundary condition on the array which removes these divergences.The resulting mobility tensor^'^ have the same numerical values as those derived by Beenakker.16 The form derived by Beenakker is rather different from that derived in this work because he started with a representation of the Oseen and quadratic tensors eqn (7) different from that used in ref (15). Beenakker has used the Poisson summation formula to derive his results. In doing this he considered a diverging term (the k = 0 term in the Fourier space sum from the Poisson formula) to be wholly proportional to the net force on all the particles and so zero by Newton’s third law. This is incorrect: there are other less rapidly diverging terms in the k = 0 contribution. The Poisson summation formula is incapable of calculating them unless radically modified.Beenakker’s system is unphysical, since it contains both linearized hydrodynamics and diverging fluid velocities. His estimate of the lattice sum for the PBC mobility tensor is incorrect as he leaves out the finite part of the k = 0 term. The physically consistent system considered in this work introduces a physical mechanism for keeping the k = 0 term finite. This is the rigid wall boundary condition at the edge of the array. In fact it suppresses the k = O term entirely. Thus Beenakker’s results are the correct ones for a system he did not study, but not correct for the system he did study. One defect of the lattice sum representations given in ref.(15) for the mobility tensor in ‘PBC with hard wall at infinity’, is that they are probably too costly in computing time to use with any but a few bench-mark simulations. Some other way must be found and work on a reaction field type hydrodynamic boundary condition is presently in progress. An important question is how much does the ‘hard wall at infinity’ condition affect the properties measured in the simulation? As yet no answer to this can be given. However, the structure of this analysis is fundamentally similar to the scalar analogue of it all which arises in the theory of simulation of systems with dipolar electrostatic interactions. There Poisson’s equation replaces eqn ( 1 ) . In the electrostatic case, the main result of defining the simulation system properly is to change completely the formula for calculating zero-frequency dielectric constants from the mean-square polariz- ation fluctuations.” Accordingly, the main place to expect marked changes in theory due to the ‘hard wall at infinity’ boundary condition is in the formula connecting diffusion constant estimates with the zero-frequency effective viscosity of the system.Work on this question is currently in progress. References 1 K. J. Gaylor, 1. K. Snook, W. van Megen and R. 0. Watts, J. Chem. SOC., Faraday Trans. 2, 1980, 76, 1067.198 Colloid Simulation Boundary Conditions 2 K. J. Gaylor, I. K. Snook, W. van Megen and R. 0. Watts, J. Phys. A, 1980, 13, 2513. 3 K. J. Gaylor, I. K. Snook and W. van Megen, J. Chem. Phys., 1981, 75, 1682. 4 C. T. Havens, Ph. D. Thesis (Case Western Reserve University, Cleveland, Ohio, U.S.A., 1979). 5 J. Rotne and S. Prager, J. Chem. Phys., 1969, 50, 4831. 6 P. Mazur and D. Bedeaux, Physica, 1974, 76, 235. 7 P. Mazur and W. van Saarloos, Physicu A, 1982, 115, 21. 8 B. U. Felderhof, Physica A, 1976, 84, 557. 9 C. W. J. Beenakker and P. Mazur, Physica A, 1985, 131, 311. 10 D. L. Ermak, personal communication to I. K. Snook, who passed it to the author. 11 I . K. Snook and W. van Megen, personal communication. 12 J. Bacon, E. Dickinson and R. Parker, Furaday Discuss. Chem. SOC., 1983, 76, 165; 173. 13 C. W. J. Beenakker, personal communication to W. van Megen and thence to author. 14 S. W. DeLeeuw, J. W. Perram and E. R. Smith, Proc. R. SOC. London, Ser. A, 1980, 373, 27. 15 E. R. Smith, I. K. Snook and W. van Megen, Physica A, in press. 16 C. W. J. Beenakker, J. Chem. Phys., 1986, 85, 1581. Received 2nd December, 1986

ISSN:0301-7249    

DOI:10.1039/DC9878300193

出版商:RSC    

年代:1987

数据来源: RSC

摘要:

Faraday Discuss. Chem. SOC., 1987, 83, 199-211 Brownian Dynamics of Chain Polymers Marshall Fixman Department of Chemistry, Colorado State University, Fort Collins, Colorado 80523, U.S.A. The effect of constrained bond lengths and angles on the Brownian dynamics of wormlike polymer chains is calculated analytically and compared with dynamical simulations. An analytical correspondence between constraints and internal friction is developed and illustrated. The approximate numerical evaluation of the constraint matrix and its correlations, characteristic of previous work, is replaced by an equilibrium simulation. Good agreement is found between the constraint formalism and dynamical simulations. With fluctuating hydrodynamic interaction included, the constraints generate a coupling between local motions and translational diffusion which causes a small decrease in the translational diffusion constant.The effects are similar to those found in a previous phenomenological study of the effects of internal friction on translational diffusion and do not vanish with increasing chain length. It is well known that the classical Gaussian model192 of polymer chain dynamics is inadequate for short, stiff chains or for local motions that affect the moderate and high-frequency properties of flexible chains. The necessary refinements in the model seem fairly obvious for the systems considered here: dilute solutions at the 0 point, Instead of the very flexible bond lengths and bond angles allowed in the Gaussian model, only modest vibrations should be permitted and motion in the dihedral angle subspace should be restricted by rotational barriers.Progress with these improvements has been slow in each of the two approaches that promise a wide generality of application. On the one hand, a fundamental approach involving the representation of the diffusion operator on a complete basis ~ e t , ~ ' ~ although successful for very short gets bogged down in numerical difficulties for chains with more than a few bonds. An obvious analytic approximation involves the pre-average of an inverse constraint If bond angles as well as lengths are constrained, this approximation fails qualitatively for the local motions of long chains. On the other hand, the same improve- ments in predictive power which are the goal of fundamental approaches are easily accomplished by the addition of internal friction to Gaussian models." The manoeuvre can be formally justified as an approximation to memory integrals in Zwanzig-Mori projection but the actual applications depart considerably from this solid foundation.An attempt is made here to supply some of the missing structure.I3 Our confidence that a useful link can be developed between formal theory and internal friction goes back to Brownian simulations on chains with constrained bond lengths and angles and barriers to internal r0tati0n.l~ Virtually all of the normal mode structure of Rouse-Zimm theory seemed to be preserved by the more realistic model, and the similarity to the predictions of internal friction models was noted. But the results were sketchy and the connection was not made in a systematic way.A more extensive treatment has now been started. The objective is to determine first whether the analytic treatment of constraints is successful and practical for the treatment of freely rotating chains with arbitrary stiffness and length, initially without rotational barriers, and whether there is a natural correspondence to the treatment of such chains by the addition of internal friction to Gaussian models. It will be seen that these questions can be answered affirmatively. 199200 Brownian Dynamics of Chain Polymers To some extent the technical difficulties are circumvented rather than confronted. The required averages and correlations involving the constraint matrix are obtained from equilibrium simulations rather than by clever approximations.In retrospect, however, the averages are seen to have simple properties; they converge rapidly with increasing chain length to intensive parameters characteristic of the backbone structure and friction constants. Basic Equations A very brief summary of the constraint formalism will be given. Earlier work4 has been modifiedI3 by the addition of internal friction forces to the basic Langevin equation. The purpose of this addition is to allow a translation of constraints into the language of internal friction. At no point in the actual calculations is internal friction included as an independent contribution. The Langevin equation with neglect of inertia has the symbolic form m d V / d t = F H + F 1 + F P + F R --+ 0 ( 1 ) where FH is the frictional drag force due to the relative motion of beads and solvent, F' is the frictional drag due to the relative motion of beads with respect to each other, FP is a potential force and FR is a random force.These symbols all stand for column matrices of order 3( N + 1). The bead positions R,, i = 0, . . . , N, and velocities V, are likewise written as 3( N + 1 ) dimensional arrays R and V, respectively. where 1I' is the distribution function of R. In the simulation algorithm" which will be referred to below, FR is a stochastic variable. The frictional drag force is F H = -Po( V- v"), where Vo is the unperturbed velocity of the solvent and Po is a friction matrix. Po= (Pi1 + T ) - ' , where Pd is a diagonal array of bead friction constants and T is the matrix of hydrodynamic interaction HI (e.g.the Oseen interaction, although in practice the Rotne-Prager interactionI6 has always been used because it keeps Po positive definite for all configurations). The internal frictional force F' = - Qo V and the coefficient matrix Qo will be taken independent of R in order to make simple contact with memory integrals. The constraining forces have known direction (determined by the form of the backbone vibrational potential) and form one part of the total potential force F'. The unknown magnitudes are determined by the criterion that bond lengths and bond angles remain constant in time. The details of the calculation may be consulted The result is that In analytical work FR= -kTV In v = G ( p , V + F " + F R ) (2) where Fs is the remainder of FP after constraint forces are subtracted, and G= G-GMAG, G=(p,+Qo)-' and MA=AMAT. The matrix M is identical to that constructed b e f ~ r e .~ MAG is a projection matrix that prevents the bead motion from violating the constraints; the inverse of M, not M itself, is known explicitly. A,,, i, j = 0, . . . , N, stands for the elements of a constant matrix. It arises from the transforma- tion to bond vector coordinates b = ATR ( b , = R, - except for bo which is the centre of mass); AT is the transpose of A. G is a diffusion matrix in eqn (2). In the absence of constraints or internal friction the diffusion matrix would be H = P i ' . It is seen that constraints and internal friction affect the diffusion matrix in a superficially different way.If only one of these effects is present then the former leads to a subtraction from H, and the latter leads to an addition to H - ' . The subtraction is a major one if bond lengths and angles are both constrained, and approximations to M can lead to unphysical results such as negative eigenvalues in the diffusion matrix. The effect of internal friction is saved from that danger by the easily satisfied requirement that Qo is non-negative, i.e. that it has non-negative eigenvalues.M. Fixman 20 1 Introduction of the drift velocity into the conservation equation for 9 gives a diffusion equation for ZIr. This equation is usefully written in terms of @ = */ZIre, where qe is the equilibrium distribution function: With a scalar product defined as ( g , , g2) = (glg2),, i.e.as an equilibrium average, 2 is self-adjoint and for any functions g1 and g,. The explicit form of the source term S ( t ) is known in terms of V(' and external f o ~ c e s , ~ . * ~ but will not be needed here. The development of quasi-analytical approximations for transport coefficients is generally based on the representation of 3, S ( t ) , and @ ( t ) on a truncated basis set. Alternatively, if the transport coefficient has already been expressed in terms of time correlation functions of dynamical variables yk( t ) , y k ( t ) = exp ( - 2 t ) y k , then it suffices to represent 3 and the yk on the truncated basis set. We will refer here to the Zwanzig-Mori projection formalism3.12 rather than the Fixman-Kovac truncated rep- resentation~~ as a theoretical basis for this application, because the former provides an explicit remainder in terms of memory integrals.However, the memory integrals will be discarded after a brief contrast between the ways in which memory integrals and internal friction correct the truncated formulae, and the formalisms then become equivalent. The column array of y k ( t ) is designated y ( t ) . The generalized Langevin equation is X? d y ( t )/ d t = - LT, y ( t ) + K ( T ) r o l y ( t - T ) d T + A( t ) ( 5 ) 0 where I'o = ( yyT)e and eqn (4) gives Lij (yiTYj)e= kT((VTYi)G(Vyi)>e- (6) The kernel K ( T ) is a time correlation function of random forces A which are not correlated with y. Without the memory terms eqn ( 5 ) is a truncated matrix representation of the diffusion equation, as used in Fixman-K~vac.~ Eqn ( 5 ) is applicable to the caluclation of a correlation matrix I?( T) = ( y ( ~ ) y ~ ( 0 ) ) , .Improvements on dT( t)/dt =: -LI';'y( t ) + A( t ) or r( T ) = exp (- Lr&-)ro ( 7 ) which are obtained from eqn ( 5 ) after suppression of the memory integral, mist be based on an expansion of the set of yk or on an empirical treatment of the memory kernel. The imposition of constraints is one e ~ c e p t i o n , ' ~ and amounts to an evaluation of the time integral of K ( T ) over a short time interval that encompasses the vibrational relaxation times and a restriction of the yk to those variables which are essentially constant over such an interval. The use in eqn (6) of the constrained diffusion matrix G rather than the unconstrained G implies that the part of K ( T) representing the decay of vibrational forces has been assumed to have infinitely rapid decay and can be written as S ( T ) & , where KO is a constant matrix.These forces then modify the dissipative factor L in eqn (7), which is replaced by L - K O , rather than the equilibrium factor r i l , and in this respect behave similarly to internal friction. The first question to be addressed is the adequacy of eqn (7) in its application to a standard worm model of polymer chains with a minimum set of y k .202 Brownian Dynamics of Chain Polymers The Y k will be chosen as Rouse modes or as the direct products of such modes. The Rouse vector modes are defined as N qi = C Q.-b. j = O ?I J where QOo = 1, QOj = Qjo = 0 for j > 0, and An additional transformation b = ArR allows the q values to be expressed in terms of bead coordinates.For vector-mode correlations the basis set is composed of Y k = q f , 1 C k S N, where qz is the x component of qk. The tensor basis set consists of the products yk = q;qi. Time correlations of such modes suffice to describe dielectric relaxation, if the dipoles are rigidly affixed to the backbone, and stress relaxation, if HI is pre-a~eraged.~”~ However, that sufficiency does not guarantee than eqn (7) is accurate, i.e. that the basis set is large enough that memory integrals are negligible. This adequacy is judged by a comparison with dynamical simulations. The matrix elements of To are just ( ~ ~ y , ) ~ and their evaluation from equilibrium or dynamical simulations is clear in principle, once the y1 and the model are defined.For the calculation of L,J from eqn (6) it is desirable to transform from the original column matrix V y J to gradients with respect to Rouse modes, since the yJ have a simple dependence on them. The transformation to b derivatives, V = AVb, gives L I J = kT((v ’YI lTB(’ bYJ i > e (10) where B = ATGA = B - BMB; B = ATGA. (11) Except for the possible inclusion of internal friction Qo in G = (Po+ QO)-’ and the retention of the centre of mass coordinate bo, B and M are as defined previ~usly.~ For the orthogonal transformation from b to q coordinates, we use V h = QVq, where Q is the matrix of Qkl and V f = d / d q , . With B4 = QBQ, L I J = kT((V4r,)TBq(VqyJ)),’ (12) Approximations and Internal Friction We briefly summarize the major simplifications of the constraint formalism that have been tested by a comparison with results from the full formalism or with dynamical simulations.Several conclusions are anticipated. ( a ) Suppression of Constraints The suppression of constraints means that the constraint matrix M is discarded; i.e. G is replaced by G, or B is replaced by B in the expression for L, eqn (12). Suppression of constraints and internal friction restores the theory to the Bixon-Zwanzig formalism in a version that neglects memory integrals; this is very nearly a Gaussian formalism. The results are qualitatively reasonable for mode correlations but not for stress relaxatim. The results are always quantitatively poor for short, stiff chains, and are worse for short-wavelength, rapidly relaxing modes than for long-wavelength modes.The failure for rapidly relaxing modes does not diminish with increasing chain length. However, the low-frequency motion of long, flexible chains does not require retention of con- straints.M. Fixman 203 ( b ) Diagonal L Approximation In this approximation the off-diagonal elements of the matrix representation of 3 on the Rouse basis set are discarded. This will be seen to be an excellent approximation under almost all circumstances for uniform worm models. It is worse f,, non-uniform chains. When the approximation is valid the formalism and the results ma, ' - e interpreted in terms of internal friction, as will be shown below.We expect that this approximation will remain valid for chains with rotational barriers. (c) Factored Tensor L Approximation A related approximate factorization of the diagonal elements of L formed on the basis of tensor modes will be seen to be reasonably accurate. For the tensor modes eqn (12) We define the factored approximation to Lkk as A superscript v on L or Lik 3 kT( B;Ik), and coupled to ( q i ) 2 , then Lkk = Lfkk will be a good approximation. = 2 k ~ ( B 3 , ( ( ~ ; ) ~ ) , = x;,r;,, . (14) designates an evaluation on the vector basis set. Thus = (( qi)2)e. If fluctuations in the constraint matrix are weakly ( d ) Internal Friction Description This is not presented as another approximation. Our purpose is to show that the constrained diffusion matrix without internal friction can be replaced by the uncon- strained B with internal friction, and that internal friction has simple properties in Rouse coordinates if approximations ( b ) and ( c ) are valid. The internal friction matrix Qo appears in combination with the ordinary friction matrix Po= H-' in the matrix elements L,j of the time evolution operator. Specifically the diffusion matrix B in bond vector coordinates is, from its definition, where Bo= ATHA and A.is the 3 ( N + 1) x 3 ( N + 1 ) matrix A o = A-'QoAT-'. If Qo generates a frictional force on bead i which is due to beads i* 1 and which is equal to A ( V l i l - V ; ) , as in Bazua-Williams,'* then A. is A times a unit matrix in any representa- tion. More generally, eqn (15) gives where A:= QA,Q and B: = QBoQ We will assume that A:, although not generally proportional to the unit matrix, is to a sufficient approximation diagonal with elements A&.This is also a well known approximation' that we use for B,9; the diagonal elements are designated B&. Consequently, eqn (16) associates a friction constant with each Rouse mode. It remains to be demonstrated that internal friction of this form can emulate constraints, i. e. that BY with internal frictjon can approximate B4. With the substitution of BY from eqn (16) for BY in eqn (12) for L, and the use of a pre-averaged, diagonal B,9 and a diagonal A; in eqn (16), the factored approximation is always valid and we may write, in analogy to eqn (14), B = ATGA = AT( H - ' + &)'-'A = (B;' + Ao)-l (15) Bq QBQ = (Bg-' + A:)-' (16) Superscript I on L indicates use of the internal friction approximation, and superscript v on L or To indicates the vector as against the tensor basis set.It is now evident that204 Brownian Dynamics of Chain Polymers 0.0 n 0 h v 0 -1.0 > 0 % v U E - - 2.0 100 2 200 Fig. 1. Comparison of simulation and analytical tensor mode correlation functions for a 10 bead chain with b = 10, constrained bond angles equal to 160°, d s = 5, mode k = 5 , and no HI. L,/p = 0.6. ( a ) Dynamical simulation, ( b ) diagonal L approximation, (c) full analytic result, ( d ) omission of constraint matrices, ( e ) constraints and off -diagonal elements in the equilibrium correlation matrix To are omitted. the h z k can be chosen to make L& = L L k or L L k = L k k , and that both equalities can be satisfied simultaneously if the factored approximation is adequate.In summary, the adequacy of the diagonal L approximation guarantees that the effects of constraints can be imitated by an internal friction matrix which is diagonal in Rouse coordinates. The adequacy of the factored L approximation guarantees further that these diagonal elements are the same for vector and tensor relaxation properties. Results We begin with a summary of conventions. All of the graphical and tabular results are expressed in reduced units. The unit of length is a reference bond length brer, the unit of energy if kT, the unit of time is P,,,bf,,/kT, where Pref is a bead friction constant given by 3 rrqObref. Friction constants are specified through the equivalent Stokes diameter d as ratios to bref.No special notation is used to designate quantities expressed in these reduced units. Chain stiffness is measured by the ratio LJp, where L, is the contour length and p is the persistence length. Results obtained from the truncated matrix representation, eqn (7), will be described here as ‘analytical’ in order to distinguish them from the results of dynamical simulations. Comparison of Results with Dynamical Simulations Dynamical simulations were performed only for freely rotating chains without HI. The algorithm used for constrained simulations, with or without HI, was recently de~cribed.’~ Results are shown only for a nine-bond chain with L J p ==: 0.6. Quite similar comparisons were found for five-bond chains with larger and smaller values of LJp.Fig. 1 showsM. Fixman 205 Table 1. Full and approximate analytical results for a constrained bond angle chain with an angle of 163", b = 5, and d s = 3 (pre-averaged HI is included, L,/p =: 2); N = 40 ( a ) tensor mode friction constants ~ k P k 1 0.258 x lo3 0.229 x lo3 0.292 x lo2 0.254 x lo3 2 0.856 x lo2 0.645 x lo2 0.210 x lo2 0.834 x lo2 3 0.465 x lo2 0.319 x lo2 0.146 x lo2 0.457 x lo2 30 0.868 x 10 0.466 x 10 0.402 x 10 0.877 x 10 20 0.701 x 10 0.195 x 10 0.506 x 10 0.720 x 10 30 0.943 x 10 0.141 x 10 0.802 x 10 0.940 x 10 40 0.103 x lo2 0.128 x 10 0.899 x 10 0 . 1 1 0 ~ lo2 vector tensor k analytic diagonal L no constraint analytic diagonal L no constraint 1 0.213 x lop4 0.218 x lop4 0.235 x lop4 0.696 x lop4 0.699 x 0.788 x 2 0.179 x 0.182 x lop3 0.229 x 0.450 x 0.449 x 0.591 x lop3 3 0.107 x lop2 0.108 x 0.152~ 0.176 x lop2 0.178 x lop2 0.258 x 10 0.772 x lo-' 0.751 x lo-' 0.146 0.119 0.119 0.228 20 0.298 0.317 0.1 19 X 10' 0.490 0.497 0.173 x 10' 30 0.517 0.525 0.331 x 10' 0.912 0.896 0.592 x 10' 40 0.894 0.889 0.613 x 10' 0.151 x 10' 0.149 x 10' 0 .1 1 2 ~ lo2 the time correlation function C ( t ) = ( y k ( t ) Y k ( 0 ) ) e vs. t, where y k = q i q i and k = 5. Results from the various analytical approximations are shown along with results from the dynamical simulation. The value of C(0) must always be given correctly by any of the analytical approximations and is not displayed. The initial time derivative of the time correlations must always be given correctly by the full constraint formalism and the agreement seen in fig.1 for short times indicates only the (apparent) absence of programming errors. Of greater import is first that there is no statistically significant discrepancy between the constraint formalism and the dynamical simulations throughout the duration of the trajectories. The quantitative error caused by omission of constraints increases with mode number k, but it is significant even for k = 1 and for vector as well as tensor modes, unless L J p >> 1. Secondly, the mixing between modes demonstrated by the non-exponential decay is accounted for almost entirely by the equilibrium correlation To. Off-diagonal elements of L play essentially no role in the mixing of Rouse modes. According to much other evidence similar to fig.1, the agreement between analytical results and dynamical simulations seems to be satisfactory for constrained chains. We now proceed to examine various features of the analytical formalism and its consequences on the assumption that the constraint formalism is essentially exact for wormlike chains. Internal Friction It follows from the adequacy of the diagonal approximation to L, which will be further illustrated below, that the effect of constraints on the dynamics can be duplicated by an appropriate choice of a diagonal internal friction matrix A,Y in eqn (16). The formulae for Lkk (tensor modes) and L i k are given in eqn (17) and (18). The values of the required matrix elements A & are illustrated for tensor modes in table 1 and in fig.2. The analytical206 Brownian Dynamics of Chain Polymers 30 20 cz" 10 i ' " " " ' 0 0.2 0.4 0.6 0.8 k l N 3 Fig. 2. ( a ) The total mode friction constant Pk and ( b ) h X k , i.e. the contribution to Pk from constraints. Pre-averaged HI is included. Solid curves refer to 80 bond chains with L, == 48 and the dashed curves to 20 bond chains with L , / p =: 12. For each chain b = 5 , d s = 3, the angles are constrained at 122", and p = 8. results for the mode friction constants Pk = (BZk)-' are constructed from equilibrium simulations of L and To according to eqn (16)-( 18): P k B:; + h s k 2 k T r g k k / Lkj,. (19) h & P k - B;k1 2 k T r g k k / L k k - BZF' (20) The contribution to P k from constraints defines the internal friction h&: where I'; is the proper equilbrium correlation matrix of vector modes for the model, L is constructed for the required tensor modes with the appropriate constrained diffusion matrix B, and B,4 is constructed with omission of both constraints and internal friction. measures the effect of constraints in slowing down the relaxation; equivalently, A & is the amount of internal friction that has to be added to B:il, according to eqn (16), in order to make Bq reproduce the effects of constraints. Table 1 shows by the agreement between the exact P k and the approximate PL computed in the factored L approximation of eqn (14), that the factored approximation is reasonably good even for a very stiff chain.The approximation improves with increasing flexibility but it is reasonably good, as is the diagonal L approximation, even in the limit of a straight rod.P k is defined in such a way as to become proportional to the relaxation times T k for chains sufficiently long and flexible that To is diagonal as well as L. The relaxation times may be calculated generally, according to eqn (7), as the eigenvalues of ToL-', and a sampling of rates 1 / ~ ~ computed in this way is given in table 1. With the interpretation of A & just given, eqn (20) may be regarded as a calculation of the internal friction A & as a function of k, and the result can be contrasted with based on physical models or comparisons with experimental results. For long chains with elastic or constrained bond angles, h& becomes a function of k / NM.Fixman 207 as N + m . As the persistence length of the chain is increased, table 1 shows that h& for large k increases with k if bond angles are constrained (but it does not increase if the angles are elastic). For either model, elastic or constrained bond angles, h tjk increases for small k as the persistence length increases; see table 1. However, A & does not show the divergence as an inverse power of k / N , which was proposed in some of the works of Cerf, Peterlin and Allegra." Rather A & remains finite in the limit N -+ a. In this respect the results agree more closely with those of BazGa and Williams," who arrived at the simplest model in which htjk is independent of k and N. It seems probable that the rise in h& at small k has only an accidental relation to previous suggestions of this effect.Indeed, it also seems probable that a small change is required in the common picture of 'internal friction', or at least in the language used to describe it. In conventional language internal friction is introduced in order to slow down relative, internal motion in a real chain. It would follow that internal friction cannot affect the rigid-body translational or rotational motion of any chain. However, A & shows all the mathematical characteristics of internal friction and yet it is quite large for a rodlike molecule, which is capable only of rigid body motion. To dispose of this verbal paradox one must recall that internal friction is always used in practice as a correction to the motion of a 'Gaussian' model.This may be a generalized Gaussian model which makes the covariance matrix of bond vectors agree with that for a real chain or a rigid rod, but it is nevertheless a model which does not recognize dynamical constraints, i.e. the model does not accurately describe the velocities of the beads. This failure will apply also to the 'modes' of the Gaussian model, which are just linear combinations of bead coordinates. The velocity or relaxation rate of such a mode will not be calculated exactly in the Gaussian model even if the mode is a rigid body displacement or rotation. Discrepancies between velocities of the model chain and the real chain may emerge virtually instantaneously, because of constraints or vibrational forces, or after some time lag, because of hindered rotation barriers.In either case any initial displacement of the Gaussian chain will generate opposing forces different from those of the real chain and the relaxation rates of the Gaussian chain will generally be larger than for the real chain. For this reason 'internal friction' may be required to slow the relaxation (with the exception of a translational displacement without HI, for which the opposing force remains forever independent of configuration). First Cumulant and Translational Diffusion The first cumulant is defined" as R(q) = (p*Zp)J(p*p},, where N p = C exp (iq- &). k - 0 n(q) is an 'initial' time derivative of an observed structure factor, but the meaning of 'initial' depends on the degree of coarse-graining in time that is implied by the construc- tion of the diffusion operator 2, i.e.on whether or not constraints are used. We have performed the equilibrium simulations required for an evaluation of R(q) vs. q for various models with fluctuating HI. Results for constrained angle chains are similar to those found by Akcasu et al.," who considered freely jointed chains and used an analytical pre-average of M - ' . Results for a Gaussian model in which the bond vectors were generated independently, but the resultant bond lengths and angles were taken to be dynamically constrained, were litt!e different from those for a constrained 90" bond-angle chain. The unconstrained Gaussian model smooths out the peaks in the first cumulant relative to the unconstrained 90" chain, but has little effect on the location of the curve.Introduction of constraints into either the Gaussian or 90" chain smooths out the peaks, as does strong us. weak HI, depresses O ( q ) at large q, and makes the two models look quite similar.208 Brownian Dynamics of Chain Polymers Table 2. Translational diffusion constants as ratios of value with constraints to values without constraints. Constrained bond angles of 90” were used, constrained bond lengths b = 5 , and fluctuating HI. The Stokes diameters are shown as values of d s / b . The number of beads ranges from 14 to 112. The ratios have standard deviations of CCI. 0.003 N + l d S / b 14 28 56 112 ~ ~~ 0.5 0.977 0.966 0.965 0.961 1 .o 0.947 0.93 1 0.93 1 0.927 2.0 0.923 0.890 0.883 0.880 The depression in s l ( q ) for large q, owing to the addition of bond angle constraints, is rather modest for the relatively short chains that were studied.This is due to a significant contribution from the translational mode that would not be present for flexible chains of realistic size. Moreover, barriers to internal rotation would make the internal friction much larger than has been found for freely roating chains, and the high-frequency limit of is likely to be quite depressed. We estimate the depression for a simple model of internal friction. See also Allegra et al.19 We assume that N is sufficiently large that the sums over modes used to compute R(q) can be converted to integrals, that hydro- dynamic interaction is relatively unimportant for large Rouse mode numbers, and that A & = h for all k.This gives where Pd is here the friction constant of a single bead. If h = O the high-q limit of q-*sl( q ) is the value for Gaussian chains without constraints. With increasing internal friction the limiting value slowly decreases. The simulation results were typically about half that given by eqn (22) because of a substantial contribution from the translational mode k = 0 for short chains. Additional data on R(q) were gathered for rather large bead friction constants in order to assess the effect on the translational diffusion constant of coupling between the constraints and fluctuating HI. The full analytical constraint formalism was used for this purpose. It was anticipated that the constraints would behave similarly to internal friction2’ in affecting the diffusion constant. Similar results were indeed found, as is shown in table 2.The translational diffusion constant was considered to be proportional to R(O), and this amounts to using the Kirkwood formula as before2’ ( i e . the diffusion constant is taken p_roportional to a double sum of the averaged elements of the diffusion matrix), with G replacing G and with internal friction omitted. The diffusion constant so calculated is actually an upper bound to the true value.20 It seems reasonable on the basis of table 2 to conclude again that a coupling between fluctuating HI and local motion affects the translational diffusion constant even in the limit of infinite chain length. Stress Relaxation In fig. 3 and 4 the storage and loss moduli G’(o) and G”(o) are displayed as functions of the reduced frequency oR = wqo[ 77]/kT, where qo is the solvent viscosity and [ 771 is the intrinsic viscosity defined with concentration units of molecule volume-’.These modulii are designed [G’], and [G”]. in the book by Ferry.21 The more interestingM. Fixrnan 209 Fig. 3. Comparison of analytic approximations for the shear modulii us. wR. ( a ) The loss modulus G' and (6) the storage modulus G". Pre-averaged HI is included. The 40 bond chain has constrained bond angles of 163", 6 = 5 , and d s = 3. L,/p = 2. (-) Full analytic results; (- - -) constraint matrices omitted. The diagonal L approximation coincides with the full analytical result. -1.0 0.0 1.0 2 .o 3.0 log O R Fig. 4. Dependence of shear modulii on wR and chain length.Same as fig. 3 except that bond angle is 122" and only full analytic results are shown. (-) 80 bond chain with L,/p = 48; (- - -) 10 bond chain with L , / p = 6.210 Brownian Dynamics of Chain Polymers extreme approximations are shown along with the full analytical results in fig. 3. The complete neglect of constraints is obviously very unsatisfactory, particularly for the loss modulus, because the high-frequency limit of the viscosity vanishes in this approximation. The diagonal L approximation is generally quite satisfactory, although in some instances deviations are seen in G' at high frequencies, around wR= lo3. The storage modulus G'( w ) may become quite sensitive in the high-frequency region, around w 2 (1) in the reduced units, to sampling errors and approximations such as a restriction of L to its diagonal elements.In this high-frequency region, corresponding to t O( l ) , G ' ( w ) is levelling off at its asymptotic value. This value is determined by the residual flexibility at high frequencies, which is in turn very sensitive to coupling between Rouse modes under some conditions. These conditions are the ones which cause the constraints to leave very little freedom available to the local motions, as happens if the bond angles are constrained rather than elastic. Stress relaxation is the most sensitive of quantities ordinarily studied at high frequencies owing to its direct dependence on forces as well as configuration. The neglect of memory integrals in the constraint formalism is safest at the shortest times, but the numerical approximations and the physical applicability of the formalism become doubtful in an application to very short times.Ordinarily these short times will not be observed in mechanical relaxation experiments. Concluding Re marks The internal friction constant A & for mode k becomes a function of ( k / N ) rather than k and N separately, for moderate chain lengths with L,> p . This means that studies of sort chains suffice for the prediction of the properties of long flexible chains. For constrained but freely rotating chains the internal friction constants may be accurately evaluated from equilibrium simulations. In general the functional dependence of h& on k / N for constrained chains is different from previous suggestions.The function is not singular at small k / N, contrary to the behaviour of the total mode friction P k , which diverges as k-3'2 for chains with hydrodynamic interaction, and contrary to some of the suggestions of k-' or k-"2 divergence in A&. On the other hand, A & may show a sharp rise at small k / N , which is significant in the quantitative interpretation of stress relaxation in finite length chains, and a smooth rise at large k / N , which is significant to other experiments. Thus the most elementary physical models of internal friction seem inadequate. There remains the question of the effects of rotational barriers on these simple findings. The previous ~ i m u l a t i o n s ~ ~ of 90" constrained angle chains led to the pre- liminary suggestion that the Rouse-Zimm mode structure was well preseved in the presence of barriers.If this conjecture stands up to more thorough tests on stiffer chains than were studied before, as we expect, then the imitation of realistic backbone potentials with internal friction will be fully justified. However, one modification of the present treatment must be anticipated. The internal friction force will have to be written as an integral over the history of relative velocities, and the kernel will require a finite rather than an infinitesimal relaxation time. This work was supported in part by NIH GM27945, and was initiated as a result of conversations with R. Pecora and D. Roitman at a conference on the Rotational Dynamics of Small and Macromolecules in Liquids, at the University of Bielefeld, April, 1986, organized by Professors Th. Dorfmiiller and R. Pecora. References 1 P. E. Rouse Jr, J. Chem. Phys., 1953, 21, 1272; B. H. Zimm, J. Chem. Phys., 1955, 24, 269. 2 H. Yamakawa, Modern Theory of Polymer Solurions (Harper and Row, New York, 1971).M. Fixman 21 1 3 R. Zwanzig, J. Chem. Phys., 1974,60, 2717; M. Bixon and R. Zwanzig, J. Chem. Phys., 1978, 68, 1896. 4 M. Fixman and J. Kovac, J. Chem. Phys., 1974, 61, 4939. 5 D. B. Roitman and B. H. Zimm, J. Chem. Phys., 1984, 81, 6333; 6348. 6 K. Nagasaka and H. Yamakawa, J. Chem. Phys., 1985,83, 6480. 7 G. T. Evans, J. Chem. Phys., 1979, 70, 2362. 8 M. Fixman and J. Kovac, J. Chem. Phys., 1974,61,4950; M. Fixman and G. T. Evans, J. Chem. Phys., 1976, 64, 3474. 9 H. Yamakawa and T. Yoshizaki, J. Chem. Phys., 1981, 75, 1016; T. Yoshizaki and H. Yamakawa, J. Chem. Phys., 1986,84,4684; H. Yamakawa, T. Yoshizaki, and M. Fujii, J. Chem. Phys., 1986,84, 4693. 10 A. Z. Akcasu, B. Hammouda, W. H. Stockmayer and G. Tanaka, J. Chem. Phys., 1986, 85, 4734. 11 G. Allegra, J. Chem. Phys., 1986, 84, 5881. 12 Statistical Mechanics. Part B : Time -Dependent Processes, ed. B. Berne (Plenum Press, New York, 1977). 13 M. Fixman, unpublished results. 14 M. Fixman, J. Chem. Phys., 1978, 69, 1527; 1538. 15 M. Fixman, Macromolecules, 1986, 19, 1195; 1204. 16 M. Fixman, J. Chem. Phys., 1982, 76, 6124. 17 U. M. Titulaer and J. M. Deutch, J. Chem. Phys., 1975, 63, 4505. 18 E. R. Bazfia and M. C. Williams, J. Chem. Phys., 1973, 59, 2858. 19 G. Allegra, J. S. Higgins, F. Ganazzoli, E. Lucchelli and S. Bruckner, Macromolecules, 1984, 17, 1253. 20 M. Fixman, J. Chem. Phys., 1986, 84, 4085. 21 John D. Ferry, Viscoelastic Properties of’ Polymers (John Wiley, New York, 1980). Received 5th December, 1986

ISSN:0301-7249    

DOI:10.1039/DC9878300199

出版商:RSC    

年代:1987

数据来源: RSC

摘要:

Faraday Discuss. Chem. SOC., 1987, 83, 213-222 Trajectory Simulation Studies of Diff usion-controlled Reactions J. Andrew McCammon* and Russell J. Bacquet Department of Chemistry, University of Houston, University Park, Houston, Texas 77004, U. S. A. Stuart A. Allison Department of Chemistry, Georgia State University, Atlanta, Georgia 30303, U.S.A. Scott H. Northrup Department of Chemistry, Tennessee Technological University, Cookeville, Tennessee 38505, U.S.A. Computer simulations of the Brownian trajectories of reactant molecules in solution can be used to calculate the rates of diffusion-controlled reactions. This paper presents new derivations of some of the basic equations used in such calculations and the results of applications to two problems from biology. Trajectory simulation studies have provided the most detailed information available for bimolecular chemical reactions in the gas phase.In a typical calculation, a reactant molecule starts moving toward a distant target molecule and the appropriate equations of motion are solved numerically for a long enough period to insure that the encounter has been completed.' This straightforward procedure has to be modified to allow trajectory simulation studies of diffusion-controlled bimolecular reactions, because it is not generally possible to identify conditions that correspond to the completion of an encounter. A reactant molecule that is undergoing Brownian motion always has a chance of diffusing to a target molecule, even if the two molecules have drifted apart a sizeable We have developed a variant of the trajectory simulation approach that allows one to determine the kinetics and mechanisms for detailed models of diff usion-controlled reaction^."-^^ The basic idea behind the approach is a simple one.Trajectories are computed only for reactant molecules that are relatively close to a target molecule. The contribution of reactant motions farther from the target, where the effective interactions between reactant and target assume simpler forms, is handled analytically. In what follows, we first outline methods that can be used to combine trajectory results and analytic results to obtain rate constants for bimolecular diff usion-controlled reactions. We then illustrate these methods by describing recent studies of diffusion- controlled reactions in two biochemical systems.Theory Fundamentals We consider the diffusion number of dilute systems. ensemble of systems is2-'' of a reactant molecule toward a reactive target in Then the bimolecular rate constant at steady state a large for this ( 1 ) k = -c0-' [ s dAn . j ( A ) 213214 Simulation Studies of Difusion-controlled Reactions where co is the concentration of reactants far from the target, S is any surface that encloses the target or, more generally, is intersected by all reactive trajectories, A denotes an element of the surface S, n is the outward surface unit vector at A and j ( A ) is the ensemble average current density of diffusing reactant at A. The connection between this equation and trajectory calculations is made by decomposing j ( A ) into the condi- tional current jo(A) that counts any reactant reaching S at A for the first time, times the probability p(A, *) that a reactant at A will ultimately react rather than escape." Eqn ( 1 ) can then be written where g(A) = -c,'n jo(A).(3) If S has a simple shape and is located such that the interactions between reactant and target are simple for reactants outside of S, then g(A) can be determined analytically by solving the appropriate Smoluchowski equation. Trajectory calculations are used to obtain the remaining factor p(A, *) by simulating the diffusion of reactants starting at S. To deal with the difficulty described in the Introduction, it is necessary to truncate trajectories that diffuse beyond a certain distance from the target and to correct for the effects of this truncation.The necessary analysis is outlined in Appendix A. The results are especially simple if the initiation and truncation surfaces are spheres, with the respective radii b and q chosen large enough that the potential of mean force for the reactant displacements is centrosymmetric for target-reactant distances r > b. In this case k = kD(b)Pm (4) where in which U ( r ) is the potential of mean force and D( r ) is the relative diffusion coefficient; and P l-(l-p)n p m = in which p is the probability that trajectories initiated at r = b will react rather than be truncated at r = q, and = k,( b ) / kD( q ) . Staging Techniques The efficiency of the rate-constant calculations can be increased significantly by breaking the trajectory simulation calculations up into stages."," For example, in a problem that involves complicated short-range interactions between reactant and target, it may be useful first to calculate the probability for a reactant to diffuse to a surface relatively close to the target, and then to compute the probability of diffusing from this intermediate surface to the site of reaction (fig.1). This can be especially useful if one wants to study the dependence of k on local details of the reaction site, since it may then be necessary only to recompute the second probability. Because the intermediate surface in a staged calculation is generally close to the target, it will typically be necessary to retain the orientational dependence of the factorsJ. A.McCammon, R. J. Bacquet, S. A. Allison and S. H. Northrup 215 Fig. 1. Schematic illustration of the surfaces used in a trajectory calculation of the rate constant for a diff usion-controlled reaction. The kidney-shaped object is a target molecule, and cross-hatched patch is a reactive site. in eqn (2). This can be done to an arbitrarily good approximation by replacing integral in eqn (2) by a sum over discrete elements Si of the intermediate surface: k = C kD(Si)Pa,i I where and Px,i is the probability that a escape. The quantities kD(Si) can be trajectory started at Sj will ultimately react rather than obtained as kD(Si) = kD(b)Pm,i (9) by a procedure very similar to what has been described in connection with eqn (4). One initiates trajectories at a spherical surface at the relatively large radius b, and computes the probability pi of intersecting Si before intersecting any other surface element Sj (j # i) or reaching the truncation surface at radius q.Then, as shown in Appendix A, Pi p . = OCj3I l-(l-P)n The quantities Pm,i can be obtained from216 Simulation Studies of Diflusion-controlled Reactions Fig. 2. Schematic illustration of the active site region in the enzyme SOD. where Pi is the probability that a trajectory started at Si will lead to reaction before it is truncated at radius r = q, and kdS) = c kD(St). 1 Eqn (7), with given by eqn (12), is derived in Appendix B. Trajectory Calculations A number of algorithms are available for computing Brownian trajectories.' All of the reaction studies reported to date have used either the original Ermak-McCammon a l g ~ r i t h m ' ~ or recent extensions due to Dickinson et al." and to Lamm.I6 The Ermak- McCammon algorithm is quite general in that it allows one to compute trajectories of reactants of arbitrary shape and flexibility, and it allows for the inclusion of hydrody- namic interactions.For spherical reactants, the algorithm of Dickinson et al. incorporates an efficient treatment of translation-rotation ~ o u p l i n g . ' ~ A limitation of both of these algorithms is that they correspond to reactant motion in the absence of reactive or reflecting boundaries. As a result, it is necessary to take very small steps near such boundaries. The Lamm algorithm explicitly includes the effects of such boundaries, but has been developed only for motion near planar boundaries in the absence of hydrody- namic interactions.I6 Northrup et al.have shown that the explicit inclusion of reaction and reflection effects can be used as the basis of a very efficient procedure for studying the effects of varying the geometric conditions for a reaction, as long as these conditions are not too complicated." Superoxide Dismutase As an example of the application of the trajectory simulation method to a rather detailed model system, we consider some ongoing studies of the diffusion-controlled reaction between the substrate superoxide (Or) and the enzyme superoxide dismutase (SOD). Past studies of this system using the trajectory simulation method have employed a simple topographical model of a spherical enzyme with reactive patches at the poles corresponding to the active sites of the dimeric enzyme.17-19 In the present work, a significantly more realistic topographical model was used (fig.2). In this model the X-ray coordinates of the (dimeric) enzyme were su erimposed on a three-dimensional grid with 1 8 , resolution.*" Volume elements (1 x3) in the grid were excluded toJ. A. McCammon, R. J. Bacquet, S. A. Allison and S. H. Northrup 217 penetration by 0, if such penetration would bring the centre of the substrate closer than ca. 2.75 8, to the centre of any non-hydrogen atom of the protein. In the course of a particular trajectory, a dynamics step which placed substrate in an excluded volume element was disregarded and the trajectory continued from the previous step.Past studies have shown that the reaction rate of this system is enhanced by the electrostatic charge distribution of SOD and that this distribution tends to ‘steer’ substrate molecules into the active sites of the enzyme.”-23 In the original Brownian dynamics the distance of closest approach between the spherical model of SOD and the substrate was 30 A. Also, simple three- and five-charge models of the enzyme yielded reaction rate constants nearly identical to those of more realistic models, where charges were placed at the coordinates of charged groups (76 charge model) or at the coordinates of all non-hydrogen atoms (2196 charge model) of the enzyme.’’ The simple charge models approximately reproduce the monopole, dipole and quadrupole moments of SOD.Recently, the potential of mean force between 0, and the three charge, uniform dielectric ( E = 78 corresponding to water) spherical model of SOD was determined by Monte Carlo methods for a series of monovalent salt ~oncentrations.~~ The empirical energy functions derived from these potentials of mean force were used in this stud or more (outer zone). For shorter distances (inner zone), pairwise additive forces using the 76 charge model of SOD were used. Counterion effects in this zone were approxi- mated by Debye-Huckel (point ion) screening of the protein charges. A uniform dielectric constant not necessarily that of water was also assumed. To avoid the need of repetitive and time consuming computation in the inner zone, forces were computed at grid points ( 1 A resolution) and tabulated at the start of the simulation.The staging technique discussed in the previous section was used to determine k. The intermediate surface consisted of a sphere of radius 30 A which circumscribed the enzyme. This surface was subdivided into 20 latitudinal surfaces (S,) or spherical sectors of equal area. The axis defined by these spherical sectors was chosen to correspond to the symmetry axis of the three-charge model use$ to compute the forces in the outer zone. First, trajectories were initiated at b = 120 A and terminated when they reached one of the S, surfaces or diffused beyond q = 500 A. From a large number (14 000) of independent trajectories, the k,( S , ) were computed using eqn (9)-( 11). Secondly, many trajectories were initiated on each of the S, surfaces and terminated when substrate diffused to within a predefined distance ( r * ) of one of the active site Cu2+ atoms or diffused beyond q=5OOA.The rate constant was then determined using eqn (7) and (12)-( 14). The reaction distance, r*, is expected to have a significant effect on rate. From a series of independent simulations with r”=3.0, 4.0 and 5.0& we obtained k = 0.72(0.06) x lo9, 2.87(0.06) x lo’, and 3.97(0.23) x lo9 dm’ mol-’ s-l, respectively, where quantities in parentheses denote uncertainties in the rates. The following choices were made: monovalent salt concentration (C,) was 0.5 mmol dm-’, T = 298 K, E = 50/78 in the inner/outer zone, and D, = 1.238 x cm2 s-’ (hydrodynamic radii of 28.5 A and 2.05 A for SOD and O,, respectively). Experimental rate constants are ca.2.0 x lo9 dm3 mol-’ s - ’ . ~ ’ From the simulations, it was observed that if r* 5 4 A, the two Cu2+ atoms did not appear to be equally reactive or accessible and for r* = 3 A, only one of the Cu2+ atoms was accessible. A value of r* = 4 A was used in the simulations described in the following paragraph. The influence of electrostatic forces on reaction rate was studied by varying the dielectric constant in the inner (but not the outer) zone and by varying C,. Simulations with C, = 0.5 mmol dm-3 and &(inner) = 30, 50 and 78 yielded k = (6.29k0.17) x lo9, (2.87 f 0.06) x lo’ and (1.12 f 0.05) x lo9 dm’ mol-’ s-’, respectively. These results show how sensitive the rate is to short-range electrostatic forces and that these short-range to compute direct forces when 0, and the centre of the enzyme were separated by 30 K218 9.4 9.2 9.0 8.8 8.6 8.4 8.2 OQ I Simulation Studies of Difusion-controlled Reactions - - - - - - - I I 1 1 I 1 9.6 9-8 t X \ L 1 1 I 1 I -2.5 -2.0 -1.5 -1.0 -0.5 0.0 log z Fig.3. Bimolecular rate constant ( k is in units of mol s-’) for the SOD-0; reaction as a function of ionic strength, I. The crosses are experimental data.21 forces are inherently attractive. Finally, we consider the effect of salt on rate assuming Debye-Huckel (point ion) screening in the inner zone, where E = 50. For C, = 0.009, 0.104 and 0.507 mol dm-’; we obtain k = (3.17 f 0.06) x lo9, (2.18 k0.12) x lo9 and (0.79 f 0.09) x lo9 dm3 mol-’ s--’, respectively. Experimentally, k decreases with increasing salt for C, > 0.01 mol dmP3,” and the simulations are in qualitative agreement (fig.3). Cytochromes The diff usion-controlled association of the electron transport proteins cytochrome c and cytochrome c peroxidase is being studied by the Brownian trajectory method in current work by Northrup et al.2s726 Novel in these simulations is the treatment of two macromolecules undergoing rotation and having anisotropic charge distributions and reactivity. This requires the calculation of the influence of electrostatic torques on the rotational motion of each macromolecule. In an initial study, the proteins were modelled as spheres, each having a central monopole carrying the net charge and also two charges of opposite sign imbedded 1.0 A inside the protein surface.The magnitudes of the monopole and dipole moments generated by this simplified charge distribution were determined directly from the X-ray crystallographic coordinates of the proteins by assigning partial charges to each non-hydrogen atom. Solvent electrolyte effects were roughly modelled using the Debye-Hiickel screening law modified to include the effect of the large ion-exclusion sizes of macromolecules. Rate constants obtained by the simulation compare moderately well with experi- mentally determined rate constants for protein docking at physiological ionic ~trength.~’ The electrostatic torques generated by inclusion of the dipoles are of crucial importance in overcoming the severe rate-retarding effect of the orientational requirements for electron transfer.The computed ionic strength dependence of the rate constant is weaker than that observed in the experiment, presumably because of the failure of the Debye- Huckel potential to adequately treat screening effects from the solution. Additional simuiations are currently being performed with a detailed charge distribu- tion including 82 and 34 charges on cytochrome c peroxidase and cytochrome c, respectively. An accurate surface topography of both proteins is being implemented inJ. A. McCammon, R. J. Bacquet, S, A. Allison and S, H. Northrup 219 a manner similar to that employed for the SOD studies, except that both reactant species are treated as topographically rough macromolecules. Future studies will attempt a more accurate treatment of the electrostatic interactions, as discussed in the following section.Concluding Remarks Most of the testing and development of the Brownian-dynamics trajectory simulation method has necessarily utilized quite simple model systems, such as spheres with reactive patches and easily specified charge distributions. Initially, salt and dielectric effects were included only in the simplest possible fashion, i.e. charges were screened by a Debye-Huckel factor which assumes that all regions are permeable to salt, and F was set to 78 throughout the system. Although the models were parametrized to resemble interesting biomolecules, their simplicity precluded more than qualitative agreement with experiment. In the course of these investigations, the Brownian-dynamics method for determining rate constants was gradually extended and refined.Its application to quite sophisticated models is now both possible and warranted. Several steps have been taken in this direction in our SOD studies and are imminent for the inherently more complex cytochrome protein project. First, the influence of counterions and added salt was accurately determined by Monte Carlo simulation within the context of the simple shape, charge and dielectric model of SOD.24 These calculations provide data that can be used to test the basic assumptions in less computationally demanding models ( e.g. the Poisson- Boltzmann models discussed below), and also provide potentials that can be used directly in Brownian-dynamics simulations. The SOD work reported here is based on these Monte Carlo results except near the protein, where a highly realistic shape and charge distribution are used.However, in this inner region, salt and dielectric effects are again treated in the simplest approximation. While this model is superior to those used in our previous work, it is clear that a more consistent treatment of salt effects and some provision for the low polarizability of the protein interior are desirable. A promising vehicle for improved description of the electrostatic interaction is the finite difference form of the Poisson- Boltzmann eq~ation.”~‘”.’~ This numerical approach can handle a spatially varying dielectric response, detailed surfaces, complex charge distributions and the presence of counterions and added salt.Several applications of this method to biochemical systems have been reported, with very encouraging r e s ~ l t s . * ~ , ’ ~ - ~ ~ There are, however, several potential sources of error, both numerical and theoretical, in this method. In order to determine the magnitude and significance of such errors, we have applied the method to some simple systems for which analytic results or Monte Carlo simulation data, or both, are available.” For the cases studied, neglect of ion-ion correlation appears to be of more consequence than expected, although linearization of the equation has little additional impact, at least for physiological salt concentrations. There is considerable sensitivity to the spacing and length of the finite difference grid and to boundary value procedures.High accuracy demands enormous computing resources, but informed choices lead to an acceptable compromise. However, one may conclude that results with this method, even for detailed biochemical systems, are not only better than those previously available, but are likely quite good in an absolute sense. Careful electrostatic potential calculations for realistic models of bio- molecules are now in progress. They should provide data suitable for truly quantitative prediction of rate constants by Brownian dynamics trajectory techniques. In addition to the improvements in the electrostatic components of the models described above, several improvements in the dynamical details of the models seem desirable. The simulations to date have not included hydrodynamic interactions. The methods needed to include these interactions exist, although such inclusion requires220 Simulation Studies of Diflusion-controlled Reactions significantly more computer time.1431s Of particular interest will be the interplay between the electrostatic torques and the hydrodynamic coupling in the orientational docking of the cytochromes. On a more detailed level, the small value (2.75 A) required for the exclusion distance between 0, and the non-hydrogen atoms of SOD suggest that the final approach of 0, to the Cu2+ ions is facilitated by fluctuations in the internal structure of SOD. These fluctuations, and the motion of the water molecules and substrate in the channel, could be studied by molecular-dynamics simulations to provide improved reactive boundary conditions for the Brownian-dynamics simulations.9 Appendix A Let Pm,I denotz the probability of diffusing from a surface at radius b to an element S, of an interior surface S without escaping or encountering any other element Sl of S.The subscript co indicates that the medium is unbounded apart from the absorbing surface at S. Let p, denote the corresponding probability in the presence of a surrounding truncation surface at radius q, as described in the text. The probability Pm,, can be related to P, by branching arguments concerning the possible fates that would be available to trajectories if they were continued instead of being truncated at q. The argument given here is more direct than that used p r e v i o ~ s l y . ' ~ - ~ ~ Consider a trajectory started at b, and assume that there is a truncation surface at q.The probability of intersecting S, before intersecting Sl 0' # i) or q is P I . The probability of intersecting q first is 1 - P , where /3 = Cp,. Now suppose that a trajectory arriving at q were continued. The probability that the trajectory would hit S, before escaping or hitting any other element of S is k , ( S , ) / k , ( q ) . The total probability of hitting S, before any other element of S is therefore Here, the first term on the right is the probability of diffusing directly from b to Si and the second term is the probability of diffusing from b to q to S; if the trajectory is not truncated at q. Substituting the right-hand side of eqn (9) for kD(Sj) and solving for P.o,i yields eqn (10).Eqn (6) is obtained if S consists of only one element, S = Si. Appendix B Let P,,, denote the probability of diffusing from an dement S, of a surface S to an interior absorbing surface * rather than escaping. The subscript indicates that the medium is unbounded apart from the absorbing surface at *. Let P, denote the corresponding probability in the presence of a surrounding truncation surface at radius q. The probability P,,, can be related to P, by an argument similar to that given in Appendix A. Consider a trajectory started at S,, and assume that there is a truncation surface at q. The probability of hitting * rather than q is P,. The probability of hitting q rather than * is 1-P,. Now suppose that a trajectory arriving at q were continued.The probability that the trajectory would hit * rather than escaping is k / k , ( 4 ) . The total probability of hitting * is therefore Multiplying both sides by kD(Sj), summing over i and using eqn (7), one obtainsJ. A. McCammon, R. J. Bacquet, S. A. Allison and S. H. Northrup 22 1 Solving for k yields k = I The denominator can be rearranged slightly to give where we have used eqn ( 1 3 ) and ( 1 4 ) . J.A.M. is the recipient of a Camille and Henry Dreyfus Teacher-Scholar Award. S.A.A. is the recipient of an NSF Presidential Young Investigator Award and a Camille and Henry Dreyfus Grant for Newly Appointed Faculty in Chemistry. S.H.N. is the recipient of an NIH Research Career Development Award. This work has been supported in part by the NIH (at Houston and Tennessee Tech), by the Robert A.Welch Foundation and the Texas Advanced Technology Research Program (at Houston), and by the Donors for the Petroleum Research Fund as administered by the American Chemical Society (at Tennessee Tech). References 1 L. M. Raff and D. L. Thompson, in Theory of Chemical Reaction Dynamics, ed. M. Baer (CRC Press, Boca Raton, Florida, 1985), vol. 111, pp. 1-121. 2 D. F. Calef and J. M. Deutch, Annu. Rev. Phys. Chem., 1983, 34, 493. 3 0. G. Berg and P. H. von Hippel, Annu. Rev. Biophys. Biophys. Chem., 1985, 14, 131. 4 J. Keizer, Acc. Chem. Res., 1985, 18, 235. 5 E . Dickinson, Chem. SOC. Rev., 1985, 14, 421. 6 J. T. Hynes, Annu. Rev. Phys. Chem., 1985, 36, 573. 7 H. L. Friedman, A Course in Statistical Mechanics (Prentice Hall, Englewood Cliffs, NJ, 1985).8 J. A. McCammon, S. H. Northrup and S. A. Allison, J. Phys. Chem., 1986, 90, 3901. 9 J. A. McCammon and S . C. Harvey, D-ynamics of Proteins and Nucleic Acids (Cambridge University Press, Cambridge, 1987). 10 S. H. Northrup, S. A. Allison and J. A. McCammon, J. Chem. Phys., 1984, 80, 1517. 11 S. A. Allison, N. Srinivasan, J. A. McCammon and S. H. Northrup, J. Phys. Chem., 1984, 88, 6152. 12 S. A. Allison, S. H. Northrup and J. A. McCammon, J. Chem. Phys., 1985, 83, 2894. 13 S. H. Northrup, M. S. Curvin, S. A. Allison and J. A. McCammon, J. Chem. Phys., 1986, 84, 2196. 14 D. L. Ermak and J. A. McCammon, J. Chem. Phys., 1978, 69, 1352. 15 E. Dickinson, S. A. Allison and J. A. McCammon, J. Chem. Soc., Faraday Trans. 2, 1985, 81, 591. 16 G. Lamm, J. Chem. Phys., 1984,80, 2845. 17 S. A. Allison and J. A. McCammon, J. Phys. Chem., 1985, 89, 1072. 18 S. A. Allison, C. Ganti and J. A. McCammon, Biopolymers, 1985, 24, 1323. 19 G. Ganti, J. A. McCammon and S. A. Allison, J. Phys. Chem., 1985, 89, 3899. 20 J. A. Tainer, E. D. Getzoff, K. M. Beem, J. S. Richardson and D. C. Richardson, J. Mol. Biol., 1982, 21 A. Cudd and I. Fridovich, J. Bid. Chem., 1982, 257, 11443. 22 E. D. Getzoff, J. A. Tainer, P. K. Weiner, P. A. Kollman, J. S. Richardson and D. C. Richardson, 23 I. Klapper, R. Hagstrom, R. Fine, K. Sharp and B. Honig, Proteins: Structure, Function and Generics, 24 R. J. Bacquet and J. A. McCammon, Ann. N. Y. Acad. Sci., 1986, 482, 245. 25 S. H. Northrup, J. D. Smith, J. 0. Boles and J. C. L. Reynolds, J. Chem. Phys., 1986, 84, 5536. 26 S. H. Northrup, J. C. L. Reynolds, C. M. Miller, K. J. Forrest and J. 0. Boles, J. Am. Chem. Soc., 1986, 27 C. H. Kang, D. L. Brautigan, N. Osheroff and E. Margoliash, J. Biol. Chem., 1978, 253, 6502. 160, 181. Nature (London), 1983, 306, 287. 1986, 1, 47. 108, 8162.222 Simulation Studies of Diflusion-controlled Reactions 28 J. Warwicker and H. C. Watson, J. Mol. Biol., 1982, 157, 671. 29 N. K. Rogers and M. J. E. Sternberg, J. Mol. B i d , 1984, 174, 527. 30 N. K. Rogers, G. R. Moore and M. J. E. Sternberg, J. Mol. Biol., 1985, 182, 613. 31 J. Warwicker, D. Ollis, F. M. Richards and T. A. Steitz, J. Mol. Biof., 1985, 186, 645. 32 R. J. Bacquet, J. A. McCammon and S. A. Allison, to be published. Received 24th November, 1986

ISSN:0301-7249    

DOI:10.1039/DC9878300213

出版商:RSC    

年代:1987

数据来源: RSC

摘要:

GENERAL DISCUSSION Mr S. M. Pimblot (Oxford University) asked Dr Meakin: In your paper on computer simulations of diff usion-limited aggregation processes you consider the problem of the successive addition of single diffusing particles to a central aggregate. Particles are processed according to a first-come, first-served discipline. Have you considered the situation when this discipline breaks down: in other words, when several particles are competing in the aggregation process. This must eventually be the case, even with very dilute systems, since the length of the boundary of the aggregate increases without limit. Dr P. Meakin (du Pont, Wilmington, DE) replied: Most models for diffusion-limited aggregation (DLA) processes have been carried out in the zero-concentration limit in which particles are added, one at a time, to a growing cluster or aggregate of particles.This is appropriate for the simulations of random growth processes in which the growth probabilities are controlled by a scalar field obeying the Laplace equation. However, it is at best an approximation when the growth process is controlled by the diffusion of the material from which the aggregate is constructed. In principle, no more than one particle can be in the vicinity of the cluster at the same time, since the addition of one particle to the cluster while the second particle is in the vicinity does not lead (on average) to the same results as separate addition of two particles. In particular, one particle could add to the mouth of a ‘fjord’ in the structure while the second particle wiis isside the fjord and this would result in the second particle having a higher probability of being in the fjord than it would have had if the first particle had already been added to the mouth of the fjord. In practice, it would probably be possible to carry out DLA-like simulations in which a few particles were allowed to diffuse at the same time, and it is unlikely that this would change the asymptotic ( M + co where M is the cluster mass) behaviour.I am not aware of any simulations of this type. Both two-dimensional’-3 and three-dimensional3 simulations of growth from a ‘sea’ of many simultaneously diffusing particles have been carried out. It appears that under these conditions the average density of the growing aggregate (p) cannot become smaller than that of the ambient sea of particles (p,,) and there is a crossover from a fractal structure characteristic of the DLA model on short length scales to a uniform structure on long length scales.Simple arguments would indicate that the correlation length ( 6 ) associated with this crossover should depend on the particle density according to 5 x p ; ’ / a (1) where a is the codimensionality of DLA ( a = d - D, where d is the Euclidean dimension- ality and D is the fractal dimensionality). However, this is not entirely consistent with two-dimensional simulation results. Growth from a finite concentration of diffusing particles also has important implica- tions for the growth kinetics. Under conditions of low density and/or small cluster sizes (the DLA limit in which >> po) the growth of the mean cluster radius ( R ) can be described by the simple scaling relationship4” In the high-concentration, large-cluster-size limit the mean cluster radius grows linearly with time.(2) R x t ’ / ( D - d + 2 ) 223224 General Discussion 1 R. F. Voss, Phys. Rev. B, 1984, 30, 334. 2 R. F. Voss and M. Tomkeiwicz, J. Electrochem. SOC., 1985, 132, 371. 3 P. Meakin and J. M. Deutch, J. Chem. Phys., 1984, 80, 2115. 4 J. M. Deutch and P. Meakin, J. Chem. Phys., 1983, 78, 2093. 5 H . G. E. Hentschel, J. M. Deutch and P. Meakin, J. Chem. Phys., 1984, 81, 2496. Dr E. Dickinson (University of Leeds) asked Dr Meakin: The appropriate choice of trajectory starting and terminating positions is a matter which arises in many Brownian dynamics simulations of this type.Trajectories in which the two particles diffuse apart may only slowly become unbiased if the central one does not rotate or the interactions are of long-range (electrostatic, hydrodynamic etc.). In your simulations, is the termina- tion radius chosen on an ad hoc basis, or do you apply a more systematic approach? Dr P. Meakin (du Pont Wilmington, DE) replied: In the diffusion-limited aggregation (DLA) model developed by Witten and Sander’ we imagine that random walkers are launched from infinity (one at a time) and stick on contact with a growing cluster or aggregate of particles. Since the position at which a particle released from infinity first intersects a hyperspherical surface enclosing the absorbing cluster is unbiased, the random walk trajectories are started at a random position on a hyperspherical ‘launching’ surface which encloses the cluster.In the case of a lattice-based model, the random walk trajectory is started on the nearest lattice site to the randomly selected point on the launching surface. A small bias may be introduced at this point, but in practice good results are obtained if the radius of the launching surface is only slightly (a few lattice units) larger than the maximum radius of the cluster. In principle, the radius of the hyperspherical surface at which the random walk is terminated and a new random walk trajectory is started at a randomly selected position on the random walk were allowed to continue, it would eventually reach the launching surface with equal probabil- ity at all positions.In fact, if the random walk were to continue its return to the launching surface, the position of its return would be biased and the magnitude of this bias can be estimated by solving the Laplace equation (V2+ = 0) with the boundary conditions 4 = 0 on the launching surface and at infinity and 4 = 1 at a point on the termination surface.2 Although the selection of a termination radius large enough to reduce this bias to a negligible magnitude would guarantee that the termination of the random walk trajectories introduced no bias into the DLA simulation, this is not necessary since this bias is important only if there is a significant bias in the probabilities with which the random walkers reach different positions on the termination surface.There is a bias in the probabilities with which different positions on the termination surface are reached because a random walker started from a position on the launching surface which is close to a region occupied by the cluster has a smaller probability of reaching the termination surface before contacting the growing cluster than a random walker started from a position which is further from the cluster. Because of the complex geometry of DLA clusters, the magnitude of this bias is difficult to estimate. The biases discussed above are largest for two-dimensional simulations. However, in this case the effect of this bias is reduced because of the relatively low probability that a particle launched near to the cluster will contact the termination surface before contacting the cluster. In practice, reliable results can be obtained from small-scale (G 10 000 sites or particles) two-dimensional simulations using a termination radius of ca. 3 R,,, where R,,, is the maximum radius of the cluster, and even smaller termination radii seem to give satisfactory results for d a 3 3 ’ 4 where d is the Euclidean dimensionality of the embedding space or lattice.In more recent large-scale two- and three-dimensional simulation^,^-^ a termination radius of lOOR,,, has been routinely used, and in some cases termination radii of 1000R,,, have been used. However, there is no evidence for measurable bias if a termination radius of lOR,,, is used.General Discussion 225 Since these results are empirical, then each model must be treated on an individual basis.For example, Lyklema and Evertsz have recently shown that the radius of the termination surface is more critical in a diffusion-limited chain growth model closely related to DLA than in DLA itself. In this model only one site (the end of the growing chain) can adsorb random walkers, and it is consequently natural that the biases discussed above should be larger. Even for this model a termination radius of lOR,,, seems to be satisfactory . 1 T. A. Witten and L M. Sander, Phys. Rev. Lett., 1981, 47, 1400. 2 W. D. Jackson Classical Electrodynamics (John Wiley, New York, 1962). 3 P. Meakin, Phys. Rev. A, 1983 27, 604. 4 P. Meakin, Phys. Reu., 1983, 27, 1485. 5 P. Meakin, J. Phys. A, 1985, 18, L661 6 R. C. Ball and R.M. Brady, J. Phys. A, 1985, 18, L809. 7 P. Meakin, R. C. Ball, P. Ramanlal and L. M. Sander, Phys. Rev. A, in press. 8 R. C. Ball, personal communication. 9 J. W. Lyklema and C. Evertsz, J. Phys. A, 1986, 19, L895. 10 P. Meakin, unpublished work. Prof. B. J. Ackerson (R.S.R.E. Malvern) turned to Dr Meakin: The anisotropy shown in fig. 4 is produced without noise reduction or other bias than that introduced by the lattice. What is the nature of this bias on a square lattice? Is it simply due to the fact that the walker must take 40% more steps to move directly (not randomly) a given distance diagonally compared to the same distance in a horizontal or vertical direction? Yet the average rms displacement is isotropic as the number of random steps becomes large.Dr P. Meakin ( d u Pont, Wilrnington, DE) replied: The anisotropy associated with the structure of large square-lattice DLA clusters appears to be a result of the fact that growth is restricted to four directions in the two-dimensional space. This idea is supported by the observation that ordinary square-lattice DLA and a ‘semi-lattice’ model in which circular particles follow off -lattice random-walk trajectories until they contact the growing cluster but are allowed to add only along one of four directions corresponding to the axes of a square lattice give esssentially identical results’ (at least for clusters containing up to lo5 particles). In this model the off-lattice random walk is stopped at the point where the random walker first touches a particle in the cluster, and the random walker is then rotated about the centre of the contacted particle until the vector from the centre of the contacted particle to the centre of the contacting particle (random walker) is directed along the nearest lattice axis.Similarly, simulations on hexagonal lattices in which either the growth process or the steps of the random walk are restricted to three of the six possible directions (each at mutual angles of 27r/3) indicate that the restriction of the growth direction is much more important than restriction of the random-walk step directions. The idea that 4 2 more growth steps are required to grow by a given distance in the diagonal direction than to grow by the same distance along the lattice axes could provide an explanation of the formation of a diamond-like shape such as is seen in clusters containing ca.lo5 However, as is shown in fig. 3, of our paper the asymptotic shape of square-lattice DLA is not a diamond. Progress has been made towards understanding the asymptotic shape of square-lattice DLA clusters (and other DLA models) in terms of the solutions to the Laplace equation for various idealized shapes polygon^,^,^ stars6” and a hierarchy of wedges8*’) with absorbing boundary conditions at the edges of the shapes and a fixed value for the harmonic field at infinity. However, the development of a comprehensive theory for the asymptotic structure of DLA clusters is not yet complete.226 General Discussion 1 P. Meakin, Phys. Rev. A, 1986, 33, 3371. 2 P. Meakin, J. Phys. A, 1985, 18, L661. 3 R.C. Ball and R. M. Brady, J. Phys. A, 1985, 18, L809. 4 L. Turkevich and H. Scher, Phys. Rev. Lett., 1985, 55, 1026. 5 R. C. Ball, R. M. Brady, G. Rossi and B. R. Thompson, Phys. Rev. Lerr., 1985, 55, 1406. 6 R. C. Ball, Physica A, 1986, 140, 62. 7 F. Family and H. G. E. Hentschel, Faraday Discuss. Chern. SOC., 1987, 83, 139. 8 T. C. Halsey, P. Meakin and 1. Procaccia, Phys. Rev. Lerr., 1986, 56, 854. 9 T. C. Halsey, Europhys. Left., to be published. Dr J. C. Earnshaw (Queen’s University Belfast) said: Prof. Jullien describes the use of the hierarchical model as a computationally rapid tool to study cluster-cluster aggregation. This use has been justified by direct comparisons of results of the model with those of Monte-Carlo numerical simulations.’ However, these comparisons were restricted to considerations of the fractal dimensions of the clusters.In the present paper the hierarchical model is used to investigate other properties of the clusters: intrinsic anisotropy, restructuring etc. To what extent can the use of the model to study such other parameters in cluster-cluster aggregation by justified at present? 1 R. Botet, R. Jullien and M. Kolb, J. Phys. A, 1984, 17, L75. Prof. R. Jullien (Orsay, France) replied: The hierarchical model can be considered as the limit of the cluster-cluster model in a box when the size distribution of clusters becomes more and more peaked. If it is well established that this model leads to the same fractal dimension as the model in the box, it is, however, not completely certain that all the other properties are also the same.It is generally believed that all the static properties are well reproduced. In the case of the anisotropy properties of the cluster- cluster aggregates, we have given’ some simple renormalization group arguments which directly relate both the fractal dimension and the anisotropy ratio to the mean relative penetration width when two clusters collide. Consequently, the anisotropy ratio, as the fractal dimension, must not depend on cluster polydispersity. I think that the same general conclusion is true for restructuring effects. 1 R. Botet and R. Jullien, J. Phys. A, 1986, 19, L907. Dr J. C. Earnshaw and Mr D. J. Robinson (Queen’s University, Belfast) added: In our experimental observations of two-dimensional aggregation in interfacial colloidal sys- tems we have seen clear evidence for the existence of rather long-ranged interactions between clusters.These interactions affect the aggregate structures observed. Two photomicrographs illustrate these points. The experimental system is the same in both cases: 1 p m diameter polystyrene microspheres trapped at the liquid-air interface of a 0.05 mol dm-3 NaCl solution. The particle size provides an indication of the distance scale. The Debye screening length in these systems is rather small (ca. 1 nm). In plate 1 a group of small clusters can be seen within an area otherwise free of particles. This group, in which the constituent clusters are separated by a few particle diameters, remained quite stable in configuration over periods of weeks.This clearly suggests the existence of a very deep (and narrow) secondary minimum in the cluster- cluster interaction. This entire group of clusters underwent Brownian motion as an ensemble, moving within the bounds set by the clear area, The group never approached the surrounding clusters, showing the presence of a very long-ranged repulsive com- ponent to the inter-cluster interaction. The potential minimum within which the group moved appeared to be rather shallow in form, allowing quite free translational and rotational diffusion. The distance scale characterising the inter-cluster interactions seems surprisingly large (much greater than the Debye length). These interactions have consequences for the structure of the aggregates formed by cluster-cluster aggregation.Inspection of plate 1 shows several points where clustersFaraday Discuss. Chem. SOC., 1987, Vol. 83 Plate 1 Plate 1. Aggregation of 1 pm spheres at a water-air surface. The arrows indicate the stable group of clusters and some of the narrow gaps between clusters discussed in the text. Plate 2. A different stage of aggregation in a similar system to that of plate 1. The tendency for clusters to approach tip-to-tip is clear. (Facing p. 226)General Discussion 227 are separated by one or two particle diameters. These gaps seemed quite stable, persisting over periods of weeks. We did, however, see some evidence that clusters could rotate about these gaps, leading to restructuring of the overall ensemble. Even after such restructuring, the clusters did not always come into primary minimum contact, but remained separated by gaps of different configuration. Plate 2 shows a different stage of aggregation in which clusters are forming ‘superclusters’ with a larger-scale fractal appearance, but are not yet bound. Indeed, further observation showed the superclusters coming into contact and thereafter remaining relatively stable.We interpret the formation of these structures (and similar ones seen in other experiments) as being driven by the long-ranged repulsion between clusters, which makes it more probable that clusters will approach each other ‘tip to tip’. In summary, during the course of an experiment, clusters grow and become oriented so as to approach each other ‘tip to tip’.They come into secondary minimum contact and restructuring may continue. The nature of the ensemble thus formed appears to modify the potential, permitting a transition to primary minimum contact. This last stage is apparent in our observations of an ultimately fully connected network. In Prof. Jullien’s modified (tip-to-tip) simulations with the hierarchical model, has he seen evidence for such long-ranged correlations as appear in our fig. 2? Incorporation of inter-cluster interactions such as implied by the present observations (as opposed to inter-particle forces) would provide a most interesting extension of cluster-cluster aggregation simulations. Prof. Jullien answered: The tip-to-tip model is a rather oversimplifying model which tries to take care of simple cluster interactions such as the ones appearing in the case of electrical polarizability.However, in your experiments, I think that much more complex behaviour occurs which might be the consequence of some hydrodynamic effects. This produces indirect cluster-cluster interactions, via the fluid. In such case, a more specific (but certainly more complex) model must be built. Dr E. Dickinson (University of Leeds) said: In many experimental examples of colloidal aggregation, the aggregates are subject to viscous shear forces in addition to Brownian motion, but only when the aggregates exceed a certain size is their motion appreciably affected by the imposed flow field. It seems likely that the structure of aggregates formed in shear flow would be different from those formed in quiescent media.Would it be possible to check this point in a simulation of cluster-cluster aggregation ? Prof. Jullien replied: Our two-dimensional calculations with Paul Meakin,’ which show that rotations of colliding clusters about their sticking points can affect the short-range structure within aggregates, give a first rough idea of what could be the effect of viscous shear forces. It is quite obvious that such rotations will certainly be more frequent in presence of shear forces. We will soon extend these calculations to three dimensions. This will, however, give only a qualitative answer. Other calculations which take care more specifically of hydrodynamic effects should be done. 1 P. Meakin and R. Jullien, J. Phjx (Paris), 1985, 46, 1543.Dr D. S. Horne (Hannah Research Institute, Ayr) asked: While the fractal concept has been applied to many growth processes with success, am I correct in saying that any relation between the growth mechanism and the internal structure of the aggregate has still to be established? In simulations you have demonstrated that it is possible to produce a transition from cluster-cluster to particle-cluster dynamics and therefore a change in fractal dimension by varying the exponent, a, in the expression for the228 Gen era1 Discussion diffusivity of the particle cluster. In other work you have also shown a concentration dependence of the fractal dimension in that higher concentrations of particles lead to more compact aggregates. Are there any physical factors which might influence the value of the exponent, a? Would you like to summarize the present state of knowledge of how the fractal dimension is related to the physics of the aggregation process? Dr Meakin' has stated that a size-dependent sticking probability has a similar effect on the size distribution as a size-dependent diffusion coefficient.Does this also produce a crossover from cluster-cluster to particle-cluster dynamics? 1 P. Meakin, T. Vicsek and F. Family, Phys. Rev. B, 1985, 31, 564. Prof. Jullien replied: The relation between the w exponent of the Smoluchowski equation and the fractal dimension' 2w = (Y + ( d - d , ) / D where a is the kinetic exponent defined in our paper and d , the fractal dimension of the cluster trajectory ( d , = 2 in the Brownian case), gives a quantitative relation between dynamic and static properties, within the general assumptions of the cluster-cluster model.It is true that many details of specific aggregation mechanisms, which are not taken into account in the existing models, could affect this relation. For example, the introduction of a size-dependent sticking probability, via an extra exponent which must be added to the right-hand side of the equation, produces a shift in the location of the crossover from cluster-cluster to particle-cluster. 1 R. Botet and R. Jullien J. Phys. A, 1984, 17, 2517. Dr D. Weitz (Exxon, Annandale, NJ) remarked: In your paper, you discuss a simulation in which you allow the clusters to restructure themselves by forming loops. This is done at each stage of the growth, so that loops are formed on all length scales.One might also consider a somewhat different process of restructuring, in which the rigidity of a cluster depends on its size. As the cluster size increases, the chain length connecting the ends of the cluster also increases, making it much easier to bend. Thus one might expect the restructuring to occur primarily at longer length scales. We believe that this effect is observed in experiments in which fractal gold aggregates are subjected to shear. Scattering measurements suggest that the clusters are restructured only on long length scales, while the structure at shorter length scales is unchanged. This is presumably a more difficult process to simulate on the computer, but nevertheless I would encourage you to attempt it, as we believe that this restructuring process is rather common.Dr J. Gregory (University College, London) said: Computer simulations of cluster- cluster aggregation give fractal dimensions in good agreement with those found for aggregates of colloidal gold.' An alternative way of showing the fractal character of aggregates is to determine their effective density as a function of aggregate size. Such experiments were carried out in Japan in the 1960s by Tambo and his co-workers (see Tambo and Watanabe2 for a review of this work in English). The systems studied were of a very practical nature: clay suspensions flocculated with aluminium salts and even activated sludge from a sewage treatment plant. Individual aggregates ('floes') were allowed to settle in a glass tube and the settling rate was measured.Floc size was determined from photographs taken during sedimentation and, from the settling rate and an assumed form of the drag coefficient, the effective (buoyant) density of the flocs was determined. For floc diameters over the approximate range 20-200 pm, log-log plots of density against diameter gave reasonable straight lines with negative slopes, showing a marked decrease in density with size.General Discussion 229 With fully destabilized clay particles, the slope of the lines was of the order of -1.4, corresponding to a fractal dimension, D ==: 1.6, which is not greatly different from the values of D = 1.75 found for gcld aggregates. When destabilization was incomplete, rather lower slopes were found [ca.-1.0 ( D == 2.0)], indicating denser flocs. This is the trend found for aggregates of gold colloids formed under ‘reaction-limited’ condition^.^ In Tambo’s work flocs were produced in a stirred vessel (orthokinetic flocculation), and perhaps we should not attach too much significance to the apparent agreement with computer simulations and experiments where diffusion is the dominant transport mechanism. Nevertheless, results of simple density measurements on systems of great practical interest indicate that current ideas on aggregate structure may be of very wide applicability. Further experiments along these lines, perhaps using polymeric flocculants, would appear to be well worthwhile. 1 D. A. Weitz and M. Oliveria, Phys. Rev. Lert., 1984, 52, 1433.2 N. Tambo and Y. Watanabe, Water Research, 1979, 13, 409. 3 D. A. Weitz, J. S. Huang, M. Y. Lin and J. Sung, Phys. Rev. Lett., 1985, 54, 1416. Prof. B. U. Felderhof (R. W. T. H., Aachen, Federal Republic of Germany) addressed Dr Adler: I would like to mention that some years ago together with J. A. McCammon and J. M. Deutch I made a calculation of the friction coefficient of a spherical shell uniformly covered with beads of friction.’ We found that upon random removal of up to 70% of the beads the friction coefficient hardly changes. This shows that in such a structure the inner flow field is very effectively annihilated and that the friction coefficient is determined mainly by the outer layer of points of friction. The same is probably true for the fractal flocs considered by you.1 J. A. McCammon, J. M. Deutch and B. U. Felderhof, Biupolymers, 1975, 14, 2613. Dr P. R.I. Adler (Meudon, France) replied: Thank you for bringing this paper to my attention; the results that I present here are in total agreement with it. It appears that the hydrodynamic drag is insensitive to the precise structure of the floc because of the existence of dead fluid regions which fill the gaps. This is true unless one deals with extreme cases where these motionless regions are expected to vanish; this might happen for values of the fractal dimension close to 1 that I suggested should be investigated more thoroughly at the end of my paper. Dr A. J. Kuin (Unilever, Vlaardingen, The Netherlands) said: (1) In fig. 2 of your paper the flocs shown are small and there is hardly any self-similarity, geometrically speaking, over various length scales.Despite this you neglect the floc size in your derivation of eqn. (1 1). How can one explain this? (2) The fractal dimension of set X is defined as In n(&j D = -1im - with D The Hausdorff dimension and n ( s ) the number of sets of diameter E to cover x.’ That is to say it is a limiting process. However, if there are only few particles the limit does not exist, as the addition of particles cannot be approximated by a continuous function. How did you determine the fractal dimension? F+O In& 1 H. 0. Peitgen and P. H. Richter, The Beauty of Fractals (Springer Verlag, Berlin, 1986). Dr P. M. Adler (Meudon, France) replied: In answer to the first question, I draw attention to fig.1 of my paper, which indicates statistical self-similarity; note that even for small N the slopes of the solid lines are very close to the values accepted in the230 General Discussion literature. Moreover, fig. 4 shows that the hydrodynamic results relative to ZI" only depend on the ratio R,/h, hence eqn (11); it may be useful to recall that the points in fig. 4 are obtained for various channel widths (ranging from 9 to 34 in arbitrary units.). In answer to the second question, eqn (2) is the usual rule according to which the weight of an object of radius R scales as RDf [cJ my ref. (1) p. 401, where Df is the dimension of the object. The definition that you are giving is one of the possible ways to introduce the fractal dimension (but definitely distinct from the Hausdorff -Besicovitch dimension) of mathematical entities whose smallest scale is zero, which is not the case here.Dr P. N. Pusey (R.S.R.E., Malvern) (partly communicated) I wish to address Dr Weitz: For fractal aggregates one might expect the hydrodynamic radius Rh, determined from the diffusion coefficient, to be proportional to the radius of gyration R,, i.e. Rh = PR,. The value p = 1.75 which you quote for DLCA clusters seems surprisingly large. For example, for Gaussian polymer coils (fractal dimension d,- = 2) Kirkwood's theory, which assumes an Oseen hydrodynamic interaction between beads, gives p = 0.665. In the case of swollen coils (df= 1.67, close to the ca. 1.77 of diffusion-limited clusters) the theoretical value of p is smaller still.I would not expect such a large difference in /3 values for polymers and aggregates with similar fractal dimensions, even though their structures are quite different. If correct it implies that the volume of solvent trapped hydrodynamically by an aggregate is 15-20 times greater than that associated with a polymer coil of similar fractal dimension. Is it possible that there is an error in your quoted value of p ? If so, to what extent are the calculations of fig. 5 affected? Since the Discussion I discovered a recent paper by Wiltzius,' who has measured p directly for reaction-limited clusters and obtained p = 0.72. While this value may be changed when proper account is taken of the distribution of cluster size, it is substantially smaller than 1.75.1 P. Wiltzius, Phys. Rev. Lett., 1987, 58, 710. Dr D. Weitz (Exxon, Annandale, N J ) replied to Dr Pusey: The value of p = 1.75 was based on some preliminary, and unpublished, calculations by Chen, Deutch and Meakin, and was the only value of /3 available when this paper was first written. Following your comment at the Discussion, they rechecked the calculation of p, and did indeed find an error of a factor of 2. We have therefore corrected the text of the paper and used a value of @ =0.8. This value is also now consistent with a calculation of p by Hess, Frisch and Klein.' The effect of a smaller p is actually quite substantial. The magnitude of Rh relative to R, determines when rotational effects start to play an important role in the light scattering.Furthermore, since 0 - 1/ R i , the actual decay rate due to rotational diffusion is strongly dependent on p. Using this smaller value, we find much better agreement with our experimental data over a wide range of cluster sizes. Furthermore, our results using this lower value of p suggest that rotational diffusion may actually be the dominant decay rate of the measured autocorrelation function, and cannot be neglected in the correct determination of Rh at high kR,. The value of p = 0.8 is also consistent with the measurement by Wiltzius that you cite, provided his result is corrected to account properly for the effects of the distribution function. The static and dynamic scattering that Wiltzius uses to measure Rg and R h , respectively, probe different moments of the cluster mass distribution, and thus the value of p he cites refers to clusters of different masses. However, we can estimate the effects of the distribution function and correct Wiltzius' result if we assume that R, is obtained from the small-k expansion of the static light scattering, I ( k ) = K ( 1 - k'Ri), where K is a constant. If we further approximate the cluster mass distribution by a power law,General Discussion 23 1 N ( M ) K M - ' , with a sharp cutoff at a maximum mass M,, we can account for the effects of the distribution function on the measured values of the hydrodynamic radius, Rh, and radius of gyration, R,, and thus correct Wiltzius' result: Thus for Wiltzius' experiment clusters of the same mass have p = 1.23Rh/ Rg = 0.9, using 7 = 1.8 and d f = 2.1.We thank you for your comment pointing out the error in our calculation. 1 W. Hess, H. L. Frisch and R. Klein, Z. Ph-ys., Ted B, 1987, 64, 65. There is an error of a factor of 2 in their eqn (14). Prof. V. Degiorgio (Pavia, Italy) said: Concerning Dr Weitz's discussion of the effect of rotational diffusion on dynamic light-scattering measurements, I would like to ask whether there are calculations or measurements of the rotational diffusion coefficient of a fractal object. Is there any dependence on the fractal dimension? Dr. Weitz replied: To my knowledge, there is no calculation or experimental measure- ment explicitly of the rotational diffusion constant for a fractal object. We have used the translational diffusion constant to determine the rotational diffusion constant through the hydrodynamic radius, i.e.we assume that the rotational diffusion constant of a cluster of mass M is given by @ ( M ) = kT/87rqRh3, where q is the fluid viscosity and the hydrodynamic radius, R h , is determined by the translational diffusion constant, D( M ) = kT/6rr7Rh. The hydrodynamic radius has been calculated and measured for fractal aggregates, and does depend on the fractal dimension, as discussed above in the comment by Peter Pusey and my response. As such, O ( M ) will also depend on the fractal dimension. There is one additional point that should also be made: both D and 0 depend on the specific direction of the motion relative to the orientation of the cluster, since the overall shape of the cluster is assymetric.This is not included in any of our calculations, nor is it included in any discussion of the translational diffusion constant. in general, we expect an axial assymetry of roughly a factor of two at most. Thus this probably will not strongly affect any measure of the translational diffusion constant, but may have a stronger affect on 0 owing to the stronger dependence on Rh . Dr. D. S. Horne (Hannah Research Institute, Ayr) remarked: I should like to propose an alternative explanation, other than rotational effects, for the dependence you observe for the normalized first cumulant on the wavevector, k. The graphs in your fig. 4 confirm that translational effects and the So term are dominant only when kR, < 1. When kR, >> 1, as you have here rotational effects become important and the first cumulant, r, can be rewritten as a linear function of k' plus a positive intercept whose value reflects the gagnitude of the rotational contribution.Your data, replotted as r against k 2 , yield an excellent straight line passing through the origin, indicating, as you concur, the domin- ance of translational effects. The r us. k' plot emphasizes the effects of data obtained at high k values, but this is precisely what is required for detecting rotational contribu- tions. The absence of such contributions suggests a reappraisal of the data. If the effective diffusion coefficient, the first cumulant normalized by k 2 (i.e. r / k 2 ) is plotted as a function of k', then the following is observed. First a shallow decrease in Deff as k' is decreased from its highest values, followed by a sharp drop in the value of Deff at the lowest wavevectors, corresponding to angles of 50 and 20".This is similar to behaviour reported by Horne' in diffusion coefficient measurements of casein micelles and ascribed there to either the presence of a few large particles (fat globules) or to some dust contamination. Horne's calculations indicate that such particles need be no232 General Discussion larger than a factor of three greater than some size characteristic of the main scatterers, but produce gross effects, similar to those seen here, on diffusion coefficients measured at 30", while leaving wide-angle measurements virtually unaffected, this when present at number fraction levels no greater than 0.1% of the main scattering population.Calculation may also show that similar effects might obtain in a polydisperse system with a long tail extending to large aggregates. 1 D. S. Horne, J. Colloid Interface Sci., 1984, 98, 537. Dr Weitz responded to Dr Horne as follows: I do not feel that your alternative explanation for the k dependence of the effective diffusion constant is applicable to our experiments. Replotting our data as you suggest, either as Deff = r/ k * or as r against k2 changes neither the behaviour of the data nor the agreement with our theoretical calculation. We have chosen to normalize r by k2 to emphasize the relatively small differences from the behaviour that would be observed for purely translational diffusion. We have plotted it against k on a logarithmic scale to show that at low enough k, where rotational diffusion does not contribute, we do indeed obtain the k-independent behaviour expected.You are suggesting that the observed k dependence is due to the form of the distribution of scatterers rather than from the effects of rotational diffusion. To investi- gate this possibility, let us write the expression for the first cumulant in the absence of the effects of rotational diffusion. Then This is simply the average of the diffusion constants, D ( M ) over the distribution of scatterers, N ( M ) , weighted by their scattering intensities, M2SM( k ) . The cluster mass distribution produced by diffusion-limited cluster aggregation has been measured using electron microscopy.' It was found to be reasonably well approxi- mated by the Smoluchowski distribution for a constant kernel,' which gives a flat, or constant distribution out to some characteristic cluster size, followed by a relatively sharp cutoff.This is the cluster mass distribution that was actually used in the calculation shown in fig. 5. As stated in the paper, we also performed the calculation without including the contributions of rotational diffusion. In that case the behaviour agrees with the predictions obtained using the above equation, and r / k ' does have some k dependence due solely to the polydisperse distribution. However, this represents < ca. 15% of the total change in r / k 2 with k. We have replotted fig. 5 in our text to show by the bottom line the calculation without the inclusion of rotational diffusion effects, demonstrating the very small contribution of the cluster mass distribution to the measured k dependence.Consider now the effect of adding a small number of large particles to the distribution. In the paper cited, Horne considered fat globules of radius ca. 300 nm, whose scattering intmsity was determized by the Rayleigh form factor for solid stheies of radius, R * M so that S,(k) = {3[sin ( k R ) - kR cos ( k R ) ] } ' . For particles of the size con- sidered by Horne, the Rayleigh form factor is very strongly biased in the forward direction, with the scattering intensity increased by as much as three orders of magnitude at small scattering angles (20") as compared to larger angles (>90"). This accounts for the observations reported by Horne: at small angles these fat globules contribute strongly to the scattering and are heavily weighted in the distribution, effectively decreasing the measured diffusion constant.However, at larger scattering angles these larger particles are essentially invisible to the scattering, and thus do not contribute in the average, making the measured diffusion constant larger. Thus the very large effect of a very small number of particles is due solely to the particular choice of size and scattering form factor, and cannot be expected to be more generally observed.General Discussion 233 In our experiments we do not have fat globules of the appropriate size to produce a similar effect. Furthermore, our results are not influenced by dust, because the scattering from the fractal gold aggregates is sufficiently intense to dominate any parasitic scattering.Indeed, for aggregates of the size used to obtain the data reported in fig. 5 , we find the same behaviour even if the solutions are completely unfiltered for dust. Thus the only remaining possibility we need consider is that the distribution of fractal clusters itself has a small number of very large clusters. While these were not observed in the electron-microscopy experiments, we cannot rule this possibility out. However, larger fractal clusters will not produce the effect you suggest, because rather than the Rayleigh form factor for spheres, for fractal clusters, S,( k ) = 1 when kR, < 1, and S,( k ) = 1/ Mkdf when kR,> 1. These larger clusters will satisfy kR,> 1 for all k used, and so will contribute the same way to the averaging for all k. Indeed, the only cause of a k-dependence in the averaging is the change in S,( k ) at the crossover from kR, < 1 to kR,> 1, which affects the relative weighting of th; clusters in their contribution to the average.Thus adding some large clusters with kR,> 1 will not produce any additional k-dependence in the measured diffusion constant, but instead will actually cause the total variation with k to be less. Thus, to conclude, the alternative explanation you suggest does not seem to be applicable whatsoever in these experiments. Instead, one fully expects a contribution to the dynamic light scattering from rotations from the highly inhomogeneous structure of the fractal aggregates, and the data seem to support this conclusion.1 D. A. Weitz and M. Y. Lin, Phys. Rev. Lett., 1986, 57, 2031. 2 R. J. Cohen and G. Benedek, J. Phys. Cliem., 1982, 86, 3696. Prof. B. U. Felderhof (R. W.T.H. Aachen, Federal Republic of Germany) said: In regard to fig. 1 in your paper you say that the change in the optical absorption spectrum as you go from isolated gold spheres to aggregated clusters is dominated by short-range correlations and may be understood from the plasma resonance properties of dimers. I would think that actualiy many-sphere correlations are important and that one should consider plasma modes running along chains of particles in the cluster. Especially the peak in the absorption spectrum at long wavelengths should be understood in this way.The concept of plasma modes is an electrostatic one and I wonder whether in your calculation of the local field effects you have used an electrostatic approximation for the retarded interactions. If not, do you have an idea whether the electrostatic approxi- mation is a good one, like it is in the theory of the dielectric constant of polarizable liquids? Dr I). Weitz replied: You are quite correct in suggesting that a full determination of the behaviour of the plasma resonance must include the effects of many spheres. However, while the inclusion of more spheres will affect the quantitive behaviour of the plasma resonance, and in particular the frequency of the peak in the absorption, the qualitative features and the essential physical picture will he the same as fnr a dimer.Our primary goal was not to calculate the detailed behaviour of the plasma resonance; thus our reference to a dimer was merely meant by way of an example of a simple geometry which reflects the essential physical picture. As to the inclusion of retardation effects: they can enter in two ways into our picture, because of the way we have split up the calculation of the local fields. The first way is in the determination of the correct polarizability of each sphere. This includes the higher-order multipoles and accounts for the shift in frequency of the plama resonance due to the nearest neighbours. In our calculations' we used a polarizability obtained from an exact solution for electromag- netic-wave scattering from an infinite chain of gold particles (taken to be cyclindrical in shape), whose spacing was chosen so that the calculated adsorption peaks matched those measured experimentally.Since the solution is exact, retardation effects are234 General Discussion I! 0 -0.5L 7 Fig. 1. The second cumulant Q (as measured in dynamic light-scattering experiments) as a function of the exponent T, in the Rayleigh scattering approximation. (-) Exponential cut-off, Cv cc M-' exp -(M/C,M,); (- - -) hard cut-off, C , x M - ' M C2M,. The measured value of Q with some indication of the uncertainty is also shown. naturally included. The second way retardation must be included is in the self-consistent calculation of the local fields, which will reflect the effect of the fractal correlations within the cluster.Since these calculations include dipole interactions which are long- range, and since the clusters can be large compared to a wavelength, retardation must be included. This was done in our calculations. However, we do not believe that the effects of retardation are very important in determining the scattering structure factor of the cluster. We base this conclusion on the observation that the fluctuations of the local field about the mean are correlated only over a distance of order a sing'le sphere diameter which is very small compared to a wavelength. The range of these correlations is not likely to be strongly dependent on whether or not retardation is included. 1 Z. Chen, P. Sheng, D. A. Weitz, H. M. Lindsay, M. Y. Lin and P. Meakin, Optical Properfies of Aggregate Dr J.G. Rarity (R.S.R.E., Mafvern) turned to Dr Weitz: Your paper quotes a power-law cluster mass distribution with exponent T =. 1.5 for the case of reaction-limited aggregation (RLCA). We have been studying salt-induced aggregation of aqueous polystyrene colloids (20 nm in diameter) using low-angle static and dynamic light scattering, At low angles the scattered intensity (in the Rayleigh limit) is directly proportional to the weight-average aggregate mass M,, and a log-log plot of intensity against low-angle hydrodynamic radius, Rh (measured using photon correlation spectros- copy, PCS), leads to an estimate of the hydrodynamic fractal dimension.' More recently we have made detailed measurements of the second cumulant Q from low-angle PCS measurements.Using model power-law size distributions we have calcu- lated Q as a function of exponent T, and show our results graphically below. Clearly Q is very sensitive to T above t = 1.5. Although measurement error is large we can estimate r from Q (as shown in fig. l ) , r = 1.9k 0.1. Measurements are, however, confined to relatively early stages of the aggregation ( M , s 200 m,, rn, being the monomer mass) to ensure that the Rayleigh approximation holds for all aggregates. Clusters, to be published.General Discussion 235 0 " 1 1.25 1.5 1.75 2 7 Fig. 2. The second cumulant, Q, as a function of T for a sharp cut-off in cluster mass distribution, including clusters of mass M = 0 (-) and starting at M = 1 (- - -). An exponential cut-off starting at M = 1 is also shown (--.--.).Note that the large singularity in Q is removed when the correct bounds of the integral are used, making a determination of T from the measured Q very difficult. To confirm that this value of T corresponds to the asymptotic value (i.e., that scaling sets in very early on during aggregation) we have also estimated T using the methods suggested by Martin and Leyvraz' based on measurement of high-angle PCS linewidths in the later stages of the aggregation. Preliminary measurements indicate 1.8 d 7 6 2 in agreement with the value implied' by results published by Schaefer and Martin3 (in silica colloids) but disagreeing with your results for linewidth as a function of angle (fig. 5 of our paper). Does this disagreement indicate some non-universality of T in RLCA in contrast with recent theoretical studies ?4 1 P.N. Pusey and J. G. Rarity, Mof. Phys., in press; see also J. G. Rarity and P. N. Pusey, in On Growth and Form, ed. H. E. Stanley and N Ostrowski (Nijhoff, Dordrecht, Netherlands, 1986), pp. 218-221. 2 J. E. Martin and F. Leyvraz, Phys. Rev. A, 1986, 34, 2346. 3 J. E. Martin and D. W. Schaefer, Phys. Rev. Lett., 1984, 53, 2457. 4 R. C. Ball, D. A. Weitz, T. A. Witten and F. Leyvraz, Phys. Rev. Lett., 1987, 58, 274. Dr Weitz replied to D r Rarity: The value of T quoted for reaction-limited cluster aggregation was obtained from measurements on gold colloids using electron micro- scopy.' It is thus not sensitive to any ambiguities in interpretation that might potentially plague light scattering measurements, although it does suffer from a largqr statistical error, due to counting limitations, than scattering measurements.The measurement of Q at very small kR, is potentially a good way to determine T, as you point out. However, even here I should like to caution you on the calculation of the expected results. In the accompanying fig. 2 we show the results of our calculation of Q. If we use a power-law distribution, N( M ) MPT, up to a cutoff, M,, and perform all the necessary intergrations from M = 0 to M = M,, we obtain results identical to those presented in your figure. This is shown by the solid curve in our figure. However, there are no clusters with M = 0. If we perform the same integrals from M = 1 to M = M,, the results are substantially different.In particular, the strong divergence as T approaches 3 - 2/df= 2 is removed, and the results are quite sensitive to the value of M,, and therefore the full details of the experimental situation should really be considered. To23 6 General Discussion do so we have been more precise and replaced the integral by a sum over the discrete cluster masses, allowing us to include a more realistic scattering structure factor. We show two examples of the results in the figure. The dot-dash line is the calculation for a power-law distribution with a sharp cutoff, while the dashed line is for an exponential cut-off, N( M ) CK M-' exp (- M/M,). The results shown are calculated for M , = 1000, k = 0.004 59 nm-' (20" scattering angle with a HeNe laser) and a single particle radius of 7.5 nm, corresponding to the gold aggregates.With these parameters kR, =. 1, and we are at the limit of the small-kR, approximation for the scattering. As can be seen from the figure, with these calculations of Q, we are unable to obtain a value as large as you have measured for polystyrene using any value of T. Furthermore, reducing the value of M , results in lower values of Q for all T. Indeed, the only way we can achieve values of Q as high as you report is to extend our calculations to higher kR,. Thus I believe great caution should be exercised in attempting to use Q to determine T. If it is possible at all, one probably has also to measure the first cumulant and use the combination of the two in determining the shape of N ( M ) and the value of r.Despite the difficulty in actually determining T from the light scattering, there does appear to be some difference between the behaviour of the polystyrene colloids and the gold colloids. Following the discussion, we performed a similar measurement of Q for RLCA using the gold colloids. At small angles we measure Qz0.4k9.1. As can be seen from the figure, this value is quite consistent with what we would expect for the value of T = 1.5 that we previously measured. It is also consistent with the value of Q =. 0.5 that we calculate using the distribution obtained from the analytic solution of the Smoluchowski equations for the sum kernel. This distribution was found to agree with our measurements of N( M ) using electron microscopy.' However, it is quite different from the value of Q you report for the polystyrene.One must rule out any simple explanation for this difference: for example, your data for the growth of the mean cluster mass with time suggested that it was linear. We generally find exponential growth with time in the RLCA regime, and tend to see linear behaviour when we are in what we have labelled the crossover regime, where we have not quite reached the reaction limit, and diffusion effects still play some role. The resultant cluster mass in this case would not necessarily even have a power-law form, possibly accounting for your observation. However, it is also worth questioning the universality of T = 1.5. The theoretical result you cite was based on an analysis of the kinetic equations describing the aggregation and showed that a simple geometric scaling argument predicted a kernel for the Smoluchowski equations of the form K , = i*-'j, with A = 1 a n d i >> 1.The solution to the Smoluchowski equations using this kernel is highly singular, and is very susceptible to small corrections to the form of the kernel. Thus, while the simplest solution suggests T = 1.5, it is possible, although unlikely, that small corrections to the form of the kernel could change the value of a. This possibility would seem to be supported by the results of various computer simulations. Brown and Ball2 obtained a value of T = 1.5 in a simulation of RLCA, while Meakin et al.3 apparently obtain a value of T closer to 1.8. Thus the value of r may have some subtle dependence on the details of the reaction- limited aggregation process.This should certainly be investigated more fully, and the use of dynamic light scattering is one possible technique to do so, if the ambiguities are carefully sorted through. 1 D. A. Weitz and M. Y. Lin, Phys. Rev. Letr., 1986, 57, 2037 2 W. D. Brown and R. C. Ball, J. Phj~s. A, 1985, 18, L517. 3 P. Meakin, F. Family and F. Y. Viscek, to be published. Dr J. D. F. Ramsay (A.E.R.E., Harwell) said: I should like to ask Dr Weitz about the effect of shear on the scattering from the aggregates of gold particles. The results in fig. 6 of his paper presumably relate to the static situation before and after shearingGeneral Discussion 237 through a capillary, and thus suggest that some irreversible change in structure has occurred.Since the shear rates applied here (lo4 s-’) are high, might one also expect some aggregate break-up to occur? Such aggregates may reassociate in the static situation. It would therefore be revealing to perform experiments during progressively increasing shear and attempt to measure structural anisotropy from the scattering in the planes perpendicular and parallel to the flow. Perhaps Dr Weitz could comment. Dr Weitz replied: It is our belief that neither aggregate break-up during the shear nor aggregate reassociation after the shear occurs in our experiments. We base this on the following experiments. We have found that it is possible to stop the aggregation of the colloids by adding a small amount of surfactant to the solution.This presumably adsorbs to the surface of the colloidal particles, acting as a ‘bumper’ layer, preventing the particles from actually touching and sticking. The addition of the surfactant does not, however, change the static scattering measured from the clusters. Therefore, after the shear was applied, the colloid was collected in a cuvette containing some surfactant in water, preventing any further aggregation or reassociation of the clusters. The reason we do not think that any aggregate break-up occured at the shear rates applied is based on an additional experiment. When we first added the surfactant to the aggregated colloid, then applied the shear, we found that the shear had no measurable effect whatsoever. We interpret this result as being consistent with our picture for the effects of the shear, which is to cause the clusters to deform and bend until the ends touch and stick together, forming loops and strengthening the structure.If the ends are not ‘sticky’, permanent loops cannot be formed and the structure is not significantly altered. However, the fact that no change is observed would seem to rule out any significant aggregate break-up. Nevertheless, your suggestion of performing the scattering while the clusters are actually undergoing the shear, to measure the shear-induced anisotopy, is an excellent suggestion. Furthermore, one might expect aggregate break-up to occur at still larger shear rates, or at these rates with weaker clusters, and this regime would be of great interest to study. Dr J.D. F. Ramsay (A.E.R.E. Harwell) made a general comment on a group of papers. Dr Weitz in his paper and other speakers have referred to the relationship between the fractal dimension of colloidal clusters and the mechanism of the kinetic aggregation process. For example, two distinct processes are mentioned by Dr Weitz which depend on the interparticle potential. In the first (DLCA), where there is effectively no potential barrier, collisions between clusters are rapid and aggregation is solely limited by Brownian diffusion to give D = 1.8. Where the barrier is still several kT (RLCA) a more compact structure, having D ==: 2.1, is formed. Recently we have made small-angle neutron scattering measurements at the high-flux reactor in Grenoble on silica sols (diameter ca.30 nm) at high temperatures in water, where we have been able to follow the process of aggregation kinetically’ at fixed and variable temperatures. These studies follow from earlier work2 which showed that the repulsive interparticle interaction, as derived from the rescaled mean spherical approxi- mation (RMSA), decreased progressively with increasing temperature [see fig. 3(A)]. At a certain threshold temperature [ca. 530 K in fig. 3(A)] aggregation begins to occur and the form of the scattering changes markedly, as shown in fig. 3(B). Although these results will be described in detail later, we note that the power-law exponents for the clusters (corresponding to 0 ) are initially ca. 1.7 and increase progressively to ca. 2.1 as aggregation proceeds. The latter would accord with the RLCA mechanism and may also reflect a reorganisation of small clusters to give a more compact structure.In general therefore, we believe that experimental scattering data may provide invaluable insight into aggregation processes. Nevertheless, some caution may be required if specific mechanisms are ascribed solely on the basis of power-law exponents.238 General Discussion 2 h 0 v + 1 c 8 6 h 0 v 4 4 2 0 I I I I ( f 1 ...... t A-i ! I ( B ) 1 Q/ lo2 k' 2 Fig. 3. Small-angle neutron scattering from a silica sol in pressurised water (ca. 150 bar) at progressively increasing temperatures. Temperatures ( a ) - ( j ) are 293, 353, 443, 473, 503, 533, 548, 573, 598 and 613 K, respectively. Note the difference in scale for I ( Q ) in (A) and (B).1 J. D. F. Ramsay, R. G. Avery and A. F. Wright, to be published. 2 J. Bunce, J. D. F. Ramsay and J. Penfold, J. Chem. Soc., Faraday Trans. 1, 1985, 81, 2845. Dr. D . S . Horne (Hannah Research Institute, Ayr,) addressed Dr Dickinson: What may be surprising to some is the observation that your simulations give similar values for the fractal dimension whatever the form of the interparticle potential. This may be a consequence of the fact that all of your modelling is carried out at relatively high volume fractions. However, would you still not have expected some divergence depen- dent on potential when the volume fraction was decreased rather than the behaviour shown in your fig. 3. Alternatively, you may be confirming by simulation what I have been seeing experimentally in ethanoi-induced aggregation of casein micelles, where the measured value of the fractal dimension did not vary within experimental error over an ethanol concentration range that progresses the aggregation from reaction-limited to diffusion-limited dynamics. Would you like to comment on this aspect of your simula- tions? Further to this, 1 presume that you calculate the fractal dimension of yourGeneral Discussion 239 Fig.4. Aggregate of 512 particles simulated with potential B at 4 = 0.05. A plot in In R , against In rn gives D = 2.1 +0.05. aggregate at long reaction times; do you see any time dependence in the values obtained and, if so, is there any effect of reaction potential? Dr E. Dickinson (University of Leeds) replied: There are several possible reasons why the two interparticle potentials give essentially the same values of the effective fractal dimension. First, it could be that, although potential A has a significant energy barrier, which reduces the reaction rate by about an order of magnitude (see table 2 of my paper), the aggregation is still not slow enough to be in the reaction-limited regime.There is certainly evidence in the literature' that if the sticking probability is not very low the fractal dimension will be essentially identical to that for diffusion-limited aggregation. Secondly, our simulations are done at relatively high volume fraction, and it could be that clusters are inhibited in their movement by other clusters from exploring all possible configurations prior to coagulation.This would be particularly the case at the later stages of aggregation and at the higher volume fractions, and its effect could be to reduce the sensitivity of the fractal dimension to the interparticle potential. Thirdly, it could just be that the statistical uncertainty is masking a real difference in fractal dimension for the two potentials. Significantly, perhaps, the value for potential A is never lower than that for potential B at constant volume fraction. On the question of the time dependence of the fractal dimension, the statistics are not good enough to come to any reliable conclusions. Your experimental turbidity measurements on ethanol-induced aggregation do seem to suggest that casein micellar aggregates have a similar long-range structure to the final 512-particle aggregates produced in our simulations (see fig.4). If the fractal dimension is indeed relatively insensitive to the interparticle potential, it could be that the short- range structure of casein micellar aggregates, as measured by g ( r ) , is dependent on ethanol concentration, even though the long-range structure is not. 1 P. Meakin, Phys. Rev. A, 1983, 27, 604; 1495.240 General Discussion Prof. H. N. W. Lekkerkerker (Utrecht, The Netherlands) said: Felderhof and Jones have shown that the dynamics of colloid crystais is drastically affected by electric dipole interactions caused by the lag of the Debye cloud about the moving macroions. In your Brownian-dynamics simulations you use a screened Coulomb repulsion, i.e. you do not take into account the effect discussed by Felderhof and Jones.It is justified to omit this effect in Brownian-dynamics simulations of charged colloidal liquids. Dr Dickinson replied: The electrical double-layer around a charged colloidal sphere is distorted from spherical symmetry when the particle moves under the influence of an externally applied force or simply by Brownian motion. The distortion induces a microscopic electric dipole, which acts to oppose the particle motion, thereby reducing its diffusion coefficient below that of the equivalent uncharged sphere. Theory and experiment show’ that the effect of the Debye cloud on the single-particle diffusion coefficient is of the order of several percent for K a = 1 but negligible for Ka b 3. Since this latter condition is satisfied for the charged particles considered in our Brownian- dynamics simulations, we do not expect any effect of the drag of the Debye cloud on the calculated rates of coagulation and aggregate dissociation.The situation considered by Felderhof and Jones’ is that in which charged spheres form colloidal crystals with a lattice spacing much larger than the particle diameter. Crystalline ordering is caused by long-ranged repulsive interactions between the particles at low ionic strength ( K a << 1). Under these conditions the interactions between the Debye clouds on different particles are no longer pairwise additive, and so we have avoided this situation in our Brownian-dynamics simulations, since we did not wish to consider multibody electrostatic interactions as well as (multibody) hydrodynamic interactions.It should also be noted that the property considered by Felderhof and Jones, the damping rate of the transverse phonons, is a collective property of the colloidal crystal. This is a very different type of dynamic behaviour from that which occurs during Brownian aggregation of a colloidal liquid. 1 G. A. Schumacher and T. G. M. van de Ven, Faraday Discuss. Chem. Soc., 1987, 83, 75. 2 B. U. Felderhof and R. B. Jones, Furaday Discuss. Chem. Soc,, 1987, 83, 69. Dr M. La1 (Unilever Research, Port Sunlight Laboratory) addressed two questions to Dr Dickinson: ( 1 ) Data presented in table 1 of your paper suggest that a linear trimer is more stable than an equilateral trimer. Since the hydrodynamic and static forces acting on the terminal particles in the linear trimer would be smaller in magnitude than in the equilateral trimer, one would conclude on physical grounds that the latter configuration is more stable (both thermodynamically and kinetically) than the former.Would the authors comment on this apparent difficulty of reconciling their data with physical intuition? (2) In the simulation study of the Brownian coagulation process, have the authors attempted to obtain the rate constant of the process from the data on the variation of particie number density with time? If so, was the rate process found to be in accord with the Smoluchowski theory both in respect of the order and the value of the rate constant? Dr Dickinson replied: The relative dissociation times for linear and equilateral trimers depend essentially on two factors: the strength of the interparticle attractive forces and the arbitrary choice of criterion to denote ‘dissociation’.The data in table 1 refer to simulations for particles having only a shallow secondary minimum (< k T ) , and we know from earlier work’ that the dissociation kinetics of such a triplet system is determined much more by the Brownian motion than the interparticle forces. The criterion of ‘dissociation’ chosen here is that each particle should get to a distance fromGeneral Discussion 241 both the others of greater than three times the particle radius. What our results suggest is that, under conditions where dissociation is mainly diffusion-controlled, a trimer becomes three monomers more quickly from an equilateral configuration than from a linear one.The physical intuition implied in your question relates, I believe, to a trimer dissociat- ing into dimer plus monomer. Certainly, for that case, I would expect faster dissociation from a linear configuration, especially so as the secondary-minimum well depth becomes of the order of a few kT. Notice, however, that the very process of linear monomer-dimer dissociation has, on average, the tendency to reinforce the integrity of the dimer, since the middle particle of the triplet finds it impossible to move away from both end particles at once without making a triangle configuration, from which it might have been better to start in the first place! For dissociation into three monomers, whilst retaining a linear configuration, the outer pair of particles must move through a relative distance of 0.65 pm; this compares with a change in pair separation ofjust 0.325 p m for dissociation in which the equilateral triangle configuration is maintained. As far as the coagulation rate constant is concerned, we have not calculated the rate constant as such, although approximate values can be estimated from the half-lives in table 2 of our paper.One reason why it is not appropriate here to report explicitly the rate constant is that we have chosen deliberately in this simulation to neglect interparticle hydrodynamic interactions, which we know have a significant effect on the coagulation rate, but not, we believe, on the structure of the growing aggregates.Secondly, what one would like to have is the initial rate, but unfortunately this is not well defined for a concentrated system because it depends on the structure of the dispersion at time t = 0, i.e. the assumed form of interaction between the particles immediately prior to the onset of coagulation. This problem does not arise in dilute dispersions, whose structure prior to coagulation approximates to that of a perfect gas. The problem is, incidentally, even more severe experimentally than in the simulation, since for a non- dilute system it is just not possible in practice to change uniformly the interparticle potential (e.g. by adding electrolyte) before some of the particles have already started to coagulate owing to that same change in experimental conditions.With regard to the Smoluchowski theory, one would in any case not expect it to be applicable at the volume fractions studied here, and indeed table 3 of our paper shows this to be so for the distribution of aggregate sizes after one half-life. We note, however, that a comparison of simulated aggregate size distributions for potentials A and B, obtained from single runs at very low volume fraction ( # = 0.005), has shown’ them both to differ from the Smoluchowski distribution (table 3) by no more than the expected statistical uncertainty. 1 J. Bacon, E. Dickinson, R. Parker, N. Anastasiou and M. Lal, J. Chem. SOC., Faraday Trans. 2, 1983, 2 G. C. Ansell, Ph. D. Thesis (University of Leeds, 1986). 79, 91. Prof. J. A. McCammon (Houston, Texas) remarked: In your simulations of sediment formation, hydrodynamic interactions were calculated only between the moving particle and its nearest neighbour in the existing deposit.Do you think that including hydrody- namic interactions with more of the particles in the deposit would change the sediment structure? One might imagine, for example, that the current approximation could increase the density of the sediment slightly by slowing the approach of the moving particle toward peaks and increasing the probability that the particle will move toward valleys in the sediment. Dr Dickinson responded as follows: Our compromise choice of hydrodynamic approximation was based on an attempt to reproduce the essential physics as faithfully as possible but with the minimum of computational complexity.In the absence of any hydrodynamic interaction, the trajectories of particles at high field strengths have a242 General Discussion roughly ballistic character, which means that the sedimenting particle simply sticks to the first particle that it meets on its approximately linear path. This unphysical effect is eliminated by the inclusion of a hydrodynamic interaction between the moving particle and its nearest neighbour in the deposit. Although we have not done any computational tests, my feeling is that including more hydrodynamic interactions would slow down the deposition rate, but would not have a significant effect on the final sediment density or the short-range structure as described by g ( r ) . This is because I think that it is the short-range hydrodynamic repulsion between pairs which has the crucial influence on the structure of systems formed by irreversible directed deposition, although I would certainly not wish to rille out completely the effect which you intuitively envisage.I should perhaps emphasize here that, in our approximation, the closest particle in the deposit may in fact change its identity several times before the moving particle finally becomes immobile. So, to some extent, the hydrodynamic interaction between the moving particle and the ‘valley’ is included, but, admittedly, only after the hydrodynamic interaction with the neighbour- ing ‘peak’ has already been accounted for. Prof. M. Fixman (Colorado State University) turned to Dr Lal: Prof. Oppenheim pointed out that the presence of a wall poses a special problem in the formulation of a Langevin equation, if the presence of the wall is recognized by the imposition of boundary conditions.However, the wall can be represented as an external potential; this requires only a familiar modification of the Langevin equation. If the potential function rises continuously rather than discontinuously, which is certainly acceptable on physical grounds no special difficulty will arise in the computer algorithms. Can it be said whether a special treatment of boundaries is generally preferable to the use of a smooth external force in the simulation of these systems? Dr La1 replied: In reality, a particle diffusing in a hydrodynamic medium ‘feels’ the presence of a wall, located in close proximity, uia the hydrodynamic interaction.This interaction includes a repulsive force which increases continuously with the reduction in particle-wall separation, in the limit of zero separation the magnitude of the repulsion approaching infinity. Such a force arises from the variation of the particle diffusion tensor with distance from the wall, with the result that the divergence of the tensor is not zero. Thus a simulation model for Brownian dynamics which takes a proper account of the particle-wall hydrodynamic interaction will accord with Prof. Fixman’s suggestion that the wall be represented as an external potential. Prof. J. A. McCammon (Houston, Texas) then asked: How might one relax the approximation in which particle trajectories are initiated at the centre of the interfacial region in the Brownian-dynamics calculations of reaction times? Depending on what one wants to simulate, might it be more accurate to start from a distribution of positions or (in the case of diffusion through an interfacial region bounded on one side by a stirred bulk phase) from a plane near the entry surface to the interfacial region? Dr La1 replied: The Brownian-dynamics method discussed in the paper does include the provision for selecting the initial position of the particle according to a prescribed distribution of which the delta function distribution (fixed initial position) is a special case.In the present (rather preliminary) work, our objective was simply to establish the applicability of the Brownian-dynamics method to surface-confined particles and to show the validity ofthe resuits obtained by comparing these with the first-passage-time results.In the context of this objective, it is immaterial whether the initial position is located in the middle of the interfacial region or close to the virtual boundary, as long as it is the same in both Brownian-dynamics and first-passage-time calculations.General Discussion 243 However, as the calculations were directed at absorption at both the boundaries (in model I), location of the initial position in the middle of the region offered the advantage that comparable numbers of absorption events could be realized for the two boundaries, leading to estimates of (7,) and (7,) (as well as of p r and p , ) with not-too-different degrees of statistical accuracy. The authors agree with Prof.McCammon that in the calculation of the diffusive transport of the particle from the bulk sol to the surface, the initial particle position should be taken as that located very close to the virtual boundary which serves as the dividing plane between the bulk sol and the interfacial region. Prof. McCammon then turned to Prof. Smith: In a real fluid, the effective range of hydrodynamic interaction between a pair of particles must be limited by the physical properties of the medium (e.g. retardation). Can these limitations be exploited to eliminate the divergence problems that you described? Dr R. B. Jones (Queen Mary College, London) added: I wish to comment on Prof. McCammon’s suggestion that retardation may help to control the divergent long-range hydrodynamic interactions when using periodic boundary conditions in Brownian- dynamics simulations.This suggestion is supported by the treatment of the hydrody- namic interactions in a periodic colloidal crystal. In order to do the lattice dynamics it is necessary to write the momentum equations at finite frequency for both fluid and lattice. This circumstance means that in the linearized Navier-Stokes equations for the fluid there is a fluid inertial term of the form -iwpu(r, w ) which gives retardation in the response of the fluid velocity field v. The presence of this inertial term regularizes the k=O term in the lattice sum for the fluid velocity and by use of an appropriate Ewald summation technique, as explained in the paper of Hurd et al.’ and derived first by Hasimoto,’ the rest of the lattice sum can be evaluated to give finite and well behaved results.1 A. J. Hurd, N. A. Clark, R. C. Mockler and W. J. O’Sullivan, J. Fluid Mech., 1985, 153, 401. 2 H. Hasimoto, J. Fluid Mech., 1959, 5, 317. Prof. E. R. Smith (LaTrobe University, Australia) answered both comments: If we include retardation, we change the model. What we need are simulation data for a specific model system, to check against theory and experiment. ,Retardation may, in the end, give us a better behaved mobility tensor. At the moment it does not seem useful since we do not have results for a non-retarded system with which to compare it. Dr E. Dickinson (University ofLeeds) then asked: There are two separate but related problems in incorporating hydrodynamic interactions into Brownian-dynamics simula- tions of concentrated dispersions: the boundary condition and the multi-body interac- tions.You quote the observation that the Rotne-Prager tensor is not positive definite in a system with periodic geometry imposed through the use of the minimum-image distance convention. It is known, however, that for a cluster of spheres in close proximity, the Rotne- Prager tensor with pairwise additivity does not give an adequate representation of the diffusive motion, even in an unbounded system with no imposed periodicity. Is it not possible, therefore, that there might exist a satisfactory form of mobility tensor (including the effect of multi-body interactions, either directly or implicitly) which would remain positive definite under minimum-image distance conditions? In other words, does the boundary condition problem reflect a deficiency in the Rotne-Prager tensor, or is it generally applicable to all forms of the mobility tensor, however well refined? Prof. Smith replied: As far as I can understand the problem, explicit use of periodic boundary conditions means that the quasi-static linearized Navier-Stokes equation does244 General Discussion not have a solution without an array of unphysical constraints on the system.With the construction used in this work, a periodic mobility tensor results. Thus, if we want to use a periodic mobility tensor to avoid edge effects on the simulation sample, we must use the construction of this paper.These considerations will also apply to clusters. Prof. M. Fixman (Colorado State University) commented: The point is made that the Rotne-Prager interaction does not give a positive definite friction matrix if periodic boundary conditions are used. However, the problem can be formulated as the solution to a variational principle (this was in fact the method used by Rotne and Prager), and any trial functions used for the velocity fields will give a positive definite friction matrix. See the references cited.”* Will this a‘ ernative approach, which omits the explicit spherical boundary condition that you consider, be mathematically inconsistent or physically improper? Another alternative might usefully be considered, namely the treatment of all cells but one in the periodic lattice as a continuum, for the purpose of calculating hydrody- namic motion.Continuum environments been considered for dielectric and diffusion simulations [For the latter see ref. (3).] This procedure, which is certainly not free from mathematical difficulties, could be combined with the use of periodic short-ranged interactions. 1 J. Rotne and S. Prager, J. Chem. Phys., 1969, 50, 4831. 2 M. Fixman, J. Chem. Phys., 1982, 76, 6124. 3 M. Fixman, J. Chem. Phys., 1984, 81, 3666. Prof. Smith answered: When we put a boundary condition on a simulation sample, we change the model. In particular we change the model of the hydrodynamics. We could indeed use the minimum image convention if we used a mobility tensor derived from a variational principle analysed within the minimum image boundary condition.I suspect that this would do more damage to the physical hydrodynamics than using a periodic boundary condition. Prof. B. U. Felderhof (R. W.T.H., Aachen, Federal Republic of Germany) said: I appreciate that the minimum image convention in conjunction with the Rotne-Prager tensor for the hydrodynamic interactions leads to problems in a system with periodic boundary conditions. Your resolution of the problem consists of enclosing the system in a large sphere, but leads to a modification of the mobility tensor which is so unpleasant that you did not wish to write it down. You also stressed that it is important to treat these problems both from a microscopic and a macroscopic point of view. Would it not be possible to consider N particles enclosed in a sphere, but immersed in an infinite fluid, so that their hydrodynamic interactions may be approximated by the Rotne-Prager tensor? Since one knows the macroscopic properties of such a sphere it should be possible to extract the desired information on transport coefficients from a computer simulation.Prof. Smith: How large a computer simulation would you like? Make a bid. Prof. Felderhof: Say 862 particles. Prof. Smith: In that case there are severe layering effects. and then continued: It is certainly possible to construct a mobility tensor of the type suggested by Prof. Felderhof. The analysis is rather similar to that with a reaction field interaction, and leads to a short-ranged minimum image mobility tensor. Unfortunately this tensor is computationally as inconvenient as the mobility-tensor for the periodic boundary condition.This is because the reaction-field mobility tensor explicitly involves three-body terms, even when calculated to the first ‘Rotne-Prager’ level.General Discussion 245 Dr M. La1 (Unilever Research, Port Sunlight Laboratory) directed attention to Prof. Fixman’s paper and said: Would Prof. Fixman agree that Monte Carlo calculations on chains, carried out by Zimm among other workers, have settled the question of the adequacy of the pre-averaging approach which underlies the Kirkwood-Riseman theory, for taking into account the intersegmental hydrodynamic interactions in a polymer molecule? It is indeed surprising that such a complex many-body hydrodynamic problem as embodied in the dynamical behaviour of a macromolecule has a simple answer which has stood the test of time for four decades.Prof. Fixman replied: The question of adequacy is settled for many purposes, but not all. The Kirkwood pre-averaging approximation was later shown to provide an upper bound to the translational diffusion constant, and the Zimm approximation of rigid body motion was shown to provide a lower bound. Monte Carlo simulations of the two approximations generally agree rather closely, typically within 10-20%, when both are applied to the same model. The feature that seems most striking to me is not that Kirkwood’s work has survived as a useful approximation, but that experiments and simulations have advanced to the point of uncovering its errors. The Zimm approxima- tion is probably better for most models, but requires a great deal more computation. An indirect approach to the Zimm lower bound, based on the inclusion of ‘internal friction’, considerably eases the Computational burden. Dr N. J. B. Green (University of Notre Dame) said: I have two questions for Prof. McCammon. I note from your paper that there are some difficulties incorporating the reactive boundary condition in your simulations. One problem that arises from the non-zero timestep sampling of diffusive trajectories close to an absorbing boundary is that the only information simulated for the trajectory is at the discrete sampled time points. It is quite possible for trajectories to enter the reactive zone and leave it again during the course of such a timestep. Thus a trajectory that has the potential to react is not identified because all the sampled configurations lie outside the reactive region. Failure to identify such trajectories and allow them to react with the required probability can lead to a significant underestimate in both the rate and amount of reaction occurring. The difficulty can be circumvented, at the expense of computational resources, by using smaller timesteps near the boundary. However, we have recently shown how construction of a ‘bridging process’, which is the diffusion process interpolating between sampled points on the trajectory, leads to an estimate of the probability that the reactive zone is entered during a timestep, conditional on the configurations at the start and at the end of the timestep.’ This estimate converges to its true value under the same conditions as the time-discretisation you have used converges to the stochastic differential equation for the trajectory. How have you dealt with this problem in your simulations? 1 P. Clifford and N. J. €3. Green, Mol. Phps., 1986, 57, 123. Prof. McCammon replied: The use of simple Gaussian probability distributions to describe random displacements of particles in Brownian-dynamics simulations requires special consideration near either absorbing or reflecting boundaries. The small timestep method can be used even with boundaries of complicated geometry; it is computationally demanding, but this can be ameliorated by the use of multiple timestep techniques. Near boundaries with simple geometries, we have made use of modified probability distributions that explicitly allow for absorption or reflection.’ This is a more efficient approach, as you suggest, and also allows for the simultaneous analysis of different boundary models if appropriate book-keeping is done.’ The bridging methods that you describe may also prove quite useful in describing the effects of reactive boundaries in Brownian dynamics studies of diff usion-controlled reactions. 1 S. H. Northrup, M. S. Curvin, S. A. Allison and J. A. McCammon, J. Chem. Phys., 1986, 84, 2196.246 General Discussion Dr Green then continued: The stochastic differential equation from which you simulate is equivalent to a diffusion equation for the diffusion-reaction process. In the diffusion approximation no account is taken of particle velocities, indeed they do not exist (they diverge). Is this formalism appropriate for the systems you are modelling, and can you justify the assumptions implicit in such a formalism? Prof. McCammon again replied: Very close to a reactive boundary, the Smoluchowski description does not provide an accurate description of the particle dynamics. For the systems we have modelled, however, the momentum relaxation times are short enough that boundary layer effects can be neglected.

ISSN:0301-7249    

DOI:10.1039/DC9878300223

出版商:RSC    

年代:1987

数据来源: RSC

摘要:

Furuday Discuss. Chem. SOC., 1987, 83, 247-258 Maximum Entropy Analysis of Dynamic Parameters via the Laplace Transform Alastair K. Livesey* MRC Laboratory of Molecular Biology, Hills Road, Cambridge and Department of Applied Mathematics and Theoretical Physics, Silver Street, Cambridge Mireille De1aye-F and Pedro Licinio Laboratoire de Physique de Solides, Universite' de Paris Sud, B6timent 51 0, 91 405 Orsa-y, France Jean-Claude Brochon LURE, CNRS, MEN, CEA, Universite'de Paris, Sud, B6timent 209d, 91405 Orsay, France Research in the field of Brownian motion often raises the problem of inverting the Laplace transform. Thus, the measured multi-exponential relaxations are the Laplace transform of the spectrum of decay times. However, the inversion of the Laplace transform is very ill-conditioned.The maximum entropy (maxent) method is shown to be a very powerful technique for such data analysis, capable of handling single peaks, multiple peaks and broad distributions in the spectra. After a theoretical demonstration of why entropy chooses the optimal reconstruction, the power of the technique is illustrated with two types of experimental examples related to Brownian dynamics. In the first example, photon correlation spectroscopy (PCS) is applied to concentrated colloidal dispersions. The results, analysed by maxent, present a detailed and coherent picture of the dynamics of a concentrated colloidal solution, leading to a test of the many-body hydrodynamic theories. In the second example, the maxent analysis of the time-resolved fluorescent decay of L-tryptophan is presented.Maxent provides some evidence that one of the features is broadened, suggestive of translative motion in the structure and/or heterogeneity of the fluorophore's environment. 1. Introduction The Laplace transform often occurs at the kernel of equations describing Brownian dynamics. For the two examples considered here [photon correlation spectroscopy (PCS) and pulsed fluorimetry], the data are described, respectively, by the equations developed in ref. (1) and (2): D( tkjca'" = C g ( ~ ~ ) exp (- tk/Tj) + B ( 1 ) J D(tdCalc= E a(Tj) exp (-tk/7j)+B' J where DCalc( t k ) are the calculated data at a delay time t k , g( T) and a ( T ) are the spectra of delay times we wish to reconstruct from the two experiments, B and B' are constant backgrounds, E is the shape of the excitation flash in the second example and @ denotes a convollcltion product.t Deceased. 247248 Maximum Entropy Analysis Laplace Transform unphysical I solutions Lapiace 4 transform . I I. I data D ( t ) Fig. 1. Diagram showing the set of all g ( 7 ) or a(.) spectra. (- - -) Boundary of the set of g(T) or a(.) which agree with data D( t ) within experimental accuracy. (- - -) Boundary of the set of physically allowable g( T) or a(.). The intersection of these two sets is the feasible set (hatched). Inside the feasible set, maxent will choose a 'preferred' solution. Such systems of equations are very ill-conditioned (because of the Laplace transform) or ill-posed because the experimenter is unable to measure sufficient data points (par- ticularly at very short or very long measurement times with respect to the characteristic time constant of the phenomenon under study).Many varied data-analysis techniques have been suggested to alleviate these problems [ eigenvalue decomposition, least- squares, linear programming, determination of a digitised g ( . ) over a finite grid, see ref. ( 1 ) and (2)]. However, they all have a restricted range of validity or usefulness and cannot satisfactorily treat both uni-model and multi-model spectra nor handle equally well data with very diiferent noise levels.' Maximum entropy (maxent) is a very powerful data-analysis tool which has been demonstrated to give excellent results over a wide range of experimental conditions for other, much better conditioned, applications which include radio-astronomical imaging, medical imaging (CAT, MRT n.m.r.spectroscopy, optical deconvolution, electron energy loss spectroscopy, Raman spectroscopy etc. (for a review see Gull and Skilling'). In this paper we report that maxent still provides excellent results when the problem is very ill-conditioned and is capable of handling a wide range of complicated data from both photon-correlation spectroscopy (PCS also known as quasi-elastic light scattering) and pulsed fluorimetry. In the following section, we present the theory of maxent and discuss its general implementation for these two cases. In section 3 , PCS data from concentrated colloidal dispersions are analysed by maxent and the tests which these data provide of many-body hydrodynamic theories are discussed.Section 4 presents a similar discussion and results for pulse fluorimetry. 2. Maximum Entropy We can view the ill-conditioning of the Laplace transform, which leads to a multiplicity of allowable solutions, in the following way.' Consider the set of all possible shapes of the decay spectra a ( T) or g ( 7) displayed as a rectangle in fig. 1. We can calculate 'mock' data sets from each curve in turn [using eqn ( 1 ) or (2)] and test whether it agrees with the noisy data set. All those spectra which agree with the data within the experi- mental precision are bounded by a dot-dashed line. Many of these spectra, however,A. K. Livesey et al. Table 1. 249 ( a ) general solution for 0 s x s 1 ( b ) maximum negative correlation 1 3 0 ( c ) uncorrelated result 0 2 3 ( d ) maximum positive correlation contain unphysical features such as negative regions.We reject the unphysical spectra by drawing a dashed line. The remaining subset of spectra (shown shaded) we call the feasible set.3 Every member of this set agrees with the data and is physically allowable. A full specification of the feasible set is the complete answer to the problem. However, it is, for any non-trivial data analysis problem, much too large (strictly infinitely large) to display, comprehend or use. We are forced to choose one (or at most a few) spectra for display and subsequent analysis. Since we are forced to choose one member of the feasible set, we do so directly by maximising some function F [ a ( T)] or F [ g ( T)] of the spectrum.The function F should be chosen to introduce the fewest artefacts into out chosen distribution. In order to understand better the principles of MEM we consider a specific and simple example originally given by Gull and S k i l l i ~ ~ g . ~ This was deliberately chosen for clarity of exposition rather than practical realism. Information: one third of kangaroos are left-handed one third of kangaroos are blue-eyed. Question: on the basis of this iuformation alone, estimate the proportion of kangaroos which are both blue-eyed and left-handed? The feasible set in this example is sufficiently small to display it as a 2 x 2 contingency table (table 1). We also display three specific choices from this set; namely those with maximal positive and negative correlations and the uncorrelated result.We believe the only rational choice of a single result from such a problem is to choose the uncorrelated result and we should like to choose that function, F ( p , ) , which gives this result. It is of course quite possible that there are some (genetic?) correlations between left-handedness and eye-colour, but we have no knowledge whether this will give rise to a positive or negative correlation in our results. We certainly do not wish our mathematical optimisation function to introduce such correlations arbitrarily when the data contained no information about these correlations. Using such an example, it has been proved by Livesey and Skilling’ that only the Shannon-Jaynes entropy,6 S, will give the uncorrelated solution.This function is defined as: S = c pi-m,-pilog- Pj I mi (3) where pj = a ( ~ ~ ) or g ( T j ) as appropriate and mi is the model which encodes our prior knowledge about the system before the experiment (see below).250 Maximum Entrop-v Analysis Laplace Transform Table 2. regularisation function P1 correlations Shannon entropy - C pi log pi $=O.lll uncorrelated least squares -C p : = 0.0833 negative Burg entropy C log pi positive intermediate form pf” 0.12176 positive (d17 - 1)/24 = 0.13013 A more complete and very general proof has been given by Shore and Johnson.’ Any other ‘regularising function’ (Tichonov and Arsenin’) will impose extra structure in our spectrum without justification. Thus table 2 shows the results for some ‘popularly’ used regularisation functions.As we have stated, all the functions give either positive or negative correlations, with the exception of Shannon-Jaynes entropy. It is clear that the fundamental property of this example is our ability, given the data, to form independent sub-classes and is not dependent on the actual physics of the measuring device. The results from our example are thus very general and can be applied to any ill-conditioned or ill-posed problem when we are forced to choose a solution. Practical measurements are also noisy. It would be wrong to fit the data exactly, because the noise would then be interpreted as if it were true signal. We chose to bound the feasible set by a chi-squared statistic: where M is the number of (independent) observations times t and (T, is the estimated standard deviation of the observation at delay time t which, assuming Poisson statistics to be valid, was taken as the square-root of the actual number of counts detected.At values of C > M the calculated and observed data are not in agreement and for C < M we start to fit the noise in D ( t ) O b ’ . The model m ( ~ ) gives our best a priori guess for the structure of pJ. Jaynes [ref. ( 6 ) , p. 1141 showed that if one has no knowledge of the values of the decay constants expected then one must choose: m, = 1/73. ( 5 ) A brief resumi of his proof has been given in the appendices of two of our previous which further show that, by transforming variables and working with equally spaced points in log T space, the model rn, becomes constant.We note the MEM solution also displays many secondary, generally useful properties. ( 1 ) Because of the log term in the entropy, the spectra automatically have the same sign as the model. (Usually this is positive.) (2) The spectra are smooth. (3) It only shows features (in particular the resolution of close peaks) if demanded by the data. (4) Provided the digitisation is sufficiently fine to show all the relevant detail, the shape of the spectrum is independent of the number points used to display it. ( 5 ) It is robust to noise. ( 6 ) For data which are linear functions of the spectrum (this includes the Laplace transform and convolutions), the maxent solution is unique. Finally, we note that, in general, the maxent solution does not have a statistical interpretation.It is no more likely than any other member of the feasible set, but remains the best choice we can make for any given data set. For further discussion see Livesey and Skilling.’ Computationally, we used the Cambridge maximum entropy suite of subroutines (MEMSYS) to find the reconstructions of a ( ~ ) or g ( T ) which have the maximal entropyA. K . Livesey et al. 25 1 S [eqn (3)] subject to the constraint that its calculated data gave a chi-squared value C [eqn (4)] less than or equal to one. A full description of the original Cambridge algorithm has been given by Skilling and Bryan.' 3. Photon Correlation Spectroscopy of Concentrated Colloids In this first example, maxent was used to analyse photon correlation spectroscopy (PCS) data from concentrated colloidal dispersions.Although Brownian motion is correctly understood in dilute colloids (provided the direct interactions between colloidal particles are known), the situation is considerably more complicated in concentrated dispersions where the hydrodynamic interactions involve many-body effects. Systematic theoretical treatments of such many-body effects have recently been proposed,1031' and need to be tested experimentally. PCS is a very powerful tool to probe Brownian dynamics. However, to simplify the interpretation of the results, PCS experiments should use colloidal particles whose shape stays invariant with the concentration c of the dispersion and whose direct interactions are known. In this example we used a colloid of eye lens specific proteins: the a-crystallins.These proteins are quasi-spherical oligomers with a mean hydrodynamic radius Rh = 95& a charge 2 =: 50e at pH 6.8 and an intrinsic polydispersity (ca. 15% in size)."," Stable dispersions of these proteins can be prepared, at various ionic strengths (IS), up to high concentrations (c > 0.300 g ~ m - . ~ , i.e. volume fractions qb 3 0.5). The structures of these dispersions have been carefully investigated by small-angle X-ray scattering and are correctly described by hard-sphere interactions at high ionic strengths and hard-sphere plus a screened Coulombic repulsion at low ionic strength." We used a classical PCS set-up, but we accumulated correlation functions % ( t > = (i( t')i( t'+ t ) ) of the scattered intensity i( t ) over a time range spanning up to 7 decades by matching successive pieces of % ( t ) measured at different sample times.Further details of the experiment and the data analysis can be found in ref. ( 1 ) and (14). Typical correlation curves, from a-recrystallin dispersions at IS = 17 meq, are shown as a semi-log plot in fig. 2 and 3. Fig. 2 corresponds to four different concentrations measured at the same scattering angle 8 = 90" [the wavevector q = (47m sin 8 / 2 ) / h ; n is the refractive index of the solution and A the wavelength of the light]. Fig. 3 corresponds to scattering at four different angles from the same concentration, c = 0.152 g ~ m - ~ . The correlation function at t = t k was assumed to be adequately described by: Initially we turned this into a linear problem by taking the (positive) square root of D( t ) - B*: ~ ( t ) = t / / ~ ( l ) - ~ * ) s i g n [D(t)-B*] (7) where D( tk)Obs is the observed value of the correlation function at t k and B" is our best estimate of the background.This helped MEMSYS to converge in what is a very well conditioned problem, but an iterative procedure is needed to determine the background. Later, a more powerful version of MEMSYS (Skilling, personal communication) allowed us to solve the non-hear equation (5) directly. It also suppressed the iteration needed to find B and even gave a four-fold increase in speed (typically 450 data points can be analysed to give a 100 point spectrum in 6 min on an HP1080 or VAX 11-780 computer). Fig. 2 and 3 show, on the same semi-log plot as the o ( f k ) measurements, the distribution of delay times g( T,) obtained by maxent.At low concentrations [fig. 2( a ) ] , g(7) shows a single peak corresponding to translational diffusion of the a-crystallin252 Maximum Entropy Analysis Laplace Transform 10 103 lo5 1 0' r / p s and t/p5 Fig. 2. Intensity correlation functions % ( t ) as measured from PCS (. - .), and relaxation times distribution g( 7) as recovered by maxent (-) from a-crystallin dispersions at 8 = 90" scattering angle and various concentrations c. ( a ) c = 0.008 g cmP3, ( b ) c = 0.152 g cm-3, ( c ) c = 0.226 g cmP3, ( d ) c = 0.280 g ~ m - ~ . particle. However, as c increases, three peaks are clearly resolved [fig. 2 ( b ) - ( d ) ] . The wavevector analysis of fig. 3 indicates that the characteristic mean time of each of the three peaks shifts with q.More quantitatively, T, varies as qP2, showing that three modes of relaxation are diffusive. We also note that: ( a ) at a given c the relative amplitudes of the first two peaks are independent of q (fig. 3); at a given q the amplitude of the first peak decreases with respect to the second when c increases, while the splitting between these two peaks increases (fig. 2). ( b ) At a given c the relative amplitude of the third peak increases at low q (fig. 3); at a given q the third peak is seen to shift both toward longer times and to increase in relative amplitude when c increases (fig. 2 ) . ( c ) Further structure is detected at low q, long times and high c, whose q dependence is more difficult to determine.The interpretation of the first peak is straighforward. In the q range explored in these experiments [ q << qm, the maximum of the static structure factor S ( q)"], the fastestA. K. Livesey et al. 253 10 1 o3 lo5 10' 'TIPS and t / p s Fig. 3. Intensity correlation function %(t) as measured from PCS (. . .) and relaxation times distribution g ( T ) as recovered by maxent (-) from a-crystallin dispersions at c = 0.152 g cmP3 and various scattering angle 8. The timescale is logarithmic. ( a ) t9 = 130", ( b ) 0 = 90", ( c ) 6 = 60", ( d ) 8=40". decay is associated with mutal diffusion:I5 where D, is the mutual diffusion coefficient, F ( t, q ) is the dynamic structure factor and S ( q -+ 0) = F ( t --+ 0, q - 0) is the static structure factor as q + 0.The shift of the first peak towards shorter times when c (and 4 ) increases reflects the competition between direct and hydrodynamic interactions. Since the pair correlation function g ( r) has been measured,l3 these data should provide a test of the hydrodynamic theories. The second peak can be attributed to self-diffusion, which is detectable by PCS in polydisperse and the evolution of the self-diffusion time with c approximately follows that of the viscosity of the usp pension.'^ The amplitude of the first peak decreases254 Maximum Entropy Analysis Laplace Transform 5 4 h 23 2 1 Fig. 4. Maxent reconstructions (broken lines) and two exponential fits (symbols) to form analogous simulations of D( t ) with different noise distributions.(-1 Original (Gaussian distribution); (- - - -) and a, fit to first simulation; (- . -) and +, fit to second simulation; (- - - -) A, fit to third simulation; (- - -) and 0, fit to fourth simulation. relative to that of the second one as c increases and results from the decrease of S ( 0 ) with increasing concentration. The third diffusive peak, whose amplitude grows at low q, probably corresponds to the diffusion of large structures (of a size comparable with the wavelength of light) whose proportion and size would increase with c. These structures are reversible upon dilution and might simply be an equilibrium distributions of ‘clusters’ resulting from long-range van der Waals attractions. The peaks at longer times are not yet fully characterized. A more complete interpretation of these data will be presented in a poster at the meeting.Maxent is a very reliable method for analysing any type of PCS data, whatever the noise level, as shown by various computer simulations presented in fig. 4. We started from a Gaussian distribution (in log 7) and calculated mock data sets using eqn (8) such that each channel ‘measured’ between (2.0 and 1.0) x lo7 counts. We then added four different sets of noise drawn from a pseudo-Gaussian distribution to these curves to generate four similar sets of ‘observations’. In attempting to analyse these noisy curves by a least-squares two-exponential fit, we found that the values for the two components depend strongly on the detailed noise distribution (see fig. 4). In each case the calculated data sets were statistically indistinguishable from the ‘observations’.Nevertheless, although the chosen position and amplitudes of the peaks shows consider- able variation, in each case the mean of two peaks corresponds to the mean of the true distribution. (The data could not be fitted satisfactorily by a single exponential.) In contrast, the four maxent distributions are very similar and, as we would expect from our selection procedure, all of them show slightly less structure than the originalA. K . Livesey et al. 255 and are thus a conservative reconstruction of the original. In particular, the mean position of the peak is very accurately determined (if the tails of the reconstruction are included) and the widths are fairly reproducible. A year and a half's experience with this data-analysis technique has led us to the following general conclusions for the application of maxent to PCS.( 1 ) With ca. 80 measurement points spanning two decasdes and ca. lo7 counts per channel, two delta functions can only be reliably separated when their true values are more than a factor of 2.4 apart. (2) Maxent can reliably handle mono- and poly-disperse systems that contain single or multiple peaks at good signal-to-noise ratio, giving easily interpretable results with virtually no user intervention. Maxent can also handle sparse, noisy data sets from a Laplace transform,'* but the interpretation requires a more careful and complete under- standing of the selection theory outlined in section 2. Once these concepts are mastered, however, the interpretation becomes routine and a surprisingly large amount of informa- tion can be gleaned from apparently poor data sets.(3) Maxent will not provide conclusions unwarranted by the data. In particular, it will not set the proportion of long decay times to zero unless forced to by the data. In order to prevent maxent fitting tails to the g(T) distribution to allow for the possibility of long delay times, it is essential to measure the decay well out into the noise of the background, thus ensuring no further components are decaying. Since the signal is inherently weak here, many more measurement are needed than a naive analysis might suggest. (4) Although the Laplace transform gives us very limited resolution in decay time, it can produce reconstructed spectra with very large dynamic ranges.We have found that we can reliably detect minority peaks whose areas are only 1% of the larger peaks in the spectrum.' Computationally, however, the program must be able to cope with dynamic ranges of up to 10", which can be quite a challenge to maxent algorithms. 4. Pulse Fluorimetry of L-Tryptophan Unpolarised pulse fluorescence is being used increasingly in biophysics to study the fluctuations and/ or heterogeneity in the three-dimensional structure of macromolecules which often arises from Brownian motion of flexible components. Polarised pulse fluorescence gives a direct physical measurement of the Brownian rotational constants in solution of the whole molecule, deformations and flexions of the structure and local motions of side-chains around their chemical bonds.Both these aspects can be related to biological activity. In this paper, however, we only present an example of analysing unpolarised pulse-fluorescence decays. Further details and examples can be found in Livesey and Brochon.* The fluorescent decay of L-tryptophan was measured after excitation by synchrotron radiation from the electron storage ring 'ACO' or Orsay, France. The repetition rate was 13.6 MHz. Fluorescent decays were measured by the time-correlated single-photon- counting techniquelg using the apparatus described by Brochcen.20 L-Tryptophan was recrystallized from ethanol-water 80: 20 v/v and dissolved in 0.01 mol dm-3 acetate buffer held at 10 "C and adjusted to pH 5.4. In pulse 'Huorimetry the Laplace transform of the spectrum of decay times a ( T ) is convolved by the finite temporal shape of the exciting pulse, E ( t ) .i where @ represents a convolution. Although we can use Poisson statistics to give the variance of the fluorescence data, these have to be increased to allow for the noise on the experimentally determined256 Maximum Entropy Analysis Laplace Transform lo4 1 o3 1 o2 h 'Lr v % 1 0 -4 4 0 - 4 -6 010 h k- v ti 0.0 5 h ++++ t 4 t + + t + + u t I 1 h , 4 2 10 20 30 tlns log T/ns Fig. 5. (a) Fluorescence decay curves of L-trytophan at two emission wavelengths (Ah = 4 nm); 1, 390nrn; 2, 32Qnm; +, experimental curves; (-1 calculated curves; 3 is the experimental flash profile. The excitation wavelength was 280 nm (AA = 4 nm) and temperature 10 "C.( b ) Plot of the weighted residuals corresponding to two curves ( 1 ) and (2). ( c ) Maxent reconstructed spectra of a ( ~ ) for 390 nm emission is plotted in arbitrary units in order to normalise peak maxima for longer decay times.A. K. Liuesey et al. 257 excitation peak shape E ( tj, The total variance a: is given by standard error propagation analysis to be where E,,, is the integral of E ( t ) , 2 is the calculated decay law and E is $ for k = 1 and k = t, otherwise E = 1. The measured data points at two wavelengths, together with the (wavelength-indepen- dent ) excitation pulse, are shown in fig. 5 ( a ) with a logarithmic ordinate. On this plot, the decay at 390 nm appears to be a straight line, showing that the decay could be fitted by a single exponential.The decay at 320 nm is, however, more complicated, involving at least two significant time constants. The corresponding maxent reconstructions of the distributions of decay constants are presented in fig. 5 ( d ) and the solid lines on fig. 5(a) show the fit to the data by these spectra. At 390 nm there is only evidence for one strong peak centred at 4.46 ns with a half-width of ca. 0.1 ns. At 320 nm the spectrum of decay constants has two peaks centred at 0.82 and 4.35 ns. The relative integrated areas under these peaks are 0.623 and 0.377, respectively. These results are in good agreement with several previously published r e s ~ l t s . ~ ~ - ~ ~ Note further, though, that the peak of the shorter-time decay is considerably wider than the long decay.Simulations using delta peaks at 0.82 and 4.35 ns, and similar signal-to-noise ratios, give equally narrow reconstructions to the two peaks, suggesting there is very little evidence for the shorter-decay component being a single delta function. Nevertheless, fitting these data with two delta functions gives a chi-squared value of 1.48, which is only just rejected by the experiment, so that the intrinsic width of this shorter-decay peak could well be narrow. Better data will be needed to fully resolve these models. Although we have only presented our results from analysing the total (unpolarised) fluorescence, we have developed a maxent program capable of analysing polarised data to give a two-dimensional graph of the number of fluorphors Y(T, 0) having both an intrinsic decay T and characteristic rotational constant 8.The program has been tested on computer-simulated data and our results from analysing experimental data are presented as a poster at the meeting. A. K. L. thanks EMBO and MRC for support. P. L. gratefully acknowledge CAPES and UFMG, Brazil for joint support. The PCS was supported by MRT (Interface Physique Biologie) and INSERM (Contrat externe finance par la CNAMTS). CNRS provided the financial support for the pulse fluorescence project. References 1 A. K. Livesey, P. Licinio and M. Delaye, J. Chem. Phys., 1986, 84, 113. 2 A. K. Livesey and J-C. Brochon, Biophys. J., 1987, in press. 3 A. F. Gull and J. Skilling, ZEE Proc., 1984, 131F, 646. 4 S. F. Gull and J. Skilling, in Indirect Imaging, ed.J. A. Roberts (Cambridge University Press, Cambridge, 5 A. K. Livesey and J. Skilling, Actu Crystullogr., Secf. A, 1985. 41, 45. 6 E. T. Jaynes, Collected Works. Papers on Probability, Sfatistics and Sratistical Physics, ed. R. D. 7 J. E. Shore and R. W. Johsnon, IEEE Trans., 1980, 26, 26. 8 A. Tikhonov and V. Arsenin, in Solutions of Ill-posed Problems (Wiley, London, 1977). 9 J. Skilling and R. K. Bryan, Mon. Not. R. Astron. Soc., 1984, 211, 111. 1985), pp. 267-279. Rosenkrantz (D. Riedel, Dordrecht, 1983). 10 C . W. J. Beenakker and P. Mazur, Phys. Lett., 1982, 91A, 290. 11 P. Mazur and W. Van Saarloos, Physrca, 1982, 115A, 21. 12 A. Tardieu and M. Delaye, J. Mol. Biol., 1986, 129, 711. 13 A. Tardieu and M. Delaye, J. Phys. (Paris), 1987, 48, 1207. 14 P. Licinio, M. Delaye and A. K. Livesey, J. Phys. (Paris), 1987, 48, 1217.258 Maximum Entropy Analysis Laplace Transform 15 P. N. Pusey and R. J. A. Tough, J. Phys. A, 1982, 15, 1291. 16 M. B. Weissman, J. Chem. Phys., 1978, 68, 5069. 17 P. N. Pusey and R. J. A. Tough, in Dynamic Light Scattering and Velocimetry Applications of Photon 18 R. Papoular and A. K. Livesey, Maximum Entropy Analysis oj-Neurron-Spin Echo Data, in preparation. 19 W. R. Ware, in Creation and Detection of the Excited State, ed. A. A. Larnola (Marcel Dekker, New 20 J. C. Brochon, in Protein Dynamics and Energy Transducrion, ed. S . I. Ishiwata (Taniguchi Foundation, 21 A. G. Szabo and D. M. Rayner, J. Am. Chem. Soc., 1980, 102, 554. 22 R. J. Robbins, G. R. Fleming, G. S. Beddard, G. W. Robinson, P. J. Thistiewaite and G . Wodfe, J. 23 J. P. Privat, Ph. Wahl and J. C. Brochon, Riocheme, 1985, 67, 949. Correlation Spectroscopy, ed R. Pecora (P!enum Press, New York, 1984). York, 1971), vol. 1, pp. 213-302. Japan, 1980), pp. 163-189. Am. Chem. Soc., 1980, 102, 6271. Receiued 12 rh December, 1986

ISSN:0301-7249    

DOI:10.1039/DC9878300247

出版商:RSC    

年代:1987

数据来源: RSC