Simulations and Stochastic Mean Field Theory - American Chemical

5564

J. Phys. Chem. B 2007, 111, 5564-5572

Linking Phase Behavior and Reversible Colloidal Aggregation at Low Concentrations: Simulations and Stochastic Mean Field Theory Antonio M. Puertas Grupo de Fı´sica de Fluidos Complejos, Departamento de Fı´sica Aplicada, UniVersidad de Almerıá, 04120, Andalucıá, Spain

Gerardo Odriozola* Programa de Ingenierıá Molecular, Instituto Mexicano del Petro´ leo, La´ zaro Ca´ rdenas 152, 07730 Me´ xico, D. F., Mexico ReceiVed: December 18, 2006; In Final Form: March 2, 2007

We have studied the link between the kinetics of clustering and the phase behavior of dilute colloids with short range attractions of moderate strength. This was done by means of computer simulations and a theoretical kinetic model originally developed to deal with reversible colloidal aggregation. Three different regions of the phase diagram were accessed. For weak attractions, a gas phase of small clusters in equilibrium forms in the system. For intermediate attractions, the system undergoes liquid-gas separation, which is signatured by the formation of a few large droplike aggregates, a gas phase of small clusters, and an overall kinetics where a few seeds succeed in explosively growing at long times, after a lag time. Finally, for very strong attractions, fractal unbreakable clusters form and grow following DLCA-like (diffusion limited cluster aggregation) kinetics; liquid-gas separation is prevented by the strength of the bonds, which do not allow restructuration. Good qualitative and quantitative agreement is found between the dynamic simulations and the kinetic model in all the three regions.

I. Introduction Over the past few years, phase behavior and aggregation kinetics of colloidal systems have been largely studied. Whereas the former establishes a link with atomic and molecular systems, the latter is exclusive of colloids and macromolecules. Both of them play important roles in many natural and human controlled processes, such as sediment deposition, water clarification, and emulsion stabilization/destabilization, to cite a few. However, an overall framework for both remains elusive. This is mainly because phase behavior lies under the title “statistical mechanics of equilibrium”, whereas aggregation kinetics clearly belongs to the “out of equilibrium” processes. Hence, theoretical tools employed to deal with the former are in general not appropriate to deal with the other, and vice versa. This distinction hinders access to the different aspects of the whole phenomenon. A commonly studied colloidal system consists of a dispersion of colloidal particles in a solvent and nonadsorbing polymers, where the latter induce the so-called “depletion attraction” among colloidal particles.1,2 Here, the range of the attraction is controlled by the polymer size, whereas its strength is set by the polymer concentration. Thus, this system constitutes a very friendly colloidal model which is easily implemented in simulations. For large colloidal concentrations, the phase diagram of such a system has been studied in detail.3-5 For this purpose, different experimental setups have been used, and the range of polymer-size to colloid-size ratio has been swept out. For short range attractions, it has been shown that the liquid-gas separation becomes metastable inside the crystallization region.6,7 So, upon increasing the polymer fraction, a * Corresponding author.

fluid phase is obtained for weak attractions and crystallites are obtained for stronger attractions. If crystallization is impeded (by using a polydisperse system, for instance) at even stronger attractions, a liquid-gas separation is indeed found. For very strong attractions, the system gels, forming a percolating network that usually comprises all particles of the system and prevents phase separation.8,9 Gelation has been interpreted in different ways, and its relation with the phase diagram is as yet unclear.10-17 In addition to this phenomenology, a phase of equilibrium clusters has been recently found for attractions weaker than those inducing gelation at low densities. It has been interpreted in terms of a short range attraction and a long-range repulsion, attributed to a residual surface charge of the particles.18-21 Colloidal aggregation, on other hand, despite its long history remains largely unrationalized, although some important features have been recognized such as the existence of aggregation regimes22 and the dynamic scaling of its kinetics.23 Aggregation of sticky free-diffusing particles corresponds to diffusion limited cluster aggregation (DLCA), which produces fractal aggregates with dimension df ≈ 1.8 and a linear increase of the mean cluster mass.22,24 Aggregation of free-diffusing particles that bind on contact with a given probability, P , 1, produces the so-called reaction limited cluster aggregation (RLCA) regime, characterized by df ≈ 2.1 and a faster than linear increase of the mean cluster mass.25-28 From the point of view of the kinetics, the major contribution was given by Smoluchowski,29,30 who described the time evolution of the aggregation process by means of a population balance equation. In this framework, only binary collisions are considered and a reaction kernel is defined containing the relevant physical information such as cluster

10.1021/jp068698b CCC: $37.00 © 2007 American Chemical Society Published on Web 04/27/2007

Linking Phase Behavior and Reversible Aggregation

J. Phys. Chem. B, Vol. 111, No. 20, 2007 5565

diffusion, structure, and reactivity. Further development of the population equation has been achieved by including cluster fragmentation31-33 and cluster sedimentation.34 The relation among aggregation, gelation, and phase transitions is however far from evident.35-37 In this paper, we establish a link from the kinetics point of view between phase separation and aggregation. We describe the clustering phenomena and phase behavior of a polydisperse colloidal system with short range attractions at extremely low concentrations, where three-body collisions are rare. At low attraction strength, we find equilibrium clusters, with bigger size the stronger the attraction. Liquid-gas coexistence is found at intermediate strengths, where a quasispherical liquid drop coexists with the dilute vapor phase. Increasing the attraction, irregular ramified clusters are formed, resembling DLCA for the strongest attractions studied. The kinetics of the system can be successfully described in all regions with a Smoluchowskitype stochastic model, including cluster breakup, originally developed to describe reversible aggregation.33 Single, double, and triple bonds are considered in the model, whose population is monitored. Previous studies dealing with aggregation driven by attractions of finite strength have mainly focused on the structure of the clusters or some average kinetic properties.37-41 The results shown in those works agree with our general observations, but we focus on the kinetics of clustering. The paper is organized as follows; the details of the simulations are given in section II, and the stochastic model is explained in section III. For a proper comparison between simulations and theory, different parameters have to be fixed; this is done in section IV. Results are presented and discussed in section V. Section VI summarizes the conclusions. II. Simulations Simulations were performed considering a polydisperse system composed of N ) 8000 soft-core particles, using Brownian dynamics or, more precisely, strongly damped Newtonian dynamics. Each particle experiences a Gaussian distributed white noise force, b gj, with zero mean and δ-correlated in time, a damping force proportional to the velocity, γb r˙ , and the deterministic forces due to interactions. Hence, the equation of motion for particle j is

