Stress–stress correlations reveal force chains in gels

Soft particulate gels consist of particulate matter (polymers, colloids, proteins, etc.) aggregated in a solid matrix, which is embedded in a fluid and typically sparse and porous.¹ Gels of this type are found in the tissues of the body, food, drugs, personal care products, and even construction materials, such as the cement used in concrete.^2,3 Their structures can be very soft and re-configurable, and their stress response is determined by a complex interplay between molecular cohesion or surface interactions, microstructural reorganization, and external driving. As a consequence, understanding, predicting, and designing the mechanics and rheology of these materials remains very challenging.

While particulate gels can be extremely soft because of the overall low solid content, the strength of the interactions that drive the aggregation of the colloidal units into the gel network must be large enough, with respect to k_BT, to overcome thermal fluctuations and stabilize the gel structural elements and their connections.^4,5 This implies that the elasticity is largely enthalpic in nature in these gels, contrarily to polymer networks, and that the interaction strength, together with the total volume fraction of the solid content, controls elastic moduli, viscoelastic spectra, and nonlinear response. Remarkably, the interaction strength alone does not allow one to predict the gel properties. The morphology of the gel network, the structural elements, and their connectivity may change a lot due to different kinetics at play during the gel self-assembly, leading to huge variations of the mechanical response. This is sometimes the case even for identical compositions of the initial particle solution.

The solidification processes through which particulate gels form, initiated via, for example, irreversible fractal aggregation or equilibrium phase separation, are typically sources of frozen-in stresses⁶ and help build a memory in their history dependent response:⁷ these processes determine how local stresses and mechanical heterogeneities get embedded and remain frozen inside the gel structure during gelation. The rigidity and elastic modulus attained ultimately depend on the geometry of the microstructure that self-assembles during solidification to accommodate the local and global constraints of mechanical equilibrium. Consequently, significant structural correlations, induced, for example, by attractive interactions, are naturally incorporated in the way stresses are transmitted through the rigid backbone of the final solid. This process of self-organization can provide an explanation for the possible emergence of rigidity at very low densities in colloidal gels since the structural correlations developed are precisely those that satisfy the constraints of mechanical equilibrium, as shown in Ref. 8. The fact that this self-organization occurs in overall amorphous and spatially heterogeneous structures naturally leads to stress localization that manifests itself under deformation.^9–12 

The tendency to localize strain and stresses and the signature of mechanical heterogeneities, shown under deformation, suggest that stress transmission could be very anisotropic and strongly localized even if the material is at rest,¹³ just as a result of the microstructural development during gelation. As a consequence, force chains in a granular solid, i.e., forces transmitted across a sample along lines of grains forming a sparse network,^14–18 may be relevant to the physics of soft particulate gels. Detecting force chains and identifying the stress bearing part of gel structures are, however, incredibly difficult in experiments and elusive even in numerical simulations of model materials, as they require disentangling mechanical heterogeneities from structural heterogeneities and rely upon threshold values for local stresses that are not easily justified. For amorphous solids, recent theoretical studies have clarified that just the constraints of mechanical equilibrium guarantee that the spatial correlations of the stresses are long-ranged and anisotropic.^19–22 The emergence of elasticity can then be obtained, even in granular solids where the grain contacts cannot be regarded as elastic springs, within a new theoretical framework that can be mapped onto a tensorial electromagnetism with vector charges (VCT).^23–25 In particular, this VCT framework provides analytical predictions for the stress–stress correlations, demonstrating that they are not only long-ranged but also anisotropic, and clearly identifies the presence of force chains in disordered assemblies of grains that are jammed under an external load.

Here, we use the predictions of the VCT framework to compute the stress–stress correlations in model particulate gels obtained from 2D and 3D numerical simulations. We find that the intensity patterns of these correlations in Fourier and real space share several common features with the patterns detected in experiments and simulations of granular solids, as well as those predicted by the VCT framework.^24,25 In particular, we recover the anisotropy of the stress–stress correlations, which therefore becomes an indicator of the presence of force chains in particulate gels. Moreover, we show that the stress–stress correlations distinguish floppy from rigid gels and are sensitive to the distance of the model gel from the rigidity threshold. These results pave the way to developing a further understanding of stress transmission in soft materials and provide the first basis to develop a VCT framework for soft gel materials.

