Topological origins of yielding in short-ranged weakly attractive colloidal gels

Colloidal gels are of great interest in many technological applications^1,2 due to their wide range of mechanical and rheological properties. Depending on the solid volume fraction, and the strength and range of interactions between the particles, attractive colloids can self-assemble into a system-spanning network structure. It is well understood that this network is primarily responsible for the macroscopic mechanical properties of a moderately concentrated colloidal gel.^3,4 A hallmark of colloidal gels is the emergence of yield stress below which the system behaves similar to a viscoelastic solid.⁵ When exposed to external deformation, colloidal gels at intermediate volume fraction go through various changes in their microstructure including local rearrangement (colloid–colloid) bond rupture, cage breakage, and structural anisotropy, which are referred to as gel yielding.^6–16 At the same time, the reformation of bonds between colloids gives rise to a complex dynamic that is highly time- and rate-dependent.¹⁷ Therefore, understanding the relationship between the changes within the microstructure of the particulate network and its yielding behavior is essential for controlling the overall mechanics of colloidal gels.

Many experimental^6,7,9 and computational^11–17 studies have been devoted to understanding the connection between microstructural changes and the yielding of colloidal gels. The results point to some possible origins of the yielding behavior, including structural anisotropy, stress localization, and network rupture. Scattering experiments on model 2D colloidal gels under flow showed structural anisotropy characterized as “butterfly” scattering patterns.^6,7,9,18 A subsequent study using modified Brownian dynamics simulation also predicted structural anisotropy during the yielding of colloidal gels under strong flow.¹⁴ These anisotropies during the yielding can be also quantified by the construction of a fabric tensor that represents local rearrangements of the microstructure.¹⁷ On the other hand, Colombo and Del Gado¹² using a series of detailed simulations found that yielding is accompanied by strong localization of stress and subsequent breakage of a fraction of network connections. This results in avalanche-like fluidization of the colloidal network, which involves processes similar to the relaxation of a rubber band upon snapping.¹⁹ In experiments, two distinct yielding processes were observed in oscillatory shear of strong gels with intermediate solid volume fractions, which were thought to arise from the rupture of inter-particle bonds at small strain and from cluster distortion and breakage at large strain.^11,16 A further simulation study also predicted two-step yielding processes and suggested the inclusion of multi-body hydrodynamic interactions to predict correct dynamics.¹³ Nevertheless, the question of which bonds/particles are primarily prone to these breakage events and ultimately responsible for the yielding of a gel remains unanswered.

A recent study found that at long times and in the flowing conditions, i.e., quasi steady state conditions, it is the lifetime of the colloidal bonds that controls the rate-dependent rheology of colloidal gels.²⁰ This is similar to previous findings for colloidal glasses.²¹ However, for the case of yielding in colloidal gels, the lifetime of a bond directly correlates with whether the colloids that share the bond are within the bulk of a structure or on the surface of the particulate network, and thus is not a reasonable proxy to whether a bond breaks under applied deformation. This behavior is reminiscent of microstructural coarsening of reversible colloidal gel whose interior is glassy and immobile, and mobile particles at the surface drive the coarsening.^22,23 One could argue that these so-called buried particles are unlikely to be primarily responsible for the yielding of the gel upon applying a deformation.

It is clear that colloidal bonds that are more susceptible to the accumulation of stresses/strains play a major role in the yielding of a gel and that anisotropies appear in the microstructure upon application of deformation. In this work, we seek to understand the bond characteristics that correlate with these characteristics and provide a topological origin to the yielding mechanism in colloidal gels. To do so, we employed a series of large-scale dynamic simulations of a reversible colloidal gel at an intermediate volume fraction of ϕ = 0.2 during the start-up of an applied deformation rate and study the bond dynamics during its yielding. Here, we focus on depletion gels that generally result in coarse porous structures as a typical case-study. These structures are known to be different from fractal-like structures that form at much smaller volume fractions of solid. Since the ultimate goal is to be able to predict failure points within the particulate network, it is essential to monitor the evolution of the bond dynamics with respect to their original properties.

In this study, all simulations were performed using a Core-Modified Dissipative Particle Dynamics (CM-DPD) method described in more detail by Boromand et al.^13,24–26 In this CM-DPD model, the motion of a particle particles is governed by the following equation:

where m_i and v_i are the mass and velocity of the particles. Solvent–solvent and colloid–solvent interactions occur via the sum of pairwise conservative force $FijC$⁠, dissipative force $FijD$⁠, and random force $FijR$⁠. For colloid–colloid interactions, conservative forces are excluded and three additional forces are included. One of the hallmarks of the DPD method is the fact that, because solvent particles are explicitly accounted for, hydrodynamics of the system are intrinsically preserved. However, when the distance between two colloids is smaller than the size of a single fluid particle, this hydrodynamic at short-range breaks down and needs to be revisited.^25,27 In our formalism, $FijH$ accounts for the short-range lubrication interaction between colloidal particles based on a pair drag term and the squeezing mode lubrication between two interacting spherical particles. $FijCore$ is the semi-hard potential interaction used to prevent any non-physical overlap of colloid particles during simulation. $FijAtt$ is the attractive short-range interaction to induce gelation that was modeled using a Morse potential,

where h_ij is the surface-to-surface distance between colloids and U₀ and κ⁻¹ are the strength and range of the attractive interactions, respectively. In the present study, a system of colloidal particles at ϕ = 20% volume fraction in a cubic box with periodic boundary conditions was modeled. The simulation box consisted of n_C = 10 000 monodisperse colloids of radius a = 1 and n_S = 502 654 solvent particles at a dimensionless kinetic temperature of k_BT = 0.1. U₀ and κ were set to 6k_BT and 30/a, respectively. Note that here we focused on a very weak attraction strength between the particles, but the observations are expected to hold for stronger interactions as well. The range of attraction corresponds to 0.1a to ensure that the density of interaction across the range can be ignored.

The computational studies were divided into three parts. First, an equilibrium configuration was generated by randomly placing the colloids and solvent particles in the simulation cell, followed by running the simulation for 1 × 10⁵ iterations with no attractive interactions between the colloids, to eliminate any non-physical overlap of the colloidal particles under quiescent conditions. Second, the colloidal gel structure was generated by running the simulation for another 3 × 10⁶ iterations with attractive interactions between colloids under quiescent conditions, which corresponds to 320τ_D of the gelation time where τ_D = πηa³/k_BT is the diffusion time of a bare colloidal particle [Fig. 1(a)]. This is to ensure that a true quasi steady structure is obtained and that no significant structural evolution is observed at longer times. This is also monitored through ensemble-averaged coordination number (the average number of neighbors per particle in the system), as well as the mean squared displacement of the colloidal particles in the arrested state [Figs. 1(b) and 1(c)].

In the final part of the study, a simple shear flow was imposed using the Lees–Edwards boundary condition²⁸ to inquire about its yielding behavior. The integration time step Δt was chosen in the range of 10⁻³–10⁻⁵ depending on the value of the shear deformation rate $γ̇$ to ensure thermodynamic and numerical stability throughout all simulations. All simulations were performed using HOOMD-blue,²⁹ the open source molecular dynamics simulation toolkit.

In order to understand the topological origin of yielding in colloidal gels, we monitor the dynamics of bond formation and rupture throughout the flow protocol for a range of shear deformation rates $γ̇$⁠. Interparticle bonds are dynamically created and broken during the course of the simulation due to the cumulative effect of thermal fluctuations and shear deformation.^12,30 Here, to minimize the noise in our calculation due to thermal fluctuations in bond breakage and formation rates, we use different distance cutoff criteria for the formation and breakage events. We define a bond as being created between two colloidal particles when one particle comes within the surface-to-surface distance h_ij ≤ 0.08 of another. These bonds remain active until the particles drift apart to a surface–surface separation distance of h_ij ≥ 0.10. This is because under deformation, and as a consequence of using finite time steps in the numerical scheme, a bond may fluctuate for a number of steps in this range but retain the connection to another colloid.

