Computer Simulation of Diffusing-Wave Spectroscopy of Colloidal

A computational approach has been developed to model diffusing wave-spectroscopy (DWS) behavior in colloidal systems. This model has been applied to t...
0 downloads 0 Views 264KB Size
Download PDF
5856

Langmuir 2000, 16, 5856-5863

Articles Computer Simulation of Diffusing-Wave Spectroscopy of Colloidal Dispersions and Particle Gels Christopher M. Wijmans,† David S. Horne,*,‡ Yacine Hemar,§ and Eric Dickinson| Department of Chemical Engineering, University of Amsterdam, Nieuwe Achtergracht 166, 1018 WV Amsterdam, The Netherlands, Hannah Research Institute, Ayr KA6 5HL, U.K., Institute for Food, Nutrition and Health, Massey University, Palmerston North, New Zealand, and Procter Department of Food Science, University of Leeds, Leeds LS2 9JT, U.K. Received June 15, 1999. In Final Form: March 24, 2000 A computational approach has been developed to model diffusing wave-spectroscopy (DWS) behavior in colloidal systems. This model has been applied to the study of both colloidal dispersions and particle gels. The individual particle dynamics are computed from a Brownian dynamics simulation. Subsequently, a large number of photon paths through the system are generated. These paths, in combination with the mean-square displacements of the particles obtained from the Brownian dynamics simulations, are used to calculate the temporal autocorrelation function g(1)(t). The simulations reproduce the effect of the system thickness on g(1)(t) for a stable colloidal dispersion as found experimentally, as well as the effect of the composition of a bimodal dispersion. Simulations of particle gels show how the gelation process influences the correlation function. These simulations reproduce the profound changes in the correlation function seen experimentally during the gelation process and demonstrate that such effects are primarily due to changes in the particle dynamics rather than in the large-scale topology of the developing network.

1. Introduction Diffusing-wave spectroscopy (DWS) is a relatively new technique that is extremely useful for the study of concentrated colloidal systems, which cannot easily be investigated using traditional dynamic light scattering techniques because of the strongly multiple scattering nature of such systems.1 An important application of DWS is the study of gelation processes. Because DWS measurements depend strongly on the dynamics of the scattering particles, DWS of gelling systems can provide valuable information on the extent and kinetics of the gelation process. Examples of this application of the DWS technique are studies of the gelation of milk induced by enzyme treatment either alone2 or in combination with acidification,3 and measurements of latex particle mobility in a gelatin solution as the solution is cooled through the sol-gel transition.4 In systems such as these, the gelation process slows the particle dynamics very considerably, and this profoundly affects the intensity autocorrelation function of the scattered light. We have recently developed a computer simulation * To whom correspondence should be addressed. † University of Amsterdam. ‡ Hannah Research Institute. § Massey University. | University of Leeds. (1) Weitz, D. A.; Pine, D. J. In Dynamic Light Scattering: The Method and Some Applications; Brown, W., Ed.; OUP: Oxford, 1993; p 652. (2) Horne, D. S. In New Physico-Chemical Techniques for the Characterization of Complex Food Systems; Dickinson, E., Ed.; Blackie: London, 1995; Chapter 11, p 240. (3) Dalgleish, D. G.; Horne, D. S. Milchwissenschaft 1991, 46, 417. (4) Horne, D. S.; McKinnon, I. R. Prog. Colloid Polym. Sci. 1997, 104, 1163.

model of particle gelation by Brownian dynamics.5 In this paper we use this simulation model to predict DWS correlation functions of gelling systems. Before using the model, we test the DWS simulation algorithm by applying the model to the (simpler) case of particles undergoing unrestricted Brownian motion. An approach similar to ours has been taken by Middleton and Fisher6 and by Durian7 to test approximate theories describing DWS. The main difference between these papers and our present simulations is that we explicitly take account of an ensemble of particles, which scatter the incident light, rather than assuming that photons are scattered at random points in space. The spatial positions of the particles and their mobilities are derived from a separate Brownian dynamics simulation. Our main aim is to determine the effect of gelation on the predicted dynamic light scattering in a high volume fraction and therefore multiple scattering situation. In addition, for comparison purposes, we first study the light scattering by a stable dispersion of colloidal particles. Throughout this paper we restrict ourselves to the backscattering geometry, which is commonly used in DWS experiments. Experimentally, this geometry has the advantage that it requires access to the sample only from one side. However, the disadvantage of this geometry compared with the transmission geometry lies in the theoretical analysis of the temporal autocorrelation function. Because many length scales are probed simulta(5) Whittle, M.; Dickinson, E. Mol. Phys. 1997, 90, 739; Wijmans, C. M.; Dickinson, E. J. Chem. Soc., Faraday Trans. 1998, 94, 129; Wijmans, C. M.; Whittle, M.; Dickinson, E. In Food Emulsions and Foams: Interfaces, Interactions and Stability; Dickinson, E., Rodrı´guez Patino, J. M., Eds.; Royal Society of Chemistry: Cambridge, 1999; p 342. (6) Middleton, A. A.; Fisher, D. S. Phys. Rev. B 1991, 43, 5934. (7) Durian, D. J. Phys. Rev. E 1995, 51, 3350.

