Joost
de Graaf†
*^{a},
Wilson C. K.
Poon†
^{b},
Magnus J.
Haughey†
^{b} and
Michiel
Hermes†
^{c}
^{a}Institute for Theoretical Physics, Center for Extreme Matter and Emergent Phenomena, Utrecht University, Princetonplein 5, 3584 CC Utrecht, The Netherlands. E-mail: j.degraaf@uu.nl
^{b}SUPA, School of Physics and Astronomy, The University of Edinburgh, King's Buildings, Peter Guthrie Tait Road, Edinburgh, EH9 3FD, UK
^{c}Soft Condensed Matter, Debye Institute for Nanomaterials Science, Utrecht University, Princetonplein 5, 3584 CC Utrecht, The Netherlands

Received
6th August 2018
, Accepted 23rd November 2018

First published on 27th November 2018

Colloidal particles with strong, short-ranged attractions can form a gel. We simulate this process without and with hydrodynamic interactions (HI), using the lattice-Boltzmann method to account for presence of a thermalized solvent. We show that HI speed up and slow down gelation at low and high volume fractions, respectively. The transition between these two regimes is linked to the existence of a percolating cluster shortly after quenching the system. However, when we compare gels at matched ‘structural age’, we find nearly indistinguishable structures with and without HI. Our result explains longstanding, unresolved conflicts in the literature.

Kinetics do not determine the equilibrium structure of a system. However, in colloidal gels, which are arrested states, the kinetics of aggregation crucially affect structure formation and therefore mechanical properties.^{3} Hydrodynamic interactions (HI) clearly affect kinetics and hence may lead to a lack of predictive power for numerical simulations of gelation that do not account for these. Simulating a system with HI is, however, far from trivial, as historically this has been the most computationally expensive part of such a simulation algorithm.^{4,5} Concerted effort over the last two decades comparing simulations with and without HI has nevertheless suggested a plethora of effects.

The earliest forays into this topic considered the clustering of a few particles into “fractal” aggregates. Tanaka and Araki^{6} showed that without HI, two-dimensional (2D) Lennard-Jones-like (LJ-like) particle aggregation leads to more compact clusters. Whitmer and Luijten^{7} extended this to three-dimensional (3D) clustering for particles interacting via an Asakura-Oosawa-like (AO-like) potential, finding the diffusion-limited cluster aggregation (DLCA) result at high interaction strength.^{8}

Turning from clusters to gels, Yamamoto et al. clarified the hydrodynamic stabilization of 2D LJ-type gels formed by quenching to zero temperature.^{9} However, in 3D (colloid volume fractions ϕ = 0.173 and 0.307) they found only a minor quantitative effect of HI, observing that structural evolution does not cease in 3D as it does in 2D.^{9} Furakawa and Tanaka^{10} investigated the time evolution of 3D AO-like gels with ϕ = 0.13 and ΔU = 4.2k_{B}T and 8.4k_{B}T quenches, with ΔU the contact attraction strength, k_{B} Boltzmann's constant, and T the temperature. HI shifted the gelation line and changed the time dependence of the correlation length in the gel, ξ ∼ t^{x}, from x = 1/3 to 1/2.

Recently, Varga and collaborators performed large-scale simulations of colloidal gelation using a highly accelerated hydrodynamics algorithm.^{11} They contrasted systems with and without HI over a wide range of ϕ and ΔU, with both AO-like potentials^{12} and a combination of short-ranged attraction and long-ranged repulsion,^{13} finding that HI could lower the gelation line,^{12} slow down compaction and speed up coagulation.^{13} Royall et al. reported a similar lowering of the gelation line with the inclusion of hydrodynamic interactions, which they related to a change in the local structure of colloidal gels with the inclusion of HI.^{14}