This paper is organized as follows: Section II provides the details of the calculations performed: In Subsection II A, we briefly review the basic ingredients of the VCT framework and its predictions for the stress–stress correlations in granular solids, while Subsection II B contains information on the numerical models and simulations used. In Subsection II C, the calculations of the stress correlations are described. We then discuss the results in Sec. III and provide a summary and outlook in Sec. IV.

At the heart of the mechanical response of jammed amorphous solids is the nature of the force-bearing networks that obey the constraints of mechanical equilibrium. Implemented locally: each particle satisfies the constraints of force and torque balance. Here, bold symbols with hats denote tensors, bold symbols without hats denote vectors, and nonbold symbols denote components of tensors and vectors.

The VCT framework makes explicit predictions about the stress–stress correlations and the stress response in terms of the emergent elastic modulus tensor, $K̂$⁠. Here, we summarize crucial features of the stress–stress correlations. (i) A hallmark of the theory is the appearance of a characteristic pinch-point in the stress–stress correlations, C_ijkl(q) = ⟨σ_ij(q)σ_kl(−q)⟩. A function of the wave vector ∣q∣ has a pinch-point singularity in ∣q∣ = 0 if its value at q = 0 depends only on the orientation of any line passing through q = 0. This means the dependence of the stress–stress correlations is only on the angular coordinates in Fourier space for q → 0. These predictions have been tested against experimental measurements in frictional granular materials and against model assemblies of frictionless soft grains in 2D and 3D. (ii) The chain-like structures seen, for example, in experiments using photoelastic discs^24,28 and commonly referred to as force-chains, are a visual representation of the highly anisotropic nature of C_ijkl. (iii) The q-space correlations predicted by the theory, and observed in experiments, imply a power-law decay at large lengthscales: ⟨σ_kl(r)σ_kl(0)⟩ ∝ ±1/r^D, where the plus sign appears for longitudinal correlations and the minus sign for transverse correlations. In this paper, we analyze the stress-correlations in numerical simulations of 2D and 3D particle gels using the VCT framework.

Gel configurations in two and three dimensions are obtained using models of N interacting colloidal particles that undergo gelation as described in Refs. 8, 10, and 29–31. We perform Molecular Dynamics (MD) simulations in a cubic (square in 2D) box of size L with periodic boundary conditions. To prepare the gels, we first reach thermal equilibrium at a high temperature and then slowly quench the particle configurations to different target temperatures T, i.e., with different cooling rates C_r, using a Nosé–Hoover thermostat and allowing the particles to aggregate and form gel networks.^29,32 The mechanical equilibrium in the final gel states is obtained by slowly withdrawing all kinetic energies with overdamped dynamics, following the protocol described in Refs. 8, 10, and 31. All the simulations are performed using LAMMPS³³ modified by us to include specific particle interactions.

For 2D gels, we consider N = 10 000 monodisperse particles interacting via a short range attractive potential $U(r)=Aϵad/r18−d/r16$⁠, where r is the interparticle distance, and A and a are dimensionless constants. The values of A = 6.27 and a = 0.85 are chosen to obtain a short-ranged attractive well of depth ϵ and range r_c ≈ 0.3d (the potential is cut and shifted to 0 at large distances). In 3D, we consider N ≈ 16 000 monodisperse particles that interact via the same short range attractive interactions plus an angular term that introduces bending rigidity to inter-particle bonds.^10,30,31

All the simulation quantities described in the following are expressed in reduced units; thus, d is used as the unit of length, ϵ is used as the unit of energy, and m (the particle mass) is used as the unit of mass. The reduced temperature is expressed in units ϵ/k_B, where k_B is the Boltzmann’s constant. The unit of stress is therefore ϵ/d³. When performing the quenches into aggregation from a high temperature, the cooling rate C_r is defined as ΔT/Δt, where ΔT is the distance between the initial and target temperatures in the gelation protocol described earlier and Δt is the duration of the quench, with the unit time τ₀ being $dm/ϵ$⁠. We define an approximate volume fraction $ϕ=πd2N4L2$ and $ϕ=πd3N6L3$ in two and three dimensions, respectively, where d is chosen as the particle diameter.