10.1021/la9907681 CCC: $19.00 © 2000 American Chemical Society Published on Web 06/10/2000

Computer Simulation of Diffusing-Wave Spectroscopy

Langmuir, Vol. 16, No. 14, 2000 5857

neously, it is less apparent that one can apply the diffusion approximation to describe the light transport in the backscattering geometry.1 By conducting simulations (as presented in this paper), one can indeed test the theoretical predictions made using the diffusion approximation. 2. Simulation Method In a colloidal dispersion at relatively low concentrations and short times, the mean-square particle displacement 〈∆r2(t)〉 as a function of time t follows directly from Einstein’s equation: 〈∆r2(t)〉 ) 6Dt, where D is the diffusion constant. However, in a gelled system the particle motion is strongly restricted, so we need to conduct Brownian dynamics (BD) simulations to study the particle kinetics. We consider here a system of N spherical particles in a cubic simulation box of size l 3 with periodic boundary conditions. From these simulations we can extract the mean-square displacement of the particles 〈∆r2(t)〉 as a function of time t. We then construct a slab of material, which we use for the simulated light-scattering experiments by considering a three-dimensional array of periodic images of the system used in the BD simulation. In two directions (x and y) the slab stretches out infinitely, and in the z-direction it has a finite thickness L, which is equal to an integer times the length l of the BD simulation box. (We denote the l 3-sized system used in the BD simulation as the “primary system” to distinguish it from the larger system used in the DWS simulation.) We first describe the BD algorithm, and subsequently we discuss the DWS simulation algorithm. Brownian Dynamics Simulation. We consider N spherical particles (typically, N is of the order 103-104) at a volume fraction φ. To model a stable colloidal dispersion, the interaction between two particles is represented as a steeply repulsive core potential φC, for which we use the following form:

()

φC ) C

σ rij

36

(1)

In this equation c is an energy parameter and rij ) |ri - rj| is the interparticle distance between the interacting particles i and j. The particle diameter σ is used as the basic unit of distance in the simulation model, and the energy scale is defined by the parameter c ) kBT. The time scale is defined by the relaxation time τr, which is the time taken for a particle to diffuse a distance of half the particle diameter: τr ) (3πσ3ηs)/(4 kBT), where ηs is the solvent viscosity. For a particle with diameter σ ) 1 µm in water, the value of τr is 0.58 s. The translational update algorithm is written as a function of the interparticle forces:

∆t ri(t + ∆t) - ri(t) ) Fi(t) + Ri(t,∆t) ξ

(2)

where Fi is the force due to interparticle interactions, ξ is the Stokes friction coefficient, and Ri(t,∆t) is the familiar Gaussian random displacement. A space-filling, percolating network of particles (i.e., a particle gel) is formed by the formation of bonds between the particles. We define a probability PB that a bond is formed between two particles whose intersurface separation is less than a critical distance b1. Such a bond results in an attractive force between bond nodes that are defined on the particle surfaces (at r ) σ/2) at points collinear with the particle centers. This force is described by the following bond potential φB:

(

φB ) b

)

bij - b1 b0

2

for bmax > bij > b1

(3)

where bij is the internode distance. As default values we use b ) 1.0, b1 ) 0.1, b0 ) 0.1, and bmax ) 1.0. Only one bond can be formed between any pair of particles. The initial alignment of the bond with the particle centers soon disappears because of the subsequent motion of the particles. The force due to the bond interaction must be added to the interparticle force in eq 2. The bond interaction also leads to a torque being exerted on the particles. This is accounted for by a rotational update algorithm5 that is analogous to the translational update algorithm of eq 2. In addition to the bond interaction, we sometimes introduce an extra long-range interaction potential φLR, which has a strong influence on the structure of the gel that is formed. This potential has the form

