Cyclically sheared colloidal gels: structural change and delayed failure time

Abstract

We present experiments and simulations on cyclically sheared colloidal gels, and probe their behaviour on several different length scales. For experimental gels formed by colloid–polymer mixtures, the shearing induces structural changes, which are quantified by the topological cluster classification, bond-order parameters, and the pore size distribution. These results are mirrored in computer simulations of a model gel-former: for cyclic shear with amplitudes up to 4%, local structural analysis shows that the material evolves down the energy landscape under shearing, and the average pore size increases. We also analyze mechanical responses including the stress and the dissipation rate, revealing a crossover between elastic and plastic responses as the strain amplitude is increased. Depending on the parameters, we observe both increased compliance after shearing (thixotropy), and reduced compliance (strain hardening). We simulate creeping flow under constant shear stress, for gels that were previously subject to cyclic shear, showing that strain-hardening also increases gel stability. This response depends on the orientation of the applied shear stress, revealing that the cyclic shear imprints anisotropic structural features into the gel.

---

1. Introduction

Colloidal gels are assemblies of particles consisting of interconnected networks of strands, which are kinetically arrested far from equilibrium.^1–6 They exhibit complex relaxations and responses to external forces, whose understanding is important for applications, and as a fundamental outstanding problem in soft matter.^7–11 These materials are particularly challenging because their behaviour is influenced by structural features on different length scales, including local motifs in the microscopic structure,¹² the response of individual strands to applied stress,^13–15 and the heterogeneous network.^16–20

Gel structure depends on particle interactions (strength and range of attractive forces) and on their volume fraction. In addition, the fact that gels are far from equilibrium means that their structures also depend on their mechanical and thermal history; they also experience physical aging, so their properties depend on the time elapsed since their preparation.^21–29 This feature can be exploited to engineer materials with specific properties, but the relationships between history, structure and gel properties are complex, and theoretical predictions are limited, so that formulation of gels often requires a large component of trial-and-error.

Among the material properties that one would like to control are the linear response to external stress (compliance) and the yielding behavior.^30,31 The process of strain-hardening offers a promising route for this, in that mechanical processing of an already-formulated material can be used to suppress yielding and/or reduce compliance. Strain-hardening is familiar in other amorphous solids such as molecular and metallic glasses^32–34 where the shear provides a mechanism for the system to descend in its energy landscape and form more stable structures (see also ref. 35). Shear-thinning can also occur in glasses³⁶ and the rheology of glasses and dense suspensions is a rich topic in its own right.^37,38

In colloidal gels, the underlying network structure makes the behaviour even richer than glassy or non-Brownian suspensions. The pioneering work of Gisler et al.³⁹ showed evidence for strain-hardening in fractal colloidal gels, see also ref. 40–42; disruption of the network by shearing is a mechanism for shear-thinning (thixotropy)^30,43 as also found in theoretical and computational models.^44–46 Studies on a wide range of experimental colloidal gels show that indeed mechanical history affects the subsequent rheology, with examples including colloid–polymer mixtures,^47–51 carbon black,^52,53 silica,⁵⁴ and clays.^55–57 A well-characterised example of such history dependence is the Bauschinger effect,⁴³ in which the direction of an applied shear is reversed; the subsequent (“backward”) response retains a dependence on the original (“forward”) shear strain.⁵⁸ More recent works have used history-dependence to imprint more complex memories in gel structure,⁵⁹ building on previous work on athermal systems,^60,61 see also ref. 62. The extent to which shear flow is able to break bonds in the gel structure can be quantified by the Mason number,^63,64 which determines the relative strength of the applied shear stress and the interparticle forces (see also ref. 47).

The systems considered here are gels that form by depletion in colloid–polymer mixtures, which we analyse both experimentally and by computer simulation. The experiments combine a shear stage with in situ particle-resolved imaging by 3d confocal microscopy, enabling microscopic changes in structure to be probed. Such methods have been used extensively to investigate colloidal glasses.^65–67 In the case of colloidal gels, combined imaging and shear has enabled direct observation of shear induced microstructural changes.^51,68,69 Rheological methods such as cyclic shear can also be used to identify a gel point where the elastic and viscous moduli cross.⁷⁰ The depletion interactions in these systems lead to interparticle attractions that are comparable to the thermal energy; this may be contrasted with the stronger interactions found in carbon black, or systems dominated by van der Waals interactions. Nevertheless, depletion can produce stable gels, whose properties can be further tuned by their mechanical history, including strengthening via rejuvenation.^48,50

We briefly summarise our results, before proceeding to a more detailed discussion. We first measure structural changes under cyclic shear, using both experiments and simulations. The gel structures evolve under shear by two important processes: the breaking of gel strands,^13–15 and the evolution of the microscopic structure towards lower energies.⁷¹ The experimental data reveal this through their pore-size distributions (mesoscopic structure) and through the topological cluster classification (TCC)⁷² which characterises the microscopic structure. The same behavior is also present in our computer simulations.

The second set of results shows the dynamical behavior of gels in computer simulation, including the rheological properties during cyclic shear, and measurements of the yielding under constant stress (creep). We present evidence for strain hardening, based on compliance and dissipation during the shear cycles. We also subject the resulting gels to constant stress, which leads to creeping flow and eventual yielding. These samples survive for longer times before failure, compared to gels that were not strain-hardened. Moreover, this enhanced stability depends significantly on the relative directions of the shearing motions for the cyclic and creeping protocols. In other words, the hardening can be used to enhance gel stability and to imprint anisotropic responses to subsequent shearing.

These results demonstrate specific situations in which gel properties can be selected by mechanical processing. The particle-resolved experiments reveal the microscopic and mesoscopic changes in structure; the simulations show how these can be harnessed to enhance gel stability and to tailor anisotropic responses. We discuss how our understanding of these far from equilibrium materials might be harnessed to for prediction and design of material properties more generally.

The structure of the paper is as follows: Section 2 describes models and methods, and Section 3 presents results for structural change due to cyclic shear. Section 4 discusses the dynamical behaviour, after which Section 5 concludes with a discussion of our results and directions for future work.

2. Models and methods

We consider gels formed of spherical colloidal particles with attractive interactions, leading to reversible bonds between them. These bonds may be strong enough that the suspension supports dense (colloidal liquid) and dilute (colloidal vapour) phases. In this case gels form by arrested spinodal decomposition,^3,5 whose dynamical arrest is controlled by the colloidal liquid phase where the packing is dense, particles have many bonds, and the local structure is reminiscent of colloidal glasses,^5,73 see also ref. 74. Due to the link with spinodal decomposition,³ the effective attraction at the onset of gelation is very close to its value at the colloidal liquid-vapour critical point, which may be related to the microscopic interaction potential by the Noro–Frenkel criterion.⁷⁵

We briefly summarise the system's control parameters, with further details given below. To characterise the effective interactions in the simulations and experiments we consider the ratio q between interaction range and the colloid diameter [small script l]; and the ratio ε₀ = ε/(k_BT) between effective interaction strength ε and the thermal energy. The colloid volume fraction is ϕ. For systems under cyclic shear, the shear amplitude γ₀ is a relevant dimensionless parameter, as is the product [small gamma, Greek, dot above]τ_B of the shear rate with the Brownian time τ_B = πη_s[small script l]³/(8k_BT), where η_s is the solvent viscosity. The Mason number is proportional to [small gamma, Greek, dot above]τ_B/ε₀, the relevance of this number to our cyclic shear experiments is discussed in Section 5.

2.1. Experiment

2.1.1. Colloidal particles and effective interactions. We consider sterically stablised poly-methyl methacrylate colloids in a cis-decalin and cyclohexylbromide solvent mixture which matches the density and refractive index. The particles were fluorescently labelled with 1,1-dioctadecyl-3,3,3,3-tetramethyl-indocarbocyanine perchlorate.⁷⁶

The particles are dispersed in a solution (4 mM) of tetrabutyl ammonium bromide, to screen the electrostatic interactions. Depletion attractions are provided by a nonadsorbing polymer (polystyrene, molecular weight 8.4 MDa). The (mean) particle diameter is [small script l] = 1500 nm, we estimate the polydispersity at 5%. We estimate the Brownian time as τ_B = 4.0 s (see above for the definition). The colloid volume fractions are obtained from the imaging results, see Table 1.

Table 1 Experimental parameters for shear. The rate [small gamma, Greek, dot above] is given relative to the Brownian time

ε	[small gamma, Greek, dot above] (/τB^−1)	γ max	ϕ
	3.72	0.0337	0.17
	9.32	0.0337	0.18

---