We consider that a gel configuration is rigid if the elastic component of the linear viscoelastic shear modulus at a low frequency is larger than the viscous component, which is true if there is a percolating rigid cluster. In 2D, we characterize such rigid gel configurations using the pebble game algorithm, which identifies a rigid cluster as a cluster where all the degrees of freedom are blocked by contacts or constraints except for the degrees of freedom of the overall translations or rotations of the whole cluster. Therefore, in a rigid gel, there could be isolated particles or other small rigid clusters that are not part of the spanning cluster; however, the floppy modes associated with these isolated particles are not relevant to the low frequency elastic response of the sample. We call a gel floppy if in the particle configuration though there is a spanning network of aggregated particles, such network is not rigid according to the pebble game algorithm.^8,34 In 3D, we take another approach. We compute the linear viscoelastic response of the gels and obtain their low frequency moduli.^29,30,35 For rigid gels, the low frequency storage modulus is larger than the loss component.

In 2D, we show data for ϕ = 0.5 at T = 0.3 (floppy gels) and T = 0.18 (rigid gels). In 3D, we analyze gels obtained using the cooling rates C_r = 9 × 10⁻⁵ϵ/(τ₀k_B), and 9 × 10⁻²ϵ/(τ₀k_B) at ϕ = 0.125, 0.075, and 0.05. All the data for the stress correlations have been averaged over 200 and 20 independently generated samples in 2D and 3D, respectively.

From the stress correlation calculations, we construct maps of the correlation intensity in Fourier space in polar (q, θ) representation for 2D gels for a range of q from 2π/L to 56π/L. For the angular variation of the correlations, along with radial averaging, the data are also averaged over a small angular bin width of $≈6°$ to reduce noise. For comparison, we also compute the correlations in real space along both x and y direction with the coarse-graining length d.

All stress correlations are scaled by the maximum value of the C_xxxx correlation function.

We now discuss the results obtained for 2D and 3D model particulate gels in mechanical equilibrium.

As discussed in the Introduction, soft particulate gels have strong structural heterogeneities that depend on the aggregation kinetics and the path to gelation. These structural heterogeneities are associated with stress heterogeneities; however, the microstructural origin of the stress heterogeneities is not easily identified. Figure 1 shows snapshots of a rigid gel close to the rigidity threshold⁸ (bottom panels) and a floppy gel (top panels) from our 2D simulations. The particles are colored according to their magnitude of xx and xy components of $σ̂p$ computed using Eq. (3). The configurations in the top and bottom panels show connectivity percolation; however, the spanning rigid cluster, which was identified using the pebble game algorithm, is present only in the configuration of the bottom panel. In the floppy gel case, the particle stresses show strong, spatially uncorrelated fluctuations. In the rigid case, we observe that the particle stresses are spatially correlated and extended. For these numerical samples, these components of the stress tensor have a magnitude varying between −10ϵ/d³ and 8ϵ/d³. If we use the whole range for σ_xxp and σ_xyp in the color maps, the gels appear mechanically homogeneous, as they are overall soft. By zooming in the stress variation down to Δσ_ij ≃ 1ϵ/d³, chain-like patterns become visible in the stress maps of rigid gels. Floppy gels feature relatively large stress fluctuations even at this large stress magnification. These, however, appear randomly distributed and uncorrelated. For rigid gels, instead, the stress fluctuations display an anisotropic pattern, which strongly suggests an anisotropic and localized stress correlation across the sample and is reminiscent of force chains in granular packings.^14–18,28

These qualitative observations can be quantitatively confirmed by analyzing the intensity maps of the stress correlations for the corresponding component of the stress tensor, as shown in Fig. 2 (left). While for floppy gels the stress–stress correlation maps do not feature any distinctive pattern, for rigid gels they do. The patterns indicate both the presence of correlations that span the whole system and that they do so in a very directional and anisotropic fashion, in spite of the particle interactions and the overall gels being isotropic. The C_xxxx map obtained from the rigid gels shows a large q cutoff in the correlations, which simply reflects the fact that well inside the gel branches the material is spatially and mechanically homogeneous. However, the correlation intensity is significant (along the q_x axis) even for the smallest q, i.e., for distances up to the whole system size. The value of C_xxxx being zero along the line q_y = 0 indicates that the forces can propagate approximately only along their own direction and very unlikely perpendicular to it. At q ≈ 0, the correlation intensity depends only on θ, which is a signature of a pinch-point singularity. This is a direct consequence of force balance, as discussed in Sec. II A, and provides a quantitative evidence to the presence of force chains. The same characteristics, i.e., the long range correlations and the strong anisotropy indicated by the angular dependence, are also present in the four-lobe pattern of the C_xyxy map for the rigid gels (Fig. 2 right). We note that for C_xyxy the signal at low q is much weaker, since shear stresses are vanishing as there is no deformation applied.