mb r¨ j ) - γb r˙ j + b g j(t) +

∑i BFij

( ) r d12

(

)[

(2)

where d12 ) a1 + a2, with a1 and a2 being the radii of the interacting particles. This potential effectively behaves as a hardcore42 and avoids the problems of dealing with noncontinuous potentials.43 The short range attraction mimics the effective component induced by the entropy driven distribution of the nonadsorbing polymers around the colloidal particles. In the form derived by Asakura and Oosawa,1 the interaction range is given by the polymer size, ξ, and the strength is given by its concentration, φp:

] } (3)

III. Stochastic Mean Field Theory Population balance equations have been used for describing the aggregation kinetics of colloids since the beginning of the last century.29,30 They are a set of infinite first-order coupled differential equations, which, in general, are numerically solved. An alternative way to generate a solution to these equations is by means of a stochastic approach.27,32,45 From this point of view, the master equation rules the process, which for reversible aggregation in its simplest version reads32,33,46

dP(N B, t)

)

dt

1 2V

kij[(Ni + 1)(Nj + 1 + δij)P(N B+ ∑ ij , t) i,j

B, t)] + Ni(Nj - δij)P(N

1

∞

∑

2Vn)2

n-1

[(Nn + 1)

fi(n-i)P(N Bi(n-i) , t) ∑ i)1 n-1

Nn

-36

]

3r r3 (ξh + 1)2 + + 4ξ 16ξ3 r 2 3ξ a1 + a2 2 (ξh + 1) 4r ξ 2ξ

for d12 + 2ξ/5 e r e d12 + 2ξ and VAO ) 0 for larger distances. Here, ξh ) (a1 + a2)/2ξ. For computational efficiency, at short distances, this potential is replaced by a steep parabola with minimum at r ) d12; this ensures that the total interaction potential retains a quadratic minimum very close to r ) d12. The resulting total interaction potential, Vtot ) Vsc + VAO, is analytical for all r and allows for straightforward integration of the equations of motion using a stochastic algorithm to treat the random forces.44 The particle radii are distributed according to a flat distribution [0.9a,1.1a], where a is the mean radius, to avoid crystallization induced by the attractions. Crystallization would hide liquid-gas separation, occurring at higher attractions at short interaction range, and would affect the local structure inside the cluster.41 The attraction range is set to 2ξ ) 0.2a, the attraction is measured in units of φp, kBT is set to kBT ) 4/3, lengths in units of the mean radius, a, and the particle mass is m ) 1. The friction coefficient, γ, is fixed at γ ) 50, which gives a Brownian time of τB ) a2γ/6kBT ) 6.25 time units. The equations of motion are solved with a time step δt ) 5 × 10-4. All simulations are performed at constant colloid density, volume fraction φc ) 10-3, and for increasing attraction strength.

(1)

where the friction coefficient γ and b gj are related by means of the fluctuation-dissipation theorem: 〈g bi(t)g bj(t′)〉 ) 6kBTγδ(t t′)δij. The interaction potential has a soft-core, given by

Vsc(r) ) kBT

{[

VAO(r) ) -kBTφp (ξh + 1)3 -

fi(n-i)P(N B, t)] ∑ i)1

(4)

where P(N B, t) is the probability for finding the system in the state N B ) (N1, N2, ..., Nn) at time t, Ni is the number of i-size clusters, kij is the aggregation kernel, and fij is the fragmentation kernel. Displaced states N+ ij and Nij are given by

N B+ ij ) N Bij )

{ {

(..., Ni + 1, ..., Nj + 1, ..., Ni+j - 1, ...) for i * j (..., Ni + 2, ..., N2i - 1, ...) for i ) j (..., Ni - 1, ..., Nj - 1, ..., Ni+j + 1, ...) for i * j (..., Ni - 2, ..., N2i + 1, ...) for i ) j

Hence, N+ ij are those states that, by following a single aggregation event, may produce the state N B. Analogously, Nij are those that may produce N B by breaking a cluster in two parts. In eq 4, the aggregation kernel, kij, establishes the mean rate at which i-mers and j-mers stick to form (i + j)-mers. On the

5566 J. Phys. Chem. B, Vol. 111, No. 20, 2007

Puertas and Odriozola

other hand, fij is the mean rate at which (i + j)-mers spontaneously break producing i-mers and j-mers. Both kernels are understood as orientational and configurational averages of all particular aggregation rates. In equilibrium, dP(N B, t)/dt equals zero and so does the righthand term of eq 4. That is, the average number of events that exit state N B equals the average number of events from all other states that produce state N B. Equilibrium may be forced by setting B, t) ) fij(Ni+j + 1)P(N BB and kijNi(Nj - δij)P(N ij , t) for all states N N Bij ; namely, by imposing the average number of events from state N B toward state N Bij equal to the average number of reverse events. This condition, the so-called detailed balance, has been very successful for developing Monte Carlo algorithms, where the goal is to obtain equilibrium configurations no matter whether the movements are naturally possible (indeed, they usually are not).47 Additionally, it has also been used to deal with reversible aggregation.48,49 However, it should be pointed out that detailed balance in aggregation is a very strong restriction which has no physical basis, and so, it may obstruct the understanding of the real behavior of the processes. Hence, our choice is to handle eq 4 as it is, which is capable of describing not only the nonrestricted equilibrium but also the out of equilibrium region. Irreversible aggregation, having dP(N B, t)/dt * 0 ∀ N B and t, can be studied by setting fij ) 0. In particular, diffusion limited cluster aggregation (DLCA) is produced by the following aggregation kernel

1 kBrown ) kSmol (i1/df + j1/df)(i-1/df + j-1/df) ij 4 11

(5)

where kBrown states for the Brownian kernel, kSmol ) 8kBT/3η is ij 11 the Smoluchowsky dimer formation rate constant, df is the cluster fractal dimension, kBT is the thermal energy, and η is the solvent viscosity.29,30 In previous works, we employed the following expression for modeling the cluster breakup.33,46

fij ) eij(1 + δij)(1 - PCij)/τ

(6)

