Hydrodynamic lubrication in colloidal gels

Abstract

Colloidal gels are elasto-plastic materials composed of an out-of-equilibrium, self-assembled network of micron-sized (solid) particles suspended in a fluid. Recent work has shown that far-field hydrodynamic interactions do not change gel structure, only the rate at which the network forms and ages. However, during gel formation, the interplay between short-ranged attractions leading to gelation and equally short-ranged hydrodynamic lubrication interactions remains poorly understood. Here, we therefore study gelation using a range of hydrodynamic descriptions: from single-body (Brownian Dynamics), to pairwise (Rotne–Prager–Yamakawa), to (non-)lubrication-corrected many-body (Stokesian Dynamics). We confirm the current understanding informed by simulations accurate in the far-field. Yet, we find that accounting for lubrication can strongly impact structure at low colloid volume fraction. Counterintuitively, strongly dissipative lubrication interactions also accelerate the aging of a gel, irrespective of colloid volume fraction. Both elements can be explained by lubrication forces facilitating collective dynamics and therefore phase-separation. Our findings indicate that despite the computational cost, lubricated hydrodynamic modeling with many-body far-field interactions is needed to accurately capture the evolution of the gel structure.

---

1 Introduction

Colloidal particles dispersed in a fluid medium can form a gel, when there are short-ranged attractions between the colloids with an interaction strength that significantly exceeds the thermal energy k_BT;^1–4 here, k_B is the Boltzmann constant and T the temperature. These interactions can be induced by, e.g., the presence of polymers⁵ (depletion) or van-der-Waals interactions^6,7 (solvation). The depth and short-ranged nature of the attractive potential well interfere with thermal rearrangement of clustered colloids. This arrests the system's natural tendency to fully phase separate.^8,9 Diffusion-based aggregation together with this arrested cluster dynamics leads to the formation of an open, space-spanning network structure.^6,10,11 The resulting material possesses useful properties, e.g., it can support the gel's buoyant weight against gravity for a finite time,^12–19 behaving solid-like during this period, yet flow under moderate applied stresses,²⁰ behaving liquid-like. This ability to yield easily and reform²¹ has led to the widespread use of particle gels in industrial, medical, and academic settings, for example, care products, printing inks, foodstuffs, crop protection, and pharmaceutical suspension formulations.^22–25

Gels are intrinsically out of equilibrium and therefore age after they form.^26–28 A gel's structure and its mechanical response are intimately linked.^20,29,30 This necessitates understanding of formation and aging in order to begin to predict and ultimately control the properties of colloidal gels. Both processes have been suggested to depend (strongly) on the nature of the fluid medium, in which the colloids are suspended.^31–36 That is, unlike equilibrium systems, for which the particle dynamics are known not to affect the average arrangement of colloids, hydrodynamic interactions (HIs) can affect gel structure by impacting the kinetics of formation. HIs result from momentum conservation in the fluid, which gives rise to long-ranged forces—decaying with the inverse of the center-to-center separation 1/r—between moving suspended particles.^37,38 The forces inducing HIs can be externally applied, stem from the interaction potential, or are thermal in nature. The latter two play a prominent role in gel formation (without gravity), thus it is relevant to quantify their effect on structure when accounting for the suspending fluid.

Over the past two decades, various research groups have performed numerical studies of colloidal gels that account for HIs.^{18,19,31–36,39–45} Introducing HIs into particle simulations poses a significant challenge. In the friction-dominated regime appropriate to colloidal hydrodynamics, a Newtonian fluid satisfies the Stokes equations. Exactly solving the full set of Stokes equations, for which the particles act as boundary conditions, is practically infeasible for all but the simplest systems.^37,38 Additionally, HIs decay as 1/r, which makes the dynamics in a particle suspension an intrinsically many-body problem.^46,47 A number of approximating techniques have been proposed over the years to make progress.^48–50 These techniques weight the accuracy, by which HIs are approximated, against computational efficiency. Prominent examples of methods that have been used in studying colloidal gels include: Fluid Particle Dynamics (FPD), Lattice-Boltzmann (LB), and Stokesian Dynamics (SD). In view of the approximating nature of gel simulations with HIs, it is not surprising that many different (seemingly conflicting) results on the impact of HIs have been reported.

De Graaf et al.³⁶ recently used lattice-Boltzmann (LB) method^48,51 to account for the many-body, far-field interactions between gelling colloids. These authors, revealed two distinct regimes in their simulations with and without HIs, and a crossover at a colloid volume fraction of ϕ ≈ 0.16. For ϕ < 0.16, HIs were found to accelerate gelation, while for greater ϕ, the effect was seen to be reversed. The apparent substantial differences in the gel structure, when comparing gels with and without HIs at equivalent time, particularly at low ϕ, were found to become negligible, when comparing the systems at equivalent “structural times”.³⁶ The findings of De Graaf et al. appear to reconcile the contradictory observations reported in previous simulation studies with HIs, requiring only limited reference to the various approximations made by other authors.

However, similar to the vast majority of simulation studies, the authors of ref. 36 ignored hydrodynamic lubrication forces (HLFs). These forces arise when colloidal particles are separated by less than ≈10% of their particle diameter and are highly dissipative.⁵² Motion parallel to the line connecting the centers of two spherical colloids is associated with a hydrodynamic friction that diverges as 1/h, where h measures the surface-to-surface separation. This divergence results from the force required to squeeze the fluid out of the small gap between the colloids. Such an interaction has the potential to interfere with gelation, as the attractions induced by depletion interactions also typically fall in this separation regime. Analogous divergences exist for the other modes of relative displacement, e.g., there is a logarithmic divergence for motion orthogonal to this connecting line.

In the context of colloidal gels, the work of Bybee and Higdon⁵³ is often referenced to justify neglecting HLFs. These authors approximated HIs by considering exclusively HLFs between nearly touching particles. Their study revealed that the percolation line and microstructure of the gel did not exhibit noticeable changes, when compared to the results obtained using non-lubricated Brownian Dynamics (BD). Thus, these authors argued that near-field HIs contribute only minimally on the time scale of gelation. However, the recent work by Townsend et al.⁵⁴ on HIs in colloidal suspensions, demonstrated that neglecting far-field HIs can lead to inaccuracies, even when colloidal particles are in close proximity to each other. In light of this study, it is worthwhile to investigate the combined effect of HLFs and far-field HIs on colloidal gelation.

