Pinch-point singularities in stress–stress correlations reveal rigidity in colloidal gels

Abstract

We demonstrate that the spatial correlations of microscopic stresses in 2D model colloidal gels obtained in computer simulations can be quantitatively described by the predictions of a theory for emergent elasticity of pre-stressed solids (vector charge theory). By combining a rigidity analysis with the characterization provided by the stress correlations, we show that the theoretical predictions are able to distinguish rigid from floppy gels, and quantify that distinction in terms of the size of a pinch-point singularity emerging at large length scales, which, in the theory, directly derives from the constraints imposed by mechanical equilibrium on the internal forces. We also use the theoretical predictions to investigate the coupling between stress–transmission and rigidity, and we explore the possibility of a Debye-like screening mechanism that would modify the theory predictions below a characteristic length scale.

---

1. Introduction

Soft particulate gels are found in a wide range of materials from consumer products to bio-manufacturing. They are made of particulate matter (polymers, colloids, proteins, etc.) aggregated in a solid matrix, which is embedded in a fluid and typically both sparse and porous.¹ The soft and adaptable nature of these gel structures make them materials that can be easily manipulated in multiple ways from squeezing, stretching, spreading on surfaces and 3D printing.^2–4 However, process-control is very challenging because of our lack of understanding of the mechanistic link between structure and function of soft gels. A complex interplay between particle cohesion and external driving determine their stress response.^5–9 Additionally, the solidification process, which occurs in non-equilibrium conditions, can lead to frozen-in stresses and history dependent stress response.^10–15 All these factors make it particularly challenging to develop theoretical frameworks that can capture the collective behavior of the constituent particles giving rise to its macroscopic emergent behavior.

Within the efforts to elucidate the rigidity of a broad range of amorphous solids, including soft particulate gels, rigidity percolation has emerged in recent years as the main theoretical framework,^16–21 based on the idea that locally rigid structures (due to mechanical constraints such as contacts, chemical bonds or steric repulsion) percolate through the material. This concept has been pivotal in clarifying the role played by microscopic properties in giving rise to macroscopic mechanical response. The issue is identifying which parts of the microstructure, self-assembled during the solidification process, may satisfy, and how, the conditions corresponding to mechanical equilibrium and Maxwell's criterion for rigidity.

In the context of granular solids, it has become increasingly evident that stress correlations emerging from the constraints of mechanical equilibrium are anisotropic, sub-dimensional and long-ranged.^16,18,22 The framework developed by Nampoothiri et al.^22,23 has demonstrated a rigorous mapping of a tensor gauge theory – a Vector Charge Theory similar to field theories for quantum spin liquids²⁴ – to a theory of the mechanics of amorphous solids (VCTG), formulated in terms of mechanical constraints. The mapping leads to an “emergent elasticity” and an accurate description of stress localization in granular solids, predicting quantitatively the outcomes of experiments with photoelastic grains and computer simulations of model granular solids. VCTG demonstrates that these specific features of the stress correlations, i.e. being anisotropic, sub-dimensional, and long-ranged, emerge purely from the constraints, and dictate mechanical response and elasticity. This aspect could not be captured by the rigidity percolation framework and has relevance well-beyond granular solids, since it applies broadly to athermal amorphous solids. Soft particulate gels are a class of materials where thermal fluctuations are usually important drivers of self-assembly and solidification, however high interaction strengths are required to form rigid, stress-bearing structures that control both dynamics and elasticity.^{15,19,25–29} Interestingly, it has been shown that spatial correlations of stress fluctuations measured in model particulate gels through computer simulations have the same characteristics predicted by the VCTG mapping,³⁰ indicating that also in these very soft materials the constraints imposed by mechanical equilibrium manifest themselves in the same way. These findings suggest that there may be a connection between the rigidity percolation framework and the VCTG approach. Establishing such a connection and clarifying which aspects of soft-particulate-gels rigidity and mechanics can be quantified by the VCTG framework could allow for a deeper understanding of these two different theoretical frameworks. Further, it can clarify similarities or differences across complex amorphous materials.

Here, we wish to contribute to some of these questions and explore the implications of the VCTG paradigm for soft particulate gels. We use the predictive framework of VCTG to interpret the stress–stress correlations measured in model colloidal gels through 2D computer simulations and show consistency between the rigidity percolation framework and VCTG framework in inferring the rigidity onset. From the stress correlations, we also estimate the emergent elastic modulus of gels as predicted by VCTG, which reveals correlations between normal stresses and mechanical strength in this class of materials. Finally, the comparison between the VCTG predictions and the simulations motivates us to extend the VCTG framework by incorporating a length scale, associated with stress fluctuations, which is absent in the case of granular solids.

The paper is organized as follows: Section 2 briefly reviews the basic ingredients of the VCTG framework and its predictions for the stress–stress correlations in rigid gels. Section 3 contains information on the numerical model and simulations used. We then discuss the results in Section 4 and provide a summary and outlook in Section 5.

2. Theoretical framework

2.1. Stress-only theories of elasticity