Our experiments at these volume fractions show that the onset of gelation occurs at concentration . Here we present data for polymer concentration c_p in the range 3.42–4.10 g L⁻¹, specifically and . The sample dimensions are 1 cm × 1 cm × 100 μm. We extract the positions of the particles from confocal microscopy images using the trackpy package. We consider the particles inside a cubic visualisation region (“box”) of linear size L ≈ 33[small script l]. This results in co-ordinates of N ≈ 1.5 × 10⁴ particles, which show that the particles form a percolating network, characteristic of the gel state.

To estimate more precisely the effective interactions, we follow previous work⁷³ by inferring the parameters of an one-component Asakura-Oosawa (AO) model,⁷⁷ which is known to be a reasonably accurate description of colloid–polymer mixtures.^78–80 In addition to the particle diameter [small script l], its parameters are the polymer radius of gyration R_g and the polymer fugacity z_p (proportional to c_p, which is the “reservoir” concentration⁷³).

As discussed in ref. 73, estimating R_g directly from properties of the polymer and solvent leads to relative errors around 10%, which are significant because the interaction strength varies as R_g³.⁷⁷ Instead we infer AO parameters consistent with the experiment by assuming that the onset of gelation at coincides with the polymer concentration at the liquid–vapor critical point of the AO model; this can be related to the underlying interaction potential by the Noro–Frenkel extended law of corresponding states.⁷⁵ (This approach for inferring interactions is effective in practice;⁷³ it may justified by the fact that gelation in colloid–polymer mixtures proceeds by arrested spinodal decomposition,^3,5 note also that the onset of gelation is weakly dependent on the colloid volume fraction ϕ in the range considered here.)

This approach yields a polymer radius of gyration R_g = 155 nm, which is a consistent estimate for this polystyrene polymer, in a good solvent; the depth of the attractive potential well at the critical point is ε* = 3.16k_BT. Hence the polymer–colloid size ratio is q = (2R_g/[small script l]) = 0.207. Note that for this size ratio, the colloidal liquid–vapour critical point is metastable to fluid-crystal phase coexistence.⁶ Since ε ∝ c_p the two polymer concentrations considered here then correspond to dimensionless interaction strengths ε₀ = ε/(k_BT) = 3.95 and 4.74.

2.1.2. Shear cell and shearing protocol. We use a shear stage to apply cyclic (simple) shear with motion in the y direction and velocity gradient in the z direction. The shear stage comprised an aluminium flexure with a coverslip window, mounted on a plexiglass baseplate and separated by a rubber spacer, which served to both dictate the gap thickness and to prevent solvent evaporation. Both the top coverslip and bottom glass surface were coated with a sintered particle layer to allow the gel to adhere to both surfaces. Displacement of the top window was driven by a Thorlabs TPZ001 T-cube piezo driver connected to a Thorlabs piezo stack with free stroke displacement of 25.5 microns, controlled by a MATLAB script.

The shear is a triangular wave with constant rate [small gamma, Greek, dot above], that is, the strain is increased linearly in time from γ = 0 to γ = γ_max and then decreased to γ = −γ_max before finally increasing back to γ = 0. We alternate shearing and imaging, the acquisition time for an 3d image of the gel is approximately 30 s (7.5τ_B); there is negligible aging during this time. The location of the visualisation region does not stay constant with time, this ensures that data is not too much affected by photobleaching that would occur if all imaging occurred in the same place. The waiting time t_w between sample loading and the start of shearing is approximately 30 minutes (450τ_B).

We consider two state points that differ in the polymer concentration and the parameters of the applied shear (details given above). The parameters for the shear experiments are reported in Table 1.

2.2. Computational

2.2.1. Model. For numerical simulations, we use an established model of a size-polydisperse colloid–polymer mixture^{9,12,24,80–82} which has been accurately mapped to the standard two-component Asakura-Oosawa model.^77,80 The depletion interaction between colloidal particles is modeled by a truncated and shifted Morse potential

U(r) = ε[e^{−2α(r−[small script l]_ij)} − 2e^{−α(r−[small script l]_ij)} + c_sh], r ≤ r_c (1)

where [small script l]_ij is the average diameter of particles i and j; the interaction strength is ε and α sets the range; the cutoff parameter is r_c = 1.4[small script l]_ij, and we truncate so that U(r) = 0 for r > r_c; the constant shift c_sh is chosen so that U is continuous at r = r_c. We consider N particles in a cubic box of volume L³, with periodic boundaries.

The overdamped colloid motion is modeled through Langevin dynamics with a large friction constant (see below for further details). Particles have mass m and they evolve with Langevin dynamics so the position r_i of particle i obeys dr_i/dt = v_i, and

where is the total potential energy, λ is the friction constant, u^aff_i is the local velocity of the (implicit) solvent, and ξ is a standard Gaussian white noise. (In the absence of shear flow u^aff_i = 0, the sheared case is discussed below.) The velocity damping time is τ_d = m/λ. All simulations are performed in LAMMPS.⁸³

To avoid crystallization we consider a size polydisperse system. We have taken 7 types of particles with diameters equally spaced between and , with relative concentrations [0.0062,0.0606,0.2417,0.3829,0.2417,0.0606,0.0062] to mimic a Gaussian distribution of diameters with 8% polydispersity.

We work with non-dimensionalized parameters throughout. As in ref. 84, we estimate the critical interaction strength for spinodal decomposition as ε* ≈ 3.1k_BT, based on the Noro–Frenkel criterion.⁷⁵ We report interaction strengths relative to this boundary, we focus on two state points with ε = 4.5k_BT = 1.45ε* and ε = 10k_BT = 3.22ε*. Our unit of distance is the mean particle diameter and the unit of time is . The colloid volume fraction is . The dimensionless parameter governing the interaction range is , we set α₀ = 33 throughout. The non-dimensionalised friction is λ₀ = λτ_sim/m. We take λ₀ = 10 to mimic overdamped dynamics, larger values of λ₀ would approximate more accurately colloidal motion in the solvent but the computational cost is higher and we have found previously that it changes very little the emergent behaviour.⁸¹ The integration time step is Δt = 0.001τ_sim. A natural time scale for colloid motion is the Brownian time which is the typical time for an overdamped free particle to diffuse its radius. For the parameters chosen here τ_B = 0.417τ_sim.

This simulation method is computationally efficient and has been shown to capture the essential features of depletion gels,^9,24 including quantitative comparisons with experiment.^12,81 As such, it complements modelling approaches that focus on the gels' network structure⁴⁴ but do not resolve the internal structure of the gel strands. On the other hand, our model neglects hydrodynamic interactions, which do affect some aspects of gel structure.^85–88

In the following, we have not attempted a quantitative matching of the parameters of simulation and experiment, because the implementation of the cyclic shear is necessarily different (not least that the simulation uses periodic boundary conditions, and an idealized thermostat that has to absorb the dissipated heat during the shear but has not been parameterised for the experiment). Despite these differences, we will show that the simulations and experiments exhibit the same qualitative behaviour.

2.2.2. Simulations.
2.2.2.1. Initialisation. Simulations are initialised in random configurations at volume fraction ϕ = 0.2. Gels are formed by simulating for a time t_w, during which spinodal decomposition occurs. Unless otherwise stated we simulate t_w = 3 × 10⁴τ_sim (≈7 × 10⁴τ_B) and system size of N = 10⁴ particles; all results are averaged over many independent runs (typically 50), to enable statistically robust conclusions.
2.2.2.2. Cyclic shear. After gel preparation we perform strain-controlled cyclic shear at finite temperatures and strain rates. The strain γ_xy = γ(t) of the system is varied cyclically with a triangular wave with amplitude γ_max and rate |[small gamma, Greek, dot above]| = 0.01τ_sim⁻¹. This cycle is repeated many times. The simulated shear rates are significantly slower than those of the experiment, to ensure that the thermostat can easily absorb the energy that is injected by this external work. We would expect qualitatively similar behaviour for larger shear rates. We compute observable quantities at the stroboscopic configuration (γ = 0) after each cycle of strain and study the evolution of these quantities as a function of number of cycles n_cyc.

We use Lees–Edwards boundary conditions. The solvent velocity u^aff_i in (2) points in the x direction, with magnitude y_i[small gamma, Greek, dot above], where y_i denotes the y-coordinate of particle i.

2.2.2.3. Creep simulation. Gels are non-equilibrium states and the cyclic shear affects their structure. After shearing, we performed stress-controlled creep dynamics to probe gels' yielding behaviour. The procedure follows that of ref. 84, a non-dimensionalized stress is maintained through a feedback control scheme that is implemented^89,90 as

∂_t[small gamma, Greek, dot above] = B[σ₀ − σ_xy(t)] (3)