An additional argument used to neglect HLFs is that in many experimental systems it is not clear that the fluid is Newtonian. That is, macromolecules (polymers) may have been added to induce gelation via depletion, which could affect the medium's hydrodynamic response. Moreover, even very well synthesized colloidal spheres are not perfectly spherical or hard.⁵⁵ This has led to debate on the relevance of contact interactions.^56–59 Both of these effects limit the applicability of analytic lubrication expressions, which are derived for perfect spheres in a Newtonian solvent. Nonetheless, it is important to fundamentally understand what effect HLFs can have on gelation, gels structure, and aging, before discounting these in favor of other mechanisms. This leads us to revisit colloidal gelation in this work.

It should be noted that LB—and many other approximations for HIs—possesses a near-field increase of the hydrodynamic friction between approaching particles.³⁶ In the case of embedded particles that couple to the lattice fluid on a sub-grid scale, the increase is non-divergent and a mix between far-field flows and compressibility/approximation-level artifacts. For Ladd-type boundaries, there are also issues, but these can be explicitly lubrication corrected,⁶⁰i.e., the far-field contribution is subtracted and replaced by accurate analytical lubrication approximations. However, it is not clear how to effectively do this when studying the dynamics of suspended sub-grid particles experiencing thermal fluctuations. Lubrication corrected Ladd-type particles typically require too many lattice sites to effectively study gelation, as the number of particles is limited by the computational requirements. A different approach is thus needed.

Here, we thus consider the effect of truncating the HIs at four distinct levels, which allows us to chart the specific effect of each approximation. The four levels of approximation are: (i) Free-draining spheres (FD)—effectively captured by BD simulations—will serve as our reference point. This model considers forces between the fluid and colloids only at the one-body level, i.e., the colloids experience Stokes drag but have no long-ranged or short-ranged HIs. (ii) Rotne–Prager–Yamakawa (RPY) hydrodynamics,^61–65 which approximate HIs using far-field Green's functions at a pair level. That is, RPY HIs ignore the intrinsic many-body effects present between suspended spheres, which is typically considered reasonable for low ϕ. Here, we do not make use of lubrication corrections. (iii) Stokesian Dynamics (SD)^66–71 without lubrication corrections (SD^ff). This approach accounts for the many-body effects, but has limited accuracy for spheres that approach each other closely, as is the case in colloidal gels. Of the approaches, it is most closely comparable to the work by De Graaf et al.³⁶—the LB method is many body in nature and we will contrast ref. 36's results to those obtained here. SD is based on matrix inversion and typically more costly than LB, limiting the number of particles that we studied to about 7000. (iv) Lastly, we consider lubrication-corrected SD (simply SD throughout). This is the most accurate of all techniques in describing the interactions between spheres suspended in a Stokes fluid. However, this comes at the price of being the most computationally expensive. Note that even full SD is approximating in the mid-field, i.e., when transitioning from lubrication to the far-field expressions.

Using the above four methods, we find that the dynamics of the gel are highly sensitive to the type of hydrodynamic model used, in line with previous studies. However, the (quasi-)steady-state structure remains relatively unaffected between these, provided HLFs are not present. Interestingly, when HLFs are introduced, significant differences are observed, particularly for colloid volume fractions ϕ < 0.138. In addition to altering the structure of the gel, resulting in smaller voids and clusters at the percolation point, lubrication also accelerates gel aging. This may seem counterintuitive, as lubrication is associated with the strongest hydrodynamic dissipation. The explanation for both these effects is that lubrication interactions facilitate phase separation, by interfering with bonding and preventing the rupture of clusters once they have formed. Additionally, HLFs suppress non-collective displacement and deformation modes in the gel arms, leading to more vigorous dynamics.

The remainder of this paper is organized as follows. In Section 2, we first provide an overview of the numerical methods used in this work. Section 3 contains specific details of the hydrodynamic models considered, while in Section 4 we present the quantities used to characterize our gel systems. Next, in Section 5, we show the results of our simulations and describe the system's dynamics and structure. Here, we also provide the first insights into the effects of HIs and HLFs. In Section 6, we connect our results with the existing literature and provide an explanation for the deviations observed in SD simulations. We also discuss some of the advantages and limitations of the SD method. Finally, Section 7 concludes with a summary of our findings and an outlook on future directions for particle-based simulations of colloidal gels.

2 Simulation setup

We want to study the influence of hydrodynamic interactions on the structural formation and evolution of a colloidal gel. We do so by performing particle-based simulations of neutrally buoyant spherical particles of diameter σ, suspended in a fluid with viscosity η. We simulate only the colloids and account for the presence of the polymers that cause depletion attraction via a generalized “high-exponent” LJ potentialwhere ε the interaction strength is set to 10k_BT. This is a smooth approximation to the well-known Asakura-Oosawa-Vrij depletion attraction,⁷² when this is combined with a steric repulsion that prevents overlap.

We can write the equations of motion for the entire system of N colloids using three 6N-dimensional generalised force vectors

which group forces and torques together. The first term in the right-hand side of eqn (2), , is a 6N-dimensional vector with the first 3N components containing the inter-particle forces and the last 3N entries are zero (our interaction potentials are torque free). The second term, contains hydrodynamic drag forces. In a steady Stokes' flow, this term can be written as the product of a hydrodynamic resistance matrix, which depends only on the relative distances between particles, and a 6N-dimensional vector containing translational- and rotational-velocities:The last term in eqn (2), contains forces and torques that account for Brownian motion. These are instantaneously correlated through the HIs and their variance is given by the fluctuation–dissipation theorem.^73,74 In a simulation with time-step size Δt, we can express this correlation as:^75,76
where δ_ij represents the Kronecker delta and the angled brackets indicate a time average.

To compute the Brownian contributions, we must specify by which convention we evaluate them.⁷⁷ For numerical studies, a common choice is the Itô convention. At each (discrete) time t_i, we compute R_FU using the relative particle positions r(t_i). Following this convention, the Brownian forces acquire a non-zero average

which is usually called “Brownian drift”. Physically, its presence is required to ensure stationarity under the Gibbs–Boltzmann distribution and generate particle configurations with the correct statistics at equilibrium.⁶⁴

Combining the above, the equation for particle displacement reads

