Disentangling microstructural elements of shear thickening suspensions via computer simulations of a minimal model

Abstract

We use a minimal model for a dense suspension undergoing thickening and thinning to investigate microstructural changes in 2d simulations. Our simulations show that in steady flow the contact network contains distinct building blocks which are clearly signaled by sharp peaks in the radial distribution function, similar to what is observed in granular jamming. These structures deform during thinning. Non-Gaussian stress fluctuations that only emerge during thickening are associated to power law tails in the distribution of local contact forces, which tend to emerge when the flow-induced building blocks form large spanning assemblies. The subset of the contact network characterized by strong contact forces and connectivity large enough to be rigid or over-constrained is increasingly likely to percolate as the system starts to thicken, and to percolate over larger strain windows during thickening. The tendency of these structures to span the sample and to persist is dramatically reduced during thinning, where instead their deformation allows for a more homogeneous spatial redistribution of contact forces, significantly reducing the fluctuations of the macroscopic stress over time.

---

1 Introduction

The proliferation of frictional and complex surface contacts between particles has been recognized as an essential ingredient in the build up of shear stress in dense suspensions under shear.^1–11 Shear thickening of particle suspensions is also often accompanied by large-scale spatial inhomogeneities of the shear stress, under an imposed rate, which also features anomalous fluctuations in time, detected in experiments and simulations.^2,12–21 Large-scale spatio-temporal heterogeneity of the stress-bearing part of the contact network has been hypothesized to explain these observations. Further insight into the growth and the morphology of the frictional contact network during shear thickening, therefore, can shed new light on the nature of the shear thickening transition and on how stresses can be controlled during flow.

In dry granular assemblies undergoing jamming or shear jamming, several studies have highlighted and clarified the distinct roles that force and contact networks play, as minimal changes in the contact network correspond to large-scale reorganization of the force network,²² and over-constrained portions of the contact network have been shown to percolate in correspondence with the shear jamming transition.^23–28 Both the jamming and the shear jamming transitions have been connected to the scaling properties of shear thickening transitions in dense suspensions,^3,29 therefore suggesting that similar roles of force and contact networks may emerge also there. As a consequence, there has been significant interest in characterizing how the network of frictional contacts changes in size and topology during thickening^30–34 and, in particular, which subset of the frictional contact network is most responsible for the changes in the suspension rheology when it thickens.^21,34

To address the many open questions on the growth and morphology of contact networks in shear thickening suspensions, computational studies in 2d and 3d have analyzed subsets of the frictional contact network and their evolution during thickening. One interesting subset has been identified through groups of particles that all have at least k frictional contacts with members of the same group (k-cores). In 3d, the presence of cores with k ≤ 3 was found to display significant correlation with the increase in the average stress with increasing the shear rate.^31,34 Locally rigid clusters of particles with at least 4 frictional contacts in 3d (4-neighbor clusters) have been shown to undergo a percolation transition during thickening, following a finite-size-scaling ansatz similar to equilibrium critical phenomena.²¹ As 4 contacts in 3d provides enough constraints, in presence of friction, to guarantee isostaticity within a mean-field Maxwell criterion,³⁵ these clusters are expected to be locally rigid or overconstrained when the contacts are not lubricated. We note that the 4-neighbor clusters are different from the 4-cores, for example, in that not all of the particle contacts need to have at least 4 contacts as well, so they tend to identify the backbone of 4-cores. Close to the percolation transition of these clusters, their percolation/de-percolation during steady flow seems to coincide with giant fluctuations of the stress,²¹ potentially providing a microscopic explanation of several aspects of the shear thickening phenomenology reported in experiments.^{8,11,16,17,29}

In other recent studies, two-dimensional simulations have shown that fully shear thickened flow states at high stress always feature system spanning clusters that are rigid according to the pebble game algorithm, which identifies clusters of particles with enough constraints to satisfy the isostaticity condition.³³ It has also been shown that the growth of such clusters is always favored when large enough stresses are imposed, and close to jamming features several characteristics typical of critical phenomena.^36,37 Collectively, all these studies support the idea that flow at high rates induces some level of self-organization of particle contacts, which then, in turn, changes stress transmission and rheological responses. More investigation is however needed to fully characterize specific microstructures that control stress transmission and their statistical properties. Recent studies of activity-driven shear thickening in model active suspensions have also highlighted how system spanning networks of constraints may appear even in model active suspensions, akin to strongly far-from equilibrium, self-propelled living systems.³⁸ In sum, gaining further insight into the self-organization underpinning shear thickening would not only open the path to truly control and design the flow of a wide range of particulate dense suspensions, but would also bring new knowledge and understanding of non-equilibrium, strongly driven, many body physical systems in broad sense.

One of the first outstanding questions is how certain microscopic structures developed under flow are specific of shear thickening, as self-organization under flow in dense suspensions is a very broad phenomenon that does not only occur in those conditions. However, separating and disentangling microstructural motifs that emerge in different conditions is very difficult even in numerical simulations that have full access to microstructures, because of the microstructure complexity, especially in three dimensions.

In this study, we resort to a two-dimensional minimal model for dense suspensions, previously introduced by Grob, Heussinger, and Zippelius¹⁵ and also studied by Maiti and Heussinger,³⁹ to specifically identify how subsets of the frictional contact networks formed under shear flow may be substantially different during thickening and thinning, and how they may be related to giant stress fluctuations akin to those detected in experiments. The simplicity of the model chosen, and of the analysis of two-dimensional simulations, allow us to draw interesting general conclusions that may apply to a broad range of systems. First, our simulations show that in steady flow the contact networks of a thickening suspension contain distinct building blocks consisting of small groups of disks. These are clearly signaled by sharp peaks in the radial distribution function which has similarities to what is observed in granular jamming. These structures cannot remain stable and deform during thinning. Second, we show that giant, non-Gaussian stress fluctuations that only emerge during thickening are associated to power law tails in the distribution of local contact forces, which tend to emerge when the flow-induced building blocks form large spanning assemblies. The subset of the contact network characterized by strong contact forces (i.e., in the tail of the distribution) and connectivity large enough to be locally rigid or over-constrained (i.e., at least 3 contacts in 2d with frictional interactions) is increasingly likely to percolate as the system starts to thicken, and to percolate over larger strain windows during thickening. We call these structures 3-SFN. Their tendency to span the sample and to persist is dramatically reduced during thinning, where instead their deformation allows for a more homogeneous spatial redistribution of contact forces, significantly reducing the fluctuations of the macroscopic stress over time.