(rc - rij) for rij < rc φLR ) LR (rc - σ)

(4)

and causes a constant force, which is attractive or repulsive (depending on the sign of LR), between two particles at a separation smaller than the long-range cutoff distance rc (we use a default value of rc ) 2.5σ). Whittle and Dickinson5 have demonstrated that the model outlined above generates gel-like structures with fractal properties. The nature of these structures depends strongly on the parameters used, specifically on the bonding probability and the long-range interaction. An attractive long-range interaction leads to a coarse gel structure with a low fractal dimension, whereas a far finer structure with a higher fractal dimension is formed for a repulsive interaction. For example, for PB ) 10-3 the fractal dimension is 1.85 for a long-range interaction parameter LR ) -0.23 and 2.05 for LR ) +0.23. The distribution of pore sizes was also shown to be sensitive to the simulation parameters used. As the value of LR increases, the occurrence of large pores becomes more likely. The parameter values used in this paper have been chosen for consistency with earlier work. Especially, the parameters that determine the bond potential (eq 3) may be expected to have a large influence on the mobility of an individual particle in a gel. We will further address this issue in the Results section. DWS Simulation Algorithm. A photon enters the simulation box at z ) 0. The photon then moves in a straight line, and at regular sampling intervals we test whether the photon has “hit” a particle (that is, come within a distance of σ/2 from a particle center). If the photon does not hit any particle, it continues moving along a straight line, leaves the system at z ) L, and is discarded. If a photon hits a particle, then we have to decide its scattering angle (θ) according to a probability distribution, P(x), defined by the Rayleigh-Gans-Debye formulation for the particle-scattering form factor:

P(x) )

9 (sin x - x cos x)2 6 x

(5)

where x ) 2πσ/λ sin(θ/2) and λ is the wavelength of the light. Unlike the formulation of Middleton and Fisher6 which considered point scatterers and isotropic scattering, our approach introduces scattering-angle probabilities dependent on the finite physical size of the particles. Such probability functions, when integrated over all scattering angles, would yield the total scattering cross-section for an individual particle.

5858

Langmuir, Vol. 16, No. 14, 2000

Wijmans et al.

After being scattered, the photon continues to move in a straight line in its new direction. At regular sampling intervals we test to determine whether the photon has hit another particle, and if so, it is scattered again according to the distribution of eq 5. A photon that reaches the plane z ) L is discarded (forward scattering). Only photons that return to the plane z ) 0 (backscattering) are considered in our calculations. Such photons may reach this plane by many different paths, and the electric field at the detector is taken as the sum of the fields due to each of these paths. The scattering particles, however, are in constant motion, and this causes the length of each path to vary over time. Because the phase of the scattering field that arrives at the detector due to a particular path depends on the length of that path, this phase evolves in time and is the mechanism giving rise to the time dependence of the autocorrelation function. Following the formulation of Middleton and Fisher,6 let us consider a single path R of a photon with a wave vector of magnitude k ) 2π/λ, collected at z ) 0 after having undergone n scattering encounters with n particles. We define kjR - kj-1R as the transfer wave vector of the jth encounter, with kjR as the wave vector after the jth event, 1 e j e n, |kjR| ) k. Again, following Middleton and Fisher,6 we note that the probability distribution for the change in path length during time t is Gaussian, with mean zero and variance

time in the suspension. This fact justifies the principle of time reversal. Comparison with Theory. Pine and co-workers8 have derived an expression for the autocorrelation function within the diffusion approximation (i.e., assuming that the light transport through the sample can be describedas a diffusion process). They used an initial condition given by an instantaneous planar light source at a distance z ) z0 from the illuminated surface. The boundary conditions were specified by requiring that (for t > 0) the net flux of diffusing light into the sample be zero. This leads to the following equation for the autocorrelation function g(1)(t):

{ [x (

)] x [x ( )]}/ ) [ x ] x [ x ]}

z0 6t L + τ l* l* z0 2 6t 6t L cosh 3 τ τ l* l* L 6t L 8t 4 6t cosh 1+ sinh + 3τ l* τ 3 τ l*

g(1)(t) ) sinh

{(

6t τ

(8)

where l * is the transport mean free path of the photons, and τ ) (k2D)-1 for diffusive motion. For infinite slab thicknesses and in the limit t , τ, eq 8 reduces to

[ x]

g(1)(t) ) exp -

z0 l*

6t τ

(9)

n

yR )

