We study the linear response rheology, structure, and dynamics of colloidal gels formed by arrested phase separation with a combination of experiments and dynamic simulation, with a view toward understanding the influence of bond strength, volume fraction, and network morphology on the viscoelastic moduli. A rescaling of the data to remove the direct, equilibrium hydrodynamic, and entropic effects enables a direct comparison of experiment and simulation; the strong agreement shows that attractive forces and Brownian motion dominate relaxation, where hydrodynamic interactions play a simpler role that can be scaled out. Morphology transitions from thick, blobby strands with large solvent pores to particle-size solvent pores surrounded by concave surfaces when volume fraction increases. Unsurprisingly, generalized Stokes–Einstein relations make a poor predictor of rheology from particle dynamics. Models that connect bond dynamics to elasticity or that connect cluster dynamics to elasticity show better agreement. Prediction of age-stiffening requires a model that includes the effects of age-coarsening; surprisingly, age-stiffening is set primarily at high frequency by the dominant network length scale. A Rouse-like theory that connects dominant network length scale to elasticity provides good predictions for all volume fractions and ages, although there is an interplay between volume fraction and structural length scales. The linear viscoelastic response of the experimental system is thus well represented in a simpler computational model in which wall effects, hydrodynamics, explicit depletant molecules, and rejuvenation were neglected, suggesting that the connections between morphology, dynamics, and rheology are encoded primarily by volume fraction, attraction strength, Brownian motion, and age.

## I. INTRODUCTION

Complex fluids encompass a wide range of soft matter, many of which are composed of a multiphase structure—a microstructure—that deforms and relaxes over observable time scales. The dynamics of this internal microstructure impart a spectrum of relaxation time scales which in turn influence mechanical properties, including behaviors that are a hallmark of non-Newtonian fluids: They display both liquidlike (viscous) and solidlike (elastic) behavior, depending on the rate with which they are perturbed relative to the relaxation time scales of the microstructure. Linear viscoelastic properties can be studied by imposing a small-amplitude oscillatory shear (SAOS) flow on the material, $\gamma (t)=\gamma 0sin(\omega t)$, with amplitude $\gamma 0$ and oscillation frequency $\omega $. In this linear-response regime, the resultant shear stress $\sigma (t)$ is linear in the imposed strain $\gamma (t)$ and strain rate $\gamma \u02d9(t)$, with coefficients that form the real and imaginary parts of a complex modulus, $G\u2032(\omega )$ and $G\u2032\u2032(\omega )$. The relative strengths of hydrodynamic, entropic, and other interparticle forces play a central role in the relaxation behavior, and this evolving dominance manifests in the frequency-dependent viscoelastic moduli. This viscoelastic behavior is found even in systems as simple as a colloidal dispersion of purely repulsive hard spheres. Such interrogation decouples the roles played by hydrodynamic and entropic forces [1–3]. The diffusive and advective transport that underlies such relaxation is sensitive to changes in particle volume fraction or interparticle forces; these in turn can play a qualitative or quantitative role in frequency-dependent material properties, either directly by hindering relative particle flux or indirectly by producing new relaxation length scales in the material structure.

Colloidal gels are an important example of soft matter that exhibit such effects [4]. Details of the interparticle interaction potential, such as attraction strength and range, influence the gel formation and its morphology during gelation and structural evolution that continues post-gelation. Two main scenarios are thus foreseen, one leading to percolation gels and one triggered by phase separation [2,5–7].

In some model gels, the presence of attractive interparticle forces promotes phase separation but simultaneously hinders the diffusive motion required for phase separation to proceed. This mechanism leads arrested phase separation [2,5,6,8], often encountered in colloid-polymer gels [9–11] such as the ones studied in the present work. In this case, a bicontinuous network of thick, glassy strands—where average coordination number of particles is six or higher—is formed where a hierarchy of length scales frozen into the network produces a spectrum of relaxation time scales [8,12,13].

When the interparticle bond strength is only several $kT$, where $k$ is Boltzmann’s constant and $T$ the absolute temperature, diffusive transport continues to restructure the network post-gelation, but particle migration is slowed down by both attractive forces and by steric hindrance due to the concentration of colloids [6–8,12,13]. Ongoing structural coarsening takes place, characterized by increasing particle coordination number and dominant network length scale [12,14] as well as slowing particle dynamics. Together, these have been shown to produce age-stiffening of the linear viscoelastic moduli [12,15].

In contrast, under the appropriate conditions, “percolation gels” may be formed usually with a more open network structure stabilized by strong attractive interactions in a regime of the phase diagram where no thermodynamic phase separation is expected [16–18]. Such gelation is favored when the interparticle potential $V0$ is strong, $V0\u226bkT$, and the initial volume fraction $\varphi =4\pi a3n/3$, is low, $\varphi \u226a1$, where $a$ is the particle radius and $n$ is the number density. This results in the formation of an elastic network with a morphology and structure that can be fractal or stringy, where bonds are essentially permanent, the average coordination number is low (two or three) and no aging or structural coarsening is observed [7,19–21]. Post-gelation relaxation processes are primarily elastic, with the elastic modulus $G\u2032(\omega )$ exhibiting a greater magnitude than $G\u2032\u2032(\omega )$ over intermediate and high frequencies. Theoretical and computational models have been devised to connect structure to rheology in cases where surface-bonding potentials restrict tangential particle pair motion [22–24] or utilizing three-body potentials that provide “angular rigidity” [20] and are favored for their amenity of modeling biopolymer, protein, or casein gels, as well as for their relative computational simplicity in modeling gel formation.

Despite these interesting results, there has been relatively little study devoted to the detailed interrogation of the linear-response rheology of “reversible” colloidal gels, where bond strength is less than about $10kT$, perhaps owing in part to the only recently emergent recognition of the dramatic changes in morphology that accompany seemingly small changes in interparticle potential and volume fraction. A recent review by Zaccarelli provides a lucid recognition of the richness in colloidal gel morphology [5]. Moreover, on the experimental side, as has been seen in both model colloidal gels [15,25,26] and more industrially relevant systems [27,28], preshear history affects significantly the structure and, therefore, the linear viscoelastic response of such out-of-equilibrium states.

The motivation for deeper understanding of the structure–property relationships has also emerged from numerous recent works focused on the nonlinear response rheology of dense, reversible gels. This includes yield-stress behavior under fixed stress [29–31], fixed strain rate [15,25,32], large-amplitude oscillatory shear (LAOS) [15,26,30,33], and delayed gravitational collapse [34,35]. For example, examination of the effect of colloid concentration and interparticle bond strength on start-up flow and yield behavior has revealed that both can act to alter the yield stress, yield strain, and stress overshoots [15,30,32]. However, decoupling the two effects has proven challenging, owing in part to their hierarchy of length scales and dense morphologies and the difficulty of performing detailed structural measurements in dense, bicontinuous gels, via confocal microscopy or scattering techniques, especially under shear where several technical challenges involving slip or shear banding, thixotropy, and shear-induced gravitational collapse have to be addressed [25,30,36]. In contrast, linear-response (small-amplitude) oscillatory interrogation provides a well-characterized technique for decoupling rate processes in such materials.

Experimental studies of depletion gel rheology have revealed that the elastic moduli increase with increasing depletion attraction strength over a wide range of frequencies [13,30,37]. Moreover, similarly to other soft yield stress systems with solidlike response (such as hard sphere glasses, jammed soft particles, and emulsions [38–40]), the elastic modulus is weakly frequency dependent and higher than the viscous modulus. The latter may exhibit a minimum at intermediate frequencies [13,30,41] related to the transition from some local short-time motion to an ultra long-time structural relaxation. The first of these effects is amenable to phenomenological models connecting the moduli to the attractive potential [30,42]. Changes in depletant concentration and size that induce variation in bond strength and range have been incorporated via free volume theories [43]. The effect of the strength and range of attraction on the linear viscoelasticity was studied in concentrated depletion gels [15]. The elastic and viscous moduli both increase with increasing volume fraction of colloids as described by the combined mode-coupling theory (MCT)/polymer reference interaction site model (PRISM) predictions by Chen and Schweizer [37], if renormalized by a value corresponding to the average number of particles in a cluster. Nevertheless, detailed structural measurements to support this view remain elusive. Advances in imaging techniques for the study of particle-scale structure that can interrogate the gel in quiescent conditions and under shear show great promise in answering these questions.

To this end, rheo-confocal measurements combining confocal microscopy and rheometry [25] of an intermediate volume fraction ($\varphi =0.44$) colloid-polymer depletion gel under steady shear have revealed the effect of preshear history on the directly measured microstructure and the time evolution of the linear viscoelasticity. While extremely strong shear flow breaks down gel structure completely to single particles resulting in the formation of less heterogeneous and stronger gels upon flow cessation, weak preshear causes aggregation of clusters, leading eventually to heterogeneous and weaker gels upon flow cessation. These experimental findings are congruent with Brownian dynamics simulations performed under the same preshear conditions [25]. Similar, yet even stronger, effects are detected when such colloidal gels are presheared by LAOS; intermediate strain amplitudes are more efficient in creating weaker and more heterogeneous gels [26]. Still, the detailed relation between the structure and the mechanical properties, linear and nonlinear, of colloidal gels is not fully decoded, although several aspects have indeed been widely explored.

Simulation studies readily provide the particle- and network-scale structure and dynamics enabling creation of structure–property relationships. One such study by Zia and co-workers [12] of reversibly bonded gel of volume fraction $\varphi =0.20$ provided detailed measurements of particle-scale dynamics which showed that such reversible colloidal gels coarsen under quiescent conditions via surface migration of weakly to moderately bonded particles which gain more bonds with time. The resulting network length scale is the single measurement needed to predict the change in magnitude of the age-stiffening of the high-frequency elastic and viscous moduli [12]. In the high frequency limit, the largest length scale dominates the relaxation response in polymeric systems [44,45], which Zia and co-workers extended to colloidal gels where these large length scales are formed by particles bonded into the network structure with a characteristic size. This idea will play an important role in the present study. While this study was conducted with a freely draining system, hydrodynamic interactions are expected only to slow the rate at which particles migrate and slow the overall coarsening of the network. The simulation studies of Zia and co-workers demonstrate the utility of pairing detailed structural measurements with rheology, but direct comparison to experimental measurement and extension to higher volume fraction gels are necessary to validate the presented model.

In the present work, we conduct a comprehensive experimental and computational study of the linear-response rheology of reversible colloidal gels, with a focus on the role played by particle volume fraction, morphology, and interparticle bond strength on the linear viscoelasticity and the relaxation response of the gel. In experiments, gels are formed via depletion interactions by the addition of a nonadsorbing polymer. The computational model comprises of freely draining Brownian hard spheres interacting via an interparticle potential matched to the strength and range of particle attractions realized in experiments. A SAOS is imposed over a range of frequencies, and the stress monitored throughout the oscillation. From it, the elastic and viscous moduli are extracted for both experiment and simulation. We compare experimental to simulation results for the effects of volume fraction, gel age, and bond strength on the moduli. Toward explaining this behavior, the detailed particle-scale information available in simulation is leveraged to discover the frequency-dependent relaxation processes in reversible gels and how these evolve with the parameter space studied.

The remainder of this paper is organized as follows: In Sec. II, we present the experimental and simulation systems and techniques, including a review of the effects of preparation protocol of the experimental system. In Sec. III, we present results from experiments and simulations and provide a direct comparison of findings as a function of volume fraction, gel age, and attraction strength. Section IV presents the use of a modified generalized Stokes–Einstein relation to approximate the linear viscoelastic response from the microscopic particle dynamics, a scaling argument to model aging effects, and finally a critical assessment of existing models for the volume fraction and attraction strength dependence of the elastic modulus. In Sec. V, we present the conclusions of this work.