The fact that the model features both thickening and thinning regimes allows us to analyze microstructures that correspond to the same (or very similar) relative viscosity, for the same solid fraction and contact laws, while being respectively in the thickening vs. thinning regime. As a consequence, our work highlights that the relative viscosity (i.e. the average stress) in itself does not contain all the information about how the microstructure affects the rheology of a suspension. More specifically, we obtain new mechanistic insight into how the distribution of angles between contact bond vectors is distinctly different between the thickening and the thinning regimes. By performing simulations with a system size significantly larger than what utilized in previous literature,^2,5,15,32 we show that stress fluctuations developed under equivalent conditions (specific contact laws, particle volume fraction, and shear rate) are qualitatively different between small samples, which display shear stresses oscillating between two values,^5,15 and large enough system sizes, which allow us instead to sample, throughout a comparable strain window, a much broader distribution, similar to experimental observations.^16,17,19 Small system sizes also do not allow for the statistics required to resolve the correlation of the percolation of the 3-SFN clusters with the rheology, which instead can be captured in sufficiently large samples. The 3-SFN, moreover, are a new and complementary diagnostic for characterizing the microstructure, which points to the microscopic origin of the stress fluctuations detected in experiments,^16–18 akin to the k-neighbor clusters introduced in the 3d numerical studies.²¹ While supporting the picture of a rigidity percolation of k-neighbor clusters, we also discuss the connection between 3-SFN and the rigid clusters identified by the pebble game algorithm in 2d simulations.^33,36

The paper is organized as follows: in Section 2 we describe the model and the numerical simulations, and Sections 3–6 discuss the rheological properties measured and the analysis of the structure and contact networks, whereas Section 8 contains some concluding remarks and an outlook.

2 Methods

We consider a two-dimensional (2D) assembly of non-Brownian soft disks sheared at a constant shear rate [small gamma, Greek, dot above]. The disks are immersed in an implicit solvent of viscosity η_sol, from which they experience a drag force proportional to their velocity relative to the imposed flow. To avoid crystallization, we consider a mixture of two sizes, big (B) and small (S) with diameters d_B and d_S and a diameter ratio d_B/d_S = 1.4. The packing fraction ϕ of the disks is fixed at 0.79, which is significantly below the isotropic jamming density of ϕ_J ≈ 0.845.⁴⁰

The disks interact via contact forces described by the Cundall–Strack model⁴¹ as shown in eqn (1), where contact forces between disks i and j are calculated using a combination of viscous and elastic terms along both the normal and tangential directions.

F⃑^ij_c = (κ_nδ − m_effζ_nvⁿ_ij)n̂_ij − (k_tΔs + m_effζ_tv^t_ij)t̂_ij (1)

Here the unit vector n̂_ij points from particle i to j and vⁿ_ij is the component of the relative velocity along that direction. v^t_ij is found by calculating the relative velocity of the disks v⃑_ij = v⃑_i − v⃑_j and taking the component along the direction t̂_ij perpendicular to n̂_ij. κ_n and κ_t are spring constants acting along the normal (n) and tangential (t) directions respectively while ζ_n and ζ_t are the related damping coefficients. The effective mass is given by m_eff = m_im_j/(m_i + m_j). The normal overlap is calculated simply as δ = R_ij − r_ij where R_ij is the sum of the radii equal to (d_i + d_j)/2 and r_ij is the inter-particle center of mass separation. The frictional (tangential) force is a function of how far disks have slid past each other since establishing contact. This requires integrating the tangential velocity from the time of contact τ₀:

The maximum value of the tangential force F⃑^ij_t is given by the Coulomb criterion F^ij_t ≤ μF^ij_n, where F^ij_n is the normal force, and μ is the friction coefficient. We choose μ = 1.0 as in other studies.^6,39

The translational and rotational degrees of freedom of the disks are described by the following equations of motion:

Here v⃑_i,f is the velocity of disk i (v⃑_i) relative to the solvent flow at its location (v⃑_f), I_i, is the moment of inertia about its center of mass, and F⃑^ij_c is the contact force due to particle j whose tangential component is given by F^ij_c,t. The second term in eqn (3) represents a drag force which is proportional to the difference between the particle velocity v⃑_i and the flow velocity v⃑_f(r_i) = ê_x[small gamma, Greek, dot above]y_i. The drag coefficient ζ_sol represents the viscosity of the surrounding fluid, ζ_sol ∼ η_sol. The angular velocities are damped by the drag coefficient .

We note that the 2D nature of the model, the oversimplified description of the solvent, especially the lack of lubrication forces, and the lack of explicit boundaries in the simulations, limit the possibilities of direct and quantitative comparison with experiments. Nevertheless, as discussed in the following, the study performed did allow us to develop new explanations and to provide new insights potentially relevant to experiments.

In our simulations, lengths are reported as multiples of the small disk diameters d_S while the energy unit ε is equal to the elastic energy of two small disks undergoing full overlap (κ_nd_S²/2). All disks are given the same mass m. Based on these fundamental mass (m), energy (ε) and length (d_S) units, the time unit is calculated as , which is approximately equal to the period disks oscillate with when pushed together. The input parameters in our simulations are κ_n = κ_t = 2, , and ζ_t = ζ_n/2. The solvent damping coefficients (ζ_sol = 0.1m/τ for translations and for rotations) are chosen such that the inertial forces are not important for the small shear rates considered here. The inertial number (where P is the system pressure and ρ is the particle density) is generally smaller than O(10⁻³).