Previous studies have suggested that rate-dependent rheology of gels can be discussed with regard to a dimensionless Mason number, $Mn=6πηγ̇a3/U0$⁠, as the ratio of convective shearing forces that promote bond breakage and the attractive forces that promote bond formation.^{13,20,31–33} Thus, we examine four different values of Mason numbers (dimensionless deformation rates) in the two regimes of shear-induced structuration and fluidization. Figure 2(a) shows the number of net bonds scaled by the total number of particles as a function of the total strain for different applied Mason numbers Mn. We observe that no appreciable net bond loss occurs before the small deformation γ ≲ 0.1 for all applied Mason numbers. Upon increasing deformation, the number of net bonds exhibits a gradual downturn, followed by a sharp decrease for γ > 0.1, indicating a yield point before or near γ ∼ 0.1. Expectedly, the net number of bonds decays more rapidly with increasing applied deformation rate, consistent with previous studies.^13,17 The initial unchanged number of total bonds can be misleading, as one may deduce that no bond is broken at these small deformations. However, during this time, new bonds are also formed, owing to the additional motion of particles due to the shearing motion. To further characterize these dynamics, we compute the number of newly formed and ruptured bonds separately as a function of deformation, and the results are presented in Fig. 2(b). While the number of both formed and ruptured bonds is negligible prior to small deformation γ < 0.1, it is quite clear that in sum they balance each other out, resulting in a rather steady total number of bonds. This is also consistent with previous findings of Colombo and Del Gado.¹² For larger deformation γ > 0.1, however, the number of ruptured and formed bonds depends strictly on the magnitude of the applied deformation rate. The lower the applied deformation rate, the higher the number of rupture and formation events over a fixed value of strain. For the lowest Mason number Mn = 0.35, the number of both formed and ruptured bonds is very large simply because the overall shear time for a given deformation at this rate is very long, and Brownian forces dominate the bond dynamics. As such, each colloid at these deformation rates and during a single total strain, on average, loses two bonds and forms two new ones. As the deformation rate $γ̇$ is increased, lesser total shearing time is available for these events, and fewer bonds are formed/ruptured. In these flow regimes, each particle on average loses more bonds than it is able to form, resulting in a sharper loss of total bonds in the system. This is expected, as increase in the Mason number also indicates a more prominent role for shearing forces in these regimes compared with the attractive forces. For the highest rate examined, $γ̇=0.5$⁠, however, the number of ruptured bonds is large, and far fewer bonds are formed as the bond dynamics are now controlled by the shear deformation predominantly. This is more evident by looking at the rate of bond loss vs bond gain, which is plotted in Fig. 2(c). For very small deformation γ < 0.1, although the rate of bond formation is greater than the bond rupture rate, it shows fewer variations compared with the rapid rate of bond breakage. Interestingly, the critical strain at which the highest bond rupture rate is observed is rather rate-independent and at γ ∼ 0.1 coincident with the strains at which the gel shows departure from the linear viscoelastic regime, before monotonically decreasing back to the bond formation rate. Note that, at long times, these two values will converge to a quasi-steady state. We also compute the stress vs strain dependence in Fig. 2(b) inset to show the connection between mechanical response and bond dynamics. We can observe that the critical strain (γ ∼ 0.1) for stress overshoot approximately coincides with the strain corresponding to the highest bond rupture rate. Nonetheless, note that due to the explicit representation of the solvent particles and as a consequence of imposing the flow profile through initialization of the solvent particle velocities, stress measurements at the very beginning of the simulations may include numerical artifacts.

These results are complemented by analyzing the non-affine distortion, which quantifies the resistance of the gel to an external deformation, plotted in Fig. 3. Non-affine deformations are defined as the deviation of real displacement from the imposed deformation averaged over all colloidal particles.¹² We notice that the non-affine displacement transitions to a different slope near γ ∼ 0.1 showing a signature of the yield point that is consistent with the results of mechanical response and bond dynamics. In addition, as expected, non-affine displacements grow slowly at high shear rates due to excess net bond loss, which results in the affine displacement of the particles along the velocity direction. At low imposed shear rates, however, less net bonds are broken and particles move rather stochastically due to dominant Brownian motion leading to a faster increase in non-affine displacements. It is also qualitatively in agreement with the previous simulation study by Colombo and Del Gado.¹² 

Since goal of this work is to find the topological origins of yielding in a gel, we limit the rest of our study to the bonds that exist at the beginning of the simulation (excluding the newly formed bonds) and break during a period of one strain unit of deformation. There are several attributes to each bond within a colloidal network, including bond position, bond orientation, contact number of bonded particles, bond life, etc. Among these attributes, the spatial position, contact-number of bonded particles, and lifetime are known to not be indicative of whether the bond will break under an applied deformation. This suggests that an overall network-level description of the bonds and their importance to the connectivity of the space-spanning structure is required instead. Thus, we borrow well-established concepts from network science to depict a quantitative measure of gel microstructure. Network science generally seeks to describe the emergent phenomena in different complex systems, such as social or economic networks, by focusing on the patterns of connections that construct the system, rather than an isolated look at the individual components of a system.^34,35 Thus, network science has been recognized as an increasingly transformative tool for studying complex systems in different areas, from biology, medicine, and neuroscience to epidemiology, ecology, and social sciences.^36–38 Commonly, a simplistic network of a complex system can be constructed by the representation of each constituent as a node, and their connections as an edge connecting two nodes. In the context of particulate and granular systems, different network science analyses have been developed and used to better understand the flow dynamics in such systems.³⁹ In a recent study, these tools have been used to identifying force clusters within dense colloidal suspensions exhibiting shear thickening behavior.⁴⁰ For attractive colloidal gels, Whitaker et al.⁴ showed that network analysis can be used to identify glassy local clusters of particles as the origin of elastic response. Here, having the spatiotemporal configuration of all particles within the simulation box, we construct a similar network representation of the initial gel. In this representation, Edge Betweenness Centrality (EBC), B_ij, can be rigorously defined as an attribute of each edge/bond within the network, and as a quantitative measure of its significance to the entire network. EBC ranks a bond based on the fraction of the shortest paths that pass through the bond between all pairs of nodes in the network. The EBC value of a bond e is computed as