(kjR - kj-1R)2〈∆r2(t)〉/6 ∑ j)1

(6)

where 〈∆r2(t)〉 is the mean-square displacement of the particles. Because the scattering particles move slowly compared to the speed of light, the absolute value of the wave vector k ) 2π/λ does not change when a photon is scattered. For paths involving only a few scatterers, this Gaussian form relies on the approximation that the scatterers execute independent random walks. For paths involving many scatterers which undergo more correlated but isotropic random motion, Middleton and Fisher6 argue that the same form will be obtained as a result of the summation over many independent light paths, as long as the correlations are of finite range. Hence for our particle gels we use the value of the mean-square particle displacement averaged over all particles that have the same number of bonds. The procedure described above is repeated until we have collected Nba backscattered photon paths (and in the meantime discarded Nfo forward-scattered photon paths). We then calculate the temporal autocorrelation function g(1)(t) from a summation over all backscattered paths:

g(1)(t) )

〈E(t)E*(0)〉 E(0)2

Nba

)

∑ exp(-yR)

(7)

R)1

In the simulation approach outlined above, the photons enter the system under an angle of 90° over a limited area (a square with length of the order of 101 particle diameters), and they can exit the system at any point. This is the time-reversed equivalent of an experiment with a planar light source and a point detector. A point detector will in reality have a certain physical size, which is reflected by the size of the primary simulation box. The fairly strict constraints on the acceptance angle of the optical fiber to the photomultiplier are reflected by the perpendicular angle of incidence in the simulation. Because the photons travel at the speed of light, they have a negligible transit

3. Results and Discussion We start by considering a system consisting of monodisperse (unbonded) particles at a volume fraction φ ) 0.05. The mean-square displacement of the particles follows directly from the equation 〈∆r2(t)〉 ) 6Dt. We first investigate such a system, rather than a more complicated particle gel, as a test for the DWS simulation method. We use a primary system with N ) 104 particles (l /σ ) 47.1). By conducting independent simulation runs with different particle configurations, we confirmed that the statistical error introduced by the finite size of the primary simulation box is negligible. Figure 1A shows simulated values for g(1)(t) which were computed for a wavelength equal to the particle radius (λ ) σ/2). In this figure, the logarithm of g(1)(t) has been plotted against xt for different system sizes L: L ) 2l (curve “a”), 20l (curve “b”), and 2000l (curve “c”); in all cases, l /σ ) 47.1. For a particle diameter of ca. 1 µm, these values of the system size in reduced units correspond to slab thicknesses of 0.1 mm, 1 mm, and 10 cm, respectively. In parts B and C of Figure 1, similar graphs are shown for a different particle concentration and wavelength. In Figure 1B the volume fraction is twice as high (φ ) 0.1; λ ) σ/2) and in Figure 1C the wavelength is twice as small (φ ) 0.05; λ ) σ/4). The value of g(1)(t) is determined by the paths of the backscattered photons only. For small values of the slab thickness L, this means that only short path lengths are accounted for, where a photon is scattered by only a relatively small number of particles before it returns to the plane of incidence at z ) 0. A large number of photons are lost because they transit the slab and exit from the face parallel to the plane of incidence. As the slab thickness increases, the number of photon-particle encounters increases. Longer path lengths, and also more of them, are taken into account and a larger proportion of the incident photons is backscattered. In the limit of infinite (8) Pine, D. J.; Weitz, D. A.; Zhu, J. X.; Herbolzheimer, E. J. Phys. (Paris) 1990, 51, 2101. (9) Horne, D. S. J. Phys. D 1989, 22, 1257.

Computer Simulation of Diffusing-Wave Spectroscopy

Langmuir, Vol. 16, No. 14, 2000 5859

Figure 2. Path-length distributions of backscattered photons for the system of Figure 1A for L ) 2l (a) and L ) 2000l (b).

Figure 1. Correlation function g(1)(t) plotted against xt/τr for a stable colloidal dispersion. In (A) φ ) 0.05 and λ ) σ/2; in (B) φ ) 0.1 and λ ) σ/2; in (C) φ ) 0.05 and λ ) σ/4. The slab thickness has been varied: L ) 2l for curves “a”; L ) 20l for curves “b”; and L ) 2000l for curves “c”. All simulations use a primary system with 10 000 particles (so that l ) 47.1σ in (A) and (C) and l ) 37.4σ in (B)). In part (A) a curve is also shown which gives g(1)(t) according to eq 9 for z0/l * ) 1.9.