As discussed by Mari et al.,² shear thickening inherently requires a time scale in the system that competes with the one corresponding to the inverse shear rate [small gamma, Greek, dot above]⁻¹. For our system, this can be provided by the time scale of the oscillations that the disks experience along the normal direction after establishing overlap. Comparable time scales can be computed from the tangential spring constant, the two damping constants ζ_n and ζ_t, and the intrinsic solvent viscosity η_sol. The absence of other microscopic time scales implies that the softness of the particles must become important in the limit of large shear rate, which leads to a transition from shear thickening to shear thinning as discussed in the following and as also observed in a previous study of the same model.³⁹

From the disks relative positions and forces, we compute the stress tensor

where r⃑_ij is the vector from the center of disk i to the point of contact with disk j, f⃑_ij is the contact force, and A is the area of the system.^24,42 Once the stress tensor is obtained, we measure the shear stress as the off-diagonal component [small sigma, Greek, circumflex]_xy, referred to subsequently as σ.

Data is collected for a system size of 20 000 particles. A constant shear rate is applied to each configuration until the system reaches a steady state. For each shear rate we use a statistically independent initial condition. Once steady state is reached, an additional Δγ = 10 is applied during which the system is sampled at even strain intervals of Δγ = 0.001 for a total of 10 000 configurations. All simulations were performed in LAMMPS⁴³ and with Lees–Edwards periodic boundary conditions.

3 Rheology and fluctuations

Fig. 1 shows the flow curve computed from the simulations once the steady state flow is reached at different shear rates. The range of rates we explore encompasses several orders of magnitude during which the system exhibits quasi-Newtonian (QN), thickening, and thinning flow respectively as the shear rate is increased. Although the shear stress in the system increases monotonically with shear rate, thinning is made possible by the softness of the disks which experience larger overlaps at faster rates. This model choice allows us to directly compare the structures of the same initial microscopic configuration undergoing thickening and thinning. Overall the rheology we report here is both qualitatively and quantitatively in good agreement with the results reported by Maiti et al.,³⁹ which used the same interaction model for similar densities. The shape of the flow curve does not change significantly by changing the system size, although we note that both the onset of thickening and the onset of thinning seem to move to slightly lower shear rates for smaller samples (see SI, Fig. S1).

Fig. 1 Relative viscosity versus shear rate for a system size of 20 000 particles (black circles). The corresponding stress fluctuations in steady state given as the standard deviation of the stress normalized by the mean (red squares). The flow curve exhibits quasi-Newtonian (QN) flow at the lowest rates [small gamma, Greek, dot above] ≲ 6 × 10⁻⁷τ⁻¹. Above this regime the system thickens until [small gamma, Greek, dot above] ≈ 6 × 10⁻⁶τ⁻¹ where thinning is observed due to particle overlaps. The stress fluctuations increase during the thickening but begin to decrease before the start of the subsequent thinning regime.

Together with the relative viscosity, which is obtained from the average shear stress, we plot also the fluctuations of the shear stress, indicating that thickening is characterized by a significant increase in the shear stress fluctuations, whereas the onset of thinning features their significant decrease. The dashed line indicates the transition from QN flow and thickening (left of the dashed line) to thinning (right of the dashed line). The location of the peak and the width of the thickening and thinning regimes depend on the friction coefficient, the viscosity of the solvent, the Young moduli of the particles, and the shear rate. Increasing the stiffness of the spring used to describe the particle contacts in the model shifts the onset of thinning to higher strain rates,¹⁵ without changing significantly the overall behavior. In the following, shear rates are reported as a fraction of the peak shear rate (the dashed line) at which the system begins to thin [small gamma, Greek, dot above]′ = [small gamma, Greek, dot above]/[small gamma, Greek, dot above]_P ([small gamma, Greek, dot above]_P = 6 × 10⁻⁶ τ⁻¹), for the specific set of parameters used here (see also Supplementary Information for a comparison of two different system sizes).

Fig. 2(a) displays the probability density distributions of the shear stress, obtained from the steady state flow, during the different flow regimes. Each data point indicates the value of σ at a given time (or strain) during steady state flow, in the different regimes of the flow curve. During thickening the distribution features a pronounced high stress tail (red squares), and shrinks back dramatically once the system transitions from thickening to thinning. The question of whether and how the large stresses in thickening stem from microscopic changes can be addressed in the simulations: the behavior of the stress fluctuations appears to originate from the fluctuations in the contact forces, as shown by the change in the contact forces distributions going into thickening and thinning flows in Fig. 2(b). Each data point here corresponds to a specific contact force present across all microscopic configurations sampled during the steady flow, computed over the same time series used in Fig. 2(a). During thickening, we detect the concomitant growth of fluctuations of microscopic quantities such as the average contact number, which also decrease during thinning (not shown). We note that the comparatively large system size used here allowed us to access the full distributions of shear stress and microscopic contact forces. Other computational studies with smaller system sizes showed instead bimodal distributions,^2,15,33,45 which we also recover in small samples (see SI, Fig. S2 and S3). Overall, these findings offered motivation to investigate the presence of specific structural motifs that could be associated to the changes in the distributions of microscopic contact forces, and in particular to the development of anomalously large contact forces, as explained in the following.

Fig. 2 The steady state distributions for the system-wide shear stress σ (a) and individual contact forces f (b). Shear rates are taken from the QN ([small gamma, Greek, dot above]′ = 5.0 × 10⁻²), thickening ([small gamma, Greek, dot above]′ = 5.0 × 10⁻¹), and thinning ([small gamma, Greek, dot above]′ = 1.7 × 10⁻¹) regimes. 〈σ〉 and f_peak are the average value of the system shear stress and the peak (most probable) value of the force distribution respectively. We find that the thickening flow is associated with an increase in fluctuations, in agreement with both experimental^12,13,44 and computational^2,15 studies.