with ψ a 6N-dimensional vector containing uniformly distributed random variable with zero mean and unit variance. We must evaluate R_FU to integrate the colloid trajectories in time. In general, computing an exact hydrodynamic resistance or mobility matrix is a task with a high computational expense, involving the approximate solution of the Stokes equations in the fluid phase surrounding the particles. Section 3 provides details on the four approximations to R_FU, which we used to perform efficient hydrodynamic simulations.

We used periodic, cubic boxes with an edge length L such that the volume fraction of the systems were ϕ ∈ [0.075, 0.138, 0.225]. The number of colloids in the boxes was chosen to be N = 1000 and ≈7000, with L suitably modified to achieve a given ϕ. Unless otherwise explicitly indicated, all simulation results presented use the larger particle number. Following,³⁶ we first equilibrated each configuration for 50τ_B using a purely repulsive inter-particle potential (r_cut = σ), where τ_B = σ²/(4D) is the Brownian time of the colloids with single-particle translational diffusion coefficient D. Subsequently, the gels were formed via an instantaneous deep quench from a purely repulsive potential to one with the aforementioned 10k_BT attraction strength, and left to form and age for [20, 50, 80, 1000]τ_B.

The FD, SD^ff, and SD simulations were performed using HOOMD-blue, a GPU-compatible Python package developed in the Glotzer Lab.⁷⁸ For the SD^ff and SD simulations we used Fiore's external plugin “Positively-split Ewald” (PSEv3) to implement the Stokesian Dynamics algorithm.⁷¹ The RPY simulations were performed using UAMMD, a GPU-accelerated software infrastructure for complex fluids simulations.⁶⁵ For each data point we ran 3 independent simulations (10 for the smaller system sizes). Each typically took several hours to several days to run on a desktop (i9-10900) with a modern graphics processing unit (GPU; NVidia RTX 2060 Super).

3 Hydrodynamic models

In this section, we provide some details on the hydrodynamic approximations used in our work, aimed at contextualizing our results. For each approximation, we provide references to help the interested reader to learn more about that approach. Fig. 1 gives a visual summary of these approaches and introduces a color coding for the methods that we use throughout.

Fig. 1 Schematic representation of the four hydrodynamic models used. The relation between the displacement of the ith particle Δr_i and the total force applied on it F^tot_i and its surrounding particles F^tot_j≠i is indicated at the top; the c_ij represent coupling between the particles. (a) In FD simulations, Δr_i depends only on the force on the particle itself, c_ij = δ_ij. (b) In RPY, particles couple pairwise, i.e., each c_ij is a function of the relative distance r_ij between the ith and jth particle, respectively. (c) In lubrication-free SD, SD^ff, there is coupling between all particles, so that each c_ij has as argument the set of all distance vectors {r_ij}. (d) In full SD, short-range lubrication interactions are incorporated, thus the coefficient c_ij is modified by an additional pair-wise interaction that causes it to vanish, when the surfaces of particles i and j approach each other (spheres with diameter σ).

FD spheres are arguably the simplest model for HIs. Here, HIs are accounted for only at a one-body level. This means that each particle experiences the same drag force (given a velocity), irrespective of the presence of other particles in the simulation volume. In this case, the resistance matrix takes the simple diagonal form R^(ij)_FU = 3πησδ_ij and the thermalization is similarly reduced, see Fig. 1a. This makes the FD approximation the fastest method at our disposal; it is fully equivalent to regular BD. The efficiency comes at the cost of neglecting HIs at any relative distance, which means that this approximation produces unrealistic dynamics for suspended colloids. This is not an issue for studying equilibrium systems, though it strongly affects the dynamics of formation^31–35,41 and collapse^18,45 in colloidal gels.

A more realistic approximation involves a direct construction of the inverse of the resistance matrix, the mobility matrix. This can be done using pairwise Greens' functions for the interactions between spheres, which are called RPY tensors, as depicted in Fig. 1b. At a two-body level, these follow from the combination of Faxén's laws for spherical particles and a multi-pole expansion of the solution to Stokes' equation.⁶¹ For two particles in an unbounded solvent, the tensor takes the form

where I is the 3-dimensional identity matrix, r = |r_i − r_j|, and r̂ is the center-to-center unit vector. It is implied that r̂r̂ represents the Kronecker product of the two vectors and acts as a 3 × 3 matrix. Invoking pair-wise addition, the tensor for a system of N colloids can be written as

In periodic domains, which is the case of our simulations, a compact form of the tensor can be computed in Fourier space, resulting in a positive-definite matrix for any particle configuration.^62–64 The downside of using a RPY-based resistance matrix, is that its pair-wise nature overestimates hydrodynamic forces, when many particles are relatively close. RPY cannot reproduce screening effects at large volume fractions,⁷⁹ which are known to impact the dynamics in dense suspensions. Gels are partially dense and partially open, which raises the question whether RPY is sufficiently accurate. The Swan group^33,35 showed that RPY is able to produce distributions of the particle contact number, that are representative of those observed in experiment. As we will discuss in Section 6, the contact number distribution may not always be a telling quantifier.

The SD method takes into account far-field, many-body HIs by computing the resistance matrix from a grand-mobility matrix via matrix inversion. The latter is constructed similarly to M_RPY, but the number of hydrodynamic moments included in the multi-pole expansion is higher than in the RPY formulation.⁷⁹ As a result, the grand-mobility matrix is an 11N × 11N matrix, and R_FU is equal to the upper left 6N × 6N block of . The operation of matrix inversion gives rise to many-body hydrodynamic interactions⁷⁹ (see Fig. 1c), which would be absent otherwise, since is computed pair-wise, as with the RPY tensor.

The HIs discussed thus far are all far-field approximations, since the multipole expansion is truncated at a finite order. As a consequence, short-range hydrodynamic forces are poorly reproduced. In particular, the divergent part of such interactions—following from lubrication—is completely neglected. In the framework of SD, such interactions are instead computed from the analytical expressions derived by Jeffrey and Onishi⁵² for two-spheres at close proximity. These expressions are used to compute a two-sphere resistance matrix accurate at small separation, which is then used as a building block to assemble the lubrication gran-resistance matrix as illustrated in Fig. 1d. However, part of the two-body resistance interactions have already been included upon the inversion of and therefore, must be subtracted from the total grand-resistance matrix. The matrix containing these interactions is found by simply inverting a two-body mobility matrix. Thus, the total grand-resistance matrix, which contains both near-field lubrication and far-field many-body interactions, is