slab thickness, all photons are backscattered and the distribution of path lengths has an infinitely long tail. We find that for L ) 2l , actually only a minority of all incident photons are backscattered: for every 10 photons that return to the plane z ) 0, 150 photons leave the system at z ) L. However, for L ) 20l , on average only 8.4 photons leave the system at z ) L for every 10 photons that return to the plane z ) 0. For L ) 2000l the system virtually represents an infinitely thick slab, and in this case for every 1000 backscattered photons there are, on average, fewer than 8 that are scattered in the forward direction. Figure 2 shows the distribution of path lengths of backscattered photons for L ) 2l and for L ) 2000l. This distribution is expressed in terms of the fraction of collected photons that have suffered the plotted number of scattering events. In the case of the thin slab with L ) 2l , this fractional distribution has a clear maximum at approximately 25 scattering events, after which it falls off relatively quickly. Less than 8% of all backscattered

photons have been scattered more than 50 times. The longest path length found in the simulation was 148 scattering events. In contrast, for L ) 2000l the distribution has a far more extended shape, with more than half of all backscattered photons being scattered more than 200 times. In this case, 8% of the photons are scattered even more than 10,000 times. The influence of the contribution from the longer path lengths (i.e., those with the higher number of scattering events), as a result of increasing slab thickness, can be invoked to explain the changes in curvature in plots of g(1)(t), as seen in Figure 1. The nonexponential character of the correlation function betrays the existence of multiple time constants, indicative of differing decay rates in the paths of various lengths. Paths involving large numbers of scattering events decay the most rapidly, because each individual particle in the chain needs to move only a short distance and takes a short time to do so. By contrast, paths made up of few scattering events decay more slowly, as the smaller number of particles must move relatively further to reach an accumulated path-length change of one wavelength. Thus, for small values of L, the data in Figure 1A show a strong curvature at low t, which disappears for larger values of L. This effect of the slab thickness is in good agreement with experimental findings.9 In Figure 1A another curve has been added, which shows g(1)(t) as predicted by eq 9. A value of 1.72 has been used for z0/l * to fit the slope of this line to that of the simulated curve. In addition, log(g(1)(t)) has been shifted upward slightly (by a value of 0.04) to get the best fit between both curves. The residual curvature in the simulated curve shows up in the small difference between this curve and the fitted line at small t. The contribution of photon paths of different length towards the correlation function is further explored in Figure 3. The curve in this figure that is labeled as “total” is the same as curve “a” in Figure 1A. This curve was obtained by summing over all backscattered photon paths (eq 7). We have also drawn curves that show the result when the summation is extended only over those paths that have 10 or fewer scattering events (“short”), and when the summation is over those paths that have more than 50 scattering events (“long”). The decay length of the total signal is only slightly less than that of the short path lengths only. However, the short path-length signal does show a very prominent concavity at short times. In Figure 1B the particle concentration has been increased by a factor of 2 to φ ) 0.1. Qualitatively, the same trends are found as for φ ) 0.05. However, g(1)(t) depends slightly less strongly on the slab thickness. As L increases, g(1)(t) decreases less strongly than in the case of φ ) 0.05. Because of the increased particle concentration,

5860

Langmuir, Vol. 16, No. 14, 2000

Figure 3. The effect of the photon path-length distribution on g(1)(t). The curve labeled as “total” is the same as curve “b” in Figure 1A, which is the correlation function for a system with φ ) 0.05 and λ ) σ/2. The “short” curve shows the correlation function as it is calculated when only photons that have been scattered fewer than 11 times are taken into account. The “long” curve is calculated using only photons that have been scattered more than 50 times.

Figure 4. Correlation functions calculated for different reflectivities of the z ) 0 boundary as defined by the parameter R (see text). Parameters: φ ) 0.05; λ ) σ/2; L ) 20l ; R ) 0 and 0.99 (the curve for R ) 0 is the same as curve b in Figure 1A).