Notably, we observe that the stress fluctuations (relative to the mean) stop increasing before the thickening is over (see also Fig. 1). As well, the thinning regime coincides with a significant reduction in the size of the fluctuations. Frictional contact network percolation is a mechanism which could explain the behavior of stress fluctuations during thickening observed in our model. Networks of frictional disks that percolate across the system are self supporting (due to the periodic boundary conditions) and should be able to carry larger loads than non-percolating networks. If frictional network percolation were a driver of thickening, then thickening could be expected to occur as long as these networks could increase the amount of time they percolated at higher rate. However, fluctuations in the shear stress should be maximized when frictional contact networks alternate the most between percolating and not-percolating.

4 Local structure evolution

We start by tracking how particle microstructure evolves with shear rate by calculating the pair correlation function g(r), which describes the number density of particles in the system at a distance r from a reference particle, and is defined in Doi⁴⁶ aswhere r⃑_mn = r⃑_m − r⃑_n is the relative position of particles m and n and N is the total number of particles. One can compute g(r) for each of the two particle species separately or together, and overall find similar features. Fig. 3 shows g(r) computed for the large diameter particles with an un-sheared state as a reference. When steady flow has developed at a shear rate in the QN regime, we note instead sharp peaks in g(r) that are not present in the un-sheared system. These peaks correspond specifically to local configurations of groups of three or four large particles or large and small particles as drawn in the figure. We note that these local structures are suggestive of building blocks of chain like structures, capable of growing to rapidly span the whole sample and efficient at transmitting stresses. We therefore hypothesize that they are building blocks of larger assemblies underpinning force chains, similar to what happens in jamming of dry grains under shear,²⁴ but here in steady flow, and therefore preferentially oriented in the compression direction. Interestingly, these structures emerge immediately under flow, i.e. even when the flow is QN, and during thickening there is little change in the g(r) profile – we find only a slight decrease in the sharp peaks as the shear rate is increased (Fig. 3). However, the comparison in Fig. 3 of the same g(r) across the different flow regimes shows that with the onset of the thinning regime the peaks significantly decrease and broaden, suggesting that major changes develop in these structural building blocks.

Fig. 3 Evolution of g(r) from an un-sheared state (black) to QN (orange), thickening (red), and thinning (blue) flow. For an un-sheared system, g(r) shows broad peaks at distances corresponding to triplet chains. When shear is applied, sharp peaks corresponding to triplet chains and quadrilaterals emerge. Little change in the microstructure is seen between the QN and thickening flow while significant change occurs when the system begins to thin. Shear rates are the same as in Fig. 2.

To understand how the local structures identified by g(r) differ during thickening and thinning flows, we have quantified and tracked their changes. In particular, if large-diameter disks i, j, k are connected such that i touches j and j touches k, with r⃑_ij = r⃑_j − r⃑_i, we can characterize the shape of the i, j, k triplet with the angle

θ = cos⁻¹(r̂_ij·r̂_jk) (7)

As shown in Fig. 4, θ = 0° corresponds to a straight chain whereas θ = 180° corresponds to an equilateral triangle. Tracking the probability density distribution of θ, we find that, similar to g(r), there is little change during thickening in the shape of the structural elements already emerged under shear, indicating that those structural elements, and their assembly, remain stable. However, the peaks of the distribution of θ evolve as thinning kicks in, indicating how the shape of these structures may be changing: the significant broadening of the peaks is consistent with disk triangles being compressed and with the deformation of initially straight triplets of large particles. If straight triplets are building blocks for force chains that carry stresses during thickening, the thinning must therefore take place through deformation and breakup of those stress bearing structures.

Fig. 4 The probability density distribution of θ at shear rates from the different flow regimes. Similar to the behavior of g(r), the θ peaks decrease and broaden significantly during the thinning regime. This change in the peaks is consist with chains deforming or triangles compressing, as shown by the enlarged picture of the θ = 120° peak in the inset. Shear rates are the same as in Fig. 2 and 3.

The emerging picture is that structural elements formed under flow grow into percolating structures that can bear and transmit large stresses during thickening, as they allow for large contact forces to develop. The structural disorder and the pronounced heterogeneities in these percolating structures probably allow for the accumulation of large contact forces on a relatively small fraction of all the contacts, leading to the long-tail distributions shown in Fig. 2(b) and hence (a). While we have evidence of distinctive local structures formed under flow and of their changes between shear thickening and thinning flow, quantities like g(r) are not adequate to capture their organization into larger scale structures, potentially spanning and underpinning force chains. In spite of several attempts, identifying these larger structure by starting from the local structural elements and following their evolution remains very difficult and inconclusive, due to the complexity of the microstructural evolution even in this relatively simple two dimensional case. We have therefore focused on identifying the large scale assemblies that control the rheological changes signaled by the flow curve by directly using the information on the large contact forces, as explained in the following.

5 Strong force network percolation

To further test the interpretation of our results, and identify the stress-bearing structures that control the rheology during thickening, we now combine the information on the tail of the contact force distributions contained in Fig. 2(b) with the idea that specific subsets of the contact network may be more relevant to the rheology, as suggested for example in Goyal et al.²¹ First, using the information contained in Fig. 2(a) and (b), we have studied groups of disks whose contacts bear forces of magnitude f which satisfies the condition

f ≥ kf_peak (8)

where f_peak is the peak of the contact force distribution for that shear rate. f_peak represents the most probable force at a given shear rate, and the parameter k can be used to adjust the force threshold in the condition above, with higher values allowing us to identify subsets of the contact network that include contact forces of higher magnitude. Here and in the following we call such a subset the strong force network (SFN), which is inspired by the approach proposed in Bi et al.²⁴ for dry granular materials under shear. These SFN's are networks of disks that are each over-constrained by large forces. In our analysis we have tested all results with several values for k, ranging from 0.75–1.25, finding that they do not depend qualitatively on the value of k. We have then considered specifically the SFN made up of particles with at least 3 contacts above the threshold defined by eqn (8). Combining the information on constraint geometry (i.e. number of contacts) with the contact forces (making sure we focus on large contact forces that emerge during thickening) appears to be central to understanding the changes in the rheological behavior across the different flow regimes.