where τ is the characteristic decay time of an exponentially shaped bond breakup probability density function and PCij is the probability for two recently produced fragments to collide and reaggregate. That is, since recently produced fragments are space correlated and since mean field theories do not explicitly account for these correlations, we must introduce a correction for the fragmentation kernel. The same correction is employed for modeling an irreversible reaction limited cluster aggregation (RLCA) regime, since collisions that do not lead to aggregation also lead to space correlated clusters. This correction, 1 - PCij, is given by (N11(ij)b)-1, where N11 is the average number of consecutive collisions per encounter of nonaggregating monomers and b is a constant exponent.28,33 This expression has been contrasted with experimental data50,51 and with detailed theoretical models.52 On the other hand, the function eij represents the average number of bonds contained in the n-sized clusters that, on breakup, lead to i- and j-sized fragments. This was obtained in ref 33 for loopless aggregates yielding the following fitted expression

eij(1 + δij) ) p1(ip2 + jp2)(ip3 + jp3)(ij)p4

(7)

where p1 ) 0.4391, p2 ) 1.006, p3 ) -1.007, and p4 ) -0.1363. For the smallest aggregates, we set the exact known values 2e11 ) e12 ) 2.

The use of expressions 5 and 6 for solving the master equation models a system where all collisions lead to bond formation and where the produced bonds live, on average, τ. In the simulations performed here, we also expect all collisions to be effective, since the imposed interparticle potential has a single minimum and no energetic barrier. On the other hand, we do not expect the clusters to break into large-large child fragments, due to the high coordination of particles inside the aggregates. Thus, for solving the master equation, we impose fij ) 0 for i > 1 and j > 1, i.e., there are no breakup events where the smallest child fragment is larger than a monomer. This represents an extreme case of nonshattering fragmentation.32,46 In addition, we also expect the particles on the surface of the aggregate to coordinate with several particles. Hence, those particles that may separate from the aggregate to form a monomer should have different characteristic decay times, which should depend on the number of particles of the aggregate that directly interact with them. This leads us to propose the following expression for the fragmentation kernel

fij )

{

[∑ ]

eij(1+δij)(1 - PCij) ET|i+j 0

k

Ek|i+j for i ) 1 or j ) 1 (8) τk otherwise

where ET|i+j ) ∑k Ek|i+j is the total number of bonds at the surface of the clusters of size i + j, which equals the total number of particles on these clusters that may follow a breakup process. In addition, subindex k refers to the strength of these bonds, i.e., k ) 1 means that the bond links two particles like a dimer does, k ) 2 occurs when a particle is linked to, for instance, a dimer forming an equilateral-like trimer, and k ) 3 refers to a bond linking an equilateral-like trimer with a particle located in such a way to obtain a tetrahedral-like tetramer. In other words, k ) 1 represents a single bond, k ) 2, a bond with double strength, and k ) 3, a bond with triple strength. Hence, τk, the characteristic decay time of a k-type bond, equals τk1, according to the Arrhenius expression and assuming an additive pair potential among particles. So, highly coordinated particles at the surface of the aggregates are much more difficult to separate from the cluster than others with lower coordination numbers. We considered k ) 3 to be the largest coordination number of a particle at the surface of a cluster. It should be noted that states are now defined not only by N B, but by N B and B E, with B E ) (..., E1|i, E2|i, E3|i, ...). Hence, the time derivative of the probability for finding the system in a given state at time t, P(N B, B E, t), is a complex function of N B, B E, Pij ) (P1|ij, P2|ij, P3|ij), being Pk|ij the probability for kij, fij, and B a surface bond formation process between i- and j-sized clusters to produce a k-type bond. It would be tedious and not instructive to report this lengthy equation here, and so, we prefer to comment on the algorithm employed for solving it. For the sake of clearness, we present in Figure 1 a flow diagram that schematizes it. Note that expression 8 is practically equal to eq 13 of previous work,53 and so it is the algorithm to produce N B and B E. Hence, the flow diagram of Figure 1 is similar to that shown in Figure 9 of the same reference. Additionally, both algorithms are slightly modified versions of that given elsewhere,46 where it is described in detail. So, we will describe some generalities of this modified version and focus on the differences. In the first step, initialization, clusters and bond populations are set to their initial values and the reaction channels45,46 are calculated accordingly. Following, the main loop is started where the time increment, dt, the cluster sizes, and the type of reaction to follow is decided according to a stochastic algorithmssee


J. Phys. Chem. B, Vol. 111, No. 20, 2007 5567 by one unit. Note that the number of surface bonds cannot increase more than one unit, i.e., e1i + e1j e e1(i+j) + 1. To choose the type of bond to be increased, another random number, ζ1, is generated and compared to B Pij, as shown in Figure 1. From this comparison, the population of primary, secondary, or tertiary bonds is increased. Aggregation processes involving monomers always lead to an increase of the total bond population, ET. When two large clusters aggregate, however, ET decreases. That is, some particles at the surfaces of both aggregating clusters turn internal particles of the new cluster. In this case, we must reduce the population of surface bonds, B E|i+j. This is done by comparing a random number ζ2 with B Ji+j ) (J1|i+j, J2|i+j, J3|i+j), being Jk|i+j ) Ek|i+j/∑lEl|i+j. In this way, k is selected and then Ek|i+j is decreased one unit. This procedure is repeated ETo|i+j - ETn|i+j times, as shown in Figure 1. Analogously, N B is updated for a breakup process, but by decreasing N1+j and increasing N1 and Nj a single unit. The bond transference is also accomplished symmetrically to that described above. That is, ETn|j+1 is given by [e1(j+1)Nj+1], and ETo|j+1 - ETn|j+1 bonds are transferred toward B Ej. Immediately after, ETn|j is approached by [e1(j)Nj] and compared to ETo|j. Note that ET cannot increase, since breakup always leads to a monomer formation. In the case that ETn|j ) ETo|j, no update of B E|j is done. Otherwise, another random number ζ3 is generated and compared with B Bj ) (B1|j, B2|j, B3|j), being Bk|j ) Ek|j/τk|j/ (∑lEl|j/τl). This is done for deciding which type of bond has to be broken. Finally, in the case of breakup and reaggregation, two extra random numbers homogeneously distributed in [0,1] are generated, to be compared with B Bj and B P1j, respectively. From this comparison, B E|j is updated. B Pij is estimated as follows53

P1|ij )

Figure 1. Flow diagram that schematizes the solving procedure of the master equation.