In this work, we consider gelation with and without HI for colloids with short-ranged attractions (ΔU = 5k_{B}T and 10k_{B}T) and 0.075 ≤ ϕ ≤ 0.225. HI are accounted for using the lattice-Boltzmann (LB) method,^{15} which enables us to simulate systems comparable or even larger than those of Varga et al.^{12,13} We use the void volume (VV)^{3,16} to study the structure of the holes in the gel. This allows precise analysis of the dynamics of gelation and identification of the onset of aging, which agrees with the one obtained from the evolution of the correlation length. We find that HI indeed speeds up gelation and slightly steepens the VV power-law scaling for low ϕ. For ϕ ≳ 0.16, however, systems with HI exhibit the same power-law exponent as their non-HI counterparts and gelation is slowed down. We relate this crossover to the presence of large clusters immediately upon quenching into the spinodal region of the phase diagram, which rearrange quickly to achieve percolation. We demonstrate clearly that aging is, however, not related to percolation and typically sets in at a much later time.

We also study the structure of gels formed with and without HI. However, in contrast with most simulations to date, we do not compare the structure at equal-and-fixed time. Instead, we compare gels at equivalent points during their evolution, in which case, surprisingly, gels formed with and without HI cannot be readily distinguished by their structure factor. More sensitive VV analysis reveals a small influence of HI on structure, but the difference is much smaller than what can be achieved by minor variations in the gel age, ϕ, or ΔU.

Note that we focus throughout on ϕ ≥ 0.075. At lower ϕ, gel formation is so slow that other effects such as gravity and convection will dominate the dynamics of any real system. Moreover, our method does not give enough statistics over long enough time to determine the effect of HI on the position of the gelation boundary. Nonetheless, our results suffice to show that HI have a large effect on the speed at which the quenched system evolves, but the path that it takes through state space is nearly unaffected. We discuss how these findings relate to and reconcile existing literature on HI in colloidal gelation.

We match the bulk diffusivity of our particles between the NH and WH simulations by setting k_{B}T = 1, the particle diameter σ = 1, and the viscosity of our suspending medium to η = 10, so that D = k_{B}T/(3πησ) = 1.1 × 10^{−2}σ^{2}τ^{−1}, with unit time and particle mass m = 1. We express all times in terms of the Brownian time t_{B} = σ^{2}/(4D) = 23.5τ for a particle to diffuse its own radius.

A generalized LJ interaction with high exponents was used to model short-ranged depletion attraction:^{23}

(1) |

We equilibrate our systems using a purely repulsive potential: ε = 10k_{B}T, c = 1, and r_{c} = σ. We then quench into the spinodal region by instantaneously switching to an attractive interaction: ε = 10k_{B}T or 5k_{B}T, c = 0, and r_{c} = 1.5σ, so that ΔU = 10k_{B}T or 5k_{B}T. The reduced second virial coefficient^{24} for U^{gen}_{LJ}(r) is

(2) |

We use a time step Δt = 0.002 and equilibrate each sample for 50t_{B} before quenching to simulate gelation and aging for another 550t_{B} for ϕ > 0.1 and 5500t_{B} for ϕ ≤ 0.1. We prepare our systems in a cubic box with edge L = 48σ (with σ = LB spacing) for 12 points in ϕ ∈ [0.075,0.225], so that we simulate between N = 15841 and 47523 particles. For each data point, we run 10 independent simulations using the Molecular Dynamics package ESPResSo^{25,26} to obtain sufficient statistics.

At low enough ϕ, small clusters form initially from particles already in close contact. These clusters diffuse, bonding on contact, leading to a fractal-like network. Such rearrangements are sped up by HI, which enhance cluster translational and rotational mobility by allowing collective effects, as found previously.^{12} At higher ϕ, particles travel only a short distance to bond. Now, compaction dominates, i.e., the approach of colloids to form a bond contributes more to the formation of a gel network than motion of (small) clusters. This approach is slowed down by the need to squeeze fluid from between particles, again agreeing with recent work.^{12}