The choice to analyze 3 contact SFN's (3-SFN) is significant for a couple reasons. First, setting a minimum contact number of 3 ensures that the disks are locally over-constrained by strong forces (in the sense of eqn (8)). Secondly, if we do not filter contacts by their force strength, we find that large networks of particles with 3 or fewer contacts always exist and percolate across the system. Meanwhile, networks of particles with 4 or more contacts never become large enough to percolate.

By analyzing the particle trajectories in steady flow, we can measure how the 3-SFN continuously percolate and break down under flow. Percolation here is defined as occurring when a specific 3-SFN spans both directions of the system. During thickening the 3-SFN experience dramatic fluctuations in their size and do not always percolate, as shown in the examples of Fig. 5. To quantify how often the 3-SFN do percolate, we calculate a percolation probability – the fraction of configurations, over the whole time series, in which there is a percolating 3-SFN – and find that during thickening this probability strongly increases with the shear rate. As seen in Fig. 6, the percolation probability grows with shear rate before peaking at the end of the thickening, to then decrease during thinning. We note that the coupling between the percolation of 3-SFN and the shear thickening/thinning becomes much less clear in small samples, where the data show instead a percolation probability very similar (and relatively large) across nearly the whole range of shear rates, without changing much between thinning and thickening. These data suggest that, if a percolation transition of 3-SFN underpins the shear thickening regime, only the large system size would provide enough resolution to accurately characterize it. While the full study of the percolation transition itself is beyond the scope of this paper, we have computed the cluster size distribution of 3-SFN from the time series in the thickening regime and found that the data for the large system size are consistent with proximity to a percolation transition, because of a power-law tail that clearly extends over more than one decade, with an exponent similar to the one predicted for random percolation in 2d⁴⁷ (see SI, Fig. S4). The data show the coexistence of several small and finite clusters with percolating ones during thickening, which is a hallmark of a percolation-type of growth.