Over the last decade, U(1) tensor gauge theories with fractonic charges²⁴ have been mapped on to stress-only elasticity theories of solids with pre-stress. In mapping to elasticity theory, kinematic constraints on the stress field are mapped to Gauss's laws^22,31 of the tensor gauge theories. The observable fields are symmetric second rank tensors, and two distinct Gauss's laws are possible, one with vector charges (VCTG) and the other with scalar charges (SCT).²⁴ The pre-stress can appear in this type of theory as the analog of the familiar polarization field present in a dielectric in standard electrostatics. Dielectric screening screens out unbound charges. In contrast to standard electrostatics, however, this polarization field is a second rank tensor. In VCTG, pre-stress emerges from screening external forces, and in SCT, they appear as a mechanism to accommodate Gaussian curvature of the metric that would have existed if all pre-stress is removed.^31–40 These theories also accommodate a Debye-type screening that occurs, as in electrostatics, via free charges and where a screening length emerges.^32,41–44 For a stress-free state in mechanical equilibrium, geometric compatibility refers to a zero curvature reference metric. For a disordered state with residual stress in mechanical equilibrium, such reference metric typically involves Gaussian curvature, i.e., being geometrically incompatible. SCT describes such incompatibility using scalar charges in 2D, and allows for a duality mapping between defects in crystalline solids and these scalar charges.³¹ This theory has also been adapted to amorphous materials.³² The mapping between Gaussian curvature and scalar charge limits the SCT to 2D materials, and the mapping of momentum conservation to Faraday's law of the corresponding tensor gauge theory makes it difficult to incorporate external body forces into SCT. In systems where external forces play a role in creating pre-stress, the kinematic constraint, ∂_iσ_ij = f_j, reflecting momentum conservation (force balance) in athermal solids, maps to the VCTG Gauss's law.²⁴ The vector “charges” are the unbound forces, f. The stress compatibility condition (vanishing Gaussian curvature) appears in VCTG as the static limit of the corresponding Faraday's law.²² Therefore, the VCTG framework can accommodate both sources of pre-stress.

In this paper, we use the VCTG framework to analyze stress–stress correlations in gels undergoing a floppy to rigid transition through numerical simulations of dynamical arrest, as deduced from a rigidity percolation analysis.¹⁹

2.2. VCTG elasticity

Using VCTG, the pre-stresses existing in an amorphous solid can be described as the polarization field in a dielectric, in the electrostatic limit of VCTG.^22,23 Crucially, VCTG emerges upon ensemble averaging over all states that are subjected to the same set of unbound forces, f, i.e., forces that are not bound up in pairs as in interaction forces. In jammed granular solids, for example, this translates to the externally imposed shear stress and/or pressure. The reasoning, which is familiar in the context of many frustrated systems,⁴⁵ is that there is no unique reference state, and that states that obey the constraint of force balance (Gauss's law) are equiprobable.

VCTG is a stress-only formulation of elasticity theory: there is no unique reference state or metric, and therefore no concept of strain. In contrast, as in an electrostatic dielectric, there is an unscreened stress field, Ê, and the screened stress field is σ_ij = E_ij + P_ij, where P̂ is the induced polarization field. In the linear dielectric version of VCTG, P̂ = χÊ, or equivalently, [small sigma, Greek, circumflex] = Λ⁻¹Ê, where Λ⁻¹ is a 4th rank tensor that appears as a coupling constant in the field theory and gives rise to the emergent elastic moduli.^22,23 The picture is that a stress-bearing structure emerges as a consequence of satisfying Gauss's law with given external forces. The compatibility relation in VCTG elasticity appears as the static limit of the analog of Faraday's law: ∇ × ∇ × E = 0. There is no such condition imposed by VCTG on the induced stress, P̂.

VCTG is a continuum field theory that predicts response and correlation functions at length scales much larger than microscopic ones such as the particle size, and it provides a rigorous basis for the emergence of stress heterogeneities such as “force-chains”.^22,23,30,44 Both the pinch-point singularity and the symmetries of the response and correlation functions are features of the Maxwell's equations of VCTG.²² In this paper, we test the applicability of the linear dielectric theory across the rigid-floppy transition. We also look for the possible emergence of a Debye-type screening length in the stress–stress correlations.

3. Numerical simulations

Gel configurations in two dimensions are obtained using a 2D model of N (=10 000) interacting monodisperse colloidal disks that undergo gelation as described in ref. 19. The particles interact via a short range attractive potential U(r) = Aε(a(d/r)¹⁸ − (d/r)¹⁶) where r is the interparticle distance, d is chosen as the particle diameter, 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). We perform molecular dynamics (MD) simulations in a square box of length L ≈ 125d and 142d, 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 a target temperatures T using a Nosé–Hoover thermostat, allowing the particles to aggregate and form gel networks.^46,47 The mechanical equilibrium in the final gel state is obtained by slowly withdrawing all kinetic energies with overdamped dynamics, following the protocol described in ref. 19, which mimics the slow aging typical of these systems.^7,48 The resulting ensemble of states satisfy the athermal momentum conservation constraint since every particle is in force balance.