## II. SYSTEMS AND METHODS

Here, we present the methods by which the gels were formed, aged, and tested.

### A. Experimental system

We used sterically stabilized polymethylmethacrylate (PMMA) model hard-sphere particles sterically stabilized by chemically grafted poly-hydroxy-stearic acid (PHSA) chains ($\u223c10nm$) suspended in octadecene. Particles have a hydrodynamic radius of $a=196nm$ with a polydispersity of 12% to suppress crystallization. The depletion attractions are implemented between the particles by adding nonadsorbing linear polybutadiene(1,4-addition) (PB) chains with the molecular weight $Mw=323300g/mol$ and a radius of gyration $Rg=19nm$ (measured by light scattering). This gives a polymer-colloid size ratio $\Delta =Rg/a=0.097$ in dilute solution. We prepared the gel at particle volume fractions of $\varphi =0.20$, $0.30$, and $0.44$ with two different attraction strengths at contact $V0=6kT$ and $15kT$ according to the modified Asakura–Oosawa (AO) model [43,46].

The Brownian time scale in experiments is calculated as $a2/D=6\pi \eta sa3/kT=0.15s$ with the solvent viscosity $\eta s=4.2mPas$, assuming diffusion in the dilute limit (with a free diffusion coefficient, $D=kT/6\pi \eta sa$). However, a Brownian time can also be computed using the viscosity of the solvent and depletant (giving an effective viscosity of the medium $\eta eff$) utilized in the $6kT$ depletion gel, $6\pi \eta effa3/kT=0.34s$, more than twice than that from the lone colloid diffusivity in solvent alone. For the $15kT$ bond strength depletion gel, a Brownian time becomes 0.90 s. Thus, we take these longer Brownian times to normalize sensitive dynamical information (e.g., frequency, the shift provided by this normalization is discussed in Sec. II C 1) and introduce $D\u2217=kT/6\pi \eta effa$ to form the time scale $a2/D\u2217$ representing the diffusion of a lone colloid in the suspending medium considering both solvent and depletant in experiment. Although the new experimental Brownian time shifts the corresponding gel age examined in the present work to younger ages, we discuss in detail the validity of such normalization where age behavior is examined in Sec. III B.

### B. Simulation system

We consider a freely draining colloidal suspension of 750 000 spherical particles of characteristic size $a$ dispersed in an implicit Newtonian solvent of density $\rho $ and viscosity $\eta $. The Reynolds number formed by the relative strength of the fluid inertia to viscous forces, $Re=\rho Ua/\eta $, where $U$ is the characteristic velocity, is thus small due to the small size of the particles, permitting the use of the Stokes equations to evolve the fluid motion. The particles are neutrally buoyant and initially dispersed at a moderate concentration with particle volume fractions $\varphi \u22614\pi a3n/3=0.2,0.3,and0.44$ to match the parallel experimental study, where $n$ is the number density of the colloidal particles. An interparticle potential provides a means of representing the hard-sphere repulsion and short ranged attraction typical of depletion gels both here and broadly utilized in experimental study [8,15,30,47,48]. Here, we model this via the Morse potential,

where the interaction between particles $i$ and $j$ is a function of the center-to-center separation $rij$ of particles size $ai$ and $aj$, respectively. We select the depth of the potential $V0$ [which is reached at $r=2a$, Fig. 1(b)] to range $6kT\u2264V0\u226415kT$, allowing the study of gels with bond strength only a few times $kT$ to far stronger gels. The range of attraction is set by $\kappa =30/a$ to model a similarly short-ranged AO [46] depletion potential of a given strength [Fig. 1(b)]. Thus, the effective range of the attraction relative to the particle size is $\Delta =0.1$ as utilized in the parallel experimental study and elsewhere in the literature [15,26,49]. Polydispersity of particle sizes of 7% is introduced to suppress the formation of crystalline regions [50], which is lower than, but comparable to, experimental systems.

A periodically replicated system of randomly distributed particles is initially dispersed. The attractive interaction is then introduced, utilizing the LAMMPS molecular dynamics package [51] to form and evolve a gel in time. The gel is aged for a prescribed time in diffusive times which we define as $a2/D$, where $D$ is the Brownian diffusion of a lone colloid in pure solvent. In simulations, there is no need to differentiate diffusivity defined based on solvent viscosity alone or viscosity of solvent and depletant because attractions are implicit and the suspending medium is pure solvent. However, for notational consistency the oscillation frequency will be normalized with $a2/D\u2217$ as for experiments and gel age or elapsed time with $a2/D$ (Sec. II A). We showed in a previous study of colloidal gels with $\varphi =0.20$ and $5\u2264V0/kT\u22646$ that the selected system size is sufficient to recover formation of a reversible colloidal gel characterized by thick, glassy strands and, further, to obtain appropriate sampling of a range of length scales for long-time quiescent aging [12]. In the present study, we select a gel at a specific age from this family of dense gels with reversible to strong bonds and examine structure, dynamics, and linear-response rheology.

### C. Experimental methods

#### 1. Rheometry

Rheological experiments were performed with an ARES-HR strain controlled rheometer with a force balance transducer using homemade cone-plate geometries of diameter $40mm$, cone angle $2.3\xb0$, and gap size $0.04mm$ with serrated surfaces to avoid wall slip [25,30,36]. The temperature was set to $T=20\xb0C$ using a standard Peltier plate, and solvent evaporation was minimized by using a solvent saturation trap. As a standard protocol before each experiment, the gel was first rejuvenated by shear at a large strain amplitude of $\gamma 0=800%$ with frequency $\omega =10rad/s$ until reaching steady state values of $G\u2032$ and $G\u2033$ under shear. This protocol has been shown to lead to a reproducible rejuvenation of the gel [26] as discussed in Sec. II C 2. Subsequently, shear was stopped and the gel was left to age for a specific time of 400–20 000 $a2/D$. Following quiescent aging, we performed dynamic frequency sweep (DFS) experiments in the linear regime to determine the viscoelastic spectrum of the sample, in a typical frequency range, $0.066\u2264\omega \u226466.7rad/s$, at different aging times. Normalization of the frequency by the free (or bare) Brownian time of a lone particle in solvent gives the dimensionless frequency $\omega a2/D$, denoted also as $Pe\omega $ in some of our previous works [15,26] that, in this study, takes values in the range $0.01\u2264\omega a2/D\u226410$. However, because particle diffusion slows down with increasing depletant concentration, it is appropriate to scale the frequency by the corresponding self-diffusion coefficient, $D\u2217$, shifting the range to $0.024\u2264\omega a2/D\u2217\u226424$ for the $6kT$ bond strength gel.

#### 2. Shear rejuvenation protocol

In simulations, the starting point for aging is a quench from an equilibrated liquid state, by switching on the attractions between particles at $t^w=0$, the waiting time $t^w=tw/(a2/D)$ since the quench. In experiments, however, the gel begins to age after a rejuvenation process, where the sample loaded into the rheometer undergoes oscillatory shear at large strain amplitudes, here at $\gamma 0=800%$ and a frequency of $\omega =10rad/s$, under the assumption that bonds break during such a process, “fluidizing” any pre-existing structure and returning it to a dispersed configuration. To examine if shear flow under such conditions is strong enough to fully break the bonds between the particles, we use a modified Péclet number (where a simple Péclet number defined on the bare diffusivity is $Pe\u2261\gamma \u02d9a2/D0$) to determine the relative strength of shear forces $Fshear$ to that of the attractive forces due to depletion interactions $Fdep$, defined as [15,26]

For values of $Pedep>1$, shear forces are stronger than attractive forces and substantial net bond loss under shear is expected [15,25,26]. Conversely, for $Pedep<1$, attractive forces are dominant and shear only serves to assist particles to find a deeper local minimum of the energy landscape, which leads to bond formation and cluster densification under shear [15,31,32]. This range of responses has been verified under steady shear by rheo-confocal experiments and Brownian dynamics simulations [25], while similar findings have also been documented under LAOS [26]. Therefore, shear rejuvenation protocols demand that $Pedep>1$. Here, all of our experiments were performed by oscillatory shear rejuvenation at strain amplitude of $\gamma 0pre=800%$ and frequency $\omega =10rad/s$. For the stronger gel with $V0=15kT$, under such conditions we have $Pedep=0.15$, using the short-time self-diffusion of a dilute suspension. However, the fact that the local volume fraction in a cluster is approximately $\varphi =0.62$ the true short-time self-diffusion is approximately $10$ times larger leading to $Pedep>1$ [25]. This means that the shear rejuvenation protocol used should be appropriate to fluidize the gel structure, allowing it to take on a dispersed configuration. However, the simulated gels are “pristine,” i.e., no preshear protocol is imposed, leading to the concern that the different histories may change initial gel morphology and influence relaxation dynamics. Thus, the comparison of presheared and pristine gels cannot be exact, but one can justify the comparison as valid while acknowledging that some differences could still exist. Validity of the comparison is supported by Petekidis and co-workers who showed that preshear [26,52] is one way of preparing a sample with reproducible gel structure and “age” in experimental measurements. These prior studies (combination of experiments and Brownian dynamics simulations) support the claim that following strong steady shear, stress relaxes to zero; and that following strong oscillatory or steady shear, one recovers the same void volume and bond distribution of a pristine gel. Broadly, one may preshear the gel at any number of flow strengths (for steady shear) or frequencies and strain amplitudes (for oscillatory shear). Selecting an adequately high flow strength, strain rate, or strain amplitude was shown to evolve a gel to a flowing state that structurally resembles a liquid. Thus, we suggest “rejuvenation” of a gel via strong shear can minimize the differences between the pristine simulation gel and the preshear depletion gel.

We first subjected a gel to LAOS for $1000a2/D$ [left of the dashed line in Fig. 2(a)], taking the first harmonic as an approximate value for $G\u2032$ and $G\u2033$. As soon as the LAOS test ended, we immediately began a new SAOS test, measuring $G\u2032$ and $G\u2033$ [right of the vertical dashed line in Fig. 2(a)]. Two pretest LAOS amplitudes were tested, $\gamma 0pre=25%$ and $\gamma 0pre=800%$. The first-harmonic approximation shows a lower value for higher amplitude shear, suggesting that the gel network has been ripped apart more efficiently during this test. After the pretest LAOS ended, the modulus from the small-amplitude shear shows a moderate increase, while that for the large-amplitude preshear shows two orders of magnitude increase, which may only reflect the transition from fluidlike flow during rejuvenation to solidlike behavior during linear-response interrogation. The higher value of the linear $G\u2032$ (determined by SAOS) for the larger preshear amplitude compared to that after a smaller preshear amplitude is due to the creation of a more homogeneous gel with more abundant small length scales in the former case, as seen before [25,26].

The resultant linear viscoelastic response from DFSs at $4000a2/D$ after flow cessation of the rejuvenation protocol at different strain amplitudes is shown in Fig. 2(b). The gel presheared at a large strain amplitude of $\gamma 0pre=800%$ has larger elastic modulus over the range of frequencies shown compared to the one presheared at intermediate strain amplitude of $\gamma 0pre=25%$.

### D. Simulation methods

#### 1. Dynamic simulation