where [small gamma, Greek, dot above] is the shear rate and σ_xy is the observed (non-dimensional) shear stress [measured from the virial]; also B is the damping parameter determining how quickly the applied stress relaxes to its imposed value. Similar to the cyclic shear, we used Lees–Edwards boundary conditions with affine flow in the x-direction. We take B = 0.01τ_sim⁻² as in ref. 84 which ensures that the stress is imposed accurately, at a reasonable computational cost (such schemes involve inertial effects including ringing, see ref. 84 and 91 for further discussion).

3. Results – gel structure (experiment and simulation)

3.1. Gel visualisation

Fig. 1(top) shows slices through the experimental system. (Specifically, we visualise particles in a region of size L × L × Z with Z = 8[small script l], for various n_cyc.) One clearly sees a percolating network of colloidal particles, characteristic of the gel state. Shearing the sample for 100–600 cycles causes structural rearrangements including some pores changing their shapes. Fig. 1(bottom) presents similar slices from numerical simulation. As noted above, we have not attempted to match the parameters precisely, in particular the waiting times before start of shear are significantly longer in the simulation (leading to thicker gel strands), and the dimensionless shear rate ([small gamma, Greek, dot above]τ_B) is much smaller in the simulations. Even so, the shearing leads to similar structural changes in both samples, in that the structure and pores are rearranged.

Fig. 1 Visualization of gel under cyclic shear (experiment and simulation). (Top) Rendering of a slice through an experimental gel, to visualize the strands and pores after different strain cycles n_cyc = 0, 100, 600 for the gel with . Arrows indicate the geometry of the (simple) shear. The region shown has volume L × L × Z with Z = 8[small script l]. (Bottom) Similar slices based on simulation data for ε₀ = 1.5ε*, γ_max = 0.04. Particles are coloured by their coordination number, according to the colorbar.

In the following subsections we analyse the results in detail for both experiment and simulation. We mainly focus on the structural change at different length scales due to cyclic shear deformation. From the simulations, we also show in Section 4 that the cyclic shear hardens the gel.

3.2. Experiment

3.2.1. Microscopic structure. To probe microscopic structural changes, we measure several local quantities, which mirror those measured in the simulation study of ref. 84. Specifically, we consider the coordination number N_b, the two-fold bond-orientation parameter q₂, and the average number of different types of clusters in which particles participate, as obtained from the TCC.⁷² Particles within a distance of 1.8[small script l] of particle i are defined as its neighbours. (This identification method for neighbours is maintained for all the experimental data analysis.) The bond orientational order parameter q_l is calculated following:⁹²where Y_lm are spherical harmonics and N_b(i) is the number of neighbours of reference particle i. We consider the case l = 2: the quantity q₂ is large when a particle's bonds are collinear, which has previously been observed to correlate with bond stretching.¹⁰

In Fig. 2, we show the evolution of these quantities as the system is sheared for two different gel samples. On increasing n_cyc, the average coordination number N_b increases [Fig. 2(a)]. This is due to a coarsening effect whereby the shearing increases the thickness of the gel strands. Fig. 2(b) shows that q₂ has a decreasing trend with n_cyc. This q₂ is associated with stretching of interparticle bonds: one tends to find larger values in thinner strands, and smaller values in thicker strands where the structure resembles the bulk. Hence the decreasing q₂ is consistent with the increasing N_b due to coarsening.

Fig. 2 Microscopic structural analysis of sheared gel (experiment). (a) Coordination number N_b as a function of the number of shear cycles n_cyc. (b) Two-fold bond orientation order parameter. (c and d) TCC analysis, showing the number of trigonal bipyramids and pentagonal bipyramids in which the particles participate. Insets show the relevant clusters. See Fig. 4 for a similar analysis of simulation data.

The data show significant statistical fluctuations (each data point is based on a single set of approximately 10⁴ particle co-ordinates). Nevertheless, the systematic trends are clear, with significant changes becoming apparent after approximately 10³–10⁴ cycles.

We also performed a TCC analysis⁷² which detects cluster motifs as a sensitive probe of local structure. (We consider here clusters which minimize the energy for the Morse potential.) This provides a characterization of the pathway that the system makes through the free energy landscape, at a local level. Amorphous states with more bonds and denser packing tend to be rich in locally favoured structures – for spheres with short-ranged attractions these include trigonal bipyramids (5-particle clusters); pentagonal bipyramids (7-particle clusters); and icosahedra (13-particle clusters). Increasing numbers of such clusters indicate structures that are lower in the energy landscape,¹² as happens for example during aging⁷¹ or annealing, see also ref. 93.

Fig. 2(c and d) show TCC results, specifically the variation of the number of trigonal bipyramids 〈n_tb〉 and pentagonal bipyramids 〈n_pb〉 in which a particle participates as a function of strain cycles. These TCC structures contain 5 and 7 particles respectively, and are sensitive to details of the packing. Both the quantities show an increasing trend with n_cyc, indicating that the system evolves downhill in the energy landscape as the shear-induced coarsening (and aging) takes place. Similar to previous measurements, significant changes become apparent after approximately 10³–10⁴ cycles.

These results complement previous analysis of gels using the TCC, which considered aging,⁷³ and the effects of hydrodynamics,⁹⁴ quench rate, and polydispersity.⁹⁵ As gelation occurs, the emergence of an elastic response is accompanied by a gel network rich in tetrahedra.¹² As time progresses, these systems pass through a sequence of progressively largely clusters with a change from tetrahedral to fivefold symmetry,^73,95 as annealing (or coarsening) takes place. The trigonal bipyramids considered in this work correspond to pairs of tetrahedra while the pentagonal bipyramids are larger clusters with five-fold symmetry. Hence the results Fig. 2(c and d) show that shearing that promotes annealing. (After annealing for very long times, one may expect the emergence of crystalline packings, because the colloid liquid is metastable;⁶ this effect is strongest in monodisperse systems and can be seen by the TCC.⁹⁵ We did not see significant crystallization in this work.)

3.2.2. Mesoscopic structure. We now consider the mesoscopic structure of the gel by measuring its pore size distribution. We measure the distribution of pore sizes in the gel as defined in ref. 96, the detailed procedure is described in ref. 84: We select a random point in the gel and we find the largest possible sphere that encompasses that point, without overlapping with any colloidal particles. We identify the pore size D_pore as the diameter of this largest sphere. By repeating this process for many random points we obtain the distribution of pore sizes.

Fig. 3(a) shows the distribution of the pore sizes D_pore, as the shear accumulates, for the gel with . The distribution shifts towards larger pores as the shear n_cyc increases. The maximal D_pore also has an increasing trend with n_cyc. The associated mean pore size [Fig. 3(b)] has a mild initial increase before growing much more suddenly. Similar to the microstructural measurements, significant changes become apparent after approximately 10³–10⁴ cycles.

Fig. 3 Mesoscopic (pore-size) analysis of sheared gel (experiment). (a) Distribution of pore diameter for different numbers of strain cycles for the gel with . (b) Average pore diameter against the number of cycles for two different gels. (c) Snapshot of a slice through gel before shearing . (d) Snapshot of the same gel after 10⁵ shear cycles. Color coding in (c and d) is the same as Fig. 1. All pore-sizes D_pore are measured in units of colloid diameter [small script l]. See Fig. 5 for a similar analysis of simulation data.

We find that pores can grow very large after extensive shearing, see Fig. 3(c and d) for n_cyc = 0, 10⁵ respectively. For the sample with , the data of Fig. 3(b) are consistent with an increasing trend for the mean pore size, but the results are subject to large fluctuations because the data contains a small number of large pores which strongly affect the estimated average. (Recall that the location of the visualisation box is not constant in time, so the large changes as a function of n_cyc do not represent large local changes.)

3.3. Numerical simulation

We now analyse results of numerical simulations, and compare with experiment. We consider two state points with different interaction strengths ε = 1.45ε*, 3.22ε*. At the microscopic level it has been shown previously for this system⁸⁴ that increasing ε during gel preparation leads to strands of reduced thickness with smaller N_b and larger q₂. Here we analyse the effect of the cyclic shear.

To study microscopic structure, particles within the interaction range r_c of particle i are defined as its neighbours for the evaluation of the co-ordination number N_b and for q₂. We present the variation of N_b and q₂ with the number of strain cycles in Fig. 4(a and b). The increasing behaviour of N_b and a decreasing trend in q₂ with n_cyc mirror the experimental observations of Fig. 2. In Fig. 4(c) and (d) we present the data for 〈n_tb〉 and 〈n_pb〉 obtained from TCC analysis. For ε = 3.22ε*,we see a significant increase in the value of 〈n_tb〉 and 〈n_pb〉 mirroring experimental data. For ε = 1.45ε*, the change in these data are negligible, despite the increasing co-ordination number. This may indicate that the microscopic structure inside the gel strands is close to that of the (metastable) colloidal liquid, due to the long annealing times for the gel, and the relatively weak bonds. Then the agitation by shearing has little effect on the structure, although there is some coarsening of the strands, which increases N_b.