The two-lobe pattern in C_xxxx and four-lobe pattern in C_xyxy for the rigid gels are identical to the correlation patterns already obtained in experiments and simulations of granular solids, confirming the pinch-point structure of stress correlations in Fourier space predicted by the elasticity theory of the VCT framework.²⁴ In Fig. 3, we also show the C_yyyy (top) and C_xxyy (bottom) maps for the 2D rigid gels in Fourier space (left) and the corresponding radially averaged angular plots (center). The stress correlation maps in the real space (right) complement, and confirm, the insight gained with the maps of the correlations in Fourier space. The remaining stress correlations for 2D gels that consider all other stress tensor components are shown in  Appendix A. These maps quantitatively establish the similarity of stress transmission in soft gels and in granular solids. They support the idea that the rigidity and elasticity of soft particulate gels, when particle interactions are sufficiently large with respect to k_BT, can be fundamentally understood as emerging properties that are the result of the local and global constraints imposed by mechanical equilibrium, rather than being just a consequence of their microscopic interactions. We note that, as briefly discussed in Sec. II A [see Eq. (2)], in the VCT framework, the angular dependence of the correlation intensity (Fig. 3 center) provides direct information on the properties of the emerging elastic tensor $K̂$⁠.

In the case of 3D gels, we have computed all 21 stress–stress correlations. We discuss some of them here in detail, while the remaining ones are shown in  Appendix B. A snapshot of one of the gel networks [at volume fraction ϕ = 0.125 and obtained with cooling rate C_r = 9 × 10⁻⁵ϵ/(τ₀k_B)] is shown in Fig. 4 (top panel). We have colored the particle bonds in red if the xx component of the particle stress is positive (tension) and blue if it is negative (compression) to highlight how, even with this information, the presence of force chains or stress localization cannot be easily detected. However, the maps of the corresponding stress correlation function clearly reveal their presence; furthermore, in this case, all stress–stress correlations are long-ranged, as determined by the mechanical equilibrium constraints;^19,22 the 2D Hammer projection of the correlation intensity for the fluctuations of the xx component of the stress tensor (Fig. 4, bottom panel) are strongly anisotropic, with the angular pattern indicating that stress transmit across the sample along specific and localized directions.

When comparing our data to the results in Ref. 24, we notice interesting quantitative differences in the corresponding angular plots. This suggests that the obvious significant differences between the overall elastic behavior of gels and granular solids can also be further investigated with this approach. In Fig. 5, we compare the intensity pattern of C_zzzz in rigid gels (top) and granular solids (bottom, from Ref. 25) to point out the striking difference in the angular dependence. As we increase the dimensionality of the space from 2D to 3D, the number of possible stress–stress correlation functions increases from 6 to 21, which provides more opportunity to identify the differences in the correlation patterns between the granular and gels case. We also note that our 3D gel model also features a bending rigidity term in the particle–particle interactions, which is not present in the 2D case. This might also account for some additional differences in the intensity patterns between 2D and 3D, compared to the granular case. The map obtained for a model granular solid shows that the stress–stress correlations are basically independent on the azimuthal angle, Φ, and only depend on θ, whereas the pattern for the gels indicates a strong dependence on both θ and Φ, which confirms that this approach is able not just to reveal similarities and common traits but also to identify differences. With this respect, we think that this striking difference in the angular dependence may be related to the fact that all structures in a granular solid can only transmit compression or shear stresses, whereas in gels tension can also be transmitted. All the remaining stress correlations from 3D gels at ϕ = 0.125 are shown in  Appendix B.