The motion of colloidal particles in the Stokes regime is dictated by Brownian forces $FB$, hydrodynamic $FH$, and interparticle forces $FP$. The Langevin equation describes how these forces act to evolve the colloidal system in time,

where $m$ is the mass or moment of inertia tensor and $U$ is the particle velocity. Brownian forces, which arise physically from the bombardment of the solvent molecules on the surface of the colloidal particle, are applied as instantaneously correlated kicks which obey Gaussian statistics,

where $I$ the identity tensor, $\delta (t)$ is the Dirac delta function, and the overbar denotes an average over times long compared to time scale of the solvent. Here, we consider a freely draining system where the hydrodynamic force is proportional to motion relative to the bulk flow $\u27e8u\u27e9$ and thus the force on particle $i$ with velocity $Ui$ is given as

The interparticle force here arises from pairwise interactions between particles, dictated by the family of Morse potentials described in Sec. II B, and is thus computed as

where $r^ij$ is a unit vector along the line of centers of particles $i$ and $j$.

The system is evolved in time utilizing the LAMMPS molecular dynamics software package [51] where the velocity Verlet integration method [53] integrates the velocity forward in time and a Langevin “thermostat” enforces Brownian statistics [54]. While particle inertia is not identically zero as in Brownian dynamics simulations [55–59], the Stokes number is minimized to reduce the influence of inertia [left hand side of Eq. (3)] and achieve Stokesian physics by setting the time step to a sufficiently small size [60], as described in our previous work [12].

A macroscopic measure of interest here is the stress tensor, which for nonhydrodynamically interacting particles reduces to the ideal osmotic pressure and the elastic stress,

where $x$ is the location of particle being acted upon by $FP$. The stress tensor may be averaged over times long enough to obtain a clear signal but short enough to not smear out the effects of ongoing aging.

A SAOS, $\gamma (t)=\gamma 0sin\u2061(\omega t)$, is applied to a selected gel, and the resulting shear stress $\sigma (t)$ is monitored for a sufficient number of cycles. To measure linear response properties, we ensure that the dimensionless displacement, $Pe/(\omega a2/D)$, is small as defined by the ratio of the dimensionless shear rate (the Péclet number, $Pe\u2261\gamma \u02d9a2/D$). We define both these quantities to distinguish the linear-response regime from the LAOS discussed above for the preshear protocol. We then extract the linear viscoelastic moduli as the shear stress $\sigma (t)$ which is linear in the strain to obtain the elastic modulus $G\u2032$ and that which is linear in the strain rate to obtain the viscous modulus $G\u2033$. The range of frequencies explored in this work in simulation is $0.1\u2264\omega a2/D\u2264100$.

#### 2. Structural and dynamical measurements

Particle positions and motion are tracked throughout simulation to enable the measurement of structure and dynamics. The hierarchy of length scales in a colloidal gel can be probed in wave space via the static structure factor $S(q)$, where one expects an abundance of longer length scales to emerge upon gelation and to continue to grow in abundance and size [12]. It is defined as a correlation in particle positions

where $q$ is the wave vector normalized by the particle size $a$, $N$ is the number of particles, and $xj\u2212xk$ is the separation between a particle pair. The structure factor can be radially averaged to give $S(q)$, where $q=|q|$. We define the most abundant length scale, which is many particle sizes, as $LS(q)\u22612\pi /(qmaxa)$, where $qmax$ is the corresponding wavenumber to the low-$q$ peak in $S(q)$ [12].

A complementary measurement in physical space is the void volume $Vvoid$ which instead computes the size of spaces devoid of particles. $Vvoid$ is defined as the volume of a sphere, located in “empty” space, with a radius equal to the nearest particle surface, as defined and used in our prior works [25,26]. To calculate $Vvoid$, the simulation box is subdivided into cubes of size $dV<a3$. Probing the empty space gives the void volume probability distribution $P(Vvoid)$ and provides complementary information on the structural heterogeneity on all length scales. With this method, a specific portion of the pore space may be counted several times, robustly mapping the heterogeneous distribution of void sizes [25]. The dominant or average pore diameter $d\u27e8VV\u27e9$ is the cube root of the average void volume,

where $\u27e8VV\u27e9$ is the average volume of voids computed from *P*($Vvoid$).

Particle dynamics are monitored by tracking particle positions as they change throughout simulation, where high-resolution (small time steps) enables measurement of mean-square and average displacements. The current particle position $r\u2032(t)$ is given by the difference between the absolute position of the particle and the ensemble average of all particle positions $r\u2032(t)\u2261r(t)\u2212\u27e8r(t)\u27e9$. The mean-square displacement is thus the ensemble average of the square of the current position minus the initial position [$(r\u2032(t)\u2212r\u2032(0))2$]. In the present work, the mean-square displacement is tracked for $100a2/D$, long enough to identify particle motion but short enough to avoid significant influence from gel aging. The mean-square displacement is presented as an average over all particles; and to interrogate the dynamics imparted by the number of contact $Nc$ (bonded neighbors) of a particle, the mean-square displacement is also binned by $Nc$. “Binning” by contact number entails identifying the contact number of each particle, then the mean-square displacement is computed for this group as it evolves. The particle may change contact number over time; while the particle remains in its original bin, we can track the temporally changing contact number for each particle.

## III. RESULTS

Here, we report both the experimental and computational measurements of linear-response rheology of reversible colloidal gels formed from well-dispersed suspensions of hard Brownian spheres at different volume fractions $\varphi =0.20$, 0.30, and 0.44. Gelation was induced *in silico* and in experiment by a short-range attraction ($\Delta =0.1a$), of two different strengths, $V0=6kT$ and $V0=15kT$ [Fig. 1(b)]. Gels were prepared as described in Sec. II. The results that follow are separated into three sections where the influence of particle volume fraction, gel age, and interparticle attraction strength is varied independently. In each section, we analyze the influence of each parameter on linear viscoelasticity, comparing experimental rheological measurements with simulations, and then attempt to elucidate the micromechanical origins via the detailed measurements of structure and particle-scale dynamics obtained by the latter.

### A. Influence of volume fraction

We first interrogate the role of initial gel volume fraction in the linear response rheology of reversible colloidal gels. Increasing the volume fraction of a dispersion of purely repulsive hard spheres directly increases the magnitude of the rheological response due to an increased contribution of entropic and hydrodynamic effects [61]. Yet, in a reversible gel, this direct relationship between volume fraction and rheological stiffening may be complicated by the presence of a gel network which evolves in time. To investigate the contributions of the gel network to the linear rheological response, experimental and simulation results for $6kT$ gels at a single gel age (as defined in Sec. II) are reported in this section and results for $15kT$ bond strength gels can be found in S3 of the supplementary material [89]. For these combinations of volume fraction and bond strength, gelation produces a bicontinuous structure: A network of amorphously structured strands with an interior volume fraction $\varphi \u22480.62$ as suggested by simulations [12], in general agreement with experimental findings [6,25,62].

Changes in particle volume fraction influence gel morphology, as shown by the structural snapshots from simulation in Fig. 3. As volume fraction is increased, visual inspection of the snapshots reveals a pronounced change in morphology: a more homogeneous spatial distribution of particles with smaller solvent pores; this makes sense if one expects a limiting condition of an attractive glass at very high volume fraction where voids are smaller than a particle size. The lower volume fraction $\varphi =0.20$ gel exhibits strands of ropey, elongated condensed regions. With less free volume, the gels of higher volume fraction appear to have strands with not only closer spacing of strand intersections, but also more spheroidal pores, formed by more rapid arrest and slower age-coarsening in denser systems, discussed below. Traces of the pores are shown on the snapshots in Fig. 3 to guide the eye. The red and white circles guide the eye to strand and pore size and are explained alongside the discussion of Fig. 5.

We want to study the multibody contribution to the complex moduli obtained by SAOS [$G\u2032$ and $G\u2033$, Figs. 4(a) and 4(b)] and how the influence of three-body and higher order interactions evolves as concentration increases. The oscillation frequency is made dimensionless by scaling it with Brownian diffusion $D\u2217$ of a lone colloid in the suspending medium. The moduli are made dimensionless on the thermal energy per characteristic particle volume and plotted as a function of frequency and reveal the expected elastic behavior over all frequencies ($G\u2032>G\u2033$) and an increase in modulus with higher volume fraction (see the supplementary material [89] for Fig. S1 and discussion of adjustment of Brownian time scale in experiment based on solvent and depletant viscosity); but this measurement smears out the competing effects of hydrodynamic interactions, steric hindrance, and bond dynamics. To disentangle these effects, we subtract from $G\u2033$ the high-frequency viscosity $\eta \u221e\u2032$ since it does not depend on structure for a dispersion of purely repulsive, hydrodynamically interacting spheres, which is related to the short-time self-diffusivity [Fig. 4(b)]. We also normalize the moduli on the corresponding difference at equilibrium, $\eta (\omega =0;\varphi )\u2212\eta \u221e\u2032(\varphi )$, because the latter encodes purely the volume fraction dependence of equilibrium diffusion and gives a measure of how higher order interactions influence relaxation. The normalizing quantity, $\eta (\omega =0;\varphi )\u2212\eta \u221e\u2032(\varphi )$, increases with increasing $\varphi $ and is a measure of as three body and higher order *thermodynamic* interactions [63]. Here, these normalizing quantities are taken from a *dispersion* of purely repulsive hard spheres. (Values utilized to form the reduced moduli were taken from Stokesian dynamics simulations [64,65] for experimental curves and from Brownian dynamics simulations [55] for simulation curves. See S1 in the supplemental material [89] for these values.) Such a normalization can then scale out direct entropic and hydrodynamic effects, leaving behind the influence of volume fraction on gel structure, bond distribution, and glassy frustration, forming an “attraction-driven modulus” revealing the indirect influence of structure on linear-response rheology relaxation in reversible colloidal gels. In a dispersion of purely repulsive hard spheres, one expects a simpler shift or a total collapse of the curves, depending on the volume fraction. This is because there is only one length scale in a mobile dispersion, even when three-body and higher order interactions influence relaxation [66] (at concentrations near the glass transition, it has been shown that the different volume fraction dependence of short and long time diffusivity drives failure of time-concentration superposition [67]). However, in a gel, an infinitude of relaxation time scales are present (and, as we will show below, volume fraction shifts the contact number and length scale distributions). Thus, one expects that the viscous dissipation and elastic stresses arising from many-particle coupling due *strictly* to the combined effects of particle bonds and glassy frustration [12] should produce a qualitative dependence on volume fraction of the frequency-dependent moduli. That is, normalizing the moduli on the concentration-dependence of a hard sphere dispersion should, to a first approximation, scale out direct excluded volume effects, revealing the influence of particle bonds and glassy frustration, which can change with volume fraction. Undoubtedly, dynamics at other length scales present in a gel could provide more detailed corrections.

The normalized moduli are plotted in Figs. 4(a) and 4(b) and reveal a strong dependence on volume fraction: The curves for the different values of $\varphi $ remain well separated. This trend is followed in both simulation and experiment, with good agreement. Surprisingly, not only do the curves fail to collapse together, there is a reversal of the trend in Figs. S1(a) and S1(b) (where moduli grew with increasing volume fraction); here, when the equilibrium hard-sphere, frequency-dependent behavior is scaled out, this trend reverses. The importance of the choice of normalization is that it decouples to a first approximation the obvious boost in modulus arising from simple steric hindrance and more no-slip surfaces; when this contribution is removed, this brings the curves closer together, depending on the oscillation frequency, with a reversal in trend at high frequency. This is surprising because one might expect that increasing the volume fraction increases the surface area available for bonding (and hence the number of bonds), in addition to changing gel morphology. This suggests that the influence of volume fraction on direct steric hindrance and hydrodynamic coupling, and thence on linear response of the gel, can be partially scaled-out (or reversed)—but its influence on gel mechanics cannot. It further suggests that changes in the number of individual bonds (bonded surface area) play a less important role than gel morphology (unbonded surface area) or that the dynamics of glassy or surface particles and dominant network length scales set the relaxation dynamics at least at higher frequencies. This is explored quantitatively below.