Once all kinetic energy has been withdrawn, the stress tensor [small sigma, Greek, circumflex] is computed for each configuration using the following equation:

where L² is the simulation box area, r_p,s is the vector connecting the center of particles p and s, and f_p,s is the force between them. The index s iterates over the neighbors of particle p within the interaction range r_p,s ≤ r_c. We obtain the pressure P from the trace of the stress tensor. Following ref. 30, we compute the stress correlations in Fourier space using,

C_ijkl(q) = 〈Δσ_ij(q)Δσ_kl(−q)〉 (2)

where q denotes the wave vector, r_p is the pth particle position vector, σ_ij p denotes the ij-component of the stress tensor of particle p, and Δσ_ij represents the fluctuation of the stress tensor with respect to the average value 〈σ_ij〉 within the corresponding individual configuration. The angular brackets 〈⋯〉 in eqn (2) denote an ensemble average over several configurations corresponding to the same target temperature T used in the gel preparation and the same particle volume (surface) fraction. For each of these samples we perform an average over all the particles before taking the ensemble averages. The stress fluctuations in Fourier space, obtained through these calculations, are then coarse-grained by imposing a cutoff at large q, corresponding to q_max = 2π/d, i.e., we do not consider any stress fluctuations occurring at length scales shorter than d. The lower q cutoff in Fourier space is set by the simulation box size q_min = 2π/L.

All the simulation quantities described in the following are expressed in reduced units; 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 constant. The unit of stress is therefore ε/d² in 2D. We define an approximate volume (surface) fraction . All the simulations are performed using LAMMPS.⁴⁹ All the data for the stress correlations are averaged by degree of rigidity percolation, sampling from datasets with more than 100 independently generated configurations.

Following ref. 19,30, for the discussion here we have selected data for the same T = 0.18 ε/k_B and particle surface fraction ϕ = 0.39 and 0.5, which were located respectively in the floppy and rigid region of the rigidity phase diagram obtained for this model. A specific gel configuration is considered to be rigid if the pebble game algorithm⁵⁰ identifies a spanning rigid cluster, i.e., a network of particles where all the degrees of freedom are blocked by contacts or constraints, allowing only for overall translations or rotations of the entire network. A gel configuration is instead considered floppy if the corresponding spanning network of aggregated particles is not rigid according to the pebble game algorithm.^18,19,50 The phase boundary reported by Zhang et al.¹⁹ for the rigidity phase diagram is based on the percolation probability measured over many independently generated samples at various T and ϕ.

4. Results

We begin by verifying the predictions of VCTG theory for gels at ϕ = 0.5, a surface fraction which corresponds to the rigid region of the rigidity phase diagram.¹⁹ At this high surface fraction, the gel configurations are most likely to be rigid, and hence their stress–stress correlations in Fourier space should exhibit a pinch-point structure: we should observe long-range correlations with correlation intensity that has only angular dependence. In Fig. 1(a)–(c), we show selected stress–stress correlations maps of gel configurations at ϕ = 0.5. In order to produce these plots, the data are re-interpolated onto a polar grid in Fourier space with an interpolation distance of q = 1/d. The maps display precisely the two- and four-fold patterns predicted by VCTG²² for C_xxxx, C_xxxy, and C_xyxy, respectively, as well as the pinch-point profile at q ≈ 0. Note that in 2D we have in total six stress–stress correlations and, while for brevity we show only three of them, all of them satisfy the theoretical predictions consistent with the results in ref. 30. The angular variation of the stress correlation intensities are shown in Fig. 1(d)–(f). The polar grid used in the maps allows for grid points to be allocated into angular averaging bins with even populations. The angular variation of all six components of C_ijkl(q) measured are perfectly fit by the predictions of VCTG shown by the lines in the angular plots in Fig. 1. The fit to the theory for the angular dependence of the stress correlations can be used to extract the elastic modulus of the gels, as discussed below.

Fig. 1 Components of ensemble-averaged stress–stress correlations for ϕ = 0.5 gels. Correlation intensity of the function C_ijkl(q_x, q_y) are shown in the top row for i_jkl = xxxx (a), xxxy (b), and xyxy (c), normalized by the maximum value of C_xxxx. The red series in (d)–(f) shows the same correlation functions averaged into angular bins in Fourier space (error bars shown correspond to standard error from averaging into angular bins), with fit lines to isotropic VCTG theoretical forms shown in black. Within the set of configurations generated with these parameters, the configurations averaged here were those that exhibited rigidity percolation in two directions.

4.1. Rigidity percolation and VCTG theory