Fig. 4 Microscopic structural analysis of sheared gel (simulation), γ_max = 0.04. (a) Coordination number N_b for increasing shear n_cyc. (b) Bond orientation order parameter q₂. (c and d) TCC analysis showing the number of trigonal bipyramids and pentagonal bipyramids in which the particles participate. (These results may be compared with Fig. 2.).

For comparison with experiment, recall that we identify with ε/ε* so that the simulation state point ε = 1.45ε* is comparable to the experiments of Section 2.1; the simulations with ε = 3.22ε* have significantly stronger interactions than the experiment. These stronger interactions promote dynamical arrest via kinetic trapping: this leads to thinner gel strands and local structures that are less deep in the energy landscape, consistent with the smaller N_b and reduced number of trigonal/pentagonal bipyramids.

Next, we investigate the mesoscopic structure of the system by analyzing the pore size of the gel. As above, we extract the distribution of the pore size D_pore and we extract its mean value. Fig. 5(a) shows results for increasing n_cyc in a gel with ε = 3.22ε*. As n_cyc increases, the distribution shifts to the right with a larger mean. In Fig. 5(b) illustrates the increasing behaviour of average poresize with n_cyc for (ε/ε*) = 1.45, 3.22. To provide a visual representation we present snapshots of slices through the gel in Fig. 5(c and d) for n_cyc = 0 and n_cyc = 1300, respectively. Initially, at n_cyc = 0, the system exhibits a smaller pore size. However, after multiple applied cycles, the larger pores tend to elongate, contributing to an overall increase in the average pore size. Some examples of elongated large pores are highlighted in Fig. 5(d): quantitative analysis of these (three-dimensional) features is challenging in general, but we intend to return to it in future work.

Fig. 5 Mesoscopic (pore-size) analysis of sheared gel (simulation), γ_max = 0.04. (a) Distribution of pore diameter for different numbers of strain cycles. (b) Average pore diameter against the number of cycles. (c) Snapshot of a slice through the gel, before shearing (interaction strength ε = 3.22ε*). (d) The same gel after n_cyc = 1300 shear cycles; three different elongated pores are highlighted with different colours. (These results may be compared with Fig. 3 where one large elongated pore is also visible.).

3.4. Summary

We find a coherent picture of the microscopic and mesoscopic structural changes in these cyclically sheared colloidal gels, both in the experiment and numerical simulations. The main effect is that shearing leads to a coarsening effect, where colloids' mean co-ordination number increases, as does the pore size of the gel. This can be rationalised in a simple way by the idea that the shearing injects energy and accelerates the random bond-breaking processes that are anyway occurring due to thermal fluctuations. In this way, it accelerates the natural aging/coarsening dynamics of the gel. However, the shearing procedure also introduces anisotropy to the system, which is apparent (for simulations) in the shape of the pores in the gel [Fig. 5].

Overall, the simulation model successfully captures the key features of structural changes across different length scales resulting from the repetitive mechanical deformation of the gel, including the microscopic structure of the arms. This complements previous results for effects of shear on the network topology, in which the internal structure of strands was not modelled in detail.^44,45 In the next section, we use numerical simulations to investigate the mechanical response of the gel in more detail.

4 Results – mechanical properties (simulation)

So far, we considered cyclic shear of fixed amplitude γ_max = 0.04 in both experiment and simulation, and we analysed structural properties of the gel as a function of the number of cycles. In this section we exploit the ability of simulation to measure the time-dependent stress during cyclic shear. We also explore the dependence on the strain amplitude, revealing an interesting regime where the gel hardens under shearing.

4.1. Strain hardening before failure

4.1.1. Hysteresis loops. We consider gels with ε = 3.22ε*, as visualised in Fig. 5. Simulation results are shown in Fig. 6 for a range of γ_max. (Note, the shear rate is held constant here and throughout at |[small gamma, Greek, dot above]| = 0.01τ_sim⁻¹, as stated in Section 2.2.) Specifically, Fig. 6(a) shows stress–strain curves, obtained from long simulations (n_cyc = 10⁴), averaged over the final 200 cycles. For very small γ_max and small shear rates [small gamma, Greek, dot above] one expects an elastic response with σ_xy ∝ γ and no hysteresis loop. We find here that some hysteresis is present even for γ_max = 0.01, which we attribute to the fact that the rate |[small gamma, Greek, dot above]| is not extremely small.

Fig. 6 Cyclic shear for a simulated gel, interaction strength ε = 3.22ε*. (a) Average stress–strain curve over a strain cycle for several values of γ_max ranging from 0.02 to 0.12. (Data are averaged over the final 200 cycles, obtained from simulations of 10⁴ cycles.) (b) Zoomed view of stress–strain curves for γ_max = 0.035 at different values of n_cyc. (c) The area enclosed by the stress–strain curve as a function of n_cyc (normalized with the value of initial cycle A₀), for three different values of γ_max. (d) Maximum stress σ^xy_max at γ_xy = γ_max in a strain cycle against n_cyc for several γ_max.

On increasing the amplitude γ_max, the area of the hysteresis loop increases; the magnitude of the shear stress also depends non-trivially on γ_max. This non-trivial dependence is associated with increasing dissipation, the onset of plastic deformation, and yielding. To characterise it, we write σ^xy_max for the maximal stress during the cycle: this has a non-monotonic dependence on γ_max.^34,71 As a simple way to locate the crossover between elastic and plastic flow, we follow³⁴ and identify the maximal value of σ^xy_max with the transition to plastic flow; the associated shear amplitude is γ^Y_max. From Fig. 6(a) we identify γ^Y_max = 0.0350 ± 0.0005. Alternative methods for characterising this crossover are discussed in ref. 97.

In Fig. 6(b) we show the evolution of the stress–strain loop as the number of shear cycles increases, for γ_max = 0.035, close to the yield strain.† Close inspection of the four curves shows that the dependence on n_cyc is non-monotonic. Fig. 6(c) plots the area of the loop (averaged over 10 consecutive cycles) as a function of n_cyc, normalised by the value in the initial cycle. For γ_max = 0.035 the non-monotonic dependence is clear. We also define σ^xy_max as the maximal stress during the cycle. This quantity is shown in Fig. 6(d), again showing non-monotonic dependence on n_cyc for γ_max = 0.035 (see also Fig. S9 of ref. 71 for an analogous effect in glasses).

These non-monotonic dependencies are associated with strain-hardening of the gel, as we now explain. Note first that for small strain amplitude γ_max = 0.01, the area and the maximal stress hardly depend on n_cyc. This indicates that the gel is mostly responding elastically and the repeated shear cycles have little effect on its structure. For larger γ_max the overall trend is that the loop area increases and the maximal stress σ^xy_max decreases. This indicates that the gel's response includes plastic rearrangements that disrupt its structure and lead to extra dissipation (larger loop area); the disruption to the structure also increases the compliance (reduces σ^xy_max). However, for γ_max = 0.035, the range labelled “strain hardening” in the figure is associated with a reduction in both the dissipation and the compliance. Similar effects are observed for other values of γ_max although the signature is most pronounced for γ_max ≈ γ^Y_max.

Fig. 7 illustrates the corresponding evolution of the gel under shear. The main observation is that the changes are subtle, although the trends are consistent with Fig. 5 (which has slightly larger γ_max = 0.04). In particular, there is some coarsening of the structure with larger pores. The amorphous structure of the gel means that these small structural changes can nevertheless cause both strain softening and hardening behaviour, as a function of n_cyc.

Fig. 7 Snapshot of the gel structure as shearing takes place. The parameters are those of Fig. 6 [ε = 3.22ε*, γ_max = 0.035] so the system is strain-hardening from n_cyc = 10²–10³ and softening from 10³–10⁴, although the structural changes are subtle. The region shown has size L × L × Z with Z = 8[small script l] (similar to Fig. 5), the colouring shows the co-ordination number (similar to Fig. 1).