To understand the surprising high-frequency trend of decreasing modulus with increasing volume fraction when equilibrium structure is scaled out, we recall prior work that showed that an increase in the network length scales with gel age produces a stiffer gel and that the largest length scale sets the response [12]. However, the length scale and volume fraction are interrelated, as discussed below, which indicates a more complex dependence on length scales. Here, a change in volume fraction also affects the distribution of length scales. As volume fraction decreases, the largest and most dominant network length scale increases. This may seem counter-intuitive unless we recognize that, because the direct influence of the volume fraction on the particle diffusion has been scaled out (at least to a first approximation using the corresponding liquid state contributions), the rest of the dependence should be a result mainly of the morphology and details of the bonded network. In contrast to dispersions, where steric hindrance and hydrodynamic entrainment set relaxation over the single-particle length scale, evidently the changes in gel morphology with changing volume fraction drive the dominant relaxation behavior over all frequencies. Unlike quiescently aging gels where both the dominant length scale and the bondedness (number of bonds) increase with gel age [12], increasing volume fraction suppresses the formation of longer length scales while allowing the gel to obtain more bonds. That is, the number of bonds plays a less predictive role than the structures which they form. Overall, these trends indicate that the influence of hydrodynamic interactions and entropic hindrance corresponding to the equivalent volume fraction [e.g., the $O(\varphi 2)$ contribution] can be to a first approximation removed with dispersion scaling and that the interplay between bond dynamics and glassy frustration that sets gel morphology ultimately drives its linear response. As discussed also below, both the number of bonds and the overall structural details (i.e., dominant length scales and local volume fractions) play a role in the determination of $G\u2032$. This agrees with the findings of Bouzid *et al.*, who showed via a simple constitutive model of the frequency-dependent linear response of a gel, that the propagation of stress by structural elements is capable of introducing an elastic response even without the presence of hydrodynamic interactions [68]. The fact that their study involved an athermal gel, i.e., very strong interparticle bonds, suggests that propagation of stress response by large-scale structure, rather than individual bond dynamics or long-range hydrodynamic interactions, is quite general. It should be noted that other studies suggest that the presence of hydrodynamic interactions shifts the gelation point [69,70] and possibly influences linear response. However, comparison of relaxation modulus in a gel to corresponding in a purely repulsive hard sphere dispersions at similar volume fraction reveals that the qualitative behavior arising from hydrodynamic interactions is identical in both systems meaning that it does not manifest differently in a gel and presumably can be scaled out as well [1,71–74]. Another view point emerges from Poon and co-workers’ recent comparison [75] of freely draining simulations to full hydrodynamics simulations similarly suggests that hydrodynamics can be scaled out, i.e., that hydrodynamics simply slow down structural evolution and, therefore, structural details should be compared at equivalent aging times. A more detailed study comparing purely repulsive, hydrodynamically interacting systems to corresponding volume fraction gels is needed to fully disentangle these questions.

A quantitative comparison of the scaled moduli from experiments and simulations is highlighted in Fig. 4(c), where the ratio of the experimental to simulation moduli, $Gexpt\u2032/Gsim\u2032$ and $Gexpt\u2033/Gsim\u2033$, is plotted over the range $0.15\u2264\omega a2/D\u2217\u226424$. The moduli are compared after the normalization presented in Figs. 4(a) and 4(b); for the comparison of the moduli prior to normalization, see Fig. S2 in the supplementary material [89]. As described above, the scaling introduced in Figs. 4(a) and 4(b) was expected to reduce or eliminate inherent differences due to the multibody and higher order terms arising simply from hydrodynamically interacting particle surfaces. This is supported by ratios close to unity for the curves shown in Fig. 4(c). However, some differences persist, such as slight differences in the frequency-dependent response. This suggests that any remaining difference in frequency-dependent behavior is introduced via the presence of hydrodynamic interactions in experiments or the difference in sample preparation history. Rapidly evolving evidence of the influence of shear protocols on gel structure and rheology suggests that the exact frequency-dependent response depends on shear history, which is observed in both similar gels to the present system [26] and model thixotropic system measured via superposition rheology [76]. It would be surprising to see closer agreement between two gels with such distinct preshear history. Improved agreement at lower volume fractions highlighted in Fig. 4 suggests that the more mobile surface particles (when Brownian motion is strong) enable the gel to better erase any remaining memory introduced by the preshear protocol, thus bringing the ratio of experimental to simulation moduli closer to unity. As further shown below, gel evolution is suppressed with higher volume fraction even in the absence of hydrodynamic interactions in our simulation study and thus we infer that a similar phenomenon erases some preshear history.

We now turn our attention to structural analysis. Above, we found that once the moduli are scaled to remove contributions due to equilibrium hydrodynamic and entropic effects, increasing volume fraction reduces the attraction-induced moduli. The remaining relaxation mechanisms span length scales from individual particle diffusion over all length scales up to the network length scale, where the qualitative volume-fraction dependence left behind suggested that changing network morphology with changing volume fraction was the underlying cause of the surprising decreasing dependence of the attraction-induced modulus on increasing volume fraction. The agreement between our freely draining simulations and fully hydrodynamically interacting experiments in the present study suggests that bond dynamics and steric glassy frustration are the primary drivers of this behavior. Here, we examine the distribution of length scales in the gel and how these evolve with volume fraction to explain this behavior. Two measures were utilized to provide quantitative analysis of these effects: The static structure factor and the distribution of void volumes (pore size) plotted in Figs. 5(a) and 5(b), respectively, and defined in Sec. II D 2. The static structure factor $S(q)$ is plotted as a function of wavenumber $q$ made dimensionless with particle size in Fig. 5(a). A peak emerges for each curve that corresponds to the most dominant length scale in the gel, occurring at $qa<1$ for all volume fractions owing to the dominance of long-range correlations in structure, i.e., large network features (multiply bonded particle strands), characteristic of bicontinuous gels created by arrested phase separation [6–8,11,49,77]. The peak widens and moves downward and to the right as $\varphi $ increases, toward a finer structure, consistent with the progression toward an attractive glass as $\varphi \u2192\varphi max$. The data indicate that higher initial volume fraction promotes formation of thinner and shorter strands with more numerous network junctions, a consequence of limited free volume (i.e., arrest dominated by steric hindrance). The probability of finding correlations in the low wavenumber region also decreases as volume fraction increases, showing that higher volume fractions inhibit development of a single pore or strand size. The curves overlap at high wavenumber because the particles are the same size for all volume fractions. The peak of $S(q)$ that corresponds to the nearest neighbors is found at higher values ($qa\u223c\pi $) not shown on this plot in order to focus on longer length scales.

A real-space interrogation of network-scale structure is made by measuring the void volume distribution P($Vvoid$), which reveals the physical size of interstitial pores by approximating voids as spheres [Fig. 5(b)]. Consistent with $S(q)$, the peak value of P($Vvoid$) moves to the left as volume fraction increases, showing that the length scale of pores (essentially, the “void network”) is smaller in denser gels. Both this trend with $\varphi $ and the corresponding trend in $S(q)$ are consistent with the idea that shrinking length scales that arise from higher initial volume fraction produce a weaker attractive gel modulus. This is consistent with our prior work [12] where we proposed that the most dominant length scale frozen into the gel by attractive forces set its linear response. That is, the indirect influence of higher initial volume fraction remains when we plot attraction-driven contribution to the moduli: Only that which arises from the gel structure, which does arrest to smaller, finer structure as $\varphi $ increases, is reported (the $\varphi $-reversal in the attraction-driven moduli). Direct influences are removed by scaling out the steric or hydrodynamic effects. Similarly to $S(q)$, the most probable void volume range widens with increasing volume fraction, suggesting that a single void “size” may not be as relevant; a more widely distributed range of length scales should appear in the high-frequency moduli as a weaker dependence on a single length scale, tested in Sec. IV B. The peak value of P($Vvoid$) also increases with increasing $\varphi $, indicating that finding small voids is more probable. One explanation for the growth in the probability of the average void size with increasing $\varphi $ is the assumption of spherical pores. The measurement of void volume by design considers the size sphere that fits into a given region of space, and a pronounced increase in void probability may indicate the presence of more spherical pores. Lower $\varphi $ concentrates the pore spaces to fewer, larger domains which may not be as spherically shaped. Increasing $\varphi $ appears to produce gel morphology that drives relaxation down to smaller and smaller length scales and create a uniform, smaller void size. This supports the claim that the moduli reveal the role of dense glassy particle strands in propagating the applied disturbance and may be able to signal a transition to dominance of attractive glasslike behavior (e.g., dependence on particle cages) rather than gel-like behavior (e.g., dependence on network length scale) as volume fraction approaches maximum packing. Overall, the shift to smaller length scales with increasing volume fraction plays a key role in the $\varphi $-reversal in attraction-driven moduli.

The dominant length scale $LS(q)$ and average pore diameter $d\u27e8VV\u27e9$ are plotted together in Fig. 5(c) to compare the frequency and real space computations; the details of measuring $LS(q)$ and $d\u27e8VV\u27e9$ are reported in Sec. II D 2. The network and pore measurements evidently provide consistent information regarding changes in relaxation length scales with increasing volume fraction. The decrease in the characteristic size of both “networks” of the gel—the network of thick particle strands and the network of pores or void spaces—agrees qualitatively. The prevalence of these dominant length scales decreases as $\varphi $ increases and this in turn influences the validity of existing models for gel linear response which rely on such measures, as discussed in Sec. IV. Neither of these length scales is totally unambiguous; $LS(q)$ comes from particle correlation distances but strands are not isotropic, and the void volume calculation forces an isotropic measurement, a sphere inserted into irregularly shaped pores. It is challenging to resolve such structure in finer detail in both simulation and experiments, but together they are consistent. The length scales are illustrated on the structure plots in Fig. 3, with a white solid circle for the void volume and a yellow dashed circle for $LS(q)$.

The evolution of gel structure discussed above is set by individual particle migration [12], driven by a combination of Brownian motion and the energy landscape where particles crawl along the network seeking curvature in which they can acquire more bonds and reduce their energy [35]. Here, we examine the change in particle dynamics with increasing initial volume fraction, to gain further insight into the corresponding change in linear response of the gel. In the simplest manifestation of the changes in structure with increasing volume fraction, colloids hinder one another sterically, forming temporary cages that restrict mobility. This can be quantified by the particle contact number, $Nc$, which ranges from zero for free diffusers to a maximum of 12. The bumpy surface of the gel network is formed by particles with one to seven contacts, $1\u2264Nc\u22647$; the glassy interior of network strands comprises particles with eight or more contacts, $8\u2264Nc\u226412$ [12]. The contact-number distribution $P(Nc)$ of a gel evolves with time as a gel age-coarsens [12] and provides a basic means to describe the structure. All particles in the simulation were “binned” by the contact number and this distribution is plotted for the three volume fractions in Fig. 6(a). The peak value of the curve always occurs at $Nc\u22656$, consistent with an isostatic elastic structure [78]. Moreover, the peak narrows and moves to higher $Nc$ with increasing $\varphi $, indicating that the surface area-to-volume ratio decreases: The $1\u2264Nc\u22647$ population (surface particles) shrinks while the $Nc=8\u221210$ population (interior particles) grows. Inspection of the tails shows that there are few free diffusers and little crystalline structure, indicating that strand interiors are disordered. The shift to the right of the distribution indicates not only do the dense regions grow and surface area shrinks, but the surfaces also become less mobile. Snapshots (Fig. 3) show that the lower volume fraction gels have “blobbier” strands and larger pore spaces while more crowded systems (higher $\varphi $) drive strands to fit around more spherical pores, not only reducing the available surface area but possibly rendering it more arrested via a flatter energy landscape. This supports our proposal that the attractive contribution to the moduli is less dependent on the bondedness (e.g., that the total number of bonds increases with $\varphi $) than on the formation of longer length scales.