a photon is more likely to be backscattered in a system with the same slab thickness. The effect of the wavelength on g(1)(t) can be seen in Figure 1C. In this figure φ ) 0.05, and the wavelength used is equal to half a particle radius (λ ) σ/4). A smaller wavelength (in relation to the particle size) means that the average angle over which a photon is scattered is smaller (eq 5). As a consequence, one has to go to larger slab thicknesses to reach the limiting behavior for g(1)(t), and an increase in the slab thickness leads to a larger decrease of g(1)(t) when the wavelength is shorter. As light with a smaller wavelength is scattered over a smaller angle, the backscattered photons have, on average, been scattered a larger number of times. This causes the faster decay of g(1)(t) as compared to the correlation function in Figure 1A. Up to now we have assumed that when a photon reaches either of the planes z ) 0 or z ) L, it subsequently exits the system. However, if the photon is reflected, that will have an effect on the photon path-length distribution. We have considered the case that there is a probability R that a photon is reflected at the z ) 0 boundary. In Figure 4 the curve for R ) 0 is the same as curve “b” in Figure 1A (φ ) 0.05, L ) 20l ). The correlation function has also been simulated for the same system with R ) 0.99. Increasing the value of R shifts the photon path length distribution toward longer path lengths. As a consequence g(1)(t)

Wijmans et al.

Figure 5. Correlation functions for a bimodal dispersion with a particle size ratio of 3. The wavelength is equal to the radius of the smaller particles. In all cases the total particle volume fraction φ ) φ1 + φ2 ) 0.1. The volume fraction of the larger particles φ2 is shown in the figure for each curve. The characteristic time τr used to define the time axis is the characteristic diffusion time of the smaller particles.

decreases. At small t the curve of log(g(1)(t)) versus xt becomes slightly convex rather than concave. At longer times the curves are roughly parallel lines for different values of R. Bimodal Colloidal Dispersions. In general it is believed that, for a polydisperse system, diffusing-wave spectroscopy will give only an average relaxational mobility for the particle size distribution. This is because the diffusing photon can encounter particles of any contributing size on any contributing pathway. In Figure 5, however, we consider one more system consisting of unbonded particles, where we explore the influence of “polydispersity” in a dispersion with two distinct particle populations, one with particle diameter 3 times the other (σ and 3σ, respectively). The total particle volume fraction, φ, which is the sum of the volume fraction of the small particles, φ1, and that of the larger ones, φ2, is kept at a constant value of 0.1. For φ1 ) 0.1 and φ2 ) 0 we have the same system as studied in Figure 1B. For the opposite case that all particles are large particles (φ1 ) 0 and φ2 ) 0.1) the function g(1)(t) decays more slowly and exactly reflects in its linear t1/2 behavior the lower diffusional mobility of those large particles. For these larger particles, however, the curvature in g(1)(t) at short decay times persists to much greater slab thicknesses. This result is consistent with our observations, plotted in Figure 1C, where, although we halved the wavelength, the same effect would have been seen by doubling particle size. In mixed systems the response of g(1)(t) to changes in the system composition is very nonlinear. If a small fraction (5%) of the larger particles is replaced by the same volume fraction of smaller ones, this has a large effect on g(1)(t). However, replacing a far larger fraction (50%) of small particles by large ones has only a very small effect. This same nonlinearity has been observed experimentally for mixed systems of polystyrene latex particles with a similar particle size ratio.9 The results can only be interpreted, however, in artificial systems such as these, where the input particle size distribution is already known. DWS yields only a single relaxation time. Attempting to extract more than one piece of information is not physically feasible. Comparison with Experiment. Finally, in Figure 6 we compare experimental and simulated data for a latex dispersion. The circles are measurements of a 0.25-mmthick sample of a 1.0 vol % dispersion of polystyrene spheres (330 nm diam) using a He-Ne laser (taken from

Computer Simulation of Diffusing-Wave Spectroscopy

Figure 6. Comparison between experimental data (circles, sample thickness 0.25 mm; squares, thickness 2 mm), taken from ref 9 and simulations (full curves) for a 1.0 vol % polystyrene latex with a particle diameter of 330 nm and a wavelength (in the simulation) of λ ) 1.44σ.

Figure 7. Average mean-square displacements 〈∆r2(t)〉 of particles in a gel (φ ) 0.05) formed with PB ) 10-3 and LR ) 0. The age of the gel is ta/τr ) 8.8 (9) and 88 (b). The line shows the mean-square displacement of a freely diffusing particle.