According to the rigidity percolation framework,^18,19 gel structures require a mechanically stable spanning cluster, which is able to transmit stresses (in addition to satisfying the requirements of connectivity percolation). The rigidity percolation analysis explicitly takes into account the gel network architecture in which mechanical constraints (local and global) operate to obtain the macroscopic emergent rigidity. Zhang et al.,¹⁹ had studied the nature of the rigidity transition in this gel model and obtained a rigidity phase diagram in the temperature-density plane. For the preparation temperature T = 0.18ε/k_B, a fraction of gel configurations formed at both ϕ = 0.39 and ϕ = 0.5 may have percolating rigid clusters. This fraction was significantly higher in the latter case, placing these two particle surface fractions on different sides of the rigidity boundary. The phase of a particular location on the rigidity phase diagram is merely defined by whether it is more likely to produce a configuration that is rigid or floppy: both rigid and floppy configurations can be generated for all choices of (ϕ, T) near the phase boundary.¹⁹

In order to understand the role of the nature of rigid clusters on gel elasticity and connect rigidity percolation theory and VCTG theory, we further split the gel configurations based on the degree of percolation of the largest rigid cluster. This can be 0, 1, or 2, if the largest rigid cluster does not percolate (0), percolates in only one direction (1), or in both directions (2). The latter case corresponds to a rigid network, and the former two are considered floppy. The fact that the two datasets chosen for this analysis are on opposite sides of the rigidity phase boundary¹⁹ implies that the particle fraction of ϕ = 0.39 produces predominantly floppy configurations, whereas ϕ = 0.5 leads to predominantly rigid configurations. In Fig. 2, we show snapshots of gel configurations at ϕ = 0.39. The subfigure adjacent to Fig. 2(a) is a zoomed-in snapshot corresponding to the blue subsection of the simulation box, which has a linear dimension of ≈25 particle diameters d and where individual particles are clearly visible, whereas the other snapshots shown in Fig. 2(a)–(c) show the full simulation box, whose linear size is L = 142d. The three snapshots respectively correspond to no percolation (a), 1-direction percolation (b) and 2-direction percolation (c). More specifically, the low-density dataset on the floppy side of the phase diagram consists of 12 configurations that percolate in 0 directions, 102 configurations that percolate in 1 direction, and 65 configurations that percolate in 2 directions, resulting in a rigidity percolation probability of 0.36. The higher-density dataset on the rigid side of the phase diagram (ϕ = 0.5) consists of 1 configuration that percolates in 0 directions, 47 configurations that percolate in 1 direction, and 76 configurations that percolation 2 directions, resulting in a rigidity percolation probability of 0.61.

Fig. 2 Snapshots of gel configurations which exhibit rigidity percolation in 0 directions (a), 1 direction (b), and 2 directions (c). The above gel configurations were generated with a packing fraction of ϕ = 0.39, which lies on the floppy side of the phase boundary reported by Zhang et al.¹⁹ The subfigure adjacent to (a) is a zoomed-in snapshot of gel particles corresponding to the section outlined in blue, which has a linear dimension of 25 particle diameters, where individual particles are visible.

Within VCTG, the presence of a pinch-point singularity in the stress correlation function is a qualitative indicator of rigidity: if a pinch-point in C_ijkl(q) is present, the system should exhibit rigidity at long length scales. We are interested in understanding how this property correlates with the nature of the rigid cluster. Fig. 3 shows results for the correlation function C_xxxx(q) ≡ 〈σ_xx(q)σ_xx(−q)〉 for gel configurations at ϕ = 0.39 with the largest rigid cluster percolating along 0, 1, and 2 directions, respectively. For each case, we average over all available configurations with the same degree of rigidity percolation. The maps (Fig. 3(a)–(c)) only show the existence of the pinch-point singularity in the case where the rigid cluster percolates in both directions. The size of the pinch-point is established by the asymptotic values of C_xxxx as |q| → 0 measured along the q_y = 0 and q_x = 0 directions. The presence/absence of a pinch point can therefore be quantified from the values of C_xxxx at the smallest q (Fig. 3(d)): a strong difference in relative value measured along the q_x = 0 and q_y = 0 cross-sections indicates that a pinch-point singularity is present. This singularity is clearly visible at length scales on the order of 50 particle diameters (see Fig. 3), which is ∼3 times smaller than L, the length of the simulation box. Based on this observation, we expect the overall results to hold also at larger system sizes.

Fig. 3 Analytical behavior of the pinch point singularity in C_ijkl(q) correlates with degree of rigidity percolation. Spectra of C_xxxx (normalized by the maximum value of C_xxxx) are ensemble averaged over configurations which exhibit rigidity percolation in 2 directions (a), 1 direction (b), and 0 directions (c). Cross sections of these spectra along q_y = 0 and q_x = 0 are shown in (d)–(f) in red and blue, respectively. The population of gels that did not percolate are shown to lack flat |q|-dependence in C_ijkl(q), indicating that these configurations are not rigid by the reckoning of VCTG.

The pinch-point becomes less obvious for the gels with rigid clusters percolating along one direction, see Fig. 3(b) and (e). For the gel systems where there is no percolating rigid cluster, Fig. 3(c), we observe no signal at q ≈ 0, see Fig. 3(f). In other words, the pinch-point is clearly absent. The plots Fig. 3(d)–(f) also show that only fully rigid gels feature, at low q, the flat q-profile of the stress correlations predicted by VCTG. This analysis demonstrates that the pebble game evaluation of the rigidity is fully consistent with the behavior of the stress–stress correlations and of the pinch-point singularity predicted by VCTG.