where the last term contains the aforementioned corrections. R_FU is then equal to the upper-left block of .

It should be noted that here we applied a 1% positive shift in the colloid diameter σ used in the inter-particle potential of eqn (1), compared to the one used to construct the resistance matrix R_FU. We did this to prevent particle–particle overlaps, which would make the numerical scheme unstable. Even though lubrication effects are already dominant at surface-to-surface separations ≈10% of the colloids diameter, the imposed shift might have weakened the contribution of lubrication to both gel dynamics and structure. Reducing this shift should make the effect of lubrication more prominent.

4 System characterization

We used various characterization techniques to understand the dynamics in our colloidal gels and quantify their structure. Here, we predominantly made use of post processing using the Python packages “freud”⁸⁰ and “NetworkX”⁸¹. In terms of identifying structures, we consider colloids to be bonded, when their center-to-center distance Δr ≤ 1.05σ, as in ref. 36.

In order to highlight the effects of HIs, we study the shape of the clusters using the relative shape anisotropy. This is defined as

where the average is performed over the total number of clusters in the system, and λ_i are the eigenvalues of each gyration tensor. Values of κ² close to zero identify spherically symmetric clusters, whereas values approaching unity occur only when all the colloids in the cluster lie into a straight line.

To quantify the degree of phase-separation in the system, we make use of the void volume (VV). Here, the VV is computed by considering the volume of a sphere centered at a point in the system, which is just in contact with the surface of the nearest colloid.^20,36,82 From such a single-point sphere volume, we can compute the average void volume 〈VV〉 via Monte-Carlo integration. The VV allows us to identify structural dissimilarities on the scale of the gel arms, rather than the level of the particles (characterized by the coordination number).

Another quantity capable of indicating the level of coarsening in the gel is the tortuosity parameter ξ. This is defined as

with N_c the number of colloids in the considered cluster, D^E and D^N the Euclidean and network distances respectively. The tortuosity parameter is equal to the relative difference between Euclidean and “network distance” between two particles in the same cluster, averaged over all the pairs. Here, we use “network distance” to mean the length of the shortest path that connect two particles that belong to the same cluster. That is, ξ ∈ [0,1] indicates how tortuous is the gel network, with low values of ξ representing systems close to completing the phase separation, and values close to one for networks with a non-trivial topology and thus, far from completing the phase-separation process. As we will show in Section 5, the large differences between the hydrodynamic models studied manifested in ξ rather than VV, with the latter showing smaller deviations.

5 Results

In this section, we provide the main results of our study. We start by showing the analogy to the LB-based work by De Graaf³⁶ to demonstrate the various methods provide similar far-field results. Here, we focus mostly on the cluster formation and percolation. Next, we turn to the aging of the system and the way it phase separates. Finally, we consider the state-space trajectories and demonstrate that lubricated interactions lead to a divergent structural pathway.

5.1 Gel network formation

We first focus on the short-time dynamics of our systems. Fig. 2 provides the fraction of particles n^max_c = N^max_c/N that belong to the largest cluster. We observe that the overall trends are comparable to those of De Graaf et al.³⁶ (magenta data). For the lowest ϕ, simulations performed with FD exhibit a slower growth than any other hydrodynamic model, while for larger ϕ, FD predicts the fastest cluster growth; a crossover regime is observed for ϕ = 0.138. For completeness, we indicate the power-law scaling that was obtained by De Graaf et al., which appears to match the trends that we find using our approaches. Note that particle incorporation in the largest cluster is substantially faster, for our lubricated SD simulations at low ϕ. We will come back to this result in the following.

Fig. 2 The fraction of particles in the largest cluster n^max_c (compared to the total number N) as a function of time t (in Brownian time units τ_B). We show 1 − n^max_c to highlight the long-time, power-law scaling for three different colloid volume fractions ϕ, labelled from small to large using circles, triangles, and squares, respectively. For completeness, we also show the LB results from de Graaf et al.³⁶ The hydrodynamic approximations are indicated using coloring, as indicated in the legend. Large open symbols indicate the point for which the system percolates. Error bars provide the standard error of the mean. Black dotted lines and numbers indicate the exponent of asymptotic power laws from ref. 36.

At these low values of ϕ, the LB result of ref. 36 matches SD^ff best, as expected, because HIs play a significant role, but LB does not incorporate lubrication corrections. Interestingly, at higher ϕ, the LB result best agrees with the FD results. This is presumably, because the near-field aspects that are differently captured by SD, SD^ff, and even RPY, which is relevant already at intermediate ϕ. Nonetheless, the trends and power-law scaling match, with the ϕ = 0.138 result switching between matching on FD or a more involved hydrodynamic model. Lastly, it should be noted that there is a substantial difference in particle numbers between our simulations and the LB results of ref. 36, meaning that a direct comparison in terms of n^max_c needs to be treated with care.

We also determined the onset of percolation as a function of ϕ. Fig. 3 shows the percolation time t_p, defined as the moment at which there exist a cluster that percolates at least one spatial dimension. In line with ref. 36, we find that in all cases, HIs speed up the percolation process at low ϕ, and slow it down in denser systems. Note that the dramatic effect of lubrication interactions at low ϕ is borne out by the percolation transition as well. Contrasting SD and SD^ff allows us to assess the role of squeeze flows. Our results show that these play a role, but that the effect at intermediate to high volume fractions (for colloidal gelation) is limited. Interestingly, for all regimes, there is limited distinction between RPY and SD^ff dynamics. This suggests that many-body hydrodynamic effects are of limited importance to gelation.

Fig. 3 The percolation time t_p (in terms of τ_B) as a function of ϕ for all of our hydrodynamic methods (legend). The dashed lines connect the points belonging to specific method's data set and serve as guides to the eye. The standard error of the mean is indicated using the error bars, but is typically smaller than the symbol size.

However, following the largest cluster may not reveal all features of the evolution of the system. We discriminate the effect of various hydrodynamic contributions in the formation of a space-spanning network, by analyzing the shape and size of the clusters that are present in the sample at a given time. Fig. 4a reports the mean relative shape anisotropy, 〈κ²〉. It is clear that this generally decreases as clusters form. This is understood as follows: dimers have κ² = 1 and, as the clusters grow, the branching structures become more isotropic, leading to a decrease of the overall value of 〈κ²〉.