ref 46 for details. The only difference is that, when a fragmentation channel is selected, a random number homogeneously distributed in [0,1] is compared to 1 - PCij to make a decision between breakup and breakup followed by reaggregation. The first choice leads to changes on N B and B E, whereas the second just may change B E. Hence, breakup followed by reaggregation can be interpreted as rearrangements on the clusters surfaces. In the case of aggregation, N B is updated, i.e., Ni and Nj are decreased one unit and Ni+j is increased one unit. Additionally, a transference of bonds from sizes i and j toward i + j is accomplished. For that purpose, ETo|i - ETn|i bonds are randomly taken from B E|i ) (E1, E2, E3)|i and added to B E|i+j. The same procedure is followed to transfer bonds from B E|j toward B E|i+j. Here, ETn|i is the total number of bonds on the surface of i-sized clusters after the aggregation event, and it is approached by [e1iNi], whereas ETo|i ) ∑kEk|i. It should be noted that this procedure keeps ET|i and ET|j close to e1iNi and e1jNj, respectively. However, ET|i+j needs readjustments, since the total number of bonds, ET ) ∑i∑kEk|i, was kept constant. In order to adjust the i + j size bond population ETn|i+j is approached by e1(i+j)Ni+j and compared to ETo|i+j ) ∑kEk|i+j. In the case that ETn|i+j equals ETo|i+j, the bond populations at the surfaces of the (i + j)-sized clusters, B E|i+j, are not changed. For ETn|i+j ) ETo|i+j + 1, we increase the population of a given type of bond

A1 A1 | | AT iAT j

P2|ij )

A2 A1 A1 A2 |i |j + |j |i AT AT AT AT

P3|ij )

A3 A1 A1 A3 |i |j + |j |i AT AT AT AT

where AT is the cluster area exposed to aggregation and A1, A2, and A3 are the areas of the clusters that may lead to single, double, and triple bonds, respectively. The areas are approached by53

AT|n ≈ n4πa2 - (n - 1)2πa2 ) (n + 1)2πa2 A2|n ≈ (n - 1)2πa2R A3|n ≈ (n - 1)2πa2β(E2 + 2E3)/ET A1|n ≈ AT|n - A2|n - A3|n

(9)

2π/3 being a the particle radius, R ) ∫2π/3-∆q sin θ dθ, and β ) ∫∆q¢ 0 sin θ dθ. We employed R ) 0.1 and β ) 0.002. These estimations are expected to produce qualitative results. We need to point out that we impose B E|1 ) 0, E1|2 ) N2, and E2|2 ) E3|2 ) E3|3 ) 0, i. e., monomers do not form bonds, dimers only form k ) 1 type bonds, and trimers cannot form k ) 3 type bonds. Finally, the algorithm produces the time evolution of the cluster size distribution, N B, and bond populations B E, given kij, N11, b, and τ1. Since τ1 is directly measured from simulations, we need to establish kij, N11, and b.


Figure 2. Bond lifetime for different states. The dashed line represents an exponential fitting, τ1 ∼ exp{U0/kBTφp}, yielding U0 ) 11.98kBT.

All curves shown in this work and produced with the stochastic algorithm correspond to the average of 100 realizations, which consider the same system size as the simulations (8000 particles). IV. Parameters of the Model The theoretical model presented above has a number of free parameters, some of them state-independent, that can be fixed from the comparison with the simulations. The aggregation kernel, kij, is mostly state-independent, since it only reports the collision rate between clusters; the fragmentation one, fij, on the other hand, is strongly affected by the internal structure of the clusters, i.e., by the coordination of the particles at the clusters surfaces (the structure of the clusters also affects its diffusion and cross section in such a way that one counterbalances the other, producing small changes in kij; we do not account for this fact). As we will show, this becomes signatured by the change of the parameters N11 and b with the state. However, the direct link between simulations and the theoretical model is established in terms of the bond lifetime, τ1. Thus, we determined τ1 from simulations, and set it in the model. The aggregation kernel is assumed to be state-independent, whereas N11 and b are fitted for every state. The bond lifetime was calculated from simulations of two single particlessa bond is formed when the surface to surface distance is smaller than the attraction range, i.e., 2ξ ) 0.2. Successive realizations of bond breaking are studied and the average lifetime is presented in Figure 2 as a function of the polymer volume fraction, φp.54 The data follow an exponential law with the attraction strength, as expected, with a parameter slightly smaller than the attraction strength ≈ 13.5φpkBT due to the thermal energy and the shape of the well. The fitting allows us to map φp onto τ1 and compare simulations and theoretical predictions for the same state. The state-independent aggregation kernel, kij, can be obtained by matching the stochastic results to the initial stages of the clustering for φp ) 1.0 while setting fij ) 0. That is, for φp ) 1.0, we obtained from simulations τ1 ≈ 28 000, and hence, the effect of fragmentation is negligible for the first 7000 time units. We also avoid the first 500 time units due to the influence of the initial configuration of the simulations. Note that the aggregation kernel is independent of the pair potential. This approximation arises from the range being much shorter than the particle radius; the aggregation rate is thus equal to the collision frequency. The free parameters of eq 5 are the dimer formation rate constant, k11, and the cluster fractal dimension, df. The fractal dimension should be close to 1.8, which is the accepted DLCA


Figure 3. Time evolution of the cluster size distribution. Symbols represent the data from simulations as obtained for φc ) 10-3 and φp ) 1.0: the circle, square, diamond, upward-pointing triangle, and left8 16 pointing triangle correspond to N1, N2, N3 + N4, ∑i)5 Ni, and ∑i)9 Ni, respectively. Lines correspond to the best-fit stochastic approach for fij ) 0 and kij ) 1.44.

value. We did not obtain a good enough agreement with any combination of k11 and df, although we should point out that the solutions of the master equation are not strongly influenced by df. Hence, we try other similar expressions for the kernel. The good match shown in Figure 3 was obtained with

1 kBrown ) k11(i1/df + j1/df)(i-(1/dH) + j-(1/dH)) ij 4

(10)