Fig. 5 Filtering the networks of connected particles by the strength of the contact forces reveals heterogeneous percolation behavior. These strong force networks (SFN's) consist of particles with at least m contacts with forces f > f_peak where f_peak is the second peak (first non-zero peak) of the force distribution. (a) and (b) show the largest cluster with m = 3 at different times during the same experiment for [small gamma, Greek, dot above]′ = 5 × 10⁻¹. All other clusters are hidden for clarity. Notably, the system shear stress in (a) is smaller than the shear stress in (b).

Fig. 6 Strong force network percolation probability versus shear rate for different system sizes, N. The probability indicates the fraction of measured configurations in which a percolating network was found. For the system size of 20 000 particles, probabilities were calculated using 10 000 configurations evenly spaced in strain during which the system was subject to a strain of Δγ = 10.

The data for the large systems in Fig. 6 clearly indicate that percolating 3-SFN are predominant during thickening and become rare in the thinning flow. The deformation of the local particle structures that we discussed in the previous section are likely responsible for interrupting, during thinning, the percolation of the structures that were instead able to efficiently transmit stresses during thickening, leading to a net increase of the viscosity of the suspension with increasing shear rate. The findings just described provide further insight when we measure the amount of strain through which a particle remains in a percolating 3-SFN, which also grows during thickening. At any given time or strain, any given disk can either belong or not to a percolating 3-SFN. This status is represented through a variable C_i(γ₀), which corresponds to the affiliation of particle i at strain γ to a percolating 3-SFN and is either 0 (belongs) or 1 (does not belong). We compute the time autocorrelation function

〈C_i(γ₀) C_i(γ₀ + δγ)〉, (9)

which, averaged over all particles and strain windows of the same size at a given shear rate, is plotted in Fig. 7(a) as a function of strain. The results show how increasing the shear rate during thickening causes disks to stay in a percolating 3-SFN for longer strains. Conversely, increasing the rate during thinning decreases the strain window spent in a 3-SFN. The data therefore suggest that percolating 3-SFN become longer and longer lived during thickening and provide further evidence that they are a specific ingredient of the thickening flow. We can use these autocorrelation functions to estimate the characteristic strain δγ_D over which a percolating 3-SFN persists on average, as the strain by which the average autocorrelation falls below ≈1/e. The error bars in the plot refer to the autocorrelation falling below 1/e ± 0.05. δγ_D is shown as a function of the shear rate in Fig. 7(b) to directly relate its rate dependence to the different regimes identified in the flow curve (Fig. 1) and discussed so far. The data clearly confirm that, during the thickening, not only do the 3-SFN's percolate more often, but also disks tend to stay in percolating networks over longer amounts of strain. During the thinning, instead, δγ_D decreases with increasing the rate, confirming that the thinning flow, where the building units of stress-bearing structures become deformed, correspond to predominantly non-spanning 3-SFN. These results start to clarify that similar contact network structures could be found in both thickening and thinning, however their persistence to growing strains can be dramatically different due to the ability of their elementary building blocks to accumulate and transmit stresses, resisting the flow.

Fig. 7 As the 3-SFN's begin to percolate more frequently particles tend to stay in a percolating network longer. (a) The ensemble averaged autocorrelation of a particle's affiliation with a percolating network as a function of a change in strain δγ. (b) The de-correlation strain (the amount of strain it takes for the auto-correlation to fall below ≈1/e) grows during the thickening regime and decreases during the thinning, consistent with percolation and de-percolation transitions. The error bars in the figure have been obtained considering the correlation to fall below 1/e ± 0.05.

We would like, at this point, to gain further insight into how the large scale assembly of the microstructural units induced by the flow may grow during thickening and control the fluctuations in the shear stress.

6 3-SFN percolation and shear stress fluctuations

The size of the 3-SFN bears significant correlations with the fluctuations of the shear stress in time, during the thickening regime. This correlation remains but grows notably weaker during the QN and thinning flow. Fig. 8(a) shows the distribution of sizes of 3-SFN obtained from microscopic configurations for the steady state flow at a shear rate in the thickening regime ([small gamma, Greek, dot above]′ = 5 × 10⁻¹). The number n_s of 3-SFN that have a number of disks s is divided by the total number of particles in the system as outlined in Stauffer.⁴⁷ This is plotted against the size s, indicating that large scale assemblies of contacts that carry large forces are present and their maximum size grows with the rate. We then split the particle configurations, obtained for the same shear rate, based on whether the shear stress in the system is above or below the peak of the stress distribution σ_peak in Fig. 2(a). Fig. 8(a) reveals that an increase in the frequency of the largest 3-SFN sizes clearly occurs in configurations where the shear stress exceeds the mean stress at that shear rate. The peak around 10⁴ particles for the σ > σ_peak distribution corresponds to 3-SFN that encompass most of the particles and are system-spanning. The difference in the shape of the two distribution seems compatible with large percolating clusters for σ > σ_peak and non-percolating finite clusters for σ < σ_peak.

Fig. 8 (a) The cluster size distribution stress states below and above the peak value of Fig. 2a during thickening, showing how larger clusters occur at higher stresses. The peak around 10⁴ for the σ > σ_peak distribution corresponds to 3-SFN that encompass most of the particles and are system-spanning. (b) The stress distribution at times when the 3-SFN does and does not percolate. The stress values are rescaled by the ensemble average of the stress (calculated over all time). The stress is higher when the 3-SFN percolates. Data for both plots taken during the thickening ([small gamma, Greek, dot above]′ = 5 × 10⁻¹).

More evidence that the growth of 3-SFN is connected to the total shear stress is given in Fig. 8(b). Here the data for the distribution of the total shear stress values in steady state of Fig. 2(a) is split based on whether on not a percolating 3-SFN can be found in that configuration. We find that the presence of a percolating 3-SFN significantly shifts this distribution towards larger stresses. The results of Fig. 8(a) and (b) suggest that the 3-SFN may be predictors of whether the shear stress in a given configuration is high or low compared to its most probable value (σ_peak). Fig. 9(a) explores this for a shear rate in the thickening regime. The status of whether the system at a given strain contains a percolating 3-SFN (or a percolating rigid network as discussed below) is plotted alongside the total shear stress. To quantify how predictive changes in the percolation of 3-SFN is of changes in the shear stress, we compute the Pearson correlation coefficient between the two:

Fig. 9 (a) Percolation status of the largest 3-SFN and pebble game cluster alongside the stress in the system. Shear rate is from the thickening regime ([small gamma, Greek, dot above]′ = 5 × 10⁻¹). For higher values of the stress the 3-SFN percolates while dips in the stress are associated with de-percolation. (b) The correlation between the percolation status of the 3-SFN and the stress versus shear rate (black data). These correlation values can be compared to the correlation between the stress with a random time series (blue data). The thickening is associated with a sustained high correlation between 3-SFN percolation and shear stress while the QN and thinning flow are associated with increasing and decreasing correlation respectively.

Here σ(γ) is the total shear stress at strain γ and S(γ) is the percolation status of the 3-SFN being considered, with 0 indicating no percolation and 1 indicating percolation. The correlation coefficient is then plotted as a function of shear rate in Fig. 9(b) which shows that during the thickening regime there is reasonably good correlation (going up to ≈0.8) between the percolation status of 3-SFN and the shear stress in the system. For reference, we also show in the plot the Pearson coefficient between the same time series of the shear stress and a randomly fluctuating variable.

While the data provided here do not analyze all possible networks or possible subsets of the contact network, combined with the microstructural analysis discussed above they support the idea that 3-SFN and their percolation are distinctive ingredients of thickening in our model suspension. The 3-SFN identify the subset of the contact network that includes large contact forces developed upon thickening because their building blocks efficiently transmit stresses across the sample when they percolate. Finally, their percolation/de-percolation creates/interrupts the build-up of large contact forces and therefore underpins the large fluctuations of the total shear stress. Once the structural elements that compose them start to deform, efficient stress transmission is compromised and percolating 3-SFN can only survive over smaller and smaller strain windows. The deformation of the building blocks, as well as the more continuous percolation/de-percolation, help redistribute contact forces more homogeneously, reducing significantly the stress fluctuations over time and allowing the system to more efficiently dissipate energy, hence allowing for shear thinning.

7 Rigid structures in the contact networks

Our data suggest that the capacity of 3-SFN to transmit stress makes their percolation important for the system rheology. That is, the sensitivity of the percolation/de-percolation of 3-SFN to σ, as shown in Fig. 9(a), must be related to the rigidity of the percolating 3-SFN. At rest, and in two dimensions, rigid structures are usually effectively identified by the pebble game algorithm,⁴⁸ which accurately accounts for all constraints acting locally and globally on the microscopic degrees of freedom of a given particle configuration. Since the algorithm allows for different particle interactions, it has been effectively utilized also for colloidal gels⁴⁹ and for sheared dry granular materials in static conditions.^26,35,50,51 Recent studies have shown that this approach can identify rigid structures relevant to shear thickening of dense suspensions,^33,36 suggesting that it can capture well enough the structural self-organization underpinning stress transmission even under flow, in spite of not explicitly including moving boundaries and the effect of continuously increasing strain. Comparing therefore the percolating 3-SFN discussed so far with the rigid percolating networks detected by the pebble game algorithm could allow us to deeper understand the growth and optimization of stress-bearing structures during thickening.

Here we have applied a (3,3) pebble game algorithm to our model suspensions, to directly compare the rigid structures identified in this way, at specific shear rates corresponding to the different rheological regimes, with those emerging from our 3-SFN analysis in Fig. 6. It should be noted that the pebble-game clusters are calculated by assuming that the system has already fully achieved rigidity, which may limit their capacity to capture the development or emergence of rigidity. We also note that the pebble-game algorithm considered here is the same used for shear jamming, i.e. it does not take into account the moving boundary conditions imposed by the steady-state flow. The networks identified by the pebble game algorithm appear to percolate at all rates studied here, independent from the specific rheological regime, suggesting that they do not necessarily capture specifically the may not be sensit structures responsible for the stress transmission. The comparison in Fig. 9(a), however, shows that during the thickening the percolation/de-percolation over time of both types of networks follows the fluctuations of the shear stress, while 3-SFN percolation may be better coupled to small stress fluctuations.

A more direct analysis (SI, Fig. S5–S7) further shows that there is a strong overlap between the two types of structures when percolating 3-SFN may occur but are not persistent (i.e., during the quasi-Newtonian and the thinning regime). However, when percolating 3-SFN become more persistent during the thickening regime, they are also heterogeneous and seem to contain a significant fraction (e.g. ∼50%) of particles that are not rigid according to the pebble game algorithm. Size seems to be a predictor of whether a percolating 3-SFN will also be rigid according to the pebble game during the thickening – 3-SFN containing around 60% of all the particles are generally mostly rigid according to the pebble game algorithm and those with less particles are not. Comparatively, the percolation of rigid clusters identified by the pebble game algorithm shows a less consistent degree of correlation (i.e. slightly lower and less consistent Pearson coefficient) with the time series of the shear stress in the thickening regime.

Overall it seems possible and consistent with the data provided here, that the 3-SFN may be better at capturing the rate dependence and the stress fluctuations because they are more likely to capture the structures responsible for the stress transmission. It would be useful in future work to perform further and extended analysis of similarities and differences of the two types of structures, sampling more rates and across different models.

8 Conclusions

Through numerical simulations of a minimal 2d model for dense suspensions under flow, we have identified specific subsets of the contact network that display significant changes depending on the different flow regimes. Flow helps the system self-organize into local structures that can rapidly grow and span the whole sample. Portions of the contact network that include disks with at least three frictional contacts (rigid or over-constrained) and where large contact forces develop tend to percolate during thickening. Their percolation/de-percolation produces large stress fluctuations that stem from the localization of large contact forces. Interestingly, these subsets of the contact network, which we call 3-SFN, have some overlap with rigid structures identified by the pebble game algorithm but seem to be more sensitive to changes in shear rate and flow regimes, and better at capturing the emergence of rigidity during flow.

The minimal model chosen here has many limitations but allowed us to effectively disentangle different concurrent effects, and demonstrate the generality of the shear thickening phenomenology, with potential implications for a broad range of complex materials. One interesting aspect, with this respect, is that, in spite of the model oversimplified description of the solvent, the results clearly demonstrate that contacts are hierarchically organized under flow, building networks that allow for localization of large contact forces during thickening.

The softness of the particles in the model leads to a thinning regime at higher shear rate, through which the system avoids (or delays) shear jamming. Thanks to this feature of the model, we were able to identify the deformation of the local structures initially created by the flow, resulting in the de-percolation of the large scale assemblies that they form. The deformation of those structural units also allows for a more homogeneous redistribution of contact forces, which reduces significantly the large stress fluctuations that characterize thickening. These findings further reinforce the idea that shear thickening emerges from the growth of large scale, spanning assemblies of specific structural elements, which allow for transmitting stresses and localizing large contact forces.

When attempting a comparison with experiments on shear thickening suspensions, physical phenomena such as significant solvent migration and the associated concentration fluctuations^8,19 indicate complexities that certainly can not be captured by our minimal computational model. However, the rapid fluctuations of percolating structures observed here could be relevant to experimental observations of dramatic and rapid fluctuations in boundary stress and particle configurations.^8,16,19 The softness of the particles in the model, which allows for significant overlaps at high shear, may be relevant to specific experimental systems where particle deformability can promote the growth of concentration fluctuations during thickening.

An interesting avenue for future work would be to explore further the percolation transition of 3-SFN during thickening. In particular, the connection to the rigidity percolation transitions discussed in the literature^21,35 and to other properties of rigid clusters³⁶ should be explored.

Author contributions

W. C. J. B. performed all simulations and analysis. V. H. A. provided some of the codes. All the authors designed the research, discussed the interpretation and wrote the article.

Conflicts of interest

There are no conflicts to declare.

Data availability

The data for all figures in this article are available at a Zenodo Open Repository (https://doi.org/10.5281/zenodo.18002353). For information on LAMMPS scripts and analysis codes please contact the corresponding authors.

The supplementary information contains additional data on the system size dependence of the flow curve, on shear stress fluctuations and cluster size distributions, and on a comparison between rigid clusters and SFN clusters across different rheological regimes. See DOI: https://doi.org/10.1039/d5sm00928f.

Acknowledgements

This work was supported by the National Science Foundation (NSF DMR-2226485 and NSF STC COMPASS – Center for Complex Particle Systems, Award No. 2243104). All the simulations discussed in this work were performed using High-Performance Computing capabilities available at Georgetown University through Google Cloud Platform.

References

1. R. Seto, R. Mari, J. F. Morris and M. M. Denn, Phys. Rev. Lett., 2013, 111, 218301.
2. R. Mari, R. Seto, J. Morris and M. M. Denn, J. Rheol., 2014, 58, 1693–1724.
3. M. Wyart and M. E. Cates, Phys. Rev. Lett., 2014, 112, 098302.
4. N. Y. C. Lin, B. M. Guy, M. Hermes, C. Ness, J. Sun, W. C. K. Poon and I. Cohen, Phys. Rev. Lett., 2015, 115, 228304.
5. A. Singh, R. Mari, M. M. Denn and J. F. Morris, J. Rheol., 2018, 62, 457–468.
6. J. F. Morris, Phys. Rev. Fluids, 2018, 3, 110508.
7. E. D. Gado and J. Morris, J. Rheol., 2020, 64, 223–225.
8. E. Moghimi, D. L. Blair and J. S. Urbach, J. Rheol., 2025, 69, 343–352.
9. E. Moghimi, J. S. Urbach and D. L. Blair, J. Colloid Interface Sci., 2025, 678, 218–225.
10. E. Moghimi, D. L. Blair and J. S. Urbach, Front. Phys., 2024, 12, 1480376.
11. J. M. Miller, D. L. Blair and J. S. Urbach, Front. Phys., 2022, 10, 991540.
12. W. H. Boersma, P. J. M. Baets, J. Laven and H. N. Stein, J. Rheol., 1991, 35, 1093–1120.
13. D. Lootens, H. Van Damme and P. Hébraud, Phys. Rev. Lett., 2003, 90, 178301.
14. D. Lootens, P. Hebraud, E. Lecolier and H. V. Damme, Oil Gas Sci. Technol., 2004, 59, 31–40.
15. C. Heussinger, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2013, 88, 050201.
16. V. Rathee, D. L. Blair and J. S. Urbach, Proc. Natl. Acad. Sci. U. S. A., 2017, 114, 8740–8745.
17. E. Y. X. Ong, M. Ramaswamy, R. Niu, N. Y. C. Lin, A. Shetty, R. N. Zia, G. H. McKinley and I. Cohen, J. Rheol., 2020, 64, 343–351.
18. V. Rathee, D. L. Blair and J. S. Urbach, J. Rheol., 2020, 64, 299–308.
19. V. Rathee, J. Miller, D. L. Blair and J. S. Urbach, Proc. Natl. Acad. Sci. U. S. A., 2022, 119, e2203795119.
20. A. Gauthier, G. Ovarlez and A. Colin, J. Colloid Interface Sci., 2023, 650(Part B), 1105–1112.
21. A. Goyal, N. S. Martys and E. D. Gado, J. Rheol., 2024, 68, 219–228.
22. K. P. Krishnaraj and P. R. Nott, Phys. Rev. Lett., 2020, 124, 198002.
23. M. E. Cates, J. P. Wittmer, J.-P. Bouchaud and P. Claudin, Phys. Rev. Lett., 1998, 81, 1841–1844.
24. D. Bi, J. Zhang, B. Chakraborty and R. P. Behringer, Nature, 2011, 480, 355–358.
25. L. Papadopoulos, M. A. Porter, K. E. Daniels and D. S. Bassett, J. Complex Networks, 2018, 6, 485–565.
26. H. A. Vinutha and S. Sastry, Phys. Rev. E, 2019, 99, 012123.
27. H. A. Vinutha and S. Sastry, Nat. Phys., 2016, 12, 578–583.
28. D. Pan, Y. Wang, H. Yoshino, J. Zhang and Y. Jin, Phys. Rep., 2023, 1038, 1–18.
29. M. Ramaswamy, I. Griniasty, D. B. Liarte, A. Shetty, E. Katifori, E. D. Gado, J. P. Sethna, B. Chakraborty and I. Cohen, J. Rheol., 2023, 67, 1189–1197.
30. M. Gameiro, A. Singh, L. Kondic, K. Mischaikow and J. F. Morris, Phys. Rev. Fluids, 2020, 5, 034307.
31. O. Sedes, H. A. Makse, B. Chakraborty and J. F. Morris, Phys. Rev. Fluids, 2022, 7, 024304.
32. M. Nabizadeh, A. Singh and S. Jamali, Phys. Rev. Lett., 2022, 129, 068001.
33. M. van der Naald, A. Singh, T. T. Eid, K. Tang, J. J. de Pablo and H. Jaeger, Nat. Phys., 2024, 20, 653–659.
34. A. D’Amico, S. Tu and A. Singh, J. Chem. Phys., 2025, 162, 214906.
35. S. Henkes, D. A. Quint, Y. Fily and J. M. Schwarz, Phys. Rev. Lett., 2016, 116, 028301.
36. A. Santra, M. Orsi, B. Chakraborty and J. F. Morris, Phys. Rev. Res., 2025, 7, 013275.
37. R. Pandare, M. Orsi, M. D. Shattuck, B. Chakraborty and J. F. Morris, J. Rheol., 2026, 70, 343–359.
38. A. G. Bayram, F. J. Schwarzendahl, H. Löwen and L. Biancofiore, Soft Matter, 2023, 19, 4571–4578.
39. M. Maiti, A. Zippelius and C. Heussinger, Europhys. Lett., 2016, 115, 54006.
40. C. S. O'Hern, L. E. Silbert, A. J. Liu and S. R. Nagel, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2003, 68, 011306.
41. P. A. Cundall and O. D. L. Strack, Géotechnique, 1979, 29, 47–65.
42. E. Azéma and F. Radja, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2012, 85, 031303.
43. A. P. Thompson, H. M. Aktulga, R. Berger, D. S. Bolintineanu, W. M. Brown, P. S. Crozier, P. J. in 't Veld, A. Kohlmeyer, S. G. Moore, T. D. Nguyen, R. Shan, M. J. Stevens, J. Tranchida, C. Trott and S. J. Plimpton, Comput. Phys. Commun., 2022, 271, 108171.
44. Q. Xu, A. Singh and H. M. Jaeger, J. Rheol., 2020, 64, 321–328.
45. R. Seto, A. Singh, B. Chakraborty, M. M. Denn and J. F. Morris, Granul. Matter, 2019, 21, 82.
46. M. Doi, Soft Matter Physics, Oxford University Press, 2013.
47. Introduction to percolation theory, ed. D. Stauffer and A. Aharony, Taylor & Francis, 2018.
48. D. J. Jacobs and B. Hendrickson, J. Comput. Phys., 1997, 137, 346–365.
49. S. Zhang, L. Zhang, M. Bouzid, D. Z. Rocklin, E. Del Gado and X. Mao, Phys. Rev. Lett., 2019, 123, 058001.
50. V. Babu, H. A. Vinutha, D. Bi and S. Sastry, Soft Matter, 2023, 19, 9399–9404.
51. H. Dashti, A. A. Saberi, S. Rahbari and J. Kurths, Sci. Adv., 2023, 9, eadh5586.

---