Fig. 4 Features of the short-time clustering in the quenched suspension. (a) The mean relative shape anisotropy 〈κ²〉 as a function of time t (in terms of the percolation time t_p; emphasized by use of the vertical dashed line). (b) The relative mean cluster size 〈n_c〉 = 〈N_c〉/N also as functions of t/t_p. The use of colors, symbols, and error bars is that of Fig. 2.

At short times, newly formed clusters have heterogeneous shapes in any hydrodynamic model considered, and, as they diffuse, these clusters will also reconfigure into more compact shapes in order to minimize their energy. The characteristic time related to the latter process is in competition with percolation, as compact clusters form space-spanning network in a less efficient way compared to more elongated ones. Systems with ϕ ≳ 0.138 show a clear minimum in 〈κ²〉 at the percolation transition. At larger times, the growth of the mean shape anisotropy could be related to the network compaction, as pointed out by Tsurusawa et al.³⁰ in their experimental study. However, the values of κ² are certainly influenced by the periodic boundaries, once a cluster forms that spans the entire simulation box. The long-time increase in 〈κ²〉 could be an artifact of this. We deem it irrelevant to our study, as we focus on short-time dynamics (t ≤ t_p) in analyzing the relative shape anisotropy.

Fig. 4b shows the mean number of clustered particles 〈n_c〉 = 〈N_c〉/N as a function of time. It is clear that for intermediate to high ϕ, the early time cluster dynamics is difficult to distinguish between the various approaches. Systematically, the values for the FD simulations appear lower, but recall that a fair comparison requires an analysis in structure space.³⁶ Interestingly, the situation is distinct for low ϕ, as all curves are initially different, though the FD, RPY, and SD^ff results start to track each other beyond t = t_p. The behavior of our SD simulations always remains separate from these curves, suggesting that lubrication has the ability to modify the structure as well as the dynamics.

In general, the effect of HIs is to both decrease the bonding rate of colloids and increase the time needed to rearrange clusters into lower energy configurations. Thus, on average, HIs impede cluster growth (see Fig. 4b). For large volume fractions, the time it takes to form a large percolated cluster, which depends on both ϕ and colloids bonding rate, is small compared to the time required by clusters to relax their shape. This becomes clear from Fig. 4a, where large ϕ curves overlap once scaled by their respective percolation times. Thus, in that regime of ϕ, the only noticeable effect of HIs on the gel formation is to decrease the bonding rate of colloids, thereby slowing down short-time gel dynamics.

In contrast, for the the lowest ϕ, the characteristic time associated with cluster reshaping is comparable with t_p, as shown in Fig. 4a. In this figure, FD systems reach a plateau in the mean relative shape anisotropy κ² before percolation occurs, while the other hydrodynamic models maintain on average more elongated shapes. This promotes the emergence of large structures composed by strand-like clusters of colloids, which speeds up space exploration,⁸³ the capturing of smaller clusters, and ultimately percolation.^31,32,34 On the other hand, clusters of particles can rapidly rearrange into more compact configuration in FD simulations. Thus, disfavouring the formation of a percolated network-like structure.

5.2 Aging and phase-separation

Next, we shift our focus to the long-time dynamics in the gelled system. Fig. 5 shows the time evolution of the tortuosity ξ for our three ϕ of interest. Note that the curves peak approximately at the respective percolation points. Any mismatch is likely explained by our limited system size and choice of defining percolation as spanning in one direction. After the systems have percolated, the resulting network-like structures will keep spinodally decomposing into colloid-rich and colloid-poor phases, though in an arrested manner. Systems with HIs exhibit a markedly faster decay in ξ for all values of ϕ examined, especially when the SD approximation is used. All of the intermediate- to high-ϕ systems show a plateauing of ξ, suggesting that these display strongly arrested dynamics. Low-ϕ systems appear to remain more dynamical over the range that we could simulate. We will argue that this effect is given by the long-range hydrodynamic interactions that propagates among the gel arms.

Fig. 5 The tortuosity ξ as a function of time t (in τ_B). The curves for the various values of ϕ are shifted (as labelled), in order to better visualize the data. The use of colors, symbols, and error bars is the same as in Fig. 2.

We also studied the mean squared displacement 〈r²〉 (MSD), averaged over all colloids in the simulation volume. Fig. 6 shows that, at short times, all MSDs display a diffusive (linear) regime, as expected for systems that have not yet significantly clustered. In the long-time limit, all systems become subdiffusive, which is representative of the caging colloids experience in the network structure.⁸⁴ For our FD system, we were able to simulate sufficiently long to characterize the trend, which appears to follow a power law: 〈r²〉 ∝ t^1/4.

Fig. 6 The mean squared displacement (MSD) 〈r²〉 averaged over all colloids, expressed as a function of time t (in τ_B). The dashed lines indicate various power-law regimes, while the large symbols indicate the percolation point. The use of colors, symbols, and error bars is the same as in Fig. 2. However, here, the dashed lines represent fitted power-law trends for the MSD.

Surprisingly, for the lowest ϕ, SD-approximated HIs give rise to a transient superdiffusive regime (〈r²〉 ∝ t^1.35), right after the system percolates. This can be related to the smaller mean cluster size 〈n_c〉 that SD has before and at percolation, see Fig. 4b.

At percolation, the size of the largest cluster is roughly the same in every hydrodynamic model, meaning that in SD simulations the percolated network is surrounded by smaller “free” clusters on average. Thermal fluctuations cause displacements of the network arms, which in turn propagate long-range hydrodynamic forces onto the remaining clusters (except in FD systems). Smaller clusters will experience larger displacements due to these interactions, increasing on average the MSD in the SD systems. This is what we alluded to at the end of the opening paragraph to this section. Lastly, all long-time MSDs with HIs are larger than those of the FD systems. The data is suggestive of a collective speeding up of the aging dynamics, which we again relate to HIs generated by the arms of the colloidal gel.

In order to facilitate future comparison of our simulations to experiment, we consider scattering parameters. Fig. 7 shows the peak positions q_m of the static structure factors, multiplied by the colloid radius, as a function of time. For all values of ϕ, SD systems systematically have a stronger decrease in q_m as a function of time. This signals accelerated aging, as the effective length scale of the gel grows more rapidly. Note that in this choice of representation, comparing the absolute values of q_m at equal times is not useful, as the systems do not share a common origin point. However, RPY and SD^ff data have a striking similarity, as expected.

