Colloidal gelation with non-sticky particles

Colloidal gels are widely applied in industry due to their rheological character—no flow takes place below the yield stress. Such property enables gels to maintain uniform distribution in practical formulations; otherwise, solid components may quickly sediment without the support of gel matrix. Compared with pure gels of sticky colloids, therefore, the composites of gel and non-sticky inclusions are more commonly encountered in reality. Through numerical simulations, we investigate the gelation process in such binary composites. We find that the non-sticky particles not only confine gelation in the form of an effective volume fraction, but also introduce another lengthscale that competes with the size of growing clusters in gel. The ratio of two key lengthscales in general controls the two effects. Using different gel models, we verify such a scenario within a wide range of parameter space, suggesting a potential universality in all classes of colloidal composites.

Sticky colloidal particles diffuse and aggregate into clusters until forming a ramified, space-spanning network, i.e., colloidal gel 1,2 . Due to the load-bearing structure, colloidal gels behave as soft solids with finite yield stress beyond which flow occurs 3,4 . This rheological signature enables gels to be widely applied in industry, ranging from foodstuffs and personal care products to pharmaceutics and biotechnology [5][6][7] . Through experiments and simulations, particulate gels of monodisperse attractive (i.e., single-component) colloids have been extensively studied in various aspects (e.g., structure 2 , dynamics 5 , and rheology 8 ), while theories have been proposed to approach a priori prediction (such as ref. 9). By contrast, realistic gel-like materials, which usually contain multiple components, remain less probed. This is partially due to the lack of proper model systems in experiments, while the intrinsic complication of polydispersity also limits the progress in fundamental understanding. While recent attention on composite systems increases [10][11][12][13][14][15] , most studies still remain on the phenomenological level of ac hoc models.
Among the diverse array of multi-component systems, the combination of gel matrix and solid fillers is prototypical in practical applications. For example, polymeric nanocomposites have remarkable mechanical properties and are ubiquitous in sensors, civil engineering, and microbial applications 16 . Progress has been made in understanding such composites, which, in the absence of strong filler-matrix interactions, can be described by conventional continuum mechanics 17 . Because of a large gap between the constituent sizes, empirical approaches (such as Krieger-Dougherty law) have been found to describe the basic behavior by assuming a continuum background 18,19 .
However, the continuum assumption no longer holds if the characteristic sizes of each component are comparable 20 . This is the case when replacing the polymer matrix in polymeric nanocomposites with a colloidal gel network, where the typical lengthscales are all micron-sized. The interplay between these lengthscales generates novelty.
Though not yet fully understood, such biphasic mixtures receive increasing attention 21,22 , and recent work reports a unique flowswitched bistability 23 , which has never been observed in regular gels. These observations suggest the essential role of filler (or inclusion 17 ) particles. While researchers have recently focused on the gelled state of these composites 17,24 , the inclusion effect on gelation dynamics is still unclear.
Using numerical simulations, we aim to shed light on the understanding of colloidal gelation with non-sticky particle fillers. The system we investigate is composed of sticky colloids, which can form a percolating gel network on their own, and non-sticky (NS) particles, which are hard spheres. As the latter stick neither to the gel colloids nor to themselves, they do not participate in gelation directly. Then the intuitive assumption seems to view NS particles as confinement to the gel part by compressing the available volume. Through extensive exploration of the parameter space, we show that the interplay between comparable lengthscales also plays a vital role and, in some cases, dominates over the confining effect and greatly impedes the gelation process. The ratio of two key lengthscales, i.e., the characteristic size of gel and the NS-particle spacing, offers a robust measure for such interplay, which controls the cluster growth during gelation. We verify this scenario over a wide range of compositions and particle sizes in different types of colloidal gel.