To take a closer look, we examine particle motion as a means to describe the level of arrest of the gel. First, we ask what we might expect. In dispersions of purely repulsive hard spheres, particle self-mobility is a monotonically decreasing function of volume fraction (owing to both entropic and hydrodynamic effects) that in turn alters bulk viscoelastic response [1]; in our recent work, we showed that the slowing of particle mobility by attractive forces also plays a role in gel modulus, at least indirectly by changing the rate at which a gel age-coarsens and thus age-stiffens [12]. Carrying forward our observation that increased volume fraction indirectly influences the attractive contribution to the moduli, we here examine how particle dynamics change at fixed attraction strength and varied volume fraction. Particle mobility, or diffusion, can be interrogated as an average over all particles in the gel by monitoring mean-square displacement over time, as shown in Fig. 6(b) for a duration of 100 Brownian times, an interval long enough to permit interrogation of particle dynamics, but short enough to avoid large-scale rearrangement of structure. Particle migration is subdiffusive at short times and nearly diffusive at long times for all volume fractions and, unsurprisingly, these dynamics slow quantitatively (differing vertical shift)—but not qualitatively (same temporal growth)—as volume fraction grows.

However, not all particles in the gel are equally mobile; in fact, particles buried deep within strands are nearly immobile, while particles on the surface can be diffusive over short times until they build a cage, or over long times as they escape a cage to coarsen the gel [12]. This behavior can be extracted from particle dynamics by binning particles according to contact number and computing the mean-square displacement as an average within each bin [Fig. 6(c)]. While the fastest particles are the least numerous, the most numerous particles, those on the surface of the gel with 4–7 contacts, are also quite mobile. Focusing on the surface-particle population, we observe a monotonic suppression of binned mean-square displacement as $\varphi $ increases. The long-time behavior of even the most mobile particles—with one to four bonds—slows down as volume fraction grows. The high-contact number particles embedded deep inside strands show the nearly frozen dynamics observed in our previous work, with higher volume fraction deepening this freeze. On the log–log plot, the effect of volume fraction is minor at the highest contact numbers, because the interior of strands never exceeds 62% regardless of starting volume fraction. The dynamics of the surface particles are surprising: Surface-particle mobility is suppressed by higher system volume fraction, implying that the surface itself is more densely packed (we have previously described it as a liquid [12]), or equivalently the energetic landscape is flatter as volume fraction increases. This indirect effect of volume fraction on particle dynamics occurs via the flattening of the attraction-induced energy landscape, where increasing volume fraction promotes smaller network length scales and a decrease in the surface-area to volume ratio, a combination that changes the curvature of strands [35]. At volume fractions even higher than those studied here, it is easy to imagine that the abundance of particle contacts creates a frustrated but low energy state where single particle-size solvent pores are surrounded by concave surfaces. In this deeply arrested state, the energy landscape cannot get any flatter without external perturbation [35].

In summary, we have found that the magnitude of the viscoelastic moduli grows with volume fraction in both experiments and simulations. Removing known entropic and hydrodynamic contributions to the moduli from hard-sphere suspensions results in the attractive contribution to the moduli which exhibits a reversal with volume fraction: Less dense gels achieve greater “attraction-driven” moduli. Higher volume fraction, with attractive interactions, suppresses the prevalence of large-scale structure and produces a wider distribution over all length scales, until at the highest volume fraction the distribution becomes quite broad with a maximum size approaching a few particle sizes. The suppression of particle dynamics on the most mobile part of the gel—its surface—reveals the energetic impact of higher volume fraction: The energy landscape that drives particle migration is flattened, suppressing the formation of large structures and thus reducing the attractive contribution to the gel moduli. The direct, equilibrium hydrodynamic, and entropic effects can be scaled out of the experimental results, aligning them closely with simulation, providing convincing support that the most prominent contributor to gel response, beyond equilibrium effects, is attractive forces.

### B. Influence of gel age

We next consider the influence of age on the viscoelastic response. The aging behavior of colloidal gels and glasses reveals the nonergodic evolution of particle microstructure and commences during the arrest process [5,49] and continues indefinitely [12,79,80]. Aging continues throughout testing, requiring careful interpretation of results. In addition, shear-rejuvenation protocols that break a gel down to a flow of individual particles [25] nonetheless encode a distorted microstructure into the rejuvenated gel. Recent work has shown that stronger rejuvenation flow leads to enhanced stiffness [26]. In contrast, in simulation there is no need for preshear. Thus, one expects an initial difference in gel stiffness between experimental and simulation that depends on the preshear history. Age-stiffening of reversible colloidal gels was studied by Zia and co-workers in the linear-response [12] and nonlinear regimes [31,32] at moderate volume fraction, $\varphi =0.20$. Here, we expand the study to higher volume fraction and compare with experimental results. Data for $\varphi =0.30$ gels are presented here and for $\varphi =0.44$ gels in S2 of the supplementary material [89].

The rheological response for the $\varphi =0.30$ gel with bond strength $6kT$ is shown in Fig. 7 for several ages. The same normalization presented in Fig. 4 is applied here, to remove the equilibrium hard-sphere and hydrodynamic effects that make well-understood contributions, and instead focus on the aging effects of attractive forces. The oscillation frequency is made dimensionless as before by scaling with the Brownian frequency for a lone colloid through the suspending medium, $D\u2217$, while gel age has been made dimensionless on the pure Brownian diffusion time $D$, i.e., the diffusivity of colloids in pure solvent. In both experiments and simulations, the elastic modulus [Fig. 7(a)] age-stiffens over all frequencies shown. For the attraction strength studied in Fig. 7, such scaling of gel age on $D$ corresponds to a factor of two difference in age between simulations and experiments. This directly affects the mobility of low-contact number particles and also reduces the frequency of bond breaking. However, the strength of attractive forces and the energy landscape formed by strand curvature are the primary forces for age coarsening. While it is tempting to seek a direct scalar adjustment, one must also recognize that the preshear history of the experimental gel has already provided a measurable difference in the initial rheology of the gel, discussed further, below.

The viscous moduli for both experiment and simulation are plotted in Fig. 7(b), which both show age-thickening at high frequency. The elastic modulus is still dominant over all frequencies for both simulation and experiment. However, at lower frequency, the two disagree. Simulations indicate a decreasing low-frequency viscous modulus as the gel ages, which also decreases with decreasing oscillation frequency. In contrast, experiments show a low-frequency plateau in the viscous modulus, more consistent with an attractive glass. This difference may arise from the preshear history of the experimental system. Strong shearing motion causes particle rearrange and presumably increases average bondedness. Less mobile particles overall may reduce the ability of the gel to dissipate energy viscously via Brownian motion of the individual particles. Such age-thinning of $G\u2033$ seen in simulations may be in agreement with experimental studies of soft colloidal glasses [81] or hard-sphere attractive glasses and gels of stronger bond strength than shown here [15]. Further discussion on the appearance of a reduced viscous dissipation at lower frequencies observed in experiments is discussed further below and in Sec. IV.

Comparison of the experimental to simulation values of $G\u2032$ [Fig. 7(c), solid curves] shows a nearly constant factor of two different for all ages, beginning with the youngest gel. That is, the disparity does not increase with age, suggesting that the shear rejuvenation protocol encodes structural changes during gelation that made it a stiffer initial gel [26]. Once the side-by-side testing commences, aging proceeds at the same rate in both simulation and experiment. This behavior has been previously ascribed to the growth of the dominant network length scale as the gel coarsens; this is discussed in Sec. IV B. Comparison of the ratio of the experimental to simulation viscous moduli [Fig. 7(c), dashed curves] highlights these frequency dependent differences, with the ratio approaching unity as frequency increases. Differences in the qualitative and quantitative behavior of $G\u2033$ with gel age could arise from the difference in structural formation, arrest, and aging due to the initial gel state (preshear or instantaneous quenching). Since the equilibrium hydrodynamic contribution has already been scaled out, any remaining $O(\varphi 2)$ hydrodynamic effect is likely to be minimal compared to preshear history.

To quantify the structural evolution in Fig. 8, measurements of the static structure factor [Fig. 9(a)], the void volume distribution [Fig. 9(b)], and length scales extracted from these measurements [Fig. 9(c)] are presented. First, the static structure factor in Fig. 9(a) shows both an increase in the overall intensity of $S(q)$ at low wavenumbers and a shift in the peak to lower wavenumbers with gel age. The wavenumbers examined here, $O(qa)\u223c0.1$, are smaller than in most experimental studies in the low $q$ regime [6,8,13], where $S(q)$ increases to finite values reflecting the magnitude of the long-wavelength fluctuations (as $q\u21920$) which, at thermal equilibrium, are linked to the osmotic compressibility of the sample. In simulation, we are able to observe the decay of $S(q)$ toward small but finite values as $q\u21920$ which is expected in a bicontinuous system [11,82].

The void volume distribution [Fig. 9(b)] shows greater pore sizes with increasing age and a lower overall probability of all pore sizes with gel age. This suggests that a consequence of a growing larger pore size [or a shift to the right of P($Vvoid$)] also reduces the likelihood of finding the smaller pore spaces. Such a shift in the void sizes could be attributed to the coarsening and smoothing of surfaces with gel age which contribute to both larger pores and fewer small features.

Comparison of the strand evolution (static structure factor) and pore evolution (void volume) is presented in Fig. 9(c). The upward trend of both with age support the visual inspection of Fig. 8 that strands thicken and pores open up with age. One limitation of our measurement method for void volumes is the approximation of pores as spherical. It is likely that the weaker growth of the peak void size with age is affected by this simplification, because pores in Fig. 8 appear to be elongated or spheroidal. The results recover the qualitative trends presented by Zia *et al.* [12] for a gel with $\varphi =0.20$ initial volume fraction, such as a similar rate of growth of the dominant length scale over the ages shown here, but a far smaller length at each age. A second peak in both $S(q)$ and P$(Vvoid)$, not seen at $\varphi =0.20$, appears here at $\varphi =0.30$; this effect is even more pronounced at $\varphi =0.44$, as shown in S2 of the supplementary material [89]. This new length scale suggests that a second strand length or diameter is also favored as volume fraction increases; static structure factor does not distinguish which feature is measured. With increasing constraints on free volume, a second length scale may arise simply due to a more rapid initial arrest which does not become more uniform with ongoing coarsening due to the reduced rate found in Sec. III A.