Fig. 7 The peak position of the dimensionless structure factor q_m as a function of time, divided by the percolation time t_p. The two higher density curves have been shifted up by the amount indicated in parentheses to prevent overlaps. The use of colors, symbols, and error bars is the same as in Fig. 2.

5.3 State-space trajectories

In the previous section, we have shown that HIs strongly affect the dynamics of the system, which was expected, barring the unexpectedly large impact of HLFs at low ϕ. Here, we investigate whether this strong change in dynamics is indicative of an actual change in gel structure. Fig. 8a shows the systems' trajectory in a “state space” formed by the average void volume 〈VV〉 and the average number of nearest neighbors 〈z〉. In this representation, time is parametric and its direction is indicated by the arrow. The results are similar to those reported by De Graaf et al.³⁶ All hydrodynamic models roughly follow the same trajectory in state space, except systems that include hydrodynamic lubrication. The SD data show deviations that become larger as the ϕ is decreased.

Fig. 8 State-space plots for the structural evolution of colloidal gels. The spaces are spanned by the average void volume 〈VV〉 and (a) average number of nearest neighbors 〈z〉 and (b) the tortuosity of the largest cluster ξ^max, respectively. The paths are followed in the direction indicated by the black arrow. The use of colors, symbols, and error bars is the same as in Fig. 2.

We highlight the differences displayed by lubricated systems in Fig. 8b, which shows the gel trajectories in a state space spanned by 〈VV〉 and the tortuosity of the largest cluster ξ^max. The qualitative picture remain roughly unchanged for large and intermediate ϕ, with deviations from the simple FD predictions becoming more significant as the hydrodynamic model considered becomes more accurate. For the lowest ϕ considered here, SD systems pursue a clearly distinct path in state space from those obtained using FD, RPY, and SD^ff. We conclude that tortuosity is a better characterizer of differences in colloidal gels than void volumes, and that the HLFs can strongly impact structure at low ϕ.

6 Discussion

We showed that far-field HIs change the dynamics of colloidal gel formation, as already found in other numerical studies of gelation.^31–36 However, unless HLFs are included, this has little impact on the structural paths that these systems explore as they form, in line with De Graaf et al.³⁶ Nonetheless, HLFs can strongly impact the dynamics and structure of a gel, in contradiction to the results of Bybee and Higdon.⁵³ The important difference here is that we include far-field effects. The strong effect of HLFs on structure are thus likely a consequence of an interplay with far-field HIs. Speculating, it could be that the fluid particle dynamics simulations by Furukawa and Tanaka³¹ picked up on this effect, as these authors used relatively low volume fractions, weak gels, and numerically refined particles.

Let us now explain the counter intuitive result obtained by combining HLFs and far-field HIs during aging. The effect of far-field-only HIs on the aging process is clear: a less static network will coarsen faster over time. Adding HFLs to this, accelerates aging further, for two reasons: (i) divergent dissipative forces hinder both the creation and rupture of colloidal bonds. That is, it is harder to push the colloids in and out of contact, meaning that once arms have formed, they are less likely to break.⁸⁵ (ii) Thermal fluctuations are more collective in nature, when including HLFs and far-field hydrodynamic interactions.

To understand observation (ii), we note that Brownian velocities can be decomposed into a near- and a far-field contribution:⁷¹. The former are fundamentally different, when lubrication is included. The effect of lubrication is to filter out relative thermal displacements and rotations of colloids that are in close proximity. Thus, when there are long-ranged interactions between the colloids, the spectrum of fluctuations is mainly populated by collective modes, which displace coherently the colloids that belong to the same gel arm. Ultimately, this facilitates the coarsening process, i.e., the flux of particles migrating from the colloid-rich to the colloid-poor regions is negatively affected by HLFs, thus, favoring phase-separation. The local compaction decreases the tortuosity parameter, as there are fewer open paths. Simultaneously, the gel network becomes defined by a larger length scale—those loops that have yet to close—which increases the associated length scale. This equivalently expresses itself in a lowering of the q_m value, at which the peak in the static structure factor is located, see Fig. 8. It would be worthwhile to revisit the work on the hopping dynamics of colloids in a gel arm et al.^86,87 in view of this collective effect.

To further understand (ii), it should also be noted that, in absence of far-field HIs, the colloid displacements are not coupled over long distances. The effect of lubrication interactions is still to narrow the spectrum of local fluctuations, but because each colloid behaves independently in terms of its thermal motion, the collective effects we observe in SD would be strongly suppressed. That is in line with the results by Bybee and Higdon.⁵³ This suggests that it might be worthwhile to combine RPY and lubrication, in order to have an efficient algorithm for simulating both lubrication and an effective far-field hydrodynamics. However, RPY and lubrication are constructed in the mobility and resistance picture, respectively. Such a combination would therefore require a highly inefficient double matrix inversion in the current implementation.

In terms of computational cost, SD simulations are clearly the most expensive. This cost stems mainly from the matrix inversion needed to reproduce many-body far-field HIs, rather than the inclusion of near-field lubrication. Fortunately, HLFs slow down bonding dynamics, allowing the use of larger Δt, which makes SD effectively more efficient than SD^ff. Nonetheless, full SD remains a costly numerical scheme in systems with strong short-ranged interactions such as colloidal gels, practically limiting the system sizes that can be studied to a few thousand particles.

7 Conclusions and outlook

Summarizing, we have studied the effect of hydrodynamic interactions in the formation and aging of colloidal gels. We considered single-body (Brownian Dynamics), pair-wise (Rotne–Prager–Yamakawa), many-body (Stokesian Dynamics), and lubrication-corrected many-body interactions. As expected, we found that the dynamics of gel formation is sensitive to the exact nature of the approximation that is used. However, the steady-state structure is relatively unaffected when considering far-field hydrodynamic interactions and appropriate rescaling of time, in line with a previous analysis based on lattice Boltzmann simulations.