The patterns of the stress correlations in gels are sensitive to the distance from the rigidity threshold, as shown by the C_xzzz maps in Fig. 6 at different ϕ = 0.125, 0.075 and 0.05 obtained for cooling rate C_r = 9 × 10⁻⁵ϵ/(τ₀k_B). As the volume fraction increases the local connectivity and gel moduli increase. At ϕ = 0.05, the networks are sparse, spatially heterogeneous, and barely rigid, as also discussed in Refs. 30 and 37, and their elasticity modulus is very small. When particle volume fraction increases, the gels obtained with the same cooling rate become stiffer, more locally connected, and move away from the rigidity threshold. The stress correlations at different densities, interestingly, have the same general pattern. However, this pattern becomes increasingly blurred as the gels approach the rigidity threshold (Fig. 6 from left to right). A similar trend is observed in the correlations of all stress components. These findings support the idea that the stress correlation intensity variation is directly related to the elastic moduli, suggesting that they may be used to detect the distance from the rigidity threshold and provide information on the marginal stability of the gels. We also confirm these findings by varying, for the same particle volume fraction, the cooling rate at which the gel is formed, since the increasing this rate leads to gels that are weaker and closer to the rigidity threshold.³⁷ By comparing, for example, the stress correlation map in Fig. 7, which shows C_xzzz for ϕ = 0.125 at C_r = 9 × 10⁻²ϵ/(τ₀k_B), with the first map on the left of Fig. 6 [corresponding to the same volume fraction and C_r = 9 × 10⁻⁵ϵ/(τ₀k_B)], we can notice how again the correlation pattern becomes more blurred for weaker and more marginal gels. Clearly, the gels with the same density but different network elasticity result in changes in the stress–stress correlations. Since these are systems that are out of equilibrium, it is not obvious that there is a relation between the correlation function and the response function. Our results show that there could be such a relation. The VCT theory relates stress correlations to the elastic moduli with a one-to-one map between the full set of stress-correlations and the components of the $K̂$ tensor.^24,25

In summary, we have demonstrated the presence of force chains in soft particulate gels by computing the stress–stress correlation functions in 2D and 3D model gels. The stress correlation patterns distinguish the rigid gels from the floppy gels and display the same general characteristics as granular solids, i.e., the strong angular dependence and the pinch point singularity predicted by a stress-only theoretical framework (the VCT framework in Ref. 24), where the elastic response is a property emerging only from the constraints of mechanical equilibrium without resorting to the constitutive relation of linear elasticity. We also demonstrate that the stress–stress correlation patterns highlight distinctive differences between gels and granular solids and that these patterns are sensitive to variation of the gel moduli due to changes in the gel topology and marginal stability, not simply to the gel density. These results are consistent with the idea that the stress–stress correlation patterns can provide further information on the tensorial properties of the emerging elasticity in these materials. Stress–stress correlation patterns have been measured in experiments on granular solids using photoelastic disks (see Ref. 24); however, this is clearly much more challenging for experiments on particulate gels, where measurement techniques to extract local stresses and their spatial distributions are still being developed.^38–41 Moreover, using the angular dependence of the stress–stress correlation functions obtained from simulations to extract the properties of the elastic tensor would provide predictions for the overall mechanical response to different deformation modes (e.g., shear, compression, and extension), which are accessible in mechanical and rheological experiments. This information, testable on a broad range of gels and experimental conditions, could then be fed back to the theory to obtain the stress–stress correlation patterns from the experiments. Hence, further developing the analysis proposed here seems a promising route to gain novel insight into the emerging elasticity of soft particulate gels and its connection to stress localization.

We thank Shang Zhang, Jishnu N. Nampoothiri, and Minaspi Bantawa for valuable discussions and configurations. This work was supported by the National Science Foundation (Grant Nos. NSF-DMR-2026842, NSF-DMR-2026825, and NSF-DMR-2026834) and the Clare Boothe Luce Program. Acknowledgment is made to the donors of the American Chemical Society Petroleum Research Fund for partial support of this research.

The authors have no conflicts to disclose.

H. A. Vinutha: Conceptualization (lead); Data curation (lead); Formal analysis (lead); Writing – original draft (lead). Fabiola Diaz Ruiz: Conceptualization (equal); Data curation (equal); Formal analysis (equal); Writing – original draft (equal). Xiaoming Mao: Conceptualization (equal); Data curation (equal); Investigation (equal); Methodology (equal); Supervision (equal); Writing – review & editing (equal). Bulbul Chakraborty: Conceptualization (equal); Data curation (equal); Investigation (equal); Methodology (equal); Supervision (equal); Writing – review & editing (equal). Emanuela Del Gado: Conceptualization (equal); Data curation (equal); Formal analysis (equal); Writing – original draft (equal).

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

In Fig. 3, we show data for C_yyyy and C_xxyy. Here, we show all the remaining stress correlations for the rigid gels. The correlation functions are computed using Eq. (5). The angular dependence of these correlation functions matches with the granular case (Fig. 8).²⁴ 

In this paper, we show C_xxxx, C_zzzz, and C_xzzz correlation functions of gels at ϕ = 0.125. In Fig. 9, the remaining 18 correlation functions are shown.