ref 9). The squares are measurements on the same latex dispersion with a sample thickness of 2 mm. The simulation used a primary system with 2000 particles (l ) 47.1σ; L ) 16l, and L ) 128l ) and a wavelength of ) 1.44σ. The data in this figure show a good agreement between experimental points and simulation (solid lines). One should bear in mind that the experimental data have always to be normalized, which leads to a (rather arbitrary) vertical shift of the data. Therefore, one should primarily compare the slopes of the curves rather than their absolute values. The slopes do indeed show a good agreement. Simulations of Particle Gel Behavior. We now turn our attention to a gelled system. We form a particle gel starting with a random dispersion of unbonded particles at a volume fraction φ ) 0.05. We use a bonding probability PB ) 10-3 for a time step ∆t/τr ) 10-6 to form a suitable network of bonded particles.5 After a time ta (the gel age), the bonding probability is set to zero, stopping further bond formation and allowing us to study the dynamical properties of the particles in the network at that level of bonding. In Figure 7 the average mean-square particle displacement 〈∆r2(t)〉 is shown for two different values of ta. The line shows the mean-square displacement of a freely diffusing particle. As particles become bonded, their mobility decreases. For ta/τr ) 88, the value of 〈∆r2(t ) 10τr)〉 is only 4% of the value for a freely diffusing particle. Of course, the actual mobility of a particle in the gel network depends on the parameters used in the bond

Langmuir, Vol. 16, No. 14, 2000 5861

Figure 8. Correlation functions for a gel formed with PB ) 10-3 and LR ) 0. The curves “a”-“f” are for increasing gel age ta. Curve “a” is for ta ) 0, which corresponds to a dispersion consisting of unbonded particles. Curve “b” ta/τr ) 2.2; curve “c”, ta/τr ) 8.8; curve “d” ta/τr ) 44; curve “e”, ta/τr ) 88; curve “f”: ta/τr ) 177. The data points for curves “b”-“f” are the simulated values of g(1)(t) computed at values for the meansquare particle displacements identical to those employed in Figure 7. The curves themselves are interpolations between these points. The simulations use a primary system with 1000 particles and φ ) 0.05 (i.e., l ) 21.8σ). The slab thickness is L ) 50l ) 1.1 × 103 σ.

potential φB of eq 3. Here we used the (relatively arbitrary) values b ) 1.0, b1 ) 0.1, and b0 ) 0.1. Using a stronger bond potential, for example b ) 5.0 and b0 ) 0.01, gives smaller values for the mean-square displacement. This effect of the bond strength is most important at relatively small times (≈10-2τr) and is relatively insignificant at long times (≈10τr). At long times the displacement of an individual particle is probably mainly determined by the movement of a cluster of surrounding particles. However, the choice of parameter values for the bond potential will certainly have an effect on the initial decay of the DWS autocorrelation function. The values shown in Figure 7 were found by averaging over all particles in the system. The average mean-square displacement of a particle decreases as the number of bonds between that particle and other particles increases. When simulating DWS experiments, we therefore use a value for 〈∆r2(t)〉 for each particle that is found by averaging over all particles with the same number of bonds. The gelling system which gave rise to the data plotted in Figure 7 has also been used to compute the data for the correlation functions, g(1)(t), derived at the different gel ages indicated and shown in Figure 8. The data for the unbonded particles (i.e., ta ) 0, curve “a”) are similar to those in Figure 1A. After gelation, the correlation function decays more slowly because of the decreased particle mobility. The initial decay of g(1)(t) is relatively unaffected by the gelation. At longer times the effect of the gelation becomes more pronounced, and as the gel age increases, the relaxation time of the correlation function also increases. We explained above that a different choice of parameter values for the bond potential will lead to different values for g(1)(t). As the gel age increases, the shear modulus also increases. Previous simulations5 showed that the high-frequency shear modulus is a linear function of the number of interparticle bonds. The statistical error in the data presented in Figure 8 is mainly due to two sources. First, the finite length of the Brownian dynamics simulations leads to an uncertainty in the values found for the mean-square particle displacements. For small times t this uncertainty is completely

5862

Langmuir, Vol. 16, No. 14, 2000

Figure 9. Correlation functions for a gel with PB ) 10-3, LR ) 1, and LR ) -1. Data are also shown for these two systems using the same values for 〈∆r2(t)〉 for all particles in both gels (× and O).