We also note that the correlation profiles in Fig. 3(d)–(f) indicate the presence of a length scale beyond which the rigid or floppy behavior of the gels is clearly distinguishable, and beyond which the VCTG predictions are satisfied for the fully rigid gels. Stress fluctuations at smaller length scales (larger q values) seem to be similar across floppy and rigid gels. Such a length scale was not detected by the pebble game analysis and may be a distinctive feature of gel structures. In granular solids, by contrast, VCTG applies all the way up to the grain scale.^22,23,44

4.2. Pressure and elastic modulus

The natural ensemble for VCTG is the generalized Edwards ensemble, where all configurations with the same global stress tensor are weighted equally.²² This would mean a constant-pressure ensemble in the case of isotropic compression. In order to better connect the predictions of the theory with properties of rigid clusters obtained instead at constant particle surface fraction, we have computed the distribution of pressure across the gel configurations produced in the simulations and investigated if and how the degree of percolation of rigid clusters can be correlated with the pressure values. Fig. 4(a) and (b) show that, regardless of the particle surface fraction, fully rigid gels whose rigid clusters percolate in both directions tend to have higher pressures. Partially floppy gels, for which a rigid cluster span the system only in one direction, tend to have lower pressures, and completely floppy gels correspond to zero pressure states. While one would expect that simply increasing the particle surface fraction may be responsible for the overall increasing pressure of rigid gels, the distinct correlations with the degree of percolation of the rigid cluster points instead to an interesting possible connection between pressure heterogeneity and rigidity percolation. The fact that these correlations may be specifically due to the rigidity (or the absence thereof) is confirmed in Fig. 4(c), showing the correlation between the pinch-point size and the degree of percolation of the rigid cluster for the less dense gels, connecting the pinch-point size to the distribution of pressures in the same sample at the same particle surface fraction.

Fig. 4 Rigidity of gel configurations at particular degree of rigidity percolation correlates with configuration pressure. Subfigures (a) and (b) depict the pressure distributions of the ϕ = 0.5 rigid gel dataset and the ϕ = 0.39 floppy dataset, with each distribution separated into groups of equivalent rigidity percolation. The distribution associated with 2-directional percolation consists of mostly higher-pressure configurations. Configurations which do not percolate display lower total pressure, and vice versa. This correlation persists across different locations of the phase diagram, even as the underlying probability distribution affecting the probability of generating a configuration which percolates fully is changing. Subfigure (c) depicts the size of the pinch-point singularity measured in the ϕ = 0.39 dataset as additional statistics are taken into account, quantifying the disappearance of the pinch point seen in non-rigid configurations. In subfigure (d), gel stiffness (K_2D) measured from fit is plotted against total configuration pressure. The square root of K_2D is measured in order to have the same units as Tr(σ). Circles are configurations from the ϕ = 0.5 dataset and diamonds are from the ϕ = 0.39 dataset.

4.2.1. Elastic moduli. As discussed above, VCTG theory provides a way to measure gel elastic moduli from the set of stress–stress correlations, C_ijkl(q). As the VCTG theoretical form of C_ijkl(q) is solely dependent upon Fourier angle, the stress correlations from each rigid gel configuration was averaged across configurations with the same particle surface fraction and with the same degree of rigidity percolation, then re-interpolated onto a polar grid in Fourier space. This allowed for C_ijkl(|q|, θ) to be averaged into angular bins for the purposes of fitting to the theoretical form of C_ijkl(q) (as shown in Fig. 1(d)–(f)). The nonlinear fitting package SymFit was used to simultaneously fit all 6 components of the stress–stress correlations, using an isotropic form of the VCTG elastic modulus tensor (presented here in Voigt notation), that appears as a coupling constant in the field theory:

With this form of Λ⁻¹, VCTG predicts the following forms for the stress correlations C_ijkl(q):²²

〈〈σ_xx(q)σ_xx(−q)〉〉 = 4K_2D sin⁴(θ),

〈〈σ_yy(q)σ_yy(−q)〉〉 = 4K_2D cos⁴(θ),

〈〈σ_xx(q)σ_yy(−q)〉〉 = 〈〈σ_xy(q)σ_xy(−q)〉〉 = 4K_2D sin²(θ)cos²(θ),

〈〈σ_xx(q)σ_xy(−q)〉〉 = 4K_2D sin³(θ)cos(θ)

〈〈σ_xy(q)σ_yy(−q)〉〉 = 4K_2D sin(θ)cos³(θ)

These expressions can also be obtained via symmetry arguments if the Airy stress correlations are known/measured.⁵¹ Although this form of elastic modulus tensor is dependent upon two Lamé coefficients λ and μ, in 2D the functions C_ijkl(q) depend only on one combined stiffness parameter, K_2D ≡ μ(λ + μ)/(λ + 2μ).²² Note that there is an overall scaling of K_2D with the shear modulus, μ. The fitting quality (r² = 0.981) confirms the strong agreement between the angular dependence of the measured stress correlations and the one predicted by theory.