Colloidal gelation
A variety of attractions lead to colloidal gelation in nature, such as van der Waals forces, depletion forces, and hydrophobic interactions 25 . In this work, we investigate gelation under strong attractions (U att ≫ k B T) within a wide range of volume fractions (0.03 ≤ ϕ ≤ 0.3). We consider two representative contact models. The first model refers to typical attractions which drive particles to aggregate with only radial forces, as shown in Fig. 1a (blue). Such conservative attractions are characterized by pairwise potentials and mostly apply to smooth particles 2 without tangential constraints, i.e., particles in contact can freely slide and rotate. By contrast, the second model constrains tangential pairwise motions (sliding, rolling, and twisting) with the presence of attraction, Fig. 1a (red). As suggested in ref. 26, we term the first contact as attraction (att) and the second model with tangential constraints as adhesion (adh).
The difference in contact interactions leads to different micromechanics, such as bending rigidity 27 and isostaticity 28 . Based on the Maxwell criteria 29 , mechanical stability for frictionless particles in 3D requires an average contact number N ≥ N c = 6. Following 30 , we consider the isostaticity condition as microstructural information in attractive and adhesive systems and determine gelation accordingly. In particular, we extract particles with N ≥ N c and check their connectivity (see the "Methods" section), Fig. 1b. The gelation point determined by this method has been verified, in both experiments 31 and simulations 24 , to agree well with that from macroscopic rheology. Therefore, we apply this gelation criterion throughout this work, with N att c = 6 for attractive gels and N adh c = 2 for adhesive gels 28 . Our simulations start from a random, homogeneous configuration and evolve following Langevin dynamics for up to 10 4 times of Brownian time τ B ≡ πηd 3 /2k B T (where η refers to fluid viscosity and d to colloid diameter). More simulation details can be found in the Methods section. For both gels, the time required for gelation t g decreases with the volume fraction ϕ in a power-law manner, Fig. 1c. Fitting of attractive gel data (blue) gives an exponent of −3.7, while that of adhesive gels (red) exhibits a lower exponent −2.1. Detailed fitting results are as follows: t adh g ≈ 0:011 × ϕ À2:1 : ð2Þ At the same volume fraction ϕ, it takes adhesive colloids less time to gel than the attractive ones. The exponents roughly agree with the values in other literature [32][33][34] , justifying our gel simulations as well as the gelation criterion we used.
To capture the structural evolution during gelation, we measure the static structure factor S(q) for both attractive (ϕ = 0.1) and adhesive (ϕ = 0.05) gels, Fig. 1d. S(q) is originally flat due to the homogeneous randomization for the initial configuration. As time increases (indicated by arrows in Fig. 1d), a peak at intermediate wavenumber q appears, grows, and shifts to a lower q. Thus, a characteristic lengthscale ξ ≡ 2π/q 0 (where q 0 refers to the peak wavenumber) increases during gelation, plausibly representing the cluster growth. Note that the low-q peak is absent at ϕ = 0.5 (Supplementary Note 1), indicating a homogeneous attractive glass (AG) state. The gel-to-glass transition explains the slight deviation from the power-law scaling of t g , Fig. 1c (blue).
While the two systems exhibit similar S(q) evolutions, the fractal dimensions d f inside clusters, which can be estimated from the slope of SðqÞ ∼ q Àd f , are different 35 . According to Fig. 1d (top), the clusters in attractive gel are rather compact (d f ≈ 3) at short range, suggesting gelation via the typical arrested-phase-separation route 2,25 . As q decreases, the fractal dimension d f drops to 2 (consistent with the value reported in ref. 7), indicating a relatively open structure at larger lengthscales. We attribute this to the emergent bending rigidity of bigger building blocks, e.g., tetrahedrons composed of attractive particles.
By contrast, the gel network in the adhesive gel is more open and ramified with a lower fractal dimension d f ≈ 1.8, as expected by diffusion-limited cluster-cluster aggregation (DLCA) 25 . Since the clusters in adhesive gels are looser than those in attractive gels, it is easier for adhesive colloids to percolate at the same volume fraction ϕ, i.e., lower gelation time t g . The above results show that the two models we used lead to two different types of colloidal gels.