Fig. 2 Structural properties for gelling systems with ΔU = 10k_{B}T, ( = NH, = WH), and for three ϕ as indicated by the labels and use of line type; the result for ΔU = 5k_{B}T is analogous. (a) Void volume (VV) probability density functions (PDFs) at the aging time t_{a}. Error bars indicate the standard error, which is typically smaller than the symbol size. One in five data points is shown to improve the presentation, but the lines that serve as guides to the eye go through every data point. (b) The time dependence of the void growth length, λ_{VV}, see eqn (3), again only one in five data points is shown. (+,×) and (■,□) = percolation (t_{p}) and aging onset times, respectively. The number in parentheses indicates the factor by which λ_{VV} is multiplied to separate the data. The thin dashed black lines and numbers indicate power-law fits and associated exponents (with error ±0.005), respectively. Only for ϕ = 0.075 is there an appreciable difference in the NH/WH initial exponent. (c) Structure factors S(q) at t_{a}, corresponding to the data in (a). We only show error bars up to qσ = 1 for every third data point to improve the presentation, the curves connect all data points; the region of the central peak is not shown qσ < 0.2. (d) The peak position of the structure factor q_{m} as a function of time. The three WH/NH data sets have been shifted up by the amount indicated in parentheses to prevent overlaps. The dashed line indicates a power law with exponent −0.04 ± 0.005. |

Our low-noise PDF data over five orders of magnitude show small but systematic differences between the WH and NH cases. At t_{a}, the WH systems at low ϕ have more small and fewer big holes than the NH systems, while the trend reverses for high ϕ. The VV PDFs contain information that is complementary to the structure factor S(q); we return to this shortly.

From the PDFs, we compute the mean, 〈VV〉, whose evolution gives a time-dependent length scale

λ_{VV}(t) ≡ (〈VV〉(t) − 〈VV〉(0))^{1/3}, | (3) |

There appears to be two power-law regimes λ_{VV} ∝ t^{x}. Initially, the WH (x ≈ 0.40) and NH (x ≈ 0.36) systems are slightly different at the lowest ϕ, but for higher ϕ their power laws are identical within errors, x ≈ 0.35 and 0.33 respectively. This regime is associated with the gelation process. The second regime, which is less well developed at the lowest ϕ, has exponent 0.05 ± 0.005, which we associate with aging.

We define the aging time, t_{a}, as the crossover between these two regimes. In practice, this is done by fitting a local power-law t^{x} to obtain a running exponent x(t) and using the criterion x = 0.08 for the onset of aging, Fig. 2b (■,□). The percolation times, t_{p}, identified earlier (+,×), are also shown. In each pair of data sets, we see that the NH and WH systems have the same λ_{VV} when they successively reach t_{p} and then t_{a}. In other words, the gels have (nearly) the same structures when they have reached a comparable evolutionary stage. Observations for ΔU = 5k_{B}T are qualitatively similar, also see Fig. 2a.

Fitting a time-dependent sequence of S(q)‡ yields q_{m}(t), Fig. 2d, which can be directly compared to small-angle scattering experiments,^{23} in which gelation is associated with the sudden slowdown in the decrease in q_{m} with time. The t_{a} obtained using the VV analysis plotted in Fig. 2d (■,□) matches this transition point well within errors, showing that it is the onset of aging rather than percolation (+,×) that identifies gelation. The negative of the power law exponents obtained from our void-volume analysis, see Fig. 2b, is consistent with the corresponding regimes in Fig. 2d within the error, see the respective dashed power-law fits.

Fig. 3 Differences in the dynamics of gelation. The times of interest, t_{p} and t_{a}, for percolation and aging, respectively, as a function of the colloid volume fraction ϕ for (a) ΔU = 5k_{B}T and (b) ΔU = 10k_{B}T. The dashed (NH) and solid (WH) curves connect the points belonging to a set of data and serve as guides to the eye. The standard error is indicated using error bars, but is typically small. The gray vertical lines indicate the ϕ values, for which there is a crossover, the light-gray field indicates the error therein. Finally, the thick black dashed lines show the diffusion-limited cluster aggregation (DLCA) prediction for the gel time, t_{DLCA}, see eqn (4), with prefactor 1 and fractal dimension d_{f} = 1.8. |