where k11 ) 1.44, df ) 1.8 for φp ) 1.0, and 1/dH ) 1/df + 1/2. The kinetic constant is very close to the Smoluchowski value without hydrodynamic interactions kS11 ) 8kBT/3η ) 1.34, which shows that the aggregation is indeed diffusion-limited. The difference between dH and df, on the other hand, means that the diffusivity of a cluster with i particles, Di, decreases faster with the cluster size than that predicted by the general accepted expression Di ) kBT/(6πηai1/df). In fact, expression 11 assumes Di ) kBT/(6πηai1/dfi1/2). Note that since df ≈ 2, this expression implies Di ≈ D0/i, i.e., that the diffusion coefficient is inversely proportional to the number of particles in the cluster, or the mass, similar to the Rouse model for polymers.55 Indeed, this model is appropriate for small, noncompact clusters, as those formed at short times and strong attractions. When comparing with experiments, obviously the simulations are not capturing the true nature of cluster diffusion. The simulation algorithm moves the internal particles of an aggregate independently of its relative position inside it, as if they were free in the fluid. Thus, the aggregate center of mass moves only due to a coordinate Brownian kick of most particles (not ALL of the particles need to move cooperatively, since our dynamics is Newtonian below the damping time). For hard spheres, an alternative method for cluster diffusion needs to be used to reproduce cluster motion.37,56 If dH is interpreted as a hydrodynamic dimension of the clusters, the low values needed to describe our results point to the unphysical aspects of their motion. Nevertheless, this approximation is expected to affect the results only at high attraction strength, when large clusters form, and even in those cases, the results can be reproduced by expression 11, which is based on many experimental results.57,58 As discussed above, the fragmentation kernel, set by N11, b, and τ1, depends on every state, and is obtained by fitting the evolution of the concentration of clusters. We show below that the states with weak and moderate attractions can be described with the same values of N11 and b (and varying τ1 according to simulations), whereas different N11 and b are needed for strong attractions.


J. Phys. Chem. B, Vol. 111, No. 20, 2007 5569 TABLE 1: Parameters of the Model Which Produced the Best Agreement for Every State with the Simulationsa φp

τ1

df

N11

b

0.25 0.40 0.42 0.50 1.00 ∞

3.51 21.18 26.91 70.15 27970 ∞

2.4 2.4 2.4 2.4 1.8 1.8

9.0 9.0 9.0 9.0 1.8 6.1

1.25 1.25 1.25 1.25 0.40 0.35

a The parameters for irreversible aggregation are quoted under φp ) ∞ from ref 33.

Figure 4. Time evolution of the cluster size distribution. The circle, square, and diamond symbols correspond to N1, N2, and N3 + N4, respectively, as obtained from simulations for φc ) 10-3 and φp ) 0.25 (upper panel) and φp ) 0.40 (lower panel). Lines correspond to the stochastic approach with τ1 ) 3.51 (φp ) 0.25) and τ1 ) 21.22 (φp ) 0.40).

Figure 5. Time evolution of the cluster size distribution. The circle, square, diamond, upward-pointing triangle, left-pointing triangle, and downward-pointing triangle symbols correspond to N1, N2, N3 + N4, 8 16 ∞ ∑i)5 Ni, ∑i)9 Ni, and ∑i)17 Ni, respectively, as obtained from simulations for φc ) 10-3 and φp ) 0.42. Lines correspond to the stochastic approach with τ1 ) 26.91 (φp ) 0.42).

V. Results

At moderate attraction strength, the system enters the binodal line of liquid-gas phase separation where liquid drops form and grow (remember that crystallization is impeded by polydispersity). The growth can be described nevertheless with our kinetic model with the same parameters as above (df, N11, and b) (Figure 5). Although the evolution of the cluster distribution at short times is very similar to the evolution of systems with lower φp, a large liquid drop, identified as the biggest cluster, grows by addition of small aggregates, as presented in Figure 6 (note that the biggest cluster for φp ) 0.25 does not grow beyond a certain size in the simulations and the state φp ) 0.40 is just in the limit). The cluster distribution is thus very polydisperse, with a few very large droplike clusters and many small ones. Even at short times, larger aggregates are formed in the case of φp ) 0.42 than for lower φp values, whereas the evolution of small clusters is very similar. These few larger aggregates are potential drop-seeds, which may either grow or redisperse. The success of a drop-seed depends on the strength of its bonds, which hinders redispersion, and on the collision frequency with other aggregates. Hence, the overall aggregation behavior shows a lag time needed to produce and age the seeds, such that they do not redisperse. During this lag time, small clusters reach a quasi-steady state, which is abruptly interrupted at longer times, due to the explosive growth of the very few successful drop-seeds. Indeed, droplike aggregates produced by this mechanism turn larger than DLCA-type aggregates after enough time (see Figure 6). We should point out that for τ1 ) 21.18 (φp ) 0.40) the theoretical model produced 82 runs where phase separation took place and 18 where the clusters could not grow (during a time range of 5 × 105). Additionally, most of these 82 runs produce a single successful seed which grew for t > 105. This also points out the stochastic nature of the process, and the important role played by population fluctua-

We study here the isochore φc ) 10-3 where only binary collisions are expected and the kinetic model is applicable, for increasing attraction strength. Since crystallization is avoided by polydispersity, only the liquid-gas transition is expected at moderate attractions, whereas DLCA-like aggregation can appear for strong attractions. Gelation does not occur in this system since we are well below the percolation line for the values of φp and φc studied. For low φp, small equilibrium clusters are observed, although the majority of the particles remain as monomers; see Figure 4. The kinetic model correctly describes the evolution and the equilibrium state of the system with N11 ) 9.0 and b ) 1.25s these values were obtained for φp ) 0.42 where larger clusters form and fragmentation has a larger contribution (see below). The fractal dimension was set to df ) 2.4, based on the compact nature of the clusters produced under low values of φp and on experiments where restructuring occurs59 (note that complete filling of space, represented by df ) 3, would be an extreme value). Small variations of df, as mentioned above, do not significantly influence the kinetics of the process. Table 1 presents the relevant parameters for every state used in the model. The values for N11 and b must be compared with those obtained for simulations of sticky particles where intracluster rearrangements are not allowed: N11 ) 6.1 and b ) 0.3528,33 (reported in Table 1 as φp ) ∞). The large differences in both values can be attributed to cluster aging, which leads to denser clusters and larger particle coordination. Rearrangements are indirectly taken into account in the model by increasing the probability of re-aggregating the fragments when a cluster breaks up, i.e., enlarging function Nij.



Figure 6. Largest aggregate size as a function of time obtained from simulations (symbols) and for φp ) 0.25 (black), 0.40 (red), 0.42 (green), 0.50 (blue), and 1.0 (orange) (φc ) 10-3 in all cases). Lines show the corresponding stochastic predictions.