Overall, higher initial volume fraction leads to slower aging and slower growth of the viscoelastic moduli. This can be understood by recognizing that aging is driven by curvature of the energy landscape [12,35]. The network strands vary in radius and thickeness, both around and along a strand. Particle migration during age coarsening is driven by particle seeking to acquire more bonds and reduce their energy; surface regions of tighter curvature offer more available bonds, while flat or concave-outward surface regions offer fewer opportunities for a migrating particle to acquire more bonds [12].

Gel aging is actually driven at the particle length scale, by diffusive migration of all mobile particles [12]. Particle mobility determines how easily particles can detach from or migrate along the network, with particles with fewer bonds (lower contact number $Nc$) than particles with many bonds (higher $Nc$). The distribution of contacts for all particles in the gel (surface, strands, pores) is presented in Fig. 10(a) for three sequential ages. The peak value of the distribution occurs at $Nc\u22656$, consistent with an isostatic microstructure [78], and similar to lower volume fraction gels with the same attraction strengths [12] (see also Sec. III A). As the gel ages, the peak moves toward larger values of $Nc$ as a result of structural coarsening and a minimization of the surface area-to-volume of the strands seen qualitatively in Fig. 8. The populations of surface particles ($Nc=1\u22127$) shrink with time while the populations of interior particles ($Nc=8\u221212$) grow. The contact-number distribution also broadens and lowers with age, recovering trends found in dense gels of lower volume fraction (also seen for the $\varphi =0.44$ $6kT$ gel, S2 of the supplementary material [89]). The growth of the high-contact number population reflects an increasing volume of interior particles, while the reduction in particles with fewer than seven contacts reflects reduction in surface area with age. Since the surface particles are most mobile, this evolution in contact number brings a slowing down of particle dynamics. While this is expected at lower volume fractions, here we examine whether the overall flatter energy landscape results in a less dramatic change in particle mobility with age. Comparison of Fig. 10(a) with Fig. 6(a) suggests that the initial volume fraction encodes the initial width of the distribution, a trend that is confirmed by examining aging behavior of other volume fractions (see supplementary material [89] for $\varphi =0.44$ and [12] for $\varphi =0.20$).

The evolution of the mean-square displacement (averaged over all particles) decreases monotonically with age as seen in Fig. 10(b) for $\varphi =0.30$. The shift in contact-number distribution [cf. Fig. 10(a)] to more bonded particles suggests that this decrease arises from growth of the volume of the dense, glassy interior of the strands and concomitant decrease of mobile surface particles. To test this idea, we compute the mean-square displacement as an average within each contact-number population.

The mean-square displacement binned by contact number is plotted for 100 Brownian times in Fig. 10(c) for the same three ages. Even if the particles acquired or lost contacts during this time, they remained in the same bin, following Zia *et al.* [12]. As expected, diffusivity decreases with age for all contact-number populations owing to the net formation of bonds. We recall that in Fig. 6(c) we saw a similar trend, owing to the increase in average contact number for denser gels; here, this slowing down with age can also be recognized as a progression toward glassy dynamics. In particular, particles with $Nc\u22658$ show subdiffusive displacements over all time scales probed here, with the maximum particle displacement smaller than the attraction range. This indicates that such highly bonded particles simply wiggle inside their attractive well while large-scale displacements are essentially frozen. This is very similar with findings for attractive glasses [83]. Furthermore, the binned particle dynamics suggest that aging proceeds in a similar, if quantitatively suppressed (that is, the magnitude is smaller), manner to lower volume fraction gels [12] suggesting the morphology encoded by initial volume fraction sets the “kinetics” of age coarsening but does not alter its qualitative features. The reduction in the viscous response of gels of advancing age in simulation at low frequencies may be related to the quantitative decrease in particle mean-square displacement since long times are related to low frequencies. The slow-down of dynamics with increasing gel age, enhanced by higher initial gel volume fraction, may be more easily detected in simulation, due to the absence of effects of preshear, depletant molecules, or hydrodynamic interactions.

In summary, colloidal gels of higher volume fraction with reversible bonds age in a qualitatively similar manner to lower volume fractions. Comparison of experimental and simulation viscoelastic moduli shows agreement for the age-stiffening of the elastic structure, but differences appear in the viscous response at low frequencies. On the log–log plot, these discrepancies are exaggerated when in reality are on the order of experimental error. Exploration of differences would focus on sources such as the preshear rejuvenation history (effects discussed in Sec. II C 2) or morphological differences produced by hydrodynamic interactions since direct effects have been scaled out. In either case, the difference is small. Microstructural measurements indicate that a broader range of length scales is generated in denser gels, indicating that more length scales may be necessary to predict viscoelastic behavior for denser gels or the increasing influence of the glassy structure on the rheological response when initial volume fraction is increased. Dynamics indicate a similar qualitative coarsening behavior to that of $\varphi =0.20$ gels, with the distinction that dynamics are quantitatively suppressed with increasing age and are initially lower than $\varphi =0.20$ gels due to a gel morphology with fewer, smoother surfaces with increasing initial volume fraction gels which slows age coarsening and thus prevents the formation of longer length scales. The connection between rheology and the structural length scales will be discussed in Sec. IV.

### C. Influence of bond strength

Lastly, we examine the influence of bond strength. The primary contribution of bond strength to the evolution of a reversible colloidal gel is the duration of a bond before Brownian motion breaks it; thus, although stronger bonds may drive more rapid initial phase separation, it subsequently slows quiescent aging [12]. Bond strength exerts both direct and indirect effects on linear response, changing the stress required to rupture bonds and altering gel morphology, respectively. These roles have been examined via experimental and theoretical approaches, reviewed in Sec. IV.

The linear response rheology for $\varphi =0.30$ gels with bond strengths of $V0=6kT$ and $15kT$ via experiment and simulation is presented in Figs. 11(a) and 11(b) at a gel age of $4000$ $a2/D$. Stronger bonds produce a higher elastic modulus, which may arise directly from bond strength or from the difference in structure and length scales that emerge quickly after gelation. The interplay between bond strength and structural length scale and how these effects set the strength of gel elasticity will be explored below via measurements of structure and dynamics and in Sec. IV via the use of existing models connecting bonds or gel structure to rheology. Stronger bonds also increase the viscous modulus [Fig. 11(b)], even in the low-frequency plateau region. Recognizing that attractions hinder diffusion and thus reduce the ability of the gel to dissipate energy viscously, we expect this to arise primarily from longer length-scale structures, discussed in Sec. IV.

The moduli found in experiment are systematically higher than that found in simulation [Fig. 11(c)], similar to all other cases discussed above. Agreement between the elastic moduli only improves with stronger attraction, which may simply arise from the preshear history, which pre-ages the experimental gel—because aging is slower in stronger gels. The most pronounced discrepancy between experiment and simulation is once again the viscous moduli at low frequency oscillation. The agreement at higher frequency (where long length scales seem to dominate the response) coupled with the quantitative discrepancy at low frequency (where a distribution of length scales affect response) suggests that remaining differences arising from sample preparation history or hydrodynamic interactions are quantitative, as discussed in Sec. II. Such differences are most likely to play a quantitative role in the response most sensitive to detailed particle structure, e.g., during weak oscillatory shear.

Bond strength encodes morphology into the gel soon after gelation, as illustrated in structural snapshots in Fig. 12. Particles are colored from red for few contacts to blue for many. Visual comparison of Figs. 12(a) and 12(b) shows that weaker bonds lead to faster evolution and coarser structure [12]. The $6kT$ gel exhibits thicker strands while the $15kT$ has finer strands that are only a few particles thick although they have reached the same age. This agrees with experimental findings for intermediate volume fraction gels where, as the attraction strength is increased, heterogeneities increase approaching the gel point, but then decrease as bond strength is increased within the gel region [13]. That is, the maximum length scale of the heterogeneous structure should appear around a bond strength that corresponds to the gel point. To comment further on the morphological changes, we next quantify the structure.

The static structure factor for $\varphi =0.30$ gels of varying attraction strength is shown in Fig. 13(a). The lower $6kT$ gel exhibits a stronger peak which also occurs at smaller wavenumbers; in contrast, the $15kT$ gel exhibits a broader peak that extends to higher wavenumbers and is more difficult to distinguish. The weaker peak for the $15kT$ gel suggests that the stringier structure observed in the snapshots is not as well characterized by a network- or strand-length scale. In contrast, the void volume distribution [Fig. 13(b)] shows a well-defined peak for both gels. The weaker $6kT$ gel exhibits a flatter distribution with two or more common pore volumes given by a peak at higher void volumes and a shoulder at lower volumes. A comparison of $\varphi =0.44$ $6kT$ and $15kT$ gels give a similar picture; plots can be found in S2 and S3 of the supplementary material [89]. Considering which network, either that of particles or that of the void spaces, sets the response of the gel, it is likely that the particle structure is responsible for the gel response; that is, gel elasticity depends on particle network length scale, an idea tested in Sec. IV. However, it is important to note that while a single long length scale may no longer be favored as attraction strength increases, the narrowing void distribution may be an indicator of an evolution toward a more homogeneous overall structure noted in prior experimental work [13].

The dominant length scale $LS(q)$ and the average pore diameter $d\u27e8VV\u27e9$ are shown in Fig. 13(c). The length scale of both the particle structure and the voids decreases as bond strength increases. The uniform pore structure and small characteristic strand length scale of the $15kT$ gel is distinct from the thick, bicontinuous structure of the $6kT$ gel. The stronger $15kT$ bonds, which promote a more “stringy” structure, dramatically slow surface-particle migration that is responsible for coarsening that in turns produces thick, glassy strands in reversible gels with weaker bonds. We examine the particle dynamics next to interrogate this further.

To examine the relative prevalence of each contact number, a strong predictor of single-particle dynamics in colloidal gels, we show the contact-number distribution in Fig. 14(a). The stronger $15kT$ gel shows a narrow distribution of $P(Nc)$, peaked near a lower value of $Nc$, indicating more surface area and less enclosed volume. This shift in $P(Nc)$ with increasingly strong bonds signals that fewer particles reside within dense strands. A bond strength of $15kT$ “freezes” the gel structure and permits little bond breakage and re-formation, dramatically slowing down the progress of ongoing phase separation. With few $Nc\u22658$ particles, the $15kT$ gel appears to have little glassy structure within strands as Brownian motion is strongly suppressed and particles are not able to form energetically favorable dense strands. That is, this stronger gel experiences arrested phase separation more rapidly, curtailing the formation of dense, thick strands.

Particle dynamics are another way to elucidate the differences between stronger and weaker gels. In Fig. 14(b), the mean-square displacement for $\varphi =0.30$ gels of varying attraction strength is shown. Both a qualitative and a quantitive change in mean-square displacement are observed: The weaker gel exhibits nearly diffusive motion while the stronger gel exhibits subdiffusive motion over the range of time shown. The average particle in the gel of $15kT$ bond strength moves about one-tenth of a particle size in 100 $a2/D$ (equal to the attraction range of $\Delta =0.1$). The weaker $6kT$ gel does exhibit a temporal or transient localization (intermediate nonergodic plateau) in the mean-square displacement, but this is consistent with the formation of cages of particles followed by cage migration at longer times [12]. This may help set practical limits on the idea of a reversible bond; while a $15kT$ bond has finite probability of breaking, single particle dynamics clearly show that this is not very probable: The average mean-square displacement shows a localization length (the maximum mean-square displacement at long times) which is of the order or smaller than the range of attraction [37]. Thus, bond strength suppresses average particle motion dramatically, hindering to a large extent particle motion responsible for surface migration and transient caging characteristic of a reversible gel.