It is clear that HI speed up the system's dynamics below a certain ϕ and slow it down above this value. The crossover occurs at ϕ = 0.095 ± 0.002 for t_{p} and this effect is pronounced, while for t_{a} we find a crossover at ϕ = 0.16 ± 0.01 and 0.16 ± 0.002 for weak and strong gels, respectively. Interestingly, percolation occurs at the same time independent of ΔU, but t_{a} is significantly increased upon lowering ΔU, with t_{a} being substantially larger for systems without HI at low ϕ.

We have also plotted in Fig. 3 the time to form a spanning cluster in ϕ → 0 diffusion-limited cluster aggregation (DLCA)^{8} approximation:

t_{DLCA} ∝ t_{B}ϕ^{−df/(3−df)}, | (4) |

〈VV〉 is mainly sensitive for the structure of the voids (and by extension the gel), while 〈NN〉 probes the arrest inside the branches of the gel. The systems start with a low average number of neighbors and a low average void volume: 〈VV_{0}〉 ≈0.024ϕ^{2} holds over the entire range of considered ϕ. As expected both parameters quickly increase as the gel forms. The path taken through state space is initially nearly insensitive to HI. After longer times, small differences between the NH and WH systems appear, predominantly for weak low-ϕ gels, with the latter having a slightly larger 〈VV〉. For high ϕ and/or ΔU, these differences are small or insignificant.

Firstly, on the most coarse-grained level, the pre-existence of such large-scale structures is the basic physics behind the observation that, despite the large difference in dynamics, systems with and without HI follow nearly identical paths through state space, Fig. 4. Secondly, it explains the poor predictive power of eqn (4), which can only be expected to (and indeed does^{7}) hold in the limit of very low ϕ, where long-range diffusion dominates, and for large ΔU, where rearrangements after initial bonding are irrelevant. Finally, the crossover from the pre-existence of nearly-percolating structures to a pre-existing spanning network at ϕ ≈ 0.1 explains why we see a change in the scaling of t_{p}(ϕ) at the same ϕ, Fig. 3. Only very small local movements of these pre-existing large clusters are needed to bring about percolation. The physics here are dominated by the need to squeeze out solvent between particles that eventually bridge such clusters, so that HI slow down the dynamics, as observed.

Once a spanning network exists, further growth of the gel structure occurs by accretion of clusters and single particles onto this backbone, depending on the volume fraction considered, until eventually the aging dynamics set in. The crossover in t_{a} occurs when ϕ ≳ 0.16, because at lower volume fractions the system reaches the aging state through exploration of pre-arrested configurations by localized cluster movements. That is, at ϕ ≈ 0.1, collective effects are still important to reach the aging state, whereas for ϕ ≳ 0.16 the squeeze-flow contributions dominate the system percolating and aging.

First, there has been persistent reports in the literature of two different modes of gelation: equilibrium and non-equilibrium, with the former attributed to percolation^{32–34} and the latter to arrested phase separation.^{1,23,35–37} There have been recent attempts to unify the two pictures by simulations^{2} and by experiments,^{38} with the latter showing a crossover for percolation to arrested phase separation at ϕ ≈ 0.2 for nano-emulsions based on rheological measurements. Our results do not directly address this issue. However, they do point to the importance of whether there is a pre-existing percolating cluster at the beginning of the aggregation process.

Second, our results reconcile the apparently contradictory findings of Yamamoto et al.^{9} and Furakawa and Tanaka.^{10} The former work reports a minor effect of HI, because their ϕ range is in the compaction regime. The latter, and more recently Varga et al.^{12} and Royall et al.,^{14} report substantial structural differences at intermediate ϕ when gels are compared at equal times. We find the same when we compare systems at equal absolute times, rather than at equal structural times, especially when we do so for t < t_{a}.