Figure 8. Time evolution of the cluster size distribution. Symbols represent the data from simulations as obtained for φc ) 10-3 and φp ) 1.0, with the same meaning as in Figure 5. Lines correspond to the stochastic approach for τ1 ) 27970 with N11 ) 9.0 and b ) 1.25 (upper panel) and N11 ) 1.8 and b ) 0.40 (lower panel).

Figure 7. Time evolution of the cluster size distribution. Symbols represent the data from simulations as obtained for φc ) 10-3 and φp ) 0.50, with the same meaning as in Figure 5. Lines correspond to the stochastic approach for τ1 ) 70.15 (φp ) 0.50).

tions, at least for the states practically over the binodal line. It is noteworthy that our mean field model allows for fluctuations in the concentration of species and bond populations due to the stochastic selection of events. It is thus possible to describe nucleation-driven phenomena, such as the phase separation discussed above, given the appropriate values of the parameters. The lines in Figure 6 represent the predictions from the stochastic model for every state. The agreement with the simulations is satisfactory and shows a clear distinction between states with low attraction strength (no aggregation) and stronger attractions. As mentioned above, the prediction for the state φp ) 0.40 is that phase separation is reached in most realizations. The curve shows the average of all runs, independently of their result. The simulation of this state indeed follows the theoretical prediction at short times, but its evolution is extremely slow; from the simulations, it is not clear if the aggregation process will continue or not. Thus, this state is probably metastable inside the binodal, close to the boundary. Figure 7 represents the evolution of the cluster distribution for φp ) 0.50, and the number of particles in the biggest cluster is plotted in Figure 6. In this case, the kinetics of clustering is very different from the state φp ) 0.42; now, many large aggregates are growing and the process resembles DLCA aggregation, i.e., the quench drives the system clearly beyond the spinodal line. However, the evolution can still be described with the same set of parameters as above. Simulations cannot elucidate whether the final state will be liquid-gas coexistence or if this clustering produces long-lived fractal aggregates (as

is typical of colloidal aggregation). The theoretical model, on the other hand, predicts for these states at very long times (t > 107) the formation of one big cluster in coexistence with a gas of monomers, which corresponds to liquid-gas coexistence. Additionally, it predicts two different mechanisms for the growing of the largest drop. These are by coalescence with other smaller ones and by the diffusion of the monomer gas phase in quasi-equilibrium with the smallest drops toward the larger one, similar to the process of Ostwald ripening observed in micelles. Hence, an overall mass transference occurs from the smaller drops toward the larger ones. The time scale of both mechanisms is much larger than that accessed by the simulations. So far, we have produced quantitatively good agreement for the time evolutions of the cluster size distribution with a single set of parameters in the model (see Table 1). However, when comparing the stochastic solutions with these parameters with the simulations for φp ) 1.0, i.e., for τ1 ) 27970, breakup is clearly underestimated. This is seen in Figure 8 (upper panel) where the evolution of monomers is well above that predicted by the stochastic model, which barely presents breakup. Hence, we conclude that large φp simulations behave differently. An inspection of the clusters obtained by the simulations reveals open fractal-like structures, and so, intracluster rearrangements, i.e., aging, should be practically gone. In other words, the proportion of reinforced bonds should be lower. Additionally, fractal open structures may favor large-large breakups, which are neglected by our model. These facts are corroborated by refitting the function Nij. We observed that N11 ) 1.8 and b ) 0.40 produce the good agreement shown in Figure 8 (lower panel)sas expected, Nij decreases for all i and j. Notwithstanding the general good agreement, there is still an overestimation of the population of the largest clusters. This suggests that large-large breakup is indeed taking place in the simulations. On the other hand, the very low value of N11 tends to compensate this fact. That is, although low values of N11 have


Figure 9. Mean number of bonds per particle as a function of cluster size, j, for different states, as labeled, from simulations (upper panel) and the theoretical model (lower panel). The bonds for the state φc ) 10-3 and φp ) 1.0 in the model has been presented for the two sets of parameters used in Figure 8: closed symbols for N11 ) 9.0 and b ) 1.25 and open symbols for N11 ) 1.8 and b ) 0.40.

been suggested,52 we think that N11 ) 1.8 is too low, and it results as a consequence of the lack of large-large breakup. The origin of the change of behavior is thus the difference in the internal structure of the clusters. This is tested in the upper panel of Figure 9, where the mean number of bonds per particle is presented as a function of its corresponding cluster size, for different states in the simulations. As expected, the average number of bonds increases with the cluster size for all states. However, the curves do not reach a plateau, except that for φp ) 1.0. This suggests that the fraction of particles at the surface, which are less coordinated, decreases with the cluster size as a consequence of the droplike shape of the clusters. On the other hand, the curve for φp ) 1.0 does reach a plateau, which is a consequence of the fractal nature of the clusters. Additionally, small clusters are denser for higher φpsin fact, at φp ) 1.0, the number of bonds indicates that clusters are as compact as possible. However, this situation is reversed for bigger clusters. There, the clusters are less dense, since rearranging implies the breaking of many individual bonds. We thus confirm that, for large attraction strengths, the clusters are more strongly hindered in their way to the spherical, most compact structure and as a result are more open, in agreement with previous findings at higher densities.37 Liquid-gas coexistence is hence impeded by the strong bonds, compatible with the interpretation of gelation as the result of arrested liquidgas separation;17 here, due to the low density of the system, we do not observe gels but individual fractal aggregates.38 The theoretical model provides information about the number of bonds per particle, although restricted to those particles at the surface. In the model, the maximum number of bonds per particle at the surface is 3, and no information is given concerning the particles “inside” the cluster. Hence, a direct comparison with the upper panel of Figure 9 is not possible.