Moreover, the contact-number binned mean-square displacement [Fig. 14(c)] further highlights the differences in particle motion between the weaker gel and the stronger gel. The lower contact-number populations ($0\u2264Nc\u22643$) of the weaker $6kT$ gel are nearly diffusive over the range of time examined, as expected for $6kT$ bonds [12], and in particular, the dynamics of the weakly bonded surface particles merge as seen in the lower volume fraction $\varphi =0.2$ gels [12] and here in more concentrated $6kT$ gels, shown in Sec. III A. The stronger $15kT$ gel shows that, first, there are no free diffusers and, second, only $2\u2264Nc\u22644$ particles may diffuse further than the length of the bond, $\u27e8r\u2032r\u2032\u27e9\u223c0.01a2$, over 100 $a2/D$. However, these weakly bonded particles do not form a single dynamical group nor do more bonded particles exhibit motion beyond their cages, suggesting that cage formation happens at very early times and, because dynamics are so frozen, the “evaporation/condensation” and “freeze/melt” equilibrium processes reported for reversible gels [12] are total suppressed soon after gelation. Such slowdown in bond breakage can also be examined utilizing theory, for instance, by calculating the Kramers escape time [30,62], where the time for a particle to escape its neighbor is nearly exponentially related to the well depth if the particles in contact are near the well minimum. Utilizing a ramp potential approximation with the same depth and range of attraction in the present work, a $6kT$ strength bond results in an escape time (time it takes to travel the length of the bond given its strength and range) of $0.44a2/D$ while a bond of strength $15kT$ is instead $580a2/D$. Hence, the stronger gel is predicted to suppress displacements of the order of the range of the attraction for times longer than shown in Figs. 6(b) and 6(c). This is in agreement with our simulation results that show that on average, particles cannot move outside of the bond within $100a2/D$, and when binned by contact number, only the most mobile particles ($2\u2264Nc\u22644$) may move such distances.

In summary, interparticle bond strength appears to play a more decisive role in distinguishing bicontinuous gels with reversible yet durable bonds from gels with nearly permanent character. Examination of the attractive contribution to the linear moduli again distinguishes the roles of network structure and bond strength in setting the rheological response. Although network length scales are reduced in size with increasing bond strength, bond strength appears to play a quantitative role in enhancing the strength of elasticity and dissipation of the gel network. Stronger bonds suppress gel aging, preventing the formation of significant regions of glassy structure, and promotes the formation of a more “homogeneous” gel structure with a narrow range of pore structures and particle contact numbers. Cage formation and migration are suppressed soon after gelation for the permanently bonded gel, concomitant with their slow coarsening. This lends strong support to the “Smoluchowski’s Ratchet” model of Zia and co-workers [12]: Without this mechanism, the gel cannot coarsen.

## IV. COMPARISON TO MODELS

A number of models have been put forth in the recent literature that connect the linear-response rheology of reversibly bonded gels to cluster dynamics [84], bond strength and volume fraction [37], network length-scale structure [12], and particle dynamics via the Stokes–Einstein relation [85]. Each of these models contributes to understanding of the origin of mechanical strength and elasticity in reversible gels. We examine the relevant connections of volume fraction, bond strength, and structure to rheology found in the present study through the lens of several of these models.

### A. Link between particle dynamics and linear rheology: A generalized Stokes–Einstein relation approach

An opportunity presented by this body of work, combining experiment and simulation, is to interrogate the ability of models which leverage single particle dynamics to predict the linear response of a colloidal gel. The relaxation modes associated with hindered transport in a colloidal gel are encoded into the frequency dependence of the viscoelastic moduli. A natural place to start, and an approach taken in previous studies [85,86], is to study viscoelastic materials utilizing passive microrheology to connect the average mean-square displacement of all particles to the linear viscoelastic properties of the sample, assuming a generalized Stokes–Einstein relation (GSER) holds. While the Stokes–Einstein relation is constructed from two primary assumptions not valid for gels (equilibrium and motion in a continuum), it has nevertheless provided at least qualitative connections between diffusion and the inverse of the complex modulus.

We thus expect the self-motion of the colloids to be influenced directly by crowding and attractions that slow Brownian diffusion, as well as indirectly by network morphology, as discussed in Sec. III. Having computed the elastic and viscous moduli utilizing bulk oscillatory shear in Sec. III, here we compare those data to $G\u2032$ and $G\u2033$ computed utilizing the GSER as [85, 86]

where $\beta (\omega a2/D)$ is the instantaneous rate of mean-square displacement, the first correction to the Stokes–Einstein mobility, and is estimated to be

The quantity in Eq. (11) is computed utilizing particle mean-square displacement tracked in simulation over time and averaged over all particles and inserted into Eq. (10). The real and imaginary parts of the complex modulus give the familiar elastic and viscous moduli, where their computation in the context of single-particle motion depends on the frequency probed and the time evolution of the mean-square displacement as $G\u2032=G\u2217cos(\pi \beta (\omega a2/D)/2)$ and $G\u2033=G\u2217sin(\pi \beta (\omega a2/D)/2)$. In Fig. 15, the moduli thus computed are plotted alongside those measured via direct rheological interrogation (cf. Fig. 11) for $\varphi =0.30$ gels with bond strengths of $V0=6kT$ and $15kT$ and age $4000$$a2/D$. The moduli computed from the GSE relation compared to linear oscillatory shear for gels of other initial volume fractions and gel ages are provided in S5 of the supplementary material [87, 89].

The best agreement between the GSER prediction and bulk oscillatory shear measurements is obtained for the $15kT$ gel, where bonds are strong and particles largely immobile. Both moduli show strong qualitative agreement, and $G\u2032$ shows also strong quantitative agreement. This suggests that the average particle mobility is predictive of the dissipative and elastic character of the gel. In contrast, the more weakly bonded 6kT gel shows substantial quantitative disagreement for both $G\u2032$ and $G\u2033$ and strong qualitative disagreement—suggesting that the average particle mobility does not reflect the appropriate length or time scale for energy storage and dissipation. The discrepancy between the GSER curves and bulk rheology may arise from the large number of more readily mobile particles, which exhibit markedly higher mobility than the average. Gels like those studied here, in which bonds are reversible over observable time scales, can have substantial surface area covered by highly mobile particles that diffuse to coarsen the gel [12], at least when volume fraction is not too high. Such particle migration contributes directly the average mean-square displacement, but prior work has established a connection of larger length scales to relaxation time scales [12]. That is, their direct contribution to elastic response may not be captured by the GSER, and instead, dynamical heterogeneity and network length scales play a more prominent role.

Dynamical heterogeneity of particles, set by a distribution of number of bonds as well as steric hindrance, produces a spectrum of relaxation time scales. We thus expect the dynamical populations (free diffusers, mobile surface particles, deeply buried particles) to exhibit qualitatively different relaxation behavior, which is smeared out by averaging mean-square displacement over all particles in the gel. While bulk oscillatory shear also averages over all particles, the most prominent contributors to relaxation change with frequency and may correspond to different dynamical populations—which can only be extracted from the mean-square displacement by evaluating qualitatively such dynamical populations separately. To test the idea that these different dynamical populations not only reside in distinct regions of the gel structure, but also form subsets of gel mechanical behavior, e.g., the elastic or viscous response, we plot the linear moduli from the GSER deduced from the binned mean-square displacement in Fig. 16.

High-contact number particles that form the solidlike strand interiors, $Nc\u22658$, exhibit a strong elastic response compared to $G\u2032$ from bulk rheology. Particles with a few contacts, $0\u2264Nc\u22644$, which carry out a fast gas-liquid exchange at strand surfaces [12], give a markedly weaker elastic relaxation. This can be understood by recognizing that they experience much less steric hindrance and have fewer bonds, in turn allowing more vigorous Brownian diffusion set their relaxation time. In between these two groups are particles with five to seven bonds—the most populous dynamical group in the gel. These particles reside on the network surface and migrate along it [12]. Unfortunately, there is no obvious single or combined group with mean-square displacement that gives a strong prediction for either the viscous or elastic modulus. This should not be surprising, owing to the nonergodic, nonequilibrium, noncontinuum nature of the gel [88]. Perhaps, the most useful conclusion that can be drawn is that glassy dynamics produce the elastic response, and the viscous response cannot be well-described by single particle motion.

In summary, the prediction of linear moduli of colloidal gels from the average particle motion shows some success in strong, nearly nonreversible gels, but in reversibly bonded gels, both average and separate subgroup particle dynamics are a poor predictor of viscoelastic response.

### B. Scaling arguments based on structural aging

In our prior work, the evolution of the bicontinuous network structure of a reversible gel (of initial $20%$ volume fraction) could be described by the growth of a single dominant length scale obtained from the static structure factor $LS(q)$ [12]. When the formation of a clear, predominant long length scale is suppressed due to the slowed age-coarsening seen in denser and stronger bond strength gels, it is natural to expect that this relationship weakens and a spectrum of length scales may matter. Here, we test the scaling theory presented in [12] with denser gels by plotting the elastic modulus from oscillatory shear scaled on $LS(q)$ for $6kT$ gels of different ages in Fig. 17.

The high-frequency elastic moduli scaled on $LS(q)$ is present in Fig. 17. This scaling collapses the curves together quite well over all frequencies, indicating that age-stiffening results from growth in the dominant length scale. Overall, the dominant length scale appears to have a strong influence on the high-frequency elastic response of reversible colloidal gels, regardless of volume fraction. This is consistent with the idea that the magnitude of the high-frequency elastic response, where $G\u2032\u223c\alpha 1/2$ due to the presence of a diffusive boundary layer, is controlled by the relaxation behavior of the dominant length scale [12]. This suggests that the dominant length scale sets the elastic response, even if the dominance is not pronounced. We also scaled the linear moduli of the stronger, $15kT$, gel at different ages with $LS(q)$ which can be found in S4 of the supplementary material [89].

In summary, the high-frequency elastic modulus of reversibly bonded gels depends on the dominant length scale and predicts age-stiffening behavior for all volume fractions studied. In agreement with our analysis of the linear elastic moduli above, the gel network is a key component of gel viscoelasticity that is distinguishable from hydroentropic contributions of increasingly more dense gels. That is, within each volume fraction, coarser gel morphology contributes to gel elasticity in a manner that can be removed by a single network length scale. Volume fraction does indeed introduce differences in morphology which may not be removed by a single length scale, as found for gels of differing attraction strength by Zia and co-workers [12].

### C. Comparison with other theoretical models: MCT and cluster models

Several other classes of models have also considered the volume fraction or attraction strength dependence on the structure or dynamics of colloidal gels as a means to predict the linear moduli. For instance, MCT provides a means to calculate the shear modulus $G$ for dense colloidal suspensions in combination with PRISM theory describing the contribution of the polymer depletant [37]. The theory yields a relation between the zero-frequency elastic modulus, the particle localization length ($rloc$), and the volume fraction $\varphi $ [37],

Here, for each volume fraction and attraction strength (polymer concentration), localization length was taken from the data of Chen and Schweizer [37] and then $G\u2032$ was calculated through Eq. (12).