Gelation with NS particles
In the presence of NS particles, sticky colloids can still diffuse and aggregate into a percolating gel network, such as  colloidal gels where ϕ directly controls gelation, the gelation in binary systems depends on the volume fractions of gel colloids ϕ g and NS particles ϕ NS . Since NS particles do not directly participate in the formation of gel network, we expect them to geometrically confine gelation and compress the free volume V free for sticky colloids. The reduction in V free then leads to an effective increase in the colloidal volume fraction. If simply consider V free by subtracting the NS particle volume V NS from the total volume V tot , we can then define an effective volume fraction ϕ eff for gel colloids as below: Though the definition above neglects the exclusion shell around NS particles 36 , it has been verified to well capture the gelation diagram with varying attractions 24 . This is probably because as clusters grow and become increasingly large and porous, the cluster-NS interpenetration 37 invalidates the application of exclusion shell. We first focus on the case of large NS particles with d NS = 8d, where d NS and d refer to the sizes of NS particles and gel colloids, respectively. According to Fig. 1b and Eq. (3), the addition of NS particles is expected to decrease the gelation time t g . For both models at high ϕ g (ϕ g ≥ 0.15 for attractive systems and ϕ g ≥ 0.07 for adhesive systems), t g decreases with the volume fraction of NS-particles ϕ NS monotonically, Fig. 2b. Plotting t g versus ϕ eff collapses these high-ϕ g data on a master curve, Fig. 2c, which converges with the gel data we have shown in Fig. 1c. In this way, regardless of contact models, the effective volume fraction ϕ eff seems to well characterize the gelation time t g of sticky-NS composites.
At low ϕ g , nevertheless, the decrease in t g becomes progressively slow as ϕ NS increases. In particular, for attractive systems at ϕ g = 0.07, the gelation time t g even increases with the addition of NS particles at high ϕ NS and becomes higher than that of the pure gel, Fig. 2b (top). This conflicts with our expectation of increasing ϕ eff . Furthermore, plotting t g versus ϕ eff shows deviation from the master curve, Fig. 2c. While the data collapse at high ϕ g justifies the definition of ϕ eff , this inconsistency indicates another physics, manifesting at low ϕ g , that delays the gelation with NS particles.