Intriguingly, introducing hydrodynamic lubrication corrections results in significant departures from the rescaled trend, especially for colloid volume fractions ϕ ≲ 0.138. Not only is the structure altered, with on average smaller voids and clusters at the percolation point, lubrication also accelerates the aging of the gel. This result seems counterintuitive, as the lubricated regime is typically associated with the strongest hydrodynamic dissipation, and is not in line with previous understanding in the literature. Both the dynamic and structure effect can, however, be explained by the fact that lubrication interactions hinder the bonding and rupturing of clusters, as well as suppress non-collective Brownian modes in the gel arms. These aspects combine to facilitate phase-separation. The key point is that both far-field hydrodynamics and lubrication forces must be present to realize this enhanced separation, which is where our study improves upon the existing literature. Studying colloidal gelation thus requires a relatively costly, but accurate representation of hydrodynamic interactions.

Moving away from the idealized Stokes-flow (approximate) solutions that we considered here, our findings strongly suggests that near-contact dynamics should be given due consideration in the modeling of experimental colloidal gels. This work provides a solid foundation for future studies that make the comparison between experiment and simulation. In particular, it will be relevant to see how the lubricated dynamics impacts gels subjected to externally imposed stimuli, such as shear.

Data availability

Open data package containing the means to reproduce the results of the simulations available at https://doi.org/10.24416/UU01-5SDJAI.

Conflicts of interest

The authors declare that there are no conflicts of interest.

Acknowledgements

The authors acknowledge NWO for funding through OCENW.KLEIN.354. We are grateful to the late Prof. James Swan for initial discussions on the possible differences between LB and SD simulations of colloidal gels. We thank Dr Andrew Fiore and Dr Madhu Majji for their help with getting PSE up and running with HOOMD-blue, as well as Dr Gwynn J. Elfring and Dr Zhouyang Ge for sharing an accelerated variant of PSE and discussing minor bugs in the resistance tensor.

References