Third, contrary to Furakawa and Tanaka, we find gels at ϕ = 0.075 in both WH and NH systems at high and low ΔU. This may reflect their protocol of comparing WH and NH systems at equal times, and/or the closeness to the gelation line of their gels. Neither do we observe a crossover from t^{−1/3} to t^{−1/2} in q_{m}(t) upon including HI. Instead, we find a weak power-law aging in the long-time regime (t > t_{a}). Furakawa and Tanaka's finding here may reflect their limited data range. We can fit a slope of −1/3 to a limited range of our q_{m} data. While our data range is still limited in the long-time limit, our finding of slow aging agrees with Yamamoto et al.'s observation of a continuous gel evolution in 3D systems.^{9}

Fourth, we analyzed the gel structure using the three-point correlation approach introduced by Royall et al.^{14} However, at equal structural times, we did not observe a difference between the local structure of systems with and without HI. The three-point correlation functions were identical within the error bar. In addition, dividing these three-point functions by the regular pair correlation function to highlight certain features, as suggested by Royall et al., proved uninformative, due to rather large errors this process induced. The pronounced effect on local structure, as indicated by the authors of ref. 14, could be due to their measurements having closer proximity to the gel line.

Finally, our short-time q_{m} data does not support the power-law dependence found by Poon et al.:^{23} our limited data is consistent with a stretched exponential fit (not shown here). Clearly the mapping between the evolution of λ_{VV} and q_{m} breaks down at small times.

The results presented here suggest several opportunities for follow-up studies. It would be worthwhile to apply this methodology to explore states closer to the gelation line and to gels with other kinds of inter-particle interactions, such as short-ranged attraction coupled with medium-ranged repulsion. Repeating our analysis using other numerical methods for hydrodynamic interactions, which make different approximations for the solution of Stokes' equations, would also be beneficial and may give closure to the discussion on the effect of hydrodynamic interactions on gelation.