Lengthscale interplay and diagrams
As previously mentioned, one notable feature of gel-particle composites is the comparable lengthscales. Here we identify two key lengthscales from each constituent and attribute the abnormal deviation from ϕ eff prediction (Fig. 2c) to their interplay. The first lengthscale is the characteristic size ξ of the gel structure derived from structure factor S(q), which evolves over time (Fig. 1c). Since our focus is the gelation time t g , we measure the structure factor S(q) at t = t g (Supplementary Note 2) and extract a time-independent lengthscale ξ g ξðt g Þ. Such lengthscale in general represents the correlation length at gelation point, Fig. 3a (inset).
For both attractive and adhesive gels, ξ g decreases with volume fraction ϕ, Fig. 3a. Namely, the more concentrated a system is, the smaller clusters it requires to assemble into a percolating network. This is consistent with the gelation at the DLCA limit 38 . Power-law fittings on the two sets of data give similar exponents, while the lengthscale in adhesive gels ξ adh g is slightly lower than that in attractive gels ξ att g at the same ϕ. Fitting results are as below: ξ adh g ≈ 0:65 × ϕ À0:90 : ð5Þ Note that the above ξ g refers to the lengthscale in pure colloidal gels and ϕ to the colloid volume fraction. In binary systems, the large NS particles can easily distort the large-scale structure so that the colloid-colloid structure factor S(q) barely exhibits a resolvable peak at low q ( Supplementary Fig. 3a). We, therefore, assume that the ϕ eff scenario also applies to the ξ g scaling, by simply replacing ϕ with ϕ eff in Eqs. (4) and (5). That is, we use ξ g in an equivalent pure gel as a proxy for that in a binary composite. By comparing the void distribution 39 in pure gels and composites, we justify such assumption in Supplementary Note 3.
Apart from geometric confinement, NS particles also generate a flexible porous medium 40 , in which colloids diffuse and aggregate into a gel network. Such medium is in general characterized by the pore size 20 . Here we use the spacing between NS particles δ to represent the porosity, Fig. 3b (inset). In particular, we simulate a collection of only NS particles at different volume fractions ϕ NS and measure the average interstice δ from Voronoi cell volume (the Methods section). In the unit of d NS , the average spacing δ, diverging at ϕ NS = 0, decreases with ϕ NS as shown in Fig. 3b. Though we measure δ in pure NS-particles systems, such quantity in binary mixtures remains almost unchanged at the same ϕ NS (see gray and black scatters in Fig. 3b). Moreover, as gelation proceeds, δ barely varies over time ( Supplementary Fig. 5). Thus, δ is time-independent and scales with NS particles' absolute, rather than relative, volume fraction.
Given other parameters in a binary system fixed, both ξ g and δ can be a priori determined by ϕ eff and ϕ NS , respectively. Then their ratio γ ≡ ξ g /δ varies as a function of ϕ g , ϕ NS , and d NS /d. Since the ξ g -scaling is different in the attractive and adhesive gels, as shown in Eqs. (4) and (5), the lengthscale ratio γ also depends on the specific contact model.
We find that the t g deviation from the master curve (Fig. 2c) correlates with γ. To quantify the degree of deviation, we use the distance between the average of measured t g and the predicted t fit g from powerlaw fitting, defined as follows: where t fit g refers to the gelation time calculated by Eq. (1) or (2) but using ϕ eff instead of ϕ. With all the data shown in Fig. 2c, we map out the ϕ g -ϕ NS diagrams for both gel models at d NS = 8d, Fig. 3c. The quantified deviation dev[t g ] is represented by the colormap.
For attractive systems, while most data show small dev[t g ], we observe two regions that present visible deviations in the diagram, Fig. 3c (left). At high ϕ g and ϕ NS , the effective volume fraction ϕ eff is so high that the system falls in the AG regime with ϕ eff > 0.4 (gray). This is consistent with the deviation at high ϕ eff in Fig. 2c, where the measured t g falls below the power-law fitting (blue dashed line).
Significant deviation also occurs at low ϕ g and high ϕ NS (upper left corner of the diagram). By drawing the iso-γ lines, we find that higher γ leads to more prominent deviation dev[t g ]. Namely, when the characteristic size in gel ξ g far exceeds the spacing δ between NS particles, gelation is greatly hindered by their interplay, which dominates over the effect of ϕ eff .
We use the same method to calculate the ratio γ as well as the deviation dev[t g ] in adhesive systems and find that this scenario appears to still work, Fig. 3c (right). As the lengthscale ratio γ increases, the deviation becomes significant at low ϕ g and high ϕ NS . For both systems, visually, the iso-γ line of γ = 2 demarcates the regions with low and high deviations. This result supports our argument that, as an important factor, the lengthscale interplay primarily affects the gelation process in binary systems at high γ.
The role of lengthscale ratio γ is further verified by varying d NS from d to 12d in both gel models, Fig. 3d, e. At the same composition, the lengthscale ratio γ decreases with the NS particle size d NS . For larger NS particles of d NS = 12d, therefore, most data points fall on the master curve when plotting t g versus ϕ eff ( Supplementary Fig. 6), and t g deviation is greatly suppressed due to the small γ. As d NS decreases, deviation becomes increasingly significant. When the colloids and NS particles have comparable sizes (i.e., d NS = d), almost all data points deviate from the master curve ( Supplementary Fig. 6). Raw data of diagrams in Fig. 3d, e can be found in Supplementary Note 5.
Remarkably, the iso-γ line at γ = 2 demarcates the low-and highdeviation regions in all cases, Fig. 3c-e. The positive correlation between deviation and lengthscale ratio γ is obvious when plotting all the data together, Fig. 3f.
Regardless of the interaction (attractive or adhesive), composition (ϕ g and ϕ NS ), and particle size ratio d NS /d, the lengthscale ratio γ, as well as the effective volume fraction ϕ eff , seems to characterize the gelation process in binary mixtures well.