negligible, but it is more important for values of t that become of a similar order of magnitude as the simulation length. Second, the structure of a simulated gel is itself subject to fluctuations. As the gelation is necessarily simulated using a finite system (103 particles), there will be an uncertainty in the average number of bonds per particle after gelation. Simulations using independently generated particle configurations show that the statistical error in the curves of Figure 8 is negligible compared to the differences between two different curves. Figures 7 and 8 were calculated using a long-range interaction parameter LR equal to zero. If the gelation takes place in a system with a different value for this parameter, the resulting gel network will have a different structure. Simulations have shown5 that negative values of LR give a coarser structure characterized by a lower fractal dimension, and positive values give a less coarse structure with a higher fractal dimension. In Figure 9 simulated data are shown for gels with LR ) -1 (coarse, low fractal dimension) and LR ) 1 (fine-structured, higher fractal dimension). In the former case, the attractive interparticle forces promote the formation of bonds, thus decreasing the average particle mobility in comparison with a gel for which LR ) 0. The result is that the relaxation time of the correlation function increases. In contrast, when LR ) 1, fewer bonds are formed and the relaxation time is shorter than for the case LR ) 0. Of course, differences in the correlation functions for different values of LR can be caused by differences in both the gel structure and in the particle dynamics. To determine the relative importance of these two effects, we conducted the following simulation. We used the spatial positions of the particles in the gels formed with LR ) -1 and LR ) 1 for a simulated DWS experiment. However, in both cases we used exactly the same values for the mean-square particle displacement. For this to be possible the gels must have different ages, and hence different extents of bonding and network structures. In fact we selected the same value for 〈∆r2(t)〉 to all particles, regardless of the number of bonds on each particle. Any differences in g(1)(t) between the two gels must now be caused by the different structures of these two gels (as expressed, for example, by differences in the radial pair correlation function). In Figure 9 one can see that the values of g(1)(t) are virtually the same in both cases. This leads to the conclusion that the differences between the gels are caused by differences in the particle dynamics and not by structural differences.

Wijmans et al.

Figure 10. Comparison of the storage shear modulus and the DWS correlation function of the gel studied in Figure 7, as a function of the gel age ta. The squares (0) show G′ at a reduced frequency fτr ) 113 (left axis), and the circles (O) show log g1 for t/τr ) 0.4 (right axis).

One note of caution should be made regarding the time scales of the gelation process and the decay of the correlation function. Not too much importance should be attached to the absolute values found in the simulations. In our Brownian dynamics simulations, the gelation kinetics are primarily determined by the parameter PB, which gives the probability for two particles to form a bond when they are at a small separation. In addition, the interaction potential will have an important influence on the gelation. We have chosen to use a simple pair potential which only accounts for the repulsion between two particles when they start to overlap. In principle, it is feasible to include an activation energy that must be overcome when two particles approach each other,10 which slows the gelation. In this way it would certainly be possible to get a better fit between the simulated gelation and the kinetics of a real system, but this was not the primary aim of our investigation. The main result of the simulations is the strong increase in the correlation function relaxation time as the gelation progresses; this increase is not affected by the precise details of the gelation kinetics. Finally, in Figure 10 we compare the shear rheology and the DWS data of a gel as a function of its age. This graph shows the storage modulus G′, computed by Brownian dynamics simulation, at a reduced frequency of fτr ) 113 (for particles with a diameter of 1 µm, this corresponds to a frequency of 200 Hz). The storage modulus increases with gel age ta. It has previously been shown5 that G′ is a linear function of the number of bonds between particles. The graph also gives the corresponding values of log(g(1)(t)) for t/τr ) 0.4. The behavior of the correlation function with gel age reflects the corresponding decrease in the particle mobility and its mirroring in the increasing elasticity of the gel. 4. Conclusions We have applied a simple model to simulate diffusingwave spectroscopy measurements of colloidal systems. We have used this model to study DWS experiments, in the backscattering mode, on stable dispersions and particle gels. (10) Mellema, M.; van Opheusden, J.; van Vliet, T. In Food Emulsions and Foams: Interfaces, Interactions and Stability; Dickinson, E., Rodriguez Patı´no, J. M., Eds.; Royal Society of Chemistry: Cambridge, 1999; p 176.

Computer Simulation of Diffusing-Wave Spectroscopy

The simulations of dispersions reproduce the characteristic effect of the slab thickness on the correlation function g(1)(t) that is found experimentally. Furthermore, the simulations show the same nonlinearity of g(1)(t) as afunction of system composition that has been found experimentally for bimodal dispersions. The simulations of particle gels show how the decay of g(1)(t) slows dramatically during the gelation process. This is a direct consequence of the self-diffusion of the gelling particles, which become increasingly more restricted in

Langmuir, Vol. 16, No. 14, 2000 5863

their movement as the gelation progresses. Differences between gels with different structures have been shown to be due primarily to differences in the individual particle dynamics rather than to structural differences. Acknowledgment. This work was supported by Contract FAIR-CT96-1216 of the European Union Framework IV Programme. LA9907681