1. H. N. Lekkerkerker, W.-K. Poon, P. N. Pusey, A. Stroobants and P. Warren, Europhys. Lett., 1992, 20, 559.
2. W. Poon, J. Phys.: Condens. Matter, 2002, 14, R859.
3. J. Bergenholtz, W. C. Poon and M. Fuchs, Langmuir, 2003, 19, 4493–4503.
4. Y.-L. Chen and K. S. Schweizer, J. Chem. Phys., 2004, 120, 7212–7222.
5. S. Asakura and F. Oosawa, J. Polym. Sci., 1958, 33, 183–192.
6. C. P. Royall, M. A. Faers, S. L. Fussell and J. E. Hallett, J. Phys.: Condens. Matter, 2021, 33, 453002.
7. P. Hartley and G. Parfitt, Langmuir, 1985, 1, 651–657.
8. G. Foffi, K. A. Dawson, S. V. Buldyrev, F. Sciortino, E. Zaccarelli and P. Tartaglia, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2002, 65, 050802.
9. E. Zaccarelli and W. C. Poon, Proc. Natl. Acad. Sci. U. S. A., 2009, 106, 15203–15208.
10. M. Carpineti and M. Giglio, Phys. Rev. Lett., 1992, 68, 3327.
11. E. Zaccarelli, J. Phys.: Condens. Matter, 2007, 19, 323101.
12. R. Buscall and L. R. White, J. Chem. Soc., Faraday Trans. 1, 1987, 83, 873–891.
13. C. Allain, M. Cloitre and M. Wafra, Phys. Rev. Lett., 1995, 74, 1478.
14. D. Senis, L. Gorre-Talini and C. Allain, Eur. Phys. J. E: Soft Matter Biol. Phys., 2001, 4, 59–68.
15. L. Starrs, W. Poon, D. Hibberd and M. Robins, J. Phys.: Condens. Matter, 2002, 14, 2485.
16. S. Manley, J. Skotheim, L. Mahadevan and D. A. Weitz, Phys. Rev. Lett., 2005, 94, 218302.
17. P. Bartlett, L. J. Teece and M. A. Faers, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2012, 85, 021404.
18. R. Harich, T. Blythe, M. Hermes, E. Zaccarelli, A. Sederman, L. F. Gladden and W. C. Poon, Soft Matter, 2016, 12, 4300–4308.
19. P. Padmanabhan and R. Zia, Soft Matter, 2018, 14, 3265–3287.
20. N. Koumakis, E. Moghimi, R. Besseling, W. C. Poon, J. F. Brady and G. Petekidis, Soft Matter, 2015, 11, 4640–4648.
21. P. Smith, G. Petekidis, S. Egelhaaf and W. Poon, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2007, 76, 041402.
22. R. G. Larson, The structure and rheology of complex fluids, Oxford University Press, New York, 1999, vol. 150.
23. A. Darras, A. K. Dasanna, T. John, G. Gompper, L. Kaestner, D. A. Fedosov and C. Wagner, Phys. Rev. Lett., 2022, 128, 088101.
24. R. Mezzenga, P. Schurtenberger, A. Burbidge and M. Michel, Nat. Mater., 2005, 4, 729–740.
25. M. A. Faers, T. H. Choudhury, B. Lau, K. McAllister and P. F. Luckham, Colloids Surf., A, 2006, 288, 170–179.
26. N. A. Verhaegh, D. Asnaghi and H. N. Lekkerkerker, Phys. A, 1999, 264, 64–74.
27. L. Cipelletti, S. Manley, R. Ball and D. Weitz, Phys. Rev. Lett., 2000, 84, 2275.
28. S. Manley, B. Davidovitch, N. R. Davies, L. Cipelletti, A. Bailey, R. J. Christianson, U. Gasser, V. Prasad, P. Segre and M. Doherty, et al. , Phys. Rev. Lett., 2005, 95, 048302.
29. S. Shah, Y. Chen, S. Ramakrishnan, K. Schweizer and C. Zukoski, J. Phys.: Condens. Matter, 2003, 15, 4751.
30. H. Tsurusawa, M. Leocmach, J. Russo and H. Tanaka, Sci. Adv., 2019, 5, eaav6090.
31. A. Furukawa and H. Tanaka, Phys. Rev. Lett., 2010, 104, 245702.
32. J. K. Whitmer and E. Luijten, J. Phys. Chem. B, 2011, 115, 7294–7300.
33. Z. Varga, G. Wang and J. Swan, Soft Matter, 2015, 11, 9009–9019.
34. C. P. Royall, J. Eggers, A. Furukawa and H. Tanaka, Phys. Rev. Lett., 2015, 114, 258302.
35. Z. Varga and J. Swan, Soft Matter, 2016, 12, 7670–7681.
36. J. de Graaf, W. Poon, M. Haughey and M. Hermes, Soft Matter, 2019, 15, 10.
37. J. Happel, H. Brenner and L. R. N. Hydrodynamics, Low reynolds number hydrodynamics. With special applications to particulate media, M. Nijhoff, Leiden, 1983.
38. S. Kim and S. J. Karrila, Microhydrodynamics: principles and selected applications, Courier Corporation, 2013.
39. H. Tanaka and T. Araki, Phys. Rev. Lett., 2000, 85, 1338.
40. R. Yamamoto, K. Kim, Y. Nakayama, K. Miyazaki and D. R. Reichman, J. Phys. Soc. Jpn., 2008, 77, 084804.
41. Z. Varga and J. W. Swan, Phys. Rev. E, 2018, 97, 012608.
42. Z. Varga, J. L. Hofmann and J. W. Swan, J. Fluid Mech., 2018, 856, 1014–1044.
43. K. A. Whitaker, Z. Varga, L. C. Hsiao, M. J. Solomon, J. W. Swan and E. M. Furst, Nat. Commun., 2019, 10, 1–8.
44. L. Turetta and M. Lattuada, Soft Matter, 2022, 18, 1715–1730.
45. J. de Graaf, K. W. Torre, W. C. Poon and M. Hermes, Phys. Rev. E, 2023, 107, 034608.
46. H. Brinkman, Flow, Turbul. Combust., 1949, 1, 27–34.
47. I. Howells, J. Fluid Mech., 1974, 64, 449–476.
48. B. Dünweg and A. J. Ladd, Advanced computer simulation approaches for soft matter sciences III, 2009, pp. 89–166.
49. E. Dickinson, Adv. Colloid Interface Sci., 2013, 199, 114–127.
50. D. S. Bolintineanu, G. S. Grest, J. B. Lechman, F. Pierce, S. J. Plimpton and P. R. Schunk, Comput. Particle Mech., 2014, 1, 321–356.
51. D. Röhm and A. Arnold, Eur. Phys. J.: Spec. Top., 2012, 210, 89–100.
52. D. Jeffrey and Y. Onishi, J. Fluid Mech., 1984, 139, 261–290.
53. M. D. Bybee, Hydrodynamic simulations of colloidal gels: Microstructure, dynamics, and rheology, University of Illinois at Urbana-Champaign, 2009.
54. A. K. Townsend and H. J. Wilson, Phys. Fluids, 2018, 30, 077103.
55. C. P. Royall, W. C. Poon and E. R. Weeks, Soft Matter, 2013, 9, 17–27.
56. J. N. Immink, J. E. Maris, P. Schurtenberger and J. Stenhammar, Langmuir, 2019, 36, 419–425.
57. G. Wang and J. W. Swan, Soft Matter, 2019, 15, 5094–5108.
58. H. T. Nguyen, A. L. Graham, P. H. Koenig and L. D. Gelb, Soft Matter, 2020, 16, 256–269.
59. B. van der Meer, T. Yanagishima and R. Dullens, arXiv, 2022, preprint arXiv:2209.12703.
60. N.-Q. Nguyen and A. Ladd, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2002, 66, 046708.
61. J. Rotne and S. Prager, J. Chem. Phys., 1969, 50, 4831–4837.
62. J. W. Swan and G. Wang, Phys. Fluids, 2016, 28, 011902.
63. A. M. Fiore, F. Balboa Usabiaga, A. Donev and J. W. Swan, J. Chem. Phys., 2017, 146, 124116.
64. A. M. Fiore and J. W. Swan, J. Chem. Phys., 2018, 148, 044114.
65. R. P. Peláez, PhD thesis, Universidad Autónoma de Madrid, 2022.
66. J. F. Brady and G. Bossis, Annu. Rev. Fluid Mech., 1988, 20, 111–157.
67. T. N. Phung, J. F. Brady and G. Bossis, J. Fluid Mech., 1996, 313, 181–207.
68. A. Sierou and J. F. Brady, J. Fluid Mech., 2001, 448, 115–146.
69. A. J. Banchio and J. F. Brady, J. Chem. Phys., 2003, 118, 10323–10332.
70. M. Wang and J. F. Brady, J. Comput. Phys., 2016, 306, 443–477.
71. A. M. Fiore and J. W. Swan, J. Fluid Mech., 2019, 878, 544–597.
72. K. Miyazaki, K. Schweizer, D. Thirumalai, R. Tuinier and E. Zaccarelli, The Asakura-Oosawa theory: Entropic forces in physics, biology, and soft matter, 2022.
73. R. Kubo, Rep. Prog. Phys., 1966, 29, 255–284.
74. J. M. Deutch and I. Oppenheim, J. Chem. Phys., 1971, 54, 3547–3555.
75. D. L. Ermak and J. A. McCammon, J. Chem. Phys., 1978, 69, 1352–1360.
76. G. Bossis and J. F. Brady, J. Chem. Phys., 1987, 87, 5437–5448.
77. N. G. Van Kampen, J. Stat. Phys., 1981, 24, 175–187.
78. J. A. Anderson, J. Glaser and S. C. Glotzer, Comput. Mater. Sci., 2020, 173, 109363.
79. L. Durlofsky, J. F. Brady and G. Bossis, J. Fluid Mech., 1987, 180, 21–49.
80. V. Ramasubramani, B. D. Dice, E. S. Harper, M. P. Spellings, J. A. Anderson and S. C. Glotzer, Comput. Phys. Commun., 2020, 254, 107275.
81. A. A. Hagberg, D. A. Schult and P. J. Swart, Proceedings of the 7th Python in Science Conference, Pasadena, CA USA, 2008, pp. 11–15.
82. M. Haw, Soft Matter, 2006, 2, 950–956.
83. P. K. Pattnayak, A. Kumar and G. Tomar, arXiv, 2023, preprint arXiv:2307.04363.
84. Q. Wu, R. Higler, T. E. Kodger and J. van der Gucht, ACS Appl. Mater. Interfaces, 2020, 12, 42041–42047.
85. C. Ness and A. Zaccone, Ind. Eng. Chem. Res., 2017, 56, 3726–3732.
86. J. M. van Doorn, J. E. Verweij, J. Sprakel and J. van der Gucht, Phys. Rev. Lett., 2018, 120, 208005.
87. J. E. Verweij, F. A. Leermakers, J. Sprakel and J. Van Der Gucht, Soft Matter, 2019, 15, 6447–6454.

---