Growth of the largest cluster
The deviation from ϕ eff scenario indicates an additional hindering effect at high γ. Such hindering results from the frustration of cluster growth. In particular, we examine the evolutions of particle fraction in the largest cluster N lc /N. For each contact model, we compare systems with different compositions and d NS but the same ϕ eff (ϕ eff = 0.2 for attraction and ϕ eff = 0.1 for adhesion), Fig. 4. The value of lengthscale ratio γ is represented by the color. Compared with the pure gel at ϕ = ϕ eff (black), binary systems with small γ exhibit similar growth, while those with large γ show a delayed increase in N lc /N. Interestingly, the cluster morphology appears to be barely affected by NS particles regardless of γ (Supplementary Note 6). These results explicitly show how the lengthscale interplay, represented by the unified parameter γ, affects colloidal gelation with NS particles.

Behavior of NS particles
While colloids aggregate into a porous network as gelation proceeds, the dynamics of NS particles varies from system to system. At dilute ϕ g and small d NS , we expect NS particles to be able to diffuse even upon gelation since the pore size of the gel matrix is much larger than the NS particles. As the size of pores shrinks, the NS diffusion starts to be confined until finally 'locked' in the matrix cage as soon as a gel network is formed.
To probe the effect of ϕ g , ϕ NS , and d NS on NS dynamics, we measure the mean squared displacement (MSD) in various composite samples. For ease of comparison between colloids and NS, we normalize MSD by diffusion coefficient D diff for each species of the particles. We find that increasing ϕ g and d NS both lead to the dynamical arrest of NS, which seems to occur simultaneously with that of colloids, Fig. 5a (left and middle). This may correspond to the case where NS size exceeds the pore size of the gel matrix, which decreases with ϕ g as in Supplementary Fig. 4. We also notice that increasing ϕ NS slows down the NS dynamics, Fig. 5a (right), probably due to the increasingly crowding surroundings 41 .
Through radial distribution function (RDF), we also study the configuration of NS particles in binary composites. Without loss of generality, we use a specific composition of ϕ g = 0.1 and ϕ NS = 0.3 with four different d NS . While a small peak appears at the second-nearest neighbor for small d NS = d and 3d, such subtle spatial correlation does not present for larger d NS , Fig. 5b. We do not identify any sign of crystallization for all cases, and there is little variation in RDF before and after gelation. These results imply that the depletion effect 42,43 (and the consequent Casimir-like attraction 44 ) between NS particles is neglectable.