Following this demonstration, K_2D was then measured for individual rigid gel configurations at different packing fractions. By performing the fitting process on individual configurations, it is then possible to examine the makeup of the ensemble-averaged prediction of K_2D. Although this type of fitting would normally be too noisy to perform on a single configuration, the process of averaging into angular bins applies a degree of smoothing which enables for reasonable fits to be carried out. We note that the value of K_2D is proportional to the magnitude of the correlation function, and hence determining a value for this constant relies upon averaging over the region of C_ijkl(q) where the magnitude of the function bears no dependence on |q| (i.e., where the measured function agrees with the theoretical prediction for the |q|-dependence of the stress correlations). Although the angular dependence of the correlation functions measured for floppy gels are indeed consistent with the angular dependence prescribed by VCTG, the altered |q|-dependence precludes them from having a well-defined magnitude to fit to VCTG theoretical forms.

The results of this process are displayed in Fig. 4(d) for all the rigid gels. The data show that the stiffness increases with pressure for rigid gels, and that there is a trend of the gel modulus increasing with pressure which becomes more clear at higher particle surface fraction. Overall the analysis performed confirms further the agreement with the VCTG theory, and how exploring the theory predictions can provide new insight into the mechanics of the gels and indicate interesting directions for future investigations. For example: some signatures of the possible coupling between stiffness and pressure indicated here can be seen in the confinement effect on the shear modulus, μ, detected in experiments on agarose gels.⁵² Combining the different results in Fig. 4 could allow one to understand how structural or pressure heterogeneities may change the gel elastic moduli, and more.

4.3. Debye-like screening

In this section, we go back to the length scale associated with the stress fluctuations which is present in both rigid and floppy gels, shown in Fig. 3(d)–(f). As discussed above, there is a characteristic length scale above which the VCTG predictions strictly apply to the gels, distinguishing clearly rigid from floppy configurations. Over distances smaller than this length scale, stress correlations behave similarly for rigid and floppy gels and the predictions of VCTG do not hold. These findings raise the question of a possible screening length for the theory. We recall that the VCTG framework corresponds to dielectric screening, due to bound charges as in electrostatics, but the charges are forces, and therefore, vectors. The polarization field P̂ represents the pairwise interaction forces between gel particles or grains where the two equal and opposite forces make up the force dipole. Screening of the type emerging from mobile charges could produce a screening length. As in electrostatics, the action of the VCTG field theory can be extended to include a length scale:where σ_ij = E_ij + P_ij, and ∇ × ∇ × E = 0, Λ_ijkl is the 4th rank tensor which gives rise to the emergent elastic modulus, and Γ_ij, g, g̃ are additional coupling constants of the theory. The screening length scale, just from dimension counting, is . This scalar length scale would manifest itself in stress–stress correlations in a complicated way depending on the coupling constants Λ and Γ. However, following ref. 32, we can identify ξ from the pressure-pressure correlations 〈P(q)P(−q)〉 in 2D, which are completely isotropic, leading to a change from q⁰ dependence to a q⁻² dependence around 1/q ≃ ξ. In Fig. 5 we show that this indeed happens for all 〈P(q)P(−q)〉 computed from the simulations for the three classes of ensembles: configurations with rigid clusters percolating in (i) 2 directions, (ii) 1 direction, and (iii) 0 directions. The single configuration at ϕ = 0.5 which did not percolate in either direction has been omitted due to lack of statistics. While only rigid gels can be seen to approach a flat |q|-dependence near q = 0, all of these correlation functions display a change in power-law behavior as predicted by the theory with screening. Given that the theory describes well the rigidity of the gels, one could expect the screening length to be sensitive to it. The crossover between the two behaviors, however, occurs at approximately the same characteristic length scale (estimated around 10–12 particle diameters) for both rigid and non rigid gels. In SCT, the screening length ξ can diverge with unbinding of defects in crystalline solids,³¹ and analogs of such transitions have been observed in the loss of rigidity in amorphous materials.³² The values of ξ obtained from the VCTG analysis are only of the order of a few particle diameters (much smaller than the system size) and do not appear to diverge with the loss of rigidity. Interestingly, this length scale is comparable to the size of the locally ordered regions seen in the zoomed-in snapshot in Fig. 2(a), suggesting possible connections between the VCTG screening length and the defects of the local order in the particle packing.⁵³

Fig. 5 Pressure–pressure correlations reveal rigidity-independent length scale. Subfigures (a) and (b) depict the value of P(q)P(−q) (normalized by the maximum value of C_xxxx) averaged over all angles for ϕ = 0.5, for gels with rigidity percolation in 2 directions and 1 direction, respectively. A flat power law is plotted near q = 0 in blue, and a power law proportional to q⁻² is plotted near qd = π in black. Subfigures (c)–(e) depict the same functions for ϕ = 0.39, for gels which percolate in 2 directions, 1 direction, and 0 directions, respectively.