4.1.2. Dynamic moduli. To probe further the strain-hardening phenomena, we investigate dynamic moduli (storage modulus G′, and loss modulus G′′). More specifically, we measure “nonlinear moduli” which incorporate the dependence of the stress response on the frequency and amplitude of an oscillatory (strain-controlled) shear. These quantities are familiar from studies of large amplitude oscillatory shear (LAOS).⁹⁸ However, the standard theory is formulated for sinusoidal shear cycles, while our numerical data is obtained with a triangular (non-sinusoidal) waveform. To estimate the moduli in this case, we follow⁹⁹ and reparameterise the time t in terms of a variable s = s(t) such that γ(t) = γ_maxsin(2πs(t)). We denote the time taken for a single shear cycle as τ_shear = γ_max/(4[small gamma, Greek, dot above]), also note that s(τ_shear) = 1 and define the corresponding angular frequency as ω = (2π/τ_shear) = (8π[small gamma, Greek, dot above]/γ_max).

The stress profile for any shear cycle can then be fitted as a function of s, as σ(s) = σ₀ sin(2πs + δ) where the amplitude σ₀ and phase lag δ are fitting parameters. An example fit is shown in Fig. 8(a). This allows nonlinear moduli G′ and G′′ to be extracted from the formula

where σ₀,δ are obtained from the fit, while the amplitude γ_max is set by the underlying shear cycle. Since the gel is hardening as the shear accumulates, we perform this fitting separately for each individual cycle. We average the resulting σ₀,δ over many independent samples.‡

Fig. 8 Estimated shear moduli for a simulated gel, based on shear protocols with different amplitudes (interaction strength ε₀ = 3.22ε*). (a) Stress σ^xy as a function of rescaled time s. Solid line through the data points is a fit with σ^xy = σ₀ sin(2πs + δ). Sample averaged storage modulus G′ and loss modulus (G′′) as a function of strain cycle for (b) γ_max = 0.01, (c) γ_max = 0.035, and (d) γ_max = 0.07.

The resulting moduli depend on (ω,γ_max,n_cyc); recall that [small gamma, Greek, dot above] is fixed throughout this study so we study G′, G′′ as functions of (n_cyc,γ_max). Results are shown in Fig. 8(b–d) for different values of strain amplitude. In each case we measure how the moduli change as the oscillatory shear accumulates (increasing n_cyc). We discuss the different amplitudes in turn.

For small amplitude γ_max = 0.01 [Fig. 8(b)], both G′ and G′′ remain almost constant. For larger amplitude γ_max = 0.035 (close to the yield strain), Fig. 8(c), shows more interesting behaviour, in that both G′ and G′′ are non-monotonic functions of the number of cycles n_cyc. An initial decrease in G′ and a corresponding increase in G′′ indicate gel softening. At intermediate cycle number, however G′ begins to rise while G′′ decreases, indicating hardening of the gel. Notably, this transition region aligns with our earlier observations – specifically, the reduction in the stress–strain loop area and the increase maximum stress as shown in Fig. 6(c and d).

For a yet larger strain amplitude γ_max = 0.07 [Fig. 8(d)], a similar hardening behavior is observed, but it is restricted to a narrower range of cycles (n_cyc = 10–100). Beyond this regime the gel exhibits liquid-like characteristics, as indicated by G′ < G′′ for n_cyc > 10³.

4.1.3. Physical interpretation. Fig. 6 and 8 demonstrate both thixotropy and strain hardening, for appropriate ranges of γ_max and n_cyc. Physically, thixotropy is expected in situations where shearing disrupts the gel network, for example by breaking strands and opening up larger pores (recall Section 3). On the other hand, strain-hardening can be rationalised because shearing tends to anneal the local structure, leading to higher co-ordination and lower energy clusters (as observed by the TCC, for example).

The competition between these effects leads to a complicated phenomenology where strain hardening is strongest for intermediate γ_max close to γ^Y_max (see also ref. 54); it is also notable that the value of γ_max affects the range of n_cyc over which strain hardening is observed. Future work should investigate the interplay between structure and dynamics in more detail, in order to predict and control whether shearing leads to more or less compliant gels.

Taken together, these results complement previous works that demonstrated strain-hardening in other gel models, which had different microscopic interactions (and where thermal fluctuations were neglected).^44,45

4.2. Delayed failure time of cyclically sheared gel

Given the strain-hardening behaviour demonstrated above, one may expect that the resulting gels are also more stable: we now demonstrate that this is indeed the case. We simulated creep dynamics at constant stress, following.⁸⁴ The initial conditions for this procedure are gel configurations that have been cyclically sheared for n_cyc strain cycles. We then impose a constant (non-dimensionalised) shear stress σ₀ in the xy plane, allowing flow along the x-direction with Lees–Edwards boundary conditions (which is the same geometry used for the cyclic shear). We take σ₀ = 0.8, this value is chosen such that the system undergoes creeping flow and then fails within our computational time window.§

We measure the strain during this creep flow: its average behaviour is shown in Fig. 9(a), for samples that have been through n_cyc = 0, 10³, 10⁴ shear cycles, with γ_max = 0.035. The time-dependent strain has an initial elastic branch followed by a plateau, creeping flow, and eventually failure. The cyclic shear increases the height of the plateau: its value is comparable with γ_max indicating that the structure of the gel has rearranged during cyclic shear so that it can accommodate this strain.¹⁰⁰

Fig. 9 Creep Simulation, interaction strength ε₀ = 3.22ε* (a) Average strain against time for different values of n_cyc = 0, 10³, 10⁴. The dashed horizontal line represents γ = 0.4, which is taken as the critical value of strain beyond which system starts to flow. Vertical dashed lines indicate the failure time τ_f for each case. (b) Failure time τ_f as a function of n_cyc for different γ_max. (c) Strain against time for shear applied to different directions and (d) the corresponding failure time as a function of n_cyc.

We define the failure time τ_f of the gel as the time for 〈γ〉 to reach a threshold γ* = 0.4, as explained in ref. 84. From Fig. 9(a), the failure time is largest for n_cyc = 10³ cycles, after which it decreases again. This is consistent with Fig. 6, which indicates that the strain-hardening period ends at around n_cyc = 2 × 10³. That is, strain hardening is also accompanied by an increased failure time.

We performed creep simulations for a range of γ_max and n_cyc, from which we computed the failure times. Results are collated in Fig. 9(b), which we compare with “pristine” gels (n_cyc = 0). For γ_max = 0.035, strain hardening persists up to n_cyc = 10³, beyond which the material begins to soften again. This reinforces the conclusion of Fig. 9(a). The same characteristic non-monotonicity is also observed for other cyclic shear amplitudes in the range 0.02 ≤ γ_max ≤ 0.05, and the ranges of n_cyc over which this is observed mirror the ranges of strain hardening in Fig. 6(d). For small amplitude, γ_max = 0.01, the failure time for n_cyc = 10 is enhanced with respect to the pristine gel. After this, τ_f is weakly affected by further shearing. This mirrors the small structural change due to such cyclic shear [recall Fig. 6]. For large amplitude γ_max = 0.07 the failure time is reduced by the cyclic shear, presumably because this significantly disrupts the gel network.

These moderate enhancements of τ_f by cyclic shear (up to a factor of 3) are affected by two competing processes. The agitation by shearing promotes coarsening of the gel structure, and activates relaxation processes that allow it to descend down the energy landscape, as occurs in physical aging. This tends to increase the failure time τ_f. At the same time, the repeated cyclic strain imprints an anisotropic memory on the gel structure: this may predispose it to flow along the direction of the shear, reducing τ_f. (This process may be understood through the framework of fatigue failure,^101,102 where damage accumulates progressively until the system ultimately fails.).

To disentangle these effects, we simulated creeping flow with constant stress applied in the xz direction (the flow is still along x but the gradient is now in z). Results for γ_max = 0.035 are shown in Fig. 9(c). Comparing with Fig. 9(a), one sees that the enhancements of τ_f are much stronger when the creep flow is not in the same plane as the cyclic shear. Our interpretation is that the increase of τ_f due to coarsening and aging effects is dominating this response due to xz-stress, because memory effects due to xy-shear are less relevant in this geometry. That is, changing the direction of the shear stress in creep simulations separates the coarsening/shearing effects from the damage accumulation due to previous cyclic shear in the xy plane. Surprisingly, the increase of τ_f with n_cyc continues even up to n_cyc = 5000, at which time Fig. 6(d) suggests that the system has entered the shear-softening regime. This result further illustrates the complex dependence of gel structure and response on its mechanical history.

Fig. 9(d) plots the failure times for creeping flow with constant xz-stress, showing that τ_f is generically increased by shearing, regardless of whether the system is in the strain-hardening or softening regime. It seems that the softening effect at large n_cyc in Fig. 6(c and d) is an anisotropic effect, as may be expected from the mechanical protocol used. Given the strong anisotropy of the response [visible as differences between Fig. 9(b and d)], it would be interesting to explore in more detail the effects of different cyclic shear protocols on failure. (For example, combining cyclic shears in more than one plane⁵⁹ might generate more isotropic structures.)