DISCUSSION
To better illustrate the effect of NS particles, here we consider two limits, Fig. 6. For infinitely-large NS particles (d NS → ∞), the colloids behave as a continuum which is geometrically confined between solidwall boundaries, i.e., the surfaces of NS particles. Within the colloidal phase, the real volume fraction is ϕ eff rather than ϕ g , and the diffusion, as well as aggregation, is purely mediated by the background solvent of viscosity η f . As the gelation time is proportional to the Brownian time τ B , we expect the scaling to be t g ∼ η f ðϕ eff Þ α , where α refers to an interaction-dependent exponent.
At another limit with d NS → 0, the NS particles form a continuum background in which sticky colloids are distributed, Fig. 6 (right). In such a case, confinement for colloids is absent so that the gelation time scales with the absolute volume fraction ϕ g rather than the effective one ϕ eff . The continuum background is essentially a hard-sphere suspension, whose viscosity η NS increases with the volume fraction of NS particles 45 . Though such suspension is not a simple Newtonian fluid of viscosity η NS in general, we may expect so because the situation is close to equilibrium 46 . In this sense, the gelation time then has the form t g ∼ η NS ðϕ g Þ α . For the smallest d NS = d we investigate, the background viscosity η NS increases with ϕ NS roughly in a Krieger-Dougherty manner 47 (Supplementary Note 7).
The above arguments can be generalized by using an effective background viscosity η eff and an effective colloidal volume fraction ϕ g eff (differing from ϕ eff in Eq. (3)) as follows: The values of η eff and ϕ g eff at the two limits are shown in Fig. 6. The collapsed gelation time t g in systems of d NS = 12d validates ϕ eff at large NS particles. For binary systems with d NS = d and the same ϕ g = 0.1, the identical structure factor S(q) (Supplementary Fig. 3) and void distribution ( Supplementary Fig. 4) suggest that the effective volume fraction ϕ g eff reduces to ϕ g at small d NS . Though this conflicts with the previous assumption (ξ g -ϕ eff scaling), using ϕ g instead of ϕ eff seems to make little difference in the iso-γ line (Supplementary Note 3).
As d NS decreases from infinity to zero, therefore, we expect transitions in both η eff (from η f to η NS ) and ϕ g eff (from ϕ eff to ϕ g ). At intermediate d NS , the interpenetration between NS particles and ramified clusters, which weakens the confinement effect and thereby decreases the effective volume fraction ϕ g eff , becomes possible. Meanwhile, the further aggregation of colloidal clusters with size comparable to the NS spacing (γ ∼ 1) requires the rearrangement of NS particles, which turns on the transition in background viscosity from η f   to η NS . In particular, pinning NS particles during gelation leads to diverging gelation time t g beyond γ = 2 (Supplementary Note 8). In this way, the t g deviation occurs as a result of both the decrease in effective volume fraction ϕ g eff and the increase of background viscosity η eff . While the above discussion merely considers d NS , the lengthscale ratio γ, including ϕ g , ϕ NS , and d NS , offers a generic, unified measure in all binary composites. The lengthscale competition then essentially represents a mechanism transition between the two limit cases as shown in Fig. 6. While using global quantities in the expression of γ, we note that local details may also play a secondary role. For example, the lengthscale ratio γ assumes uniform sizes for clusters and NS spacing, which, in practice, both have size distributions (Supplementary Note 4) and irregular morphologies (such as porosity and tortuosity 40,48 ). Meanwhile, as binary systems involve different interactions, the coupling of localized glassy dynamics may also interfere with the formation of gel network 49 . These factors, as well as the others not listed, may more or less affect the gelation process and cannot be fully captured by a sole quantity γ. This is consistent with the scattering of data points in Fig. 3f.
In summary, we use Langevin dynamics simulations to investigate colloidal gelation with non-sticky particulate fillers. Through extensive exploration in the parameter space, we find that the interplay between two lengthscales (ξ g and δ), represented by their ratio γ, matters in composite gelation. At γ < 2, the NS particles act as geometric confinement and effectively increase the colloidal concentration in the form of ϕ eff . As γ increases, gelation is progressively hindered as a result of both the decrease in effective concentration and the increasingly-viscous background. Our results not only shed light on industrial formulation and processing, but also open up a new scheme for the tunability in practical gel materials. Though precise prediction is still challenging due to the missing of microscopic details, we successfully capture the generic importance of lengthscale competition in multi-component systems. Our finding will inspire the fundamental understanding and may lead to the efficient development of colloidal composite materials.