J. Phys. Chem. B, Vol. 111, No. 20, 2007 5571 However, we expect the coordination of the surface particles to behave similarly to that of all cluster particles, and so, trends would be comparable. The average number of bonds at the surface is presented in the lower panel of Figure 9 for the same states as the simulations. There, the average number of surface bonds increases with the cluster size for all states but always reaches a plateau. This plateau is obtained at larger cluster sizes for increasing φp, suggesting that aging becomes more difficult. In addition and in agreement with simulations, the bonds strength at the surface increase with φp except at φp ) 1.0. This fact can now be explained in terms of the rates τ2/τ1 and τ3/τ1. That is, for small φp, τ2/τ1 and τ3/τ1 are relatively small, and so, the many formations and breakups that the small clusters experience lead to a pseudo-steady-state populations of bonds having a relatively large number of single ones. As φp increases, τ2/τ1 and τ3/τ1 become larger, since τk ) τk1, producing pseudosteady-state bond populations having a larger fraction of double and triple bonds. The remarkably different behavior found for φp ) 1.0 again points to a much hindered aging process. In the model, this occurs because single bonds are already very stable and conversion from single to double or triple bonds is much rarer than at lower attractions. Thus, in both the simulations and the theoretical model, clusters are denser at low attraction strength, since those with a more open structure dissolve or age to denser ones, more tightly bonded. At φp ) 1.0, on the other hand, aggregates formed with simple bonds are stable and cannot restructure internally, keeping a low number of bonds per particle. The model, however, fails to predict the different behavior of small and large clusters. The pair distribution function or the structure factor could be used to study the internal structure of the clusters. However, due to the small size of the clusters, little information can be gained about their internal structure. The pair distribution function shows a minimum at distances larger than the cluster size for φp ) 1.0, not present for φp ) 0.42, which implies a typical cluster size, and a depletion layer formed during the aggregation process; the large aggregates feed on their surrounding too fast and create a depletion layer.60 At φp ) 0.42, on the contrary, the overall growth of clusters is slow enough to let the diffusion process refill the layer. VI. Conclusions Using computer simulations and a theoretical model originally developed to deal with reversible colloidal aggregation, we have studied the kinetics of clustering phenomena in colloids with short range attractions with moderate strength. For weak attractions, equilibrium clusters form in the system, which enlarge their sizes with increasing the pair potential well. Upon increasing the attraction strength, the system undergoes liquidgas separation, which we observe as the formation of a few large compact aggregates, whereas the rest of the system is comprised of small clusters (crystallization is prevented by using a polydisperse system). For very strong attractions, irreversible clusters form and grow in a DLCA-like manner, causing the decrease of (ideally) all of the monomers and small clusters. Both the simulations and the model predict that the number of seeds decreases with decreasing attraction strength and that their growth is more explosive at long times, after a lag time. The kinetic model used in this work is based on the master equation framework. This equation, as mentioned in the body of this work, properly describes not only the nonrestricted equilibrium state but also the out of equilibrium region. Thus, it is an appropriate tool for linking the phase behavior, frequently accessed by Monte Carlo simulations, and reversible aggrega-

5572 J. Phys. Chem. B, Vol. 111, No. 20, 2007 tion, from the point of view of their kinetics. The model contains two independent kernels: one describing aggregation and the other one describing fragmentation. The latter assumes that a cluster with j particles can break up only to produce a monomer and a cluster with j - 1 particles and considers aging as a consequence of the strengthening of the surface bonds. Despite its apparent complexity, the model incorporates simple events which are present in the simulations, and thus are needed in order to reproduce the “experimental” data. We found that the kinetics of clustering for weak and moderate attractions can be correctly reproduced with a single set of parameters for all states. For strong attractions, however, the parameters of the model must be readjusted to extreme values (in particular Nij, which controls the probability of reforming a recently broken bond). This is a consequence of the restriction to only monomer fragmentation. Nevertheless, the model is able to describe the kinetics of all three regimes: stable clusters, liquid-gas separation, and fractal aggregation, indicating that most of the assumptions of the model hold. Thus, it can be used to predict the final steady state of a given system, much faster than simulations. The presence of equilibrium clusters below the liquid-gas spinodal is reminiscent of the cluster phase recently found experimentally and reproduced by computer simulations with a short range attraction and a long range repulsion.18-20 In our case, the cluster distribution also reaches equilibrium but results from the competition between bonding and thermal breakup.38 Acknowledgment. A.M.P. acknowledges the financial support provided by the Ministerio de Educacioń y Ciencia, under project MAT2004-03581.

Puertas and Odriozola (25) Lin, M. Y.; Lindsay, H. M:; Weitz, D. A.; Ball, R. C.; Klein, R.; Meakin, P. Phys. ReV. A 1990, 41, 2005. (26) Gonza´lez, A. E. Phys. ReV. Lett. 1993, 71, 2248. (27) Thorn, M.; Seesselberg, M. Phys. ReV. Lett. 1994, 72, 3622. (28) Odriozola, G.; Moncho-Jorda´, A.; Schmitt, A.; Callejas-Fernańdez, J.; Martıńez-Garcıá, R.; Hidalgo-A Ä lvarez, R. Europhys. Lett. 2001, 53, 797. (29) Smoluchowski, M. Z. Phys. 1916, 17, 557. (30) Smoluchowski, M. Z. Phys. Chem. 1917, 92, 129. (31) Sontag, H.; Strenge, K. Coagulation kinetics and Structure formation; Plenum Press: London, 1987. (32) Thorn, M.; Broide, M. L.; Seesselberg, M. Phys. ReV. E 1995, 51, 4089. (33) Odriozola, G.; Schmitt, A.; Moncho-Jorda´, A.; Callejas-Fernańdez, J.; Martıńez-Garcıá, R.; Leone, R.; Hidalgo-AÄ lvarez, R. Phys. ReV. E 2002, 65, 031405. (34) Odriozola, G.; Leone, R.; Moncho-Jorda´, A.; Schmitt, A.; R.; Hidalgo-A Ä lvarez, R. Physica A 2004, 335, 35. (35) Poulin, P.; Bibette, J.; Weitz, D. A. Eur. Phys. J. B 1999, 7, 277. (36) Sandkuhler, P.; Lattuada, M.; Wu, H.; Sefcik, J.; Morbidelli, M. AdV. Colloid Interface Sci. 2005, 113, 65. (37) Babu, S.; Rottereau, M.; Nicolai, T.; Gimel, J. C.; Durand, D. Eur. Phys. J. E 2006, 19, 203. (38) Dıéz-Orite, S.; Stoll, S.; Schurtenberger, P. Soft Matter 2005, 1, 364. (39) Rottereau, M.; Gimel, J. C.; Nicolai, T.; Durand, D. Eur. Phys. J. E 2004, 15, 133. (40) Rottereau, M.; Gimel, T.; Nicolai, J. C.; Durand, D. Eur. Phys. J. E 2004, 15, 141. (41) Charboneau, P.; Reichman, D. R. Phys. ReV. E 2007, 75, 011507. (42) Melrose, J. R. EuroPhys. Lett. 1992, 19, 51. (43) Strating, P. Phys. ReV. E 1999, 59, 2175.