where σ(a, b) is the total number of shortest paths connecting any pair of nodes within a system denoted as a and b, and σ(a, b|e) is the number of shortest paths that traverse through bond e. In practice, at the first step, the shortest path between all possible combinations of a and b are found within the network. In the next step, for each bond, the number of times that it appears in those shortest paths is calculated. Since total number of node combinations within a system is constant, bonds that most frequently appear in the list of shortest paths will have the highest EBC values. Hence, a bond with a high edge betweenness centrality value represents a bridge-like connector between two parts of a network, and thus its removal may affect the communication between many pairs of nodes through the shortest paths between them. Recent studies⁴¹ showed that potential failure in the disordered materials occurs at high centrality locations. Another bond attribute that must be considered in the yielding process is its orientation. The orientation of the bond, θ_ij, is defined as the angle between the vector of the bond and the extensional axis in the plane of the velocity–velocity gradient. This is because the bonds that are oriented near the extensional axis are expected to break more readily due to shear forces pulling them apart. Figure 4 shows the probability density distribution of bond orientation and edge betweenness centrality of the initial gel structure prior to shear deformation. As expected, bond orientation, θ_ij, shows a uniform distribution, i.e., the gel is formed isotropically at the quiescent condition and without any particular orientational preference. On the other hand, the edge betweenness centrality, B_ij, has a Gamma-like distribution. This suggests that there is a small population of bonds that are of high EBC values. In this population, ∼1% of all possible paths within a network traverse the given edge, B_ij ∼ 10⁻².

Edge betweenness centrality is related to a gel network’s connectivity, which can be better understood by conceptualizing conductivity and/or percolation of a network. For a network to be considered conductive or percolated, there essentially must be a path from each given point to another. As such, the relative significance of an edge to the number of these paths also reflects its significance to the overall conductivity or percolation of the network. Consequently, the removal of a bond with a high edge betweenness centrality will be potentially more detrimental to the loss of connectivity in the system. In the same manner that electricity is transmitted through different percolated paths in a conductive material, the transmission of stress through the gel network is essential for the system to retain its elasticity, i.e., not yield. Since bonds with higher EBC values play a larger role in providing this transmission route, we hypothesize that the loss of these bonds will be more detrimental to the overall integrity of the network and hence the yielding process. This is consistent with the idea of stress localization: the higher the number of paths that traverse a certain bond, the higher the possibility of accumulating stresses on that bond.

To test this hypothesis, we calculated the fraction of ruptured bonds within a total of γ = 1.0 deformation with respect to their initial bond orientation and edge betweenness centrality and the results are plotted in Fig. 5. For small applied Mason number Mn = 0.35, bonds are almost evenly ruptured with respect to their initial θ_ij [Fig. 5(a)]. Upon increasing shear rate, pronounced anisotropies in bond rupture emerge, and more of the bonds that are aligned along the extensional axis (θ_ij ∼ 0°) are broken. For the highest applied Mason number Mn = 17.28, the fraction of ruptured bonds aligned along the extensional axis is nearly seven (7) times greater than those aligned along the compressional direction. This is anticipated, as the lower applied deformation rates provide more time to the bonds to re-arrange and escape their breakage, as opposed to a fast breakage due to large displacements under higher Mason numbers. More importantly, bonds with lower B_ij are found to break less frequently at all applied deformation rates [Fig. 5(b)]. For all applied deformations, there is a clear growing trend of bond breakage as B_ij grows; that is a bond with higher EBC is more likely to break. This is even more significant when the chain relaxation processes are taken into account. Consider two adjacent bonds within a string of particles. Due to their topology, the EBC value of the two edges will be very close, if not identical (same paths crossing one edge must also cross the other). However, if one of the edges is ruptured due to accumulation/localization of the stress, the other one immediately relaxes upon the rupture and does not break as a result. Thus, only a fraction of these bonds is expected to break under applied deformation. The results presented here are found to be consistent for three independent simulation runs.