Simulations
We perform Langevin dynamics simulations on LAMMPS 50 . Under thermostat at k B T, our system contains two species of spherical particles which differ in size and interaction. The first species consists of sticky particles with an average diameter d (bidispersed with 0.87d and 1.13d to prevent crystallization), while the second species consists of elastic spheres of diameter d NS . To ensure that both colloids and NS particles are diffusive within the relevant time range for the gelation process (≳0.1τ B ), we set the particle mass proportional to the particle size with the constant damping time m/3πηd ≪ τ B (see Supplementary Note 9 for the details). The NS-NS and NS-g interactions are simple elastic repulsions when overlapped, modeled by a modified Hertzian model with a high modulus Ed 3 ≫ k B T. To capture the interaction between sticky colloids, we use the Derjaguin-Muller-Toporov (DMT) contact model 51 with a sufficiently strong attraction U att = 20k B T. With the same modulus E, the overlap caused by cohesion is small (≈0.01d) at force balance, ensuring shortrange attraction. Tangential constraints on sliding, rolling, and twisting (Fig. 1a) are all modeled in a modified Coulomb manner with the same spring constant and friction coefficient μ. Respectively, we set μ = 0 for attraction without constraints and μ = 1 for adhesion with constraints on all three motions.
Our simulation occurs in a cubic box of side length L = 50d with periodic boundaries, which is sufficiently large for bulk condition (Supplementary Fig. 10). In the absence of colloid-colloid attraction, an initial configuration is generated through multiple relaxations. We first randomize NS particles and wait for them to relax for 100τ B , and then relax randomly-distributed, non-attractive colloids with pinned NS particles for another 100τ B . Upon equilibrium, we unpin NS particles and allow the bulk system to relax shortly for 10τ B . The initial state generated by such a pre-relaxing protocol exhibits no overlapping and no visible aggregation caused by depletion forces. We find little dependence on the pre-relaxing duration ( Supplementary Fig. 11), suggesting the robustness of our results.
Starting from a homogeneous random configuration, each system evolves up to 10 4 τ B . We run each simulation for at least three times, with the data points and error bars shown in this work representing the average and standard deviation, respectively. Visualization and part of data analysis are carried out using OVITO 52 .

Determining gelation time
A robust criterion for gelation is crucial to accurately determine the gelation time t g . The experimental convention views the liquid-to-solid transition, typically characterized by oscillatory rheology 53 , as the gelation point. Inspired by the Maxwell criteria for stability 29 , recent work, including both simulation 24 and experiment 31 , correlates the evolution of clusters of isostatic particles (contact number N ≥ N c ) with colloidal gelation and confirms the validity of such structural indicator by comparison with macroscopic rheology.
In this work, we define the gelation time t g as the time required for the isostaticity percolation. As Fig. 1b shows, we first extract all isostatic particles with N ≥ N c (Table. 1 in ref. 28) and then examine their connectivity. Gelation is determined if there exist clusters percolating through periodic boundaries in all three directions (x, y, and z). Recent work correlates rigidity percolation with gelation boundary 54 . For adhesive systems, our method (isostatic percolation with N adh c = 2) gives the same result as rigidity percolation, since each pair constitutes a minimal rigid cluster. Yet this may not hold for attractive contacts with no tangential constraints, where 3D rigidity analysis is challenging 31 (and beyond the scope of this work). Therefore, we consistently apply the isostaticity method for attractive systems with N att c = 6.
Voronoi analysis and NS-particle spacing To measure the average spacing between NS particles δ, we perform simulations of only NS particles at different volume fractions ϕ NS . Upon Brownian relaxation, Voronoi analysis is performed, and the spacing between each pair of particles δ i is estimated as follows: where V cell, i refers to the volume of the Voronoi cell of the i-th particle.
Here we assume isotropic distribution and regard each Voronoi polyhedron as an equivalent sphere. Then the average spacing δ = 〈δ i 〉.
The distribution of δ i can be found in Supplementary Fig. 5.

Data availability
The data generated in this study are provided in the Supplementary Information/Source Data file. Source data are provided with this paper.