References and Notes

(44) Paul, W.; Yoon, D. Y. Phys. ReV. E 1995, 52, 2076.

(1) Asakura, S.; Oosawa, F. J. Chem. Phys. 1954, 22, 1255. (2) Likos, C. Phys. Rep. 2001, 348, 267. (3) Vincent, B.; Edwards, J.; Emmett, S.; Croot, R. Coll, Surf. 1988, 31, 267. (4) Lekkerkerker, H.; Poon, W. C. K.; Pusey, P. N.; Stroobants, A.; Warren, P. Europhys. Lett. 1992, 20, 559. (5) Ilett, S. M.; Orrock, A.; Poon, W. C. K.; Pusey, P. N. Phys. ReV. E 1995, 51, 1344. (6) Dijkstra, M.; van Roij, R.; Evans, R. Phys. ReV. Lett. 1999, 82, 117. (7) Dijkstra, M.; van Roij, R.; Evans, R. Phys. ReV. E 1999, 59, 5744. (8) Duda, Y. J. Colloid Interface Sci. 1999, 213, 498. (9) Trappe, V.; Sandkhuler, P. Curr. Opin. Colloid Interface Sci. 2004, 8, 494. (10) Cassin, G.; Duda, Y.; Holovco, M.; Badiali, J. P.; Pileni, M. P. J. Chem. Phys. 1997, 107, 2683. (11) Vocarin, E.; Duda, Y.; Holovco, M. J. Statist. Phys. 1997, 88, 1333. (12) Bergenholtz, J.; Fuchs, M. Phys. ReV. E 1999, 59, 5706. (13) Segre`, P. N.; Prasad, V.; Schofield, A. B.; Weitz, D. A. Phys. ReV. Lett. 2001, 86, 6042. (14) Pham, K. N.; Puertas, A. M.; Bergenholtz, J.; Egelhaaf, S. U.; Moussaı¨d, A.; Pusey, P. N.; Schofield, A. B.; Cates, M. E.; Fuchs, M.; Poon, W. C. K. Science 2002, 296, 104. (15) Kroy, K.; Cates, M. E.; Poon, W. C. K. Phys. ReV. Lett. 2004, 92, 148302. (16) Sciortino, F. et al. Comput. Phys. Commun. 2005, 169, 166. (17) Manley, S.; Wyss, H.; Miyazaki, K.; Conrad, J.; Trappe, V.; Kaufman, L. Phys. ReV. Lett. 2005, 95, 238302. (18) Stradner, A.; Sedgwick, H.; Cardinaux, F. M.; Poon, W. C. K.; Egelhaaf, S.; Schurtenberger, P. Nature 2004, 432, 492. (19) Campbell, A.; Anderson, V. M.; van Duijneveldt, J.; Bartlett, P. Phys. ReV. Lett. 2005, 94, 208301. (20) Sciortino, F.; Tartaglia, P.; Zaccarelli, E. J. Phys. Chem. B 2005, 109, 21942. (21) Lu, P.; Conrad, J.; Wyss, H.; Schofield, A. B.; Weitz, D. A. Phys. ReV. Lett. 2006, 96, 028306. (22) Lin, M. Y.; Lindsay, H. M.; Weitz, D. A.; Ball, R. C.; Klein, R.; Meakin, P. Nature 1989, 339, 360. (23) van Dongen, P.; Ernst, M. Phys. ReV. Lett. 1985, 54, 1396. (24) Weitz, D. A.; Olivera, M. Phys. ReV. Lett. 1985, 52, 1433.

(45) Gillespie, D. T. J. Phys. Chem. 1977, 81, 2340. (46) Odriozola, G.; Schmitt, A.; Callejas-Fernańdez, J.; Martıńez-Garcıá, R.; Leone, R.; Hidalgo-A Ä lvarez, R. J. Phys. Chem. B 2003, 107, 2180. (47) Frenkel, D.; Smit, B. Understanding Molecular Simulation; Academic Press: New York, 1996. (48) Family, F.; Meakin, P.; Deutch, J. M. Phys. ReV. Lett. 1986, 57, 727. (49) Babu, S.; Gimel, T.; Nicolai, J. C. J. Chem. Phys. 2006, 125, 184512. (50) Lattuada, M.; Sandkuhler, P.; Wu, H.; Sefcik, J.; Morbidelli, M. AdV. Colloid Interface Sci. 2003, 103, 33. (51) Lattuada, M.; Wu, H.; Sandkuhler, P.; Sefcik, J.; Morbidelli, M. Chem. Eng. Sci. 2004, 59, 1783. (52) Lattuada, M.; Wu, H.; Sefcik, J.; Morbidelli, M. J. Phys. Chem. B 2006, 110, 6574. (53) Odriozola, G.; Leone, R.; Schmitt, A.; Callejas-Fernańdez, J.; Martıńez-Garcıá, R.; Hidalgo-A Ä lvarez, R. J. Chem. Phys. 2004, 121, 5468. (54) Two peaks in the distribution were observed of bond lifetimes: one from real bonds, which increases with φp, and another at very small times, arising from particles that have just entered the interaction range but leave it before exploring the well. Nevertheless, the bond length was kept equal to the interaction range for simplicity. (55) Doi, M.; Edwards S. F. The theory of polymer dynamics; Clarendon Press: Oxford, 1986. (56) Rottereau, M.; Gimel, J. C.; Nicolai, T.; Durand, D. Eur. Phys. J. E 2005, 18, 15. (57) Broide, M. L.; Cohen, R. J. J. Colloid Interface Sci. 1992, 153, 493. (58) Fernańdez-Barbero, A.; Schmitt, A.; Cabrerizo-Vı´lchez, M.; HidalgoA Ä lvarez, R. Physica A 1996, 230, 53. (59) Soos, M.; Sefcik, J.; Morbidelli, M. Chem. Eng. Sci. 2006, 61, 2349. (60) Carpineti, M.; Giglio, M.; Dediordio, V. Phys. ReV. E 1995, 51, 590.

Simulations and Stochastic Mean Field Theory - American Chemical

Recommend Documents