5. Discussion

Colloidal gels are arrested far from equilibrium, so it is expected that their properties are affected by their mechanical history, including cyclic shear. However, it remains challenging to predict these effects and exploit them for material design: materials' history-dependence is complex in general and there may be competing effects. For example the shearing may break gel strands, but it also leads to coarsening and hence to thicker strands. In this situation, particle-resolved experiments and simulations offer valuable opportunities to disentangle different behaviors and to understand their microscopic mechanisms.

We have made progress in this direction, showing in both experiments and simulations that shearing leads to coarsening of the gel network and to denser microscopic packing within the strands and local structure indicating states deeper in the energy landscape. Consistent with recent work on glasses, this illustrates the general principle that shearing offers a mechanism for accelerating the progress of materials towards lower free-energy states. We also showed that these structural changes also affect material responses like compliance and stability.

Previous work on sheared gels^55,63,64 has highlighted the role of the Mason number Mn as a measure of whether the shear stress is strong enough to break interparticle bonds (see also the depletion Peclet number of ref. 47). Following⁶³ we define Mn = 3πη_s[small script l]³[small gamma, Greek, dot above]q/(2ε) = 12[small gamma, Greek, dot above]τ_Bq/ε₀ (recall that q is the interaction range measured in units of the diameter, it appears because a typical interparticle force is the well-depth ε divided by the range). As already discussed, the dimensionless shear rates for the experiments (Table 1) are much larger than those of the simulations, which is fixed throughout at [small gamma, Greek, dot above]τ_B ≈ 0.004 (Section 2.2). This leads to Mn ∼ 10⁻⁴ for the simulations, but Mn ∼ 1 for the experiments. At face value, this suggests that the two systems might be in quite different regimes, despite their similar behaviour.

However, we emphasize that we have performed cyclic shear with relatively small strains (a few percent), contrary to^55,63,64 which considered either steady flow or shear cessation, where the strains are much larger. Our interpretation is that the amplitude γ_max of the cyclic shear strain is controlling the behaviour, while the shear rate (hence also the Mason number) plays a secondary role. Physically: even if Mn is larger than 1, the system can accommodate small shear strains without bonds breaking. This interpretation is supported by Fig. 6 in which increasing γ_max causes a qualitative change in behaviour, while Mn remains constant. In this sense the regime studied here is different to that of ref. 64. For even larger Mn one may expect that the shear rate becomes relevant again but this situation is not considered here. This issue could be usefully addressed in future work.

More generally, our understanding of gel rheology will inform future work that aims to predict effects of mechanical processing, and use them to design material properties. Theoretical insights would be very valuable in this endeavour, because the multi-scale structure of gels means that many colloidal particles have to be simulated (or tracked) in order to model the behaviour of a single gel strand, while much of the important physics is happening on larger length scales, via properties of the gel network. This points towards the development of simplified meso-scale models, although there are significant couplings between microscopic structure and the responses of gel strands, which indicate that these theoretical models should not be too simple. In any case, we are optimistic that prediction of history-dependent properties can be improved beyond the current state-of-the-art, via a combined approach based on theory, experiment, and computer simulations.

Conflicts of interest

There are no conflicts to declare.

Data availability

The data in this article are available at https://doi.org/10.17863/CAM.119276. This includes the data that underlies the figures, the custom analysis scripts, and input parameters used to run molecular dynamics simulations via LAMMPS.

Acknowledgements

We thank Kris Thijssen, Abraham Mauleon-Amieva, Rui Cheng, Malcolm Faers, and Srikanth Sastry for helpful discussions. This work was supported by the EPSRC through grants EP/T031247/1 (HB and RLJ) and EP/T031077/1 (CPR and TBL). CPR acknowledges the grant DiViNew from the Agence National de Recherche. JEH and CPR also acknowledge the European Research Council (ERC consolidator grant NANOPRS, project 617266). HB thanks SUPRA (India) for support through grant SPR/2021/000382 during a part of the period when this work was performed.

References