- S. Manley, H. Wyss, K. Miyazaki, J. C. Conrad, V. Trappe, L. J. Kaufman, D. R. Reichman and D. A. Weitz, Phys. Rev. Lett., 2005, 95, 238302 CrossRef CAS.
- E. Zaccarelli, J. Phys.: Condens. Matter, 2007, 19, 323101 CrossRef.
- N. Koumakis, E. Moghimi, R. Besseling, W. Poon, J. Brady and G. Petekidis, Soft Matter, 2015, 11, 4640 RSC.
- G. Batchelor and J. Green, J. Fluid Mech., 1972, 56, 401 CrossRef.
- J. Brady and G. Bossis, Annu. Rev. Fluid Mech., 1988, 20, 111 CrossRef.
- H. Tanaka and T. Araki, Phys. Rev. Lett., 2000, 85, 1338 CrossRef CAS.
- J. Whitmer and E. Luijten, J. Phys. Chem. B, 2011, 115, 7294 CrossRef CAS.
- T. Witten Jr. and L. M. Sander, Phys. Rev. Lett., 1981, 47, 1400 CrossRef.
- R. Yamamoto, K. Kim, Y. Nakayama, K. Miyazaki and D. R. Reichman, J. Phys. Soc. Jpn., 2008, 77, 084804 CrossRef.
- A. Furukawa and H. Tanaka, Phys. Rev. Lett., 2010, 104, 245702 CrossRef PubMed.
- J. Swan and G. Wang, Phys. Fluids, 2016, 28, 011902 CrossRef.
- Z. Varga, G. Wang and J. Swan, Soft Matter, 2015, 11, 9009 RSC.
- Z. Varga and J. Swan, Soft Matter, 2016, 12, 7670 RSC.
- C. Royall, J. Eggers, A. Furukawa and H. Tanaka, Phys. Rev. Lett., 2015, 114, 258302 CrossRef PubMed.
- B. Dünweg and A. Ladd, Adv. Polym. Sci., 2009, 221, 89 Search PubMed.
- M. Haw, Soft Matter, 2006, 2, 950 RSC.
- E. J. Hinch, J. Fluid Mech., 1975, 72, 499 CrossRef.
- D. Röhm and A. Arnold, Eur. Phys. J.: Spec. Top., 2012, 210, 89 Search PubMed.
- P. Ahlrichs and B. Dünweg, J. Chem. Phys., 1999, 111, 8225 CrossRef CAS.
- D. Röhm, S. Kesselheim and A. Arnold, Soft Matter, 2014, 10, 5503 RSC.
- M. Bybee, Hydrodynamic simulations of colloidal gels: Microstructure, dynamics, and rheology, University of Illinois, Urbana-Champaign, 2009 Search PubMed.
- N. Y. C. Lin, B. M. Guy, M. Hermes, C. Ness, J. Sun, W. C. K. Poon and I. Cohen, Phys. Rev. Lett., 2015, 115, 228304 CrossRef.
- W. C. K. Poon, A. Pirie and P. L. Pusey, Faraday Discuss., 1995, 101, 65 RSC.
- M. Noro and D. Frenkel, J. Chem. Phys., 2000, 113, 2941 CrossRef CAS.
- H.-J. Limbach, A. Arnold, B. Mann and C. Holm, Comput. Phys. Commun., 2006, 174, 704 CrossRef CAS.
- A. Arnold, S. Lenz, O. Kesselheim, R. Weeber, F. Fahrenberger, D. Röhm, P. Košovan and C. Holm, Meshfree methods for partial differential equations VI, Springer, 2013, pp. 1–23 Search PubMed.
- C. Allain, M. Cloitre and M. Wafra, Phys. Rev. Lett., 1995, 74, 1478 CrossRef CAS PubMed.
- C. Grimaldi, J. Chem. Phys., 2017, 147, 074502 CrossRef PubMed.
- M. Miller and D. Frenkel, J. Chem. Phys., 2004, 121, 535 CrossRef CAS.
- S. Zhang, L. Zhang, M. Bouzid, D. Zeb Rocklin, E. Del Gado and X. Mao, arXiv, 2018, 1807.08858, 1.
- S. Hayward, D. Heermann and K. Binder, J. Stat. Phys., 1987, 49, 1053 CrossRef.
- H. Verduin and J. Dhont, J. Colloid Interface Sci., 1995, 172, 425 CrossRef CAS.
- C. de Kruif and J. van Miltenburg, J. Chem. Phys., 1990, 93, 6865 CrossRef CAS.
- M. Grant and W. Russel, Phys. Rev. E: Stat. Phys., Plasmas, Fluids, Relat. Interdiscip. Top., 1993, 47, 2606 CrossRef CAS.
- M. Carpineti and M. Giglio, Phys. Rev. Lett., 1992, 68, 3327 CrossRef CAS PubMed.
- N. Verhaegh, D. Asnaghi, H. Lekkerkerker, M. Giglio and L. Cipelletti, Phys. A, 1997, 242, 104 CrossRef CAS.
- P. Lu, E. Zaccarelli, F. Ciulla, A. Schofield, F. Sciortino and D. Weitz, Nature, 2008, 453, 499 CrossRef CAS PubMed.
- M. Helgeson, Y. Gao, S. Moran, J. Lee, M. Godfrin, A. Tripathi, A. Bose and P. Doyle, Soft Matter, 2014, 10, 3122 RSC.

## Footnotes |

† Conceptualization, J. d. G., W. C. K. P., and M. H.; methodology, J. d. G., W. C. K. P., and M. H.; software, J. d. G. and M. H.; validation, J. d. G. (lead) and M. J. H. (supporting); investigation, J. d. G. (lead) and M. J. H. (supporting); writing – original draft, J. d. G. and M. H.; writing – review & editing, W. C. K. P.; funding acquisition, W. C. K. P.; resources, W. C. K. P.; supervision, J. d. G. and W. C. K. P. |

‡ We obtained the location of the peak using simple sorting and binning, and fitted it to a local 4th order polynomial. Nevertheless, large errors remain due to the relatively large uncertainty in the measured S(q). |

This journal is © The Royal Society of Chemistry 2019 |