Since both the bond orientation angle and the edge betweenness centrality show evident correlations to the likelihood of bond breakage during the yielding process, it is logical to consider their collective effects as well. This is shown using color-coded breakage kymographs in Fig. 6, where the color intensity is mapped onto the breakage event likelihood. Overall, evidently, the bonds with a high EBC value and oriented in the extensional axis are clearly most frequently broken under all deformation rates. At small deformation rates, the EBC appears to be more significant, and by increasing the applied deformation, rate bond orientation begins to play a more significant role. One should note that at the smallest Mason numbers studied, Brownian effects are more significant and affect the dynamics of yielding over a much longer time compared with the fluidization regime.

Having identified bonds that are broken during the simulation and having the likelihood of their breakup events statistically, one can return to the initial gel and visualize those bonds as the topological origins of the yielding in colloidal gels. Snapshots of the initial gel with flagged bonds that were identified as broken during shear deformation are visualized in Fig. 7. Initial and broken bonds with B_ij in the top 10 percentile during a total of 0.1 strain of deformation in the ∇v–v plane are colored blue and red, respectively, and all existing bonds are colored gray for visual clarity. It is quite evident that the ruptured bonds are heterogeneous at both Mason numbers, such that some local structures are completely broken up by the end of the simulation and some structures remain rather intact for both gels. At the same time, the bonds with the highest values of EBC that are broken appear to be uniformly distributed throughout the structure. However, while those are more evenly orientated for the small Mason number, at the larger Mason number they appear as highly oriented along the extensional axis. Additionally, consistent with the chain relaxation micromechanical perspective described previously, the broken bonds with the highest EBCs appear as disconnected. This further confirms that upon rupture of a bond as a result of stress localization, the adjacent bonds relax and do not necessarily break. We emphasize here that across different simulation runs of the same gel structure, for a given applied Mason number, the probability that bonds with these EBC values will break remains the same. The exact bonds that break can vary, because upon the rupture of one bond, the adjacent bonds may relax in a different order. Therefore, this does not predict which specific bonds are breaking, but rather gives us insight into which characteristics lead to breakage events and the eventual yielding of the structure.

One of the most significant characteristics of attractive colloidal gels is their ability to resist relatively small stresses and exhibit critical yield stress demarcating the departure from a solid-like behavior to that of a viscoplastic fluid. Therefore, understanding the mechanisms under which a gel yields upon application of an external deformation is of utmost significance in the design and manipulation of its mechanical properties. In this work, employing a series of detailed particle-based simulations, we studied the topological origins of gel yielding, i.e., structural failure points under applied deformation. To characterize which bond rupture events are responsible for the yielding of the particulate network, we monitored the evolution of bond rupture/formation over the course of the shear simulation under a range of applied deformation rates. We found that two bond attributes, namely the bond orientation angle with respect to the extensional axis, and the edge betweenness centrality are the most informative quantities that correspond to the breakage of a bond under deformation and eventually yielding of a colloidal gel. Namely, bonds with the highest edge betweenness centrality that are aligned along the deformation field are significantly more likely to break under flow, regardless of the rate of applied deformation. Since the orientation angle is not a quantity that can be controlled in colloidal gels formed isotropically under quiescent conditions, our findings suggest that the edge betweenness centrality can be used independently to predict failure points within a gel structure a priori and only by knowing the initial topology of the particle network. As such, a reverse design of a network with a specific EBC configuration can be pursued as a viable route for the fabrication of gel structures with superior properties. Finally, our study was limited to short-range and weak attractions of the order of few thermal energies $[U0∼O(kBT)]$⁠. Future work will be required to better understand the applicability of these findings to different strengths/ranges of attractions, and perhaps develop a unified view of the yielding mechanism in colloidal gels.

This work was supported by the National Science Foundation (Grant Nos. CBET-2118962 and CBET-2104869). Computational resources were generously provided by the Massachusetts Green High-Performance Computing Center in Holyoke, MA.

The authors have no conflicts to disclose.

Deepak Mangal: Conceptualization (supporting); Data curation (equal); Formal analysis (equal); Investigation (equal); Methodology (equal); Software (equal); Validation (equal); Visualization (equal); Writing – original draft (equal). Mohammad Nabizadeh: Conceptualization (supporting); Data curation (equal); Formal analysis (equal); Investigation (equal); Methodology (equal); Software (equal). Safa Jamali: Conceptualization (lead); Funding acquisition (lead); Investigation (equal); Methodology (equal); Project administration (lead); Resources (lead); Software (supporting); Supervision (lead); Writing – review & editing (lead).

The data that support the findings of this study are available within the article.