Future studies with varying system sizes and gel models could help elucidate the origin of the screening and clarify its physical meaning in the context of VCTG theory.

5. Conclusion

We have tested quantitatively the predictions of a vector charge theory (VCTG) for amorphous solids with computer simulations of 2D model colloidal gels. The stress correlations in the model colloidal gels feature precisely the pinch-point singularity predicted by the theory as a result of the constraints imposed by the mechanical equilibrium to the internal forces. We have demonstrated that the size of the pinch-point singularity allows one to quantitatively distinguish rigid from floppy gels, and introduces a direct connection between the vector charge theory and the rigidity percolation framework previously introduced to describe colloidal gelation and rigidity transitions for a range of soft materials. The angular dependence of the stress correlations in Fourier space is also quantitatively predicted by the theory and allows for extraction of the emergent elastic modulus of the gels. The striking agreement of the results obtained through the numerical simulations of the gel model with the theory predictions suggest that the theoretical description of pre-stress may adequately capture the different sources of pre-stresses in systems like colloidal gels.

We used the analysis performed to investigate correlations between gel rigidity and distributions of internal pressures across gels obtained in similar conditions (same preparation protocol and particle density). The results suggest that, for the same densities, there is an interplay between the presence of a rigid spanning cluster and the heterogeneities of pressures in gel configurations, which could have interesting potential implications for the mechanical properties. The stress correlations computed in the model gels also reveal some deviations from the VCTG theory, indicating a length scale beyond which the theory applies and rigid/floppy states can be clearly distinguished. A Debye-like screening mechanism seems to be able to capture these deviations from the VCTG theory in terms of a screening length, suggesting that a screening of the type due to mobile charges could be limited to very short lengthscale (i.e., the grain size) in granular solids but extends over much larger lengthscales in self-assembled gels, pointing to a possible connection with geometric frustration in the gel structures. Future investigations of the screening length and of its microscopic origin across different gel models and self-assembly processes could provide more insight into these questions.

Author contributions

A. C. and H. A. V. performed all the analysis. F. D. R. and H. A. V. performed the numerical simulations of gels. A. C., H. A. V., X. M., E. D. G. and B. C. designed the research, discussed the interpretation and wrote the article.

Data availability

Data used to generate plots for this article are available at Zenodo Open Repository at https://zenodo.org/records/15058606.

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

We thank Jishnu N. Nampoothiri for valuable discussions. This work was supported by the National Science Foundation (NSF-DMR-2026842, NSF-DMR-2026825, NSF-DMR-2026834, and NSF CBET-2228681).

References