1. E. Zaccarelli, Colloidal gels: equilibrium and non-equilibrium routes, J. Phys.: Condens. Matter, 2007, 19, 323101.
2. L. Cipelletti and L. Ramos, Slow dynamics in glassy soft matter, J. Phys.: Condens. Matter, 2005, 17, R253.
3. P. J. Lu, E. Zaccarelli, F. Ciulla, A. B. Schofield, F. Sciortino and D. A. Weitz, Gelation of particles with short-range attraction, Nature, 2008, 453, 499.
4. P. J. Lu and D. A. Weitz, Colloidal particles: Crystals, glasses, and gels, Annu. Rev. Condens. Matter Phys., 2013, 4, 217.
5. C. P. Royall, M. A. Faers, S. L. Fussell and J. E. Hallett, Real space analysis of colloidal gels: Triumphs, challenges and future directions, J. Phys.: Condens. Matter, 2021, 33, 453002.
6. W. Poon, The physics of a model colloid--polymer mixture, J. Phys.: Condens. Matter, 2002, 14, R859.
7. P. Segre, V. Prasad, A. B. Schofield and D. Weitz, Glasslike kinetic arrest at the colloidal-gelation transition, Phys. Rev. Lett., 2001, 86, 6042.
8. H. Sedgwick, K. Kroy, A. Salonen, M. Robertson, S. Egelhaaf and W. Poon, Non-equilibrium behavior of sticky colloidal particles: beads, clusters and gels, Eur. Phys. J. E: Soft Matter Biol. Phys., 2005, 16, 77.
9. B. J. Landrum, W. B. Russel and R. N. Zia, Delayed yield in colloidal gels: Creep, flow, and re-entrant solid regimes, J. Rheol., 2016, 60, 783.
10. H. Tsurusawa, M. Leocmach, J. Russo and H. Tanaka, Direct link between mechanical stability in gels and percolation of isostatic particles, Sci. Adv., 2019, 5, eaav6090.
11. J. Rouwhorst, C. Ness, S. Stoyanov, A. Zaccone and P. Schall, Nonequilibrium continuous phase transition in colloidal gelation with short-range attraction, Nat. Commun., 2020, 11, 3558.
12. C. Patrick Royall, S. R. Williams, T. Ohtsuka and H. Tanaka, Direct observation of a local structural mechanism for dynamic arrest, Nat. Mater., 2008, 7, 556.
13. J. M. van Doorn, J. E. Verweij, J. Sprakel and J. van der Gucht, Strand plasticity governs fatigue in colloidal gels, Phys. Rev. Lett., 2018, 120, 208005.
14. J. E. Verweij, F. A. Leermakers, J. Sprakel and J. Van Der Gucht, Plasticity in colloidal gel strands, Soft Matter, 2019, 15, 6447.
15. K. Thijssen, T. B. Liverpool, C. P. Royall and R. L. Jack, Necking and failure of a particulate gel strand: signatures of yielding on different length scales, Soft Matter, 2023, 19, 7412.
16. A. D. Dinsmore, V. Prasad, I. Y. Wong and D. A. Weitz, Microscopic structure and elasticity of weakly aggregated colloidal gels, Phys. Rev. Lett., 2006, 96, 185502.
17. L. C. Hsiao, R. S. Newman, S. C. Glotzer and M. J. Solomon, Role of isostaticity and load-bearing microstructure in the elasticity of yielded colloidal gels, Proc. Natl. Acad. Sci. U. S. A., 2012, 109, 16029.
18. M. Bouzid, J. Colombo, L. V. Barbosa and E. Del Gado, Elastically driven intermittent microscopic dynamics in soft solids, Nat. Commun., 2017, 8, 15846.
19. M. Nabizadeh, F. Nasirian, X. Li, Y. Saraswat, R. Waheibi, L. C. Hsiao, D. Bi, B. Ravandi and S. Jamali, Network physics of attractive colloidal gels: Resilience, rigidity, and phase diagram, Proc. Natl. Acad. Sci. U. S. A., 2024, 121, e2316394121.
20. D. Mangal, M. Nabizadeh and S. Jamali, Predicting yielding in attractive colloidal gels, Phys. Rev. E, 2024, 109, 014602.
21. S. M. Fielding, P. Sollich and M. E. Cates, Aging and rheology in soft materials, J. Rheol., 2000, 44, 323.
22. V. Trappe, V. Prasad, L. Cipelletti, P. Segre and D. A. Weitz, Jamming phase diagram for attractive particles, Nature, 2001, 411, 772.
23. L. Cipelletti, S. Manley, R. C. Ball and D. A. Weitz, Universal aging features in the restructuring of fractal colloidal gels, Phys. Rev. Lett., 2000, 84, 2275.
24. R. N. Zia, B. J. Landrum and W. B. Russel, A micro-mechanical study of coarsening and rheology of colloidal gels: Cage building, cage hopping, and smoluchowski's ratchet, J. Rheol., 2014, 58, 1121.
25. M. Nabizadeh and S. Jamali, Life and death of colloidal bonds control the rate-dependent rheology of gels, Nat. Commun., 2021, 12, 4274.
26. P. Bartlett, L. J. Teece and M. A. Faers, Sudden collapse of a colloidal gel, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2012, 85, 021404.
27. S. Patinet, D. Vandembroucq and M. L. Falk, Connecting local yield stresses with plastic activity in amorphous solids, Phys. Rev. Lett., 2016, 117, 045501.
28. J. T. Parley, S. M. Fielding and P. Sollich, Aging in a mean field elastoplastic model of amorphous solids, Phys. Fluids, 2020, 32, 127104.
29. J. Pollard and S. M. Fielding, Yielding, shear banding, and brittle failure of amorphous materials, Phys. Rev. Res., 2022, 4, 043037.
30. D. Bonn, M. M. Denn, L. Berthier, T. Divoux and S. Manneville, Yield stress materials in soft condensed matter, Rev. Mod. Phys., 2017, 89, 035005.
31. S. Aime, D. Truzzolillo, D. J. Pine, L. Ramos and L. Cipelletti, A unified state diagram for the yielding transition of soft colloids, Nat. Phys., 2023, 19, 1673.
32. H.-N. Lee, K. Paeng, S. F. Swallen and M. D. Ediger, Direct measurement of molecular mobility in actively deformed polymer glasses, Science, 2009, 323, 231.
33. C. Deng and C. A. Schuh, Atomistic mechanisms of cyclic hardening in metallic glass, Appl. Phys. Lett., 2012, 100, 251909.
34. P. Leishangthem, A. D. S. Parmar and S. Sastry, The yielding transition in amorphous solids under oscillatory shear deformation, Nat. Commun., 2017, 8, 14653.
35. S. Kim, I. Svetlizky, D. A. Weitz and F. Spaepen, Work hardening in colloidal crystals, Nature, 2024, 630, 648.
36. H. Mizuno, A. Ikeda, T. Kawasaki and K. Miyazaki, Universal mechanism of shear thinning in supercooled liquids, Commun. Phys., 2024, 7, 199.
37. T. Divoux, E. Agoritsas, S. Aime, C. Barentin, J.-L. Barrat, R. Benzi, L. Berthier, D. Bi, G. Biroli, D. Bonn, P. Bourrianne, M. Bouzid, E. Del Gado, H. Delanoë-Ayari, K. Farain, S. Fielding, M. Fuchs, J. van der Gucht, S. Henkes, M. Jalaal, Y. M. Joshi, A. Lematre, R. L. Leheny, S. Manneville, K. Martens, W. C. K. Poon, M. Popovic, I. Procaccia, L. Ramos, J. A. Richards, S. Rogers, S. Rossi, M. Sbragaglia, G. Tarjus, F. Toschi, V. Trappe, J. Vermant, M. Wyart, F. Zamponi and D. Zare, Ductile-to-brittle transition and yielding in soft amorphous materials: perspectives and open questions, Soft Matter, 2024, 20, 6868.
38. J. Mewis and N. J. Wagner, Colloidal Suspension Rheology, Cambridge Series in Chemical Engineering, Secondoftwo publisher Cambridge University Press, 2011.
39. T. Gisler, R. C. Ball and D. A. Weitz, Strain hardening of fractal colloidal gels, Phys. Rev. Lett., 1999, 82, 1064.
40. M. Pouzot, T. Nicolai, L. Benyahia and D. Durand, Strain hardening and fracture of heat-set fractal globular protein gels, J. Colloid Interface Sci., 2006, 293, 376.
41. K. M. Schmoller, P. Fernández, R. C. Arevalo, D. L. Blair and A. R. Bausch, Cyclic hardening in bundled actin networks, Nat. Commun., 2010, 1, 134.
42. G. de Oliveira Reis, T. Gibaud, B. Saint-Michel, S. Manneville, M. Leocmach, L. Vaysse, F. Bonfils, C. Sanchez and P. Menut, Irreversible hardening of a colloidal gel under shear: The smart response of natural rubber latex gels, J. Colloid Interface Sci., 2019, 539, 287.
43. R. G. Larson and Y. Wei, A review of thixotropy and its rheological modeling, J. Rheol., 2019, 63, 477.
44. J. Colombo and E. Del Gado, Stress localization, stiffening, and yielding in a model colloidal gel, J. Rheol., 2014, 58, 1089.
45. M. Bouzid and E. Del Gado, Network topology in soft gels: Hardening and softening materials, Langmuir, 2018, 34, 773.
46. H. A. Lockwood, M. H. Agar and S. M. Fielding, Power law creep and delayed failure of gels and fibrous materials under stress, Soft Matter, 2024, 20, 2474.
47. N. Koumakis and G. Petekidis, Two step yielding in attractive colloids: transition from gels to attractive glasses, Soft Matter, 2011, 7, 2456.
48. N. Koumakis, E. Moghimi, R. Besseling, W. C. Poon, J. F. Brady and G. Petekidis, Tuning colloidal gels by shear, Soft Matter, 2015, 11, 4640.
49. E. Moghimi, A. R. Jacob and G. Petekidis, Residual stresses in colloidal gels, Soft Matter, 2017, 13, 7824.
50. E. Moghimi, A. R. Jacob, N. Koumakis and G. Petekidis, Colloidal gels tuned by oscillatory shear, Soft Matter, 2017, 13, 2371.
51. B. Rajaram and A. Mohraz, Microstructural response of dilute colloidal gels to nonlinear shear deformation, Soft Matter, 2010, 6, 2246.
52. C. Perge, N. Taberlet, T. Gibaud and S. Manneville, Time dependence in large amplitude oscillatory shear: A rheo- ultrasonic study of fatigue dynamics in a colloidal gel, J. Rheol., 2014, 58, 1331.
53. V. Grenard, T. Divoux, N. Taberlet and S. Manneville, Timescales in creep and yielding of attractive gels, Soft Matter, 2014, 10, 1555.
54. M. Das and G. Petekidis, Shear induced tuning and memory effects in colloidal gels of rods and spheres, J. Chem. Phys., 2022, 157, 234902.
55. I. Sudreau, S. Manneville, M. Servel and T. Divoux, Shear-induced memory effects in boehmite gels, J. Rheol., 2021, 66, 91.
56. I. Sudreau, M. Auxois, M. Servel, É. Lécolier, S. Manneville and T. Divoux, Residual stresses and shear-induced overaging in boehmite gels, Physical Review, Materials, 2022, 6, L042601.
57. I. Sudreau, M. Servel, E. Freyssingeas, F. Liénard, S. Karpati, S. Parola, X. Jaurand, P.-Y. Dugas, L. Matthews, T. Gibaud, T. Divoux and S. Manneville, Shear-induced stiffening in boehmite gels: A rheo-x-ray-scattering study, Phys. Rev. Mater., 2023, 7, 115603.
58. D. Kushnir, C. Ruscher, E. Bartsch, F. Thalmann and P. Hébraud, Stress overshoot, hysteresis, and the Bauschinger effect in sheared dense colloidal suspensions, Phys. Rev. E, 2022, 106, 034611.
59. E. M. Schwen, M. Ramaswamy, C.-M. Cheng, L. Jan and I. Cohen, Embedding orthogonal memories in a colloidal gel through oscillatory shear, Soft Matter, 2020, 16, 3746.
60. N. C. Keim and S. R. Nagel, Generic transient memory formation in disordered systems with noise, Phys. Rev. Lett., 2011, 107, 010603.
61. D. Fiocco, G. Foffi and S. Sastry, Encoding of memory in sheared amorphous solids, Phys. Rev. Lett., 2014, 112, 025702.
62. P. Edera, M. Bantawa, S. Aime, R. T. Bonnecaze and M. Cloitre, Mechanical tuning of residual stress, memory, and aging in soft glassy materials, Phys. Rev. X, 2025, 15, 011043.
63. Z. Varga and J. W. Swan, Large scale anisotropies in sheared colloidal gels, J. Rheol., 2018, 62, 405.
64. S. Jamali, R. C. Armstrong and G. H. McKinley, Time-rate-transformation framework for targeted assembly of short-range attractive colloidal suspensions, Mater. Today Adv., 2020, 5, 100026.
65. C. P. Royall, P. Charbonneau, M. Dijkstra, J. Russo, F. Smallenburg, T. Speck and C. Valeriani, Colloidal hard spheres: Triumphs, challenges and mysteries, Rev. Mod. Phys., 2024, 96, 045003.
66. R. Besseling, E. R. Weeks, A. B. Schofield and W. C. Poon, Three-dimensional imaging of colloidal glasses under steady shear, Phys. Rev. Lett., 2007, 99, 028301.
67. P. Schall, D. A. Weitz and F. Spaepen, Structural rearrangements that govern flow in colloidal glasses, Science, 2007, 318, 1895.
68. P. Ballesta, N. Koumakis, R. Besseling, W. C. K. Poon and G. Petekidis, Slip of gels in colloid--polymer mixtures under shear, Soft Matter, 2013, 9, 3237.
69. J. C. Conrad and J. A. Lewis, Structure of colloidal gels during microchannel flow, Langmuir, 2008, 24, 7628.
70. S. M. Fenton, P. Padmanabhan, B. K. Ryu, T. T. D. Nguyen, R. N. Zia and M. E. Helgeson, Minimal conditions for solidification and thermal processing of colloidal gels, Proc. Natl. Acad. Sci. U. S. A., 2023, 120, e2215922120.
71. H. Bhaumik, G. Foffi and S. Sastry, The role of annealing in determining the yielding behavior of glasses under cyclic shear deformation, Proc. Natl. Acad. Sci. U. S. A., 2021, 118, e2100227118.
72. A. Malins, S. R. Williams, J. Eggers and C. P. Royall, Identification of structure in condensed matter with the topological cluster classification, J. Chem. Phys., 2013, 139, 234506.
73. C. P. Royall, S. R. Williams and H. Tanaka, Vitrification and gelation in sticky spheres, J. Chem. Phys., 2018, 148, 044501.
74. B. Keshavarz, D. G. Rodrigues, J.-B. Champenois, M. G. Frith, J. Ilavsky, M. Geri, T. Divoux, G. H. McKinley and A. Poulesquen, Time-connectivity superposition and the gel/glass duality of weak colloidal gels, Proc. Natl. Acad. Sci. U. S. A., 2021, 118, e2022339118.
75. M. G. Noro and D. Frenkel, Extended corresponding-states behavior for particles with variable range attractions, J. Chem. Phys., 2000, 113, 2941.
76. A. I. Campbell and P. Bartlett, Fluorescent hard-sphere colloids for confocal microscopy, J. Colloid Interface Sci., 2002, 256, 325.
77. S. Asakura and F. Oosawa, On interaction between 2 bodies immersed in a solution of macromolecules, J. Chem. Phys., 1954, 22, 1255.
78. R. Verma, J. Crocker, T. Lubensky and A. Yodh, Entropic colloidal interactions in concentrated dna solutions, Phys. Rev. Lett., 1998, 81, 4004.
79. C. P. Royall, A. A. Louis and H. Tanaka, Measuring colloidal interactions with confocal microscopy, J. Chem. Phys., 2007, 127, 044507.
80. J. Taffs, A. Malins, S. R. Williams and C. P. Royall, A structural comparison of models of colloid-polymer mixtures, J. Phys.: Condens. Matter, 2010, 22, 104119.
81. A. Razali, C. J. Fullerton, F. Turci, J. E. Hallett, R. L. Jack and C. P. Royall, Effects of vertical confinement on gelation and sedimentation of colloids, Soft Matter, 2017, 13, 3230.
82. S. Griffiths, F. Turci and C. P. Royall, Local structure of percolating gels at very low volume fractions, J. Chem. Phys., 2017, 146, 014905.
83. A. P. Thompson, H. M. Aktulga, R. Berger, D. S. Bolintineanu, W. M. Brown, P. S. Crozier, P. J. Int Veld, A. Kohlmeyer, S. G. Moore, T. D. Nguyen, R. Shan, M. J. Stevens, J. Tranchida, C. Trott and S. J. Plimpton, LAMMPS - a flexible simulation tool for particle-based materials modeling at the atomic, meso, and continuum scales, Comput. Phys. Commun., 2022, 271, 108171.
84. H. Bhaumik, T. B. Liverpool, C. P. Royall and R. L. Jack, Yielding in colloidal gels: From local structure to mesoscale strand breakage and macroscopic failure, Phys. Rev. E, 2025, 111, 055412.
85. C. P. Royall, J. Eggers, A. Furukawa and H. Tanaka, Probing colloidal gels at multiple length scales: The role of hydrodynamics, Phys. Rev. Lett., 2015, 114, 258302.
86. Z. Varga, G. Wang and J. Swan, The hydrodynamics of colloidal gelation, Soft Matter, 2015, 11, 9009.
87. Z. Varga and J. Swan, Hydrodynamic interactions enhance gelation in dispersions of colloids with short-ranged attraction and long-ranged repulsion, Soft Matter, 2016, 12, 7670.
88. J. de Graaf, W. C. K. Poon, M. J. Haughey and M. Hermes, Hydrodynamics strongly affect the dynamics of colloidal gelation but not gel structure, Soft Matter, 2019, 15, 10.
89. T. A. Vezirov, S. Gerloff and S. H. L. Klapp, Manipulating shear-induced non-equilibrium transitions in colloidal films by feedback control, Soft Matter, 2015, 11, 406.
90. R. Cabriolu, J. Horbach, P. Chaudhuri and K. Martens, Precursors of fluidisation in the creep response of a soft glass, Soft Matter, 2019, 15, 415.
91. R. H. Ewoldt and G. H. McKinley, Creep ringing in rheometry, Rheol. Bull., 2007, 76, 4.
92. P. J. Steinhardt, D. R. Nelson and M. Ronchetti, Bond-orientational order in liquids and glasses, Phys. Rev. B: Condens. Matter Mater. Phys., 1983, 28, 784.
93. V. V. Vasisht, G. Roberts and E. Del Gado, Emergence and persistence of flow inhomogeneities in the yielding and fluidization of dense soft solids, Phys. Rev. E, 2020, 102, 010604.
94. C. P. Royall and S. R. Williams, The role of local structure in dynamical arrest, Phys. Rep., 2015, 560, 1.
95. C. P. Royall and A. Malins, The role of quench rate in colloidal gels, Faraday Disc., 2012, 158, 301.
96. L. D. Gelb and K. E. Gubbins, Pore size distributions in porous glasses: A computer simulation study, Langmuir, 1999, 15, 305.
97. M. Kamkar, R. Salehiyan, T. B. Goudoulas, M. Abbasi, C. Saengow, E. Erfanian, S. Sadeghi, G. Natale, S. A. Rogers, A. J. Giacomin and U. Sundararaj, Large amplitude oscillatory shear flow: Microstructural assessment of polymeric systems, Prog. Polym. Sci., 2022, 132, 101580.
98. K. Hyun, M. Wilhelm, C. Klein, K. Cho, J. Nam, K. Ahn, S. Lee, R. Ewoldt and G. McKinley, A review of nonlinear oscillatory shear tests: Analysis and application of large amplitude oscillatory shear (LAOS), Prog. Polym. Sci., 2011, 36, 1697.
99. W.-T. Yeh, M. Ozawa, K. Miyazaki, T. Kawasaki and L. Berthier, Glass stability changes the nature of yielding under oscillatory shear, Phys. Rev. Lett., 2020, 124, 225502.
100. N. C. Keim, J. D. Paulsen, Z. Zeravcic, S. Sastry and S. R. Nagel, Memory formation in matter, Rev. Mod. Phys., 2019, 91, 035002.
101. S. Maity, H. Bhaumik, S. Athani and S. Sastry, Fatigue failure in glasses under cyclic shear deformation, arXiv, 2024, preprint, arXiv:2409.17384 DOI:10.48550/arXiv.2409.17384.
102. B. P. Bhowmik, H. G. E. Hentschel and I. Procaccia, Scaling theory for wöhler plots in amorphous solids under cyclic forcing, Phys. Rev. E, 2022, 105, 015001.

---

Footnotes

† Each loop is averaged over 10 cycles, they are not perfectly symmetric under inversion through the origin because the system is evolving structurally throughout the shearing process.

‡ We verified that similar results are obtained by averaging the time-dependent stress over the samples and then fitting the average.

§ See ref. 84 for further discussion of the σ₀-dependence: we would expect yielding to be similar at smaller σ₀, if sufficiently long simulations were performed. The behaviour for larger σ₀ is different because the system yields almost immediately.

---