Figures 18(a) and 18(b) show a comparison of the elastic moduli from experiments and simulations at $\omega a2/D=0.1$ and $0.125$ (the lowest available) respectively, with theoretical predictions from MCT-PRISM theory. $G\u2032$ from simulation is smaller than that from experiment for all volume fractions and attraction strengths, as shown in supplementary material [89], when not considering equilibrium hydroentropic effects. For the weaker $6kT$ gels [Fig. 18(a)], MCT-PRISM-2 theory shows a stronger increase of elastic modulus with particle volume fraction compared to experiments and simulations. A source of the difference between this theory and the present results could arise from age or volume fraction-dependence of $rloc$. From the mean-square displacement measured in simulations, particle caging that reduces particle migration increases in higher volume fraction gels and in older gels, and the two effects are convolved (e.g., the progression of gel aging depends on $\varphi $ and $V0/kT$, thus influencing cage building and hopping). While the prior works consider reduction of the localization length with increasing $\varphi $, gel age and the rate of arrest or aging have not been considered. A smaller $rloc$ would impart a stiffer modulus, following Eq. (12) where $rloc<1$.

For the stronger $15kT$ gels [Fig. 18(b)], the model shows a similar increase of the elastic modulus with particle volume fraction compared to experiments, but it overpredicts $G\u2032$ by a factor of 30 for all volume fractions, indicating that $rloc$ far over predicts the level of arrest of particle motion. The theoretical elastic modulus can be reduced through an empirical approximation by the average number of particles participating in a cluster, $Ncluster$, as $Gexp\u2032=GMCT\u2032/Ncluster$. With $Ncluster=30$, the MCT-PRISM-2 predictions are in good agreement with experiments although Fig. 12(b), representing the corresponding simulations, does not appear to have such a number of particles, nor a length scale corresponding to $\u223c30a$. Further studies, involving the effect of aging on $rloc$ and the role of heterogeneites (or cluster) size, would necessary to provide a more detailed assessment of the applicability of this model in such colloidal gels.

Next, we compare our results to the model by Zaccone *et al.* [84] which takes into account the heterogeneous nature of the gel microstructure, based on a hierarchical arrest scheme differentiating the microscopic (primary particle-level) from the mesoscopic (cluster-level) contributions to the macroscopic elasticity. The model considers the clusters as compact spherical or quasispherical renormalized particles of diameter $L$, with the effective volume fraction $\varphi c$. The elasticity originating from the intercluster contribution $Gc$ is calculated as

where $z(\varphi c)$ is the mean coordination number of the clusters and $k$ is the bond stiffness which is measured as the second derivative of the potential at the minimum [84]. Here, we considered the dominant length scale $LS(q)$ from simulations as the average size of clusters $L$ in Eq. (13). Alternatively, Zaccone *et al.* propose elasticity can also arise due to the intracluster contribution $Gg$, also computed from Eq. (13). For $Gg$, $\varphi c$ and $z(\varphi c)$ are reduced to the particle level and $L$ becomes the diameter of particle, representing the dense strand interiors. For both $Gc$ and $Gg$, bond stiffness was computed using the Morse potential and $z(\varphi c)$ was taken from Zaccone *et al.* [84].

Figure 18(c) shows a comparison of $G\u2032$ from experiment and simulation, at $\omega a2/D=0.1$ and $0.125$, respectively, to the theoretical predictions of $Gc$ and $Gg$ [84]. The predicted intercluster elasticity $Gc$ shows a similar increase of elastic modulus with $\varphi $ as in experiments and simulations, suggesting that the increase of the elastic modulus with $\varphi $ is due to an increase of intercluster bonds and decrease of $LS(q)$, that is, a higher density of strands are connected by more particle bonds. However, this theory also overestimates $G\u2032$ for all $\varphi $, except for $\varphi =0.3$ experimental measurements. This might be due to the assumption that the bond stiffness is calculated at minimum of the potential whereas not all particles might be located at the minimum attraction well and hence they may have smaller values of bond stiffness. As in [84], $Gg$ largely over predicts $G\u2032$ and does not vary with $\varphi $ due to the assumption that the volume fraction inside the strands or clusters is constant ($\varphi c=0.62$, computed for a $\varphi =0.20$ gel in [12]) and independent of overall $\varphi $. That is, this model suggests that the elasticity of the gels is driven by the larger network structure rather than the predicted contribution from the glassy strand interiors over all volume fractions studied, emphasizing the enduring role of network structure in modeling approaches for even dense gels. Adaptation of this model to predict aging behavior remains to be done.

In summary, the comparison to several existing theories points to different physical interpretations of the primary driver of elasticity of the gels. Models based on localization lengths as a predictor of gel elasticity qualitatively fit simulation and experimental trends with volume fraction. Further study of the effects of age and heterogeneity on the localization length may prove useful. A comparison to a theoretical approach based on cluster dynamics instead suggests that the elasticity arises from the clusters themselves rather than the glassy interior of a cluster. Such an insight suggests that future avenues for theoretical treatments must still consider the network structure, and its temporal evolution, even at higher volume fractions and attraction strengths.

## V. CONCLUSIONS

Comparison of experimental study of moderately concentrated, reversible colloidal gels to *in silico* study of a freely draining model system yields two surprising results: first, the high-frequency linear-elastic response modulus is approximately linear in the dominant network length scale, although local volume fraction plays an equivalently important role in the determination of $G\u2032$. This seems to hold at any gel age, for volume fractions ranging from 20% to 44% and bond strength less than $\u223c10kT$. This is consistent with the predictions of prior models by Zia *et al.* [12]. Second, the full linear viscoelastic response of the experimental system is well represented in a simpler computational model in which wall effects, hydrodynamics, explicit depletant molecules, and rejuvenation were neglected, suggesting that the connections between morphology, dynamics, and rheology are encoded primarily by volume fraction, attraction strength, Brownian motion, and age. Overall, the comparison between presheared experimental gels and instantaneously quenched simulation gels is valid, but quantitative differences remain. Comparison of the effect of hydrodynamics on linear response in purely repulsive systems to that in gels suggests that these effects can be scaled out (notwithstanding the need for future work to interrogate these effects in detail). These conclusions are valid for colloidal gels formed from a special balance of crowding and attractive colloidal forces that together interrupt phase separation to produce a network of thick, ropey particle strands in which Brownian motion is strong enough to overcome attractive forces and steric hindrance to restructure the gel over time. Such gels exhibit a wide spectrum of length scales, from single-particle to strands to the entire network, where one might expect bond strength to set elastic response and single-particle diffusion that actively restructures the gel to set viscous response. Since dynamics and structure depend at least quantitatively on the strength of hydrodynamic interactions, we decoupled these contributions, at least to a first approximation by taking into account the dynamics in the corresponding liquid, to elucidate the attraction-induced linear response.

Removal of the well-known equilibrium hydroentropic contributions in purely repulsive hard-sphere suspensions from the gel’s viscoelastic response revealed the contributions arising directly from attractive interparticle interactions and indirectly from the network structure created by the attractive interactions, placing the experimental and simulation viscoelasticity on the same footing. Once the *volume-fraction dependent equilibrium contribution* to the elastic modulus has been scaled out, to a first approximation, the remaining response reflects the “attractive viscoelasticity,” arising from pair-level bonds and the gel network, which at high frequencies reverses the hard-sphere monotonic increase of elastic strength with steric hindrance. Attractive forces, via larger length scales, thus make a substantially higher contribution to the overall linear response of less-dense gels compared to higher density ones, although of course local volume fraction plays an important role as well. It is the most prevalent long length scale itself rather than simply the increasing relative amount of no-slip surfaces or steric hindrance that dictates the increase of the attractive contribution to the linear moduli. Overall, network-length structural relaxation is responsible for the elastic moduli of a broad range of reversible colloidal gels to an extent that was not previously recognized: Gel morphology is described by one length scale, even while it includes the influence of local volume fraction.

Particle dynamics, coarsening rate, and length-scale distribution can vary widely between gels of differing initial volume fraction, but in all cases, reversibly bonded gel morphology is characterized by a bicontinuous structure with a dominant network length scale. Though one single length scale seems to matter most, the network structure fills space in a qualitatively different manner with differing volume fraction. The ropey, blobby, interconnected network of the $\varphi =0.20$ gels forms outwardly curving surfaces, providing a rough energy landscape with numerous open energy wells for particles to sample as they diffuse along the network surface, ideal for the progression of age coarsening. Once the volume fraction is increased to $\varphi =0.44$, a different structural picture emerges: Pore spaces become round as the strands themselves form smooth, but concave surfaces. These concave surfaces have consequences for surface migration: The surface itself is no longer as “rough,” weakening the gradients in surface mobility, and slowing the growth of the dominant length scale. Thus, ongoing age-coarsening (growth of the dominant length scale) is further suppressed in dense gels, both directly via increased steric hindrance and indirectly by smoother, reduced surface area.

The most strongly bonded gels studied here, $15kT$, are more difficult to characterize using the same metrics applied to the reversible gels. First, even though a $15kT$ gel may have the same viscoelastic response as an aged, reversible gel, their morphologies may be dramatically different. Second, in contrast to the clear dominance of a single length scale in reversible gels, used to scale out age coarsening, no such dominant length scale appears in very strong gels, preventing such scaling. Thus, models such as the Rouse-like theory from our prior work do not give a good connection of structure to elastic modulus, while application of the Generalized Stokes–Einstein relation utilizing the average particle mobility provides a good prediction of gel elasticity. This suggests that elastic response in more deeply arrested gels emerges from single-particle dynamics.

Comparison to models that connect rheology to passive microrheology [85], network-length scale [12], single-particle dynamics [37], and cluster dynamics [84] show that, for reversible gels, the dominant or cluster length scale is the primary predictor of gel elasticity. Models based on averaged single-particle dynamics neglect the dynamical heterogeneity and an increasing depth of arrest with time of reversibly bonded colloidal gels. Yet, gel elasticity from bulk rheology is qualitatively recovered by a generalized Stokes–Einstein relation to convert glassy particle mobility into the elastic modulus [85], confirming the role of the glassy network strands as a driver of gel elasticity. While incorporation of the age-evolving network length scale is a key to predict gel elasticity, gel dissipation remains poorly modeled due to its complex dependence on dynamical heterogeneity in a reversible gel.

Hydrodynamic coupling between purely repulsive hard-sphere colloidal particles has long been recognized as a key factor in suspensions where no other force is present to transmit the motion of one particle to another over long distances. In a freely draining *suspension* of freely diffusing colloids, particles do not interact until they contact one another directly. In a freely draining colloidal *gel*, the network propagates information similar to many-body hydrodynamic interactions, coupling motion over long distances. This network structure is a key contributor to the ability of the gel to store and dissipate energy, consistent with prior work for strongly bonded colloidal gels [20,21,68]. Inclusion of hydrodynamic interactions in the simulation model may improve the modest quantitative disagreement between simulation and experiment here, but is not necessary to recover colloidal gel viscoelasticity measured in experiments. Future study of the nonlinear rheology of reversible colloidal gels will enable interrogation of the influence of hydrodynamic interactions on such response and determine if differences are observed beyond those known from the hard-sphere dispersion literature.

## ACKNOWLEDGMENTS

The authors thank E. Del Gado (Georgetown University) for many helpful and insightful conversations and A. B. Schofield (University of Edinburgh) for providing the PMMA particles. The authors also acknowledge helpful conversations with Benjamin J. Landrum. G.P. and E.M. acknowledge financial support by EU projects, “EUSMI” and “DiStruc.” G.P. acknowledges support from Humboldt Foundation. This work was supported, in part, for R.N.Z. by an Office of Naval Research Young Investigator Award (No. N00014-14-1-0744), a National Science Foundation CAREER Award (No. CBET-1351184), and a National Science Foundation Extreme Science and Engineering Discovery Environment (XSEDE) Research Award (No. OCI-1053575) to utilize high-performance computational resources at the Texas Advanced Computing Center (TACC) at the University of Texas at Austin.