1. E. Del Gado, D. Fiocco, G. Foffi, S. Manley, V. Trappe and A. Zaccone, in Colloidal Gelation, John Wiley & Sons, Ltd, 2016, ch. 14, pp. 279–291.
2. T. F. Tadros, Interfacial Phenomena and Colloid Stability: Industrial Applications, Walter de Gruyter GmbH & Co KG, 2015.
3. G. Petekidis and N. J. Wagner, in Rheology of Colloidal Glasses and Gels, Cambridge University Press, 2021, pp. 173–226.
4. Y. Cao and R. Mezzenga, Nat. Food, 2020, 1, 106–118.
5. S. Aime, L. Ramos and L. Cipelletti, Proc. Natl. Acad. Sci. U. S. A., 2018, 115, 3587–3592.
6. Z. Filiberti, R. Piazza and S. Buzzaccaro, Phys. Rev. E, 2019, 100, 042607.
7. J. Colombo and E. Del Gado, J. Rheol., 2014, 58, 1089–1116.
8. N. Koumakis, E. Moghimi, R. Besseling, W. C. Poon, J. F. Brady and G. Petekidis, Soft Matter, 2015, 11, 4640–4648.
9. T. Gibaud, T. Divoux and S. Manneville, Statistical and Nonlinear Physics, Springer, 2022, pp. 313–336.
10. S. Alexander, Phys. Rep., 1998, 296, 65–236.
11. V. J. Anderson and H. N. Lekkerkerker, Nature, 2002, 416, 811–815.
12. C. P. Royall, M. A. Faers, S. L. Fussell and J. E. Hallett, J. Condens. Matter Phys., 2021, 33, 453002.
13. X. Mao, P. M. Goldbart, X. Xing and A. Zippelius, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2009, 80, 031140.
14. T. Divoux, V. Grenard and S. Manneville, Phys. Rev. Lett., 2013, 110, 018304.
15. M. Bouzid, J. Colombo, L. V. Barbosa and E. Del Gado, Nat. Commun., 2017, 8, 1–8.
16. R. P. Behringer and B. Chakraborty, Rep. Prog. Phys., 2018, 82, 012601.
17. S. Henkes and J. Schwarz, Statistical and Nonlinear Physics, Springer, 2022, pp. 427–448.
18. H. A. Vinutha and S. Sastry, Phys. Rev. E, 2019, 99, 012123.
19. S. Zhang, L. Zhang, M. Bouzid, D. Z. Rocklin, E. Del Gado and X. Mao, Phys. Rev. Lett., 2019, 123, 058001.
20. H. Tong, S. Sengupta and H. Tanaka, Nat. Commun., 2020, 11, 4863.
21. H. Dashti, A. A. Saberi, S. Rahbari and J. U. R. Kurths, Sci. Adv., 2023, 9, eadh5586.
22. J. N. Nampoothiri, M. D'Eon, K. Ramola, B. Chakraborty and S. Bhattacharjee, Phys. Rev. E, 2022, 106, 065004.
23. J. N. Nampoothiri, Y. Wang, K. Ramola, J. Zhang, S. Bhattacharjee and B. Chakraborty, Phys. Rev. Lett., 2020, 125, 118002.
24. M. Pretko, Phys. Rev. B, 2017, 96, 1–26.
25. K. A. Whitaker, Z. Varga, L. C. Hsiao, M. J. Solomon, J. W. Swan and E. M. Furst, Nat. Commun., 2019, 10, 2237.
26. B. Keshavarz, D. G. Rodrigues, J.-B. Champenois, M. G. Frith, J. Ilavsky, M. Geri, T. Divoux, G. H. McKinley and A. Poulesquen, Proc. Natl. Acad. Sci. U. S. A., 2021, 118, e2022339118.
27. L. C. Hsiao, R. S. Newman, S. C. Glotzer and M. J. Solomon, Proc. Natl. Acad. Sci. U. S. A., 2012, 109, 16029–16034.
28. M. Bantawa, B. Keshavarz, M. Geri, M. Bouzid, T. Divoux, G. H. McKinley and E. Del Gado, Nat. Phys., 2023, 19, 1178–1184.
29. S. M. Fenton, P. Padmanabhan, B. K. Ryu, T. T. Nguyen, R. N. Zia and M. E. Helgeson, Proc. Natl. Acad. Sci. U. S. A., 2023, 120, e2215922120.
30. H. A. Vinutha, F. D. Diaz Ruiz, X. Mao, B. Chakraborty and E. Del Gado, J. Chem. Phys., 2023, 158, 114104.
31. M. Pretko, Z. Zhai and L. Radzihovsky, Phys. Rev. B, 2019, 100, 134113.
32. N. S. Livne, A. Schiller and M. Moshe, Phys. Rev. E, 2023, 107, 055004.
33. M. Pretko and L. Radzihovsky, Phys. Rev. Lett., 2018, 120, 195301.
34. A. Gromov and P. Surówka, SciPost Phys., 2020, 8, 065.
35. S. Kobayashi and R. Tarumi, R. Soc. Open Sci., 2025, 12, 241568.
36. A. Kumar, M. Moshe, I. Procaccia and M. Singh, Phys. Rev. E, 2022, 106, 015001.
37. M. Moshe, I. Levin, H. Aharoni, R. Kupferman and E. Sharon, Proc. Natl. Acad. Sci. U. S. A., 2015, 112, 10873–10878.
38. M. Moshe, E. Esposito, S. Shankar, B. Bircan, I. Cohen, D. R. Nelson and M. J. Bowick, Phys. Rev. E, 2019, 99, 013002.
39. Y. Bar-Sinai, G. Librandi, K. Bertoldi and M. Moshe, Proc. Natl. Acad. Sci. U. S. A., 2020, 117, 10195–10202.
40. M. Moshe, E. Esposito, S. Shankar, B. Bircan, I. Cohen, D. R. Nelson and M. J. Bowick, Phys. Rev. Lett., 2019, 122, 048001.
41. A. Lematre, C. Mondal, M. Moshe, I. Procaccia, S. Roy and K. Screiber-Re'em, Phys. Rev. E, 2021, 104, 024904.
42. A. Kumar, M. Moshe, I. Procaccia and M. Singh, Phys. Rev. E, 2022, 106, 015001.
43. Y. Fu, Y. Jin, D. Pan and I. Procaccia, Phys. Rev. Lett., 2025, 134(17), 178201.
44. S. Chakraborty, J. N. Nampoothiri, S. Bhattacharjee, B. Chakraborty and K. Ramola, In Preparation.
45. L. Balents, Nature, 2016, 540, 534–535.
46. M. G. Noro and D. Frenkel, J. Chem. Phys., 2000, 113, 2941–2944.
47. J. Colombo and E. Del Gado, Soft Matter, 2014, 10, 4003–4015.
48. M. Bantawa, W. A. Fontaine-Seiler, P. D. Olmsted and E. Del Gado, J. Phys.: Condens. Matter, 2021, 33, 414001.
49. S. Plimpton, J. Comput. Phys., 1995, 117, 1–19.
50. D. J. Jacobs and B. Hendrickson, J. Comput. Phys., 1997, 137, 346–365.
51. A. Lematre, J. Chem. Phys., 2018, 149, 104107.
52. B. Mao, T. Divoux and P. Snabre, J. Rheol., 2016, 60, 473–489.
53. R. Maharana, D. Das, P. Chaudhuri and K. Ramola, Phys. Rev. E, 2024, 109, 044903.

---