Dynamics of equilibrium-linked colloidal networks

Colloidal gels are commonly formed by quenching a dispersion with a low volume fraction, rapidly increasing the strength of short-range interparticle attractions relative to the thermal energy scale. For deep quenches, a strong thermodynamic driving force for the macroscopic phase separation into colloid-poor and colloid-rich fluid phases competes with glassy frustration and increasing structural relaxation times.^1,2 As the latter exceeds the time scale for observation, a dynamically arrested percolated colloidal network forms.^3,4 Structure and structure-dependent (e.g., rheological and optical) properties of these nonequilibrium gels depend on the colloidal interactions and the quench protocol, as does the gel aging process.^5–8 

Colloids with orientation-dependent interactions that allow bonding to only a few neighbors (e.g., patchy particles) can alternatively form equilibrium gels.^9–11 The fluid–fluid demixing transition for particles with such a microscopic valence restriction occurs at low volume fractions, opening up paths on the phase diagram between the demixing and glass transitions to directly form an equilibrium gel from a single-phase fluid. Equilibrium gels are soft and viscoelastic, exhibiting distinctive, spatially uniform dynamics that are stable in time.¹² The dynamic density correlations in an equilibrium gel decay via a characteristic two-step process.^12–14 The intermediate-time plateau between the two, which reflects a solid-like vibrational motion, has a height that is sensitive to temperature, increasing upon cooling as bonds form to extend the colloid network. The slow, secondary relaxation process, characterizing viscous flow, shows an Arrhenius temperature dependence with an apparent activation energy that is proportional to the strength of an attractive bond.^12–14 

Scalable synthesis of colloids with precisely defined valence poses formidable technical challenges, indicating a need for simpler and more generalizable strategies for assembling equilibrium gels. Linked colloidal assemblies, which restrict average colloidal valence, offer an attractive alternative that has motivated both modeling^15–25 and experimental^26–35 studies. Examples include nanocrystals or DNA nanostars with functionalized ligands and complementary linkers as well as emulsions with bridging telechelic polymers. Molecular thermodynamic theory and computer simulation of coarse-grained models have recently established that reducing the number of linkers per colloid can change the phase diagram in a way analogous to reducing the valence of equilibrium gel-forming patchy colloids.^16,18,36 Furthermore, a re-entrant percolation transition as a function of linker concentration was predicted for a model-linked colloidal network,¹⁶ mirroring the re-entrant gelation reported experimentally for colloidal dispersions with bridging species concentration.^37–39 These results suggest the need for a more complete understanding of how the dynamic properties of linked colloidal gel-forming systems relate to linker concentration, bonding motifs, and network structure.

Here, we use simulation to characterize the dynamics of the coarse-grained PolyPatch (Polymer linked, Patchy colloid) model,^18,20 in which the ends of ditopic molecules link neighboring nanoparticles into a network by attractively bonding to patches on their surfaces. Our simulations find dynamic signatures of a re-entrant equilibrium gel as a function of linker concentration, centered at stoichiometric conditions where the number of linker ends matches the number of nanoparticle surface bonding sites. The characteristic slow relaxation time of the self-intermediate scattering function scales with the colloid–colloid bond persistence time of the networks. As predicted from Wertheim theory,²¹ deviating from stoichiometry further restricts nanoparticle valence and the number of links per colloid that extend the colloidal network, producing linker-starved or linker-saturated networks. These structural changes correlate with reduced relaxation time scales and activation energies. Together, the results of this study help clarify how linked colloidal gels relate to patchy particle networks and suggest strategies for tuning the former’s properties via macroscopic control parameters, some of which may enable real-time, in situ dynamic control.

We adopt the PolyPatch model, a coarse-grained representation of dispersed colloids with surfaces decorated by six patchy sites that can attractively bond to the end beads of bifunctional linker molecules.^18,20,21,25 Here, the linkers are treated as flexible, linear chains of eight beads with diameter σ and mass m. Linkers, which also include bending stiffness or have different numbers of beads, can be chosen to tune the balance between bridging and looping in the PolyPatch model networks, and the corresponding implications of such linkers for phase behavior, structure, and structure-dependent mechanical and optical properties have been explored elsewhere.^18,20,25 Each colloid is modeled as a spherical particle with diameter σ_c = 5σ and mass m_c = 125m; its six interaction sites of diameter σ are rigidly fixed in an octahedral arrangement a distance r ≈ 3.12σ from the colloid center.

All particles in the PolyPatch model interact with the hard sphere-like repulsions of a core-shifted Weeks-Chandler-Andersen potential,⁴⁰ i.e.,

where $β=(kBT)−1$⁠, k_B is the Boltzmann constant, r is the center-to-center separation of particles of types i and j, and δ_ij = (σ_i + σ_j)/2 − σ modifies the divergence to reflect the respective diameters, σ_i and σ_j. This repulsive interaction is truncated at $σij*=21/6σ+δij$⁠.

Two bonded monomers of a linker chain also interact by a finitely extensible nonlinear elastic potential,⁴¹ 

with βk_s = 30 σ⁻² and r_s = 1.5 σ being the commonly adopted parameters that avoid unphysical chain crossings⁴² These coarse-grained interactions correspond to short polymers in a good solvent with an average bond length of 0.97σ.¹⁸ 

Attractions between colloidal interaction sites and linker ends can be expressed as follows:

where ɛ_b is the bond strength. The decay and range of the attractive interaction in Eq. (3) ensure that two linker ends cannot bond to the same colloid site while the repulsions of Eq. (1) prevent two colloid sites from bonding to the same linker end.¹⁸ Experimentally, thermoreversible linker-site bonds have been realized using DNA,^43–49 dynamic covalent bonding,³² and a variety of other reversible and orthogonal chemistries that are available for use in the linker-mediated assembly.^30,50–53

Molecular dynamics simulations of the PolyPatch model in the canonical ensemble were carried out using lammps (22 August 2018).⁵⁴ An integration time step of 0.001 τ was chosen, where $τ=βmσ2$ is a characteristic time, and the temperature was maintained using a Langevin thermostat with a friction coefficient of 0.1 m/τ applied to each linker bead and colloid. To study how linker concentration impacts static and dynamic properties, the simulations included N_c = 1000 colloids and a range of linker-to-colloid ratios Γ ≡ N_l/N_c from 1 to 6. It is helpful to recast the latter as the ratio of linker ends to colloid sites Γ* = Γ/3, which varies from 1/3 to 2. This range encompasses the behavior of linker-starved and linker-saturated state compositions, i.e., conditions below and above the stoichiometric ratio of Γ* = 1, respectively. The dimensions of the periodically replicated cubic simulation cell of volume V were chosen to investigate colloid packing fractions $ηc=Ncπσc3/(6V)$ from 0.01 to 0.15.

For each choice of η_c and Γ*, configurations were thermalized for 0.1 × 10⁵ τ in the absence of linker-site attractions (βɛ_b = 0). Next, ɛ_b was increased (at fixed temperature) for 0.1 × 10⁵ τ using a linear ramp until ɛ_b = 10. After repeating the thermalization protocol at ɛ_b = 10, ɛ_b was again increased using a linear ramp of 0.8 × 10⁻³ ɛ_b/τ. At a given βɛ_b, molecular dynamics simulations were run until the potential energy converged; the required run time varied in the range 1.5 × 10³ τ–0.1 × 10⁷ τ, depending on βɛ_b and Γ*. Finally, production simulations were performed to characterize the microstructure and the dynamics, with typical run times in the range of 0.2 × 10⁵ τ–0.2 × 10⁷ τ for the state points considered.

The behavior of the static structure factor at the near-zero wavevector was used to estimate state points where macroscopic phase separation occurs spontaneously by spinodal decomposition. The partial static structure factor of the colloids was computed as follows:⁵⁵ 

Here, k = (2π/V^1/3)n is a wavevector in the cubic box of volume V, n is a vector of integers, and r_j is the position of the jth colloid. The angle brackets indicate an ensemble average, which was calculated as a temporal average. Specifically, we averaged S_cc(k) for the 22 smallest nonzero wavevector magnitudes k = |k| and extrapolated to k = 0 by fitting the data to a Lorentzian function as follows:⁵⁶ 

where ξ_cc is a correlation length. As in the previous work, we used S_cc(0) > 10 as an approximate indicator of the spinodal region of the phase diagram.^16,18,57 Phase separation is also thermodynamically favored for metastable state points that lie between the spinodal and the binodal (i.e., equilibrium phase coexistence boundary), where it occurs via an activated process.

To characterize the colloidal dynamics, we calculated the mean squared displacement of colloidal particles $(Δr)2(t)$⁠) as follows:

where r_j(t) is the position vector of the jth colloid at time t. We also computed the colloidal self-intermediate scattering function as follows:

For sufficiently strong bonding interaction βɛ_b, F_s(k*, t) decayed via a two-step relaxation, where k* is the position of the static structure factor’s primary peak. To estimate a characteristic structural relaxation time τ_s and the nonergodicity factor A_s, we fit the slow relaxation of the self-intermediate scattering to the stretched exponential function, $Fs(k*,t)=Asexp[−(t/τs)β1]$ (Fig. S1).

We computed the following two characteristic times for the bond dynamics: continuous bond lifetime τ_b and bond persistence time τ_c.^58–62 Here, a bond refers to a linker-mediated connection between two colloids. Continuous bond lifetime was obtained by analyzing the bond survival probability S_b(t), which can be expressed as follows:

Here, i and j are colloidal indices; s_ij(t₀, t) = 1, if colloids i and j are continuously bonded in the time interval [t₀, t₀ + t], and s_ij(t₀, t) = 0, otherwise. The bond survival probability data were fit with an exponential function, S_b(t) = exp[−(t/τ_b)], to obtain an estimate for the continuous bond lifetime τ_b (Fig. S2). A linker-mediated bond connecting the same pair of colloids that effectively persists, despite intermittent breaking and rapid reformation, will not be apparent from analyzing the bond survival probability alone.

To measure longer-lasting effective bonds between two colloids, we also compute the more permissive bond persistence time, which is obtained from the bond correlation function C_b(t):

Here, h_ij(t) = 1, if colloids i and j are bonded at a given time t, and h_i,j(t₀ + t) = 0, if the colloids are not bonded at that time (i.e., independent of intermittent bond formation or breaking in the interval [t₀, t₀ + t]). The bond correlation function was fit by a stretched exponential function, $Cb(t)=exp[−(t/τc)β2]$⁠, to obtain an estimate of bond persistence time (Fig. S3).

Results for all dynamic properties were averaged over five independent molecular dynamics trajectories at each state point. Bars plotted with the simulated data points in the figures indicate standard deviations.

Building on the successful application of Wertheim theory for describing the equilibrium properties of patchy-particle fluids,^9,63–66 we use a recently introduced extension^20,21 of Wertheim’s first-order thermodynamic perturbation theory^67–70 (TPT1) to estimate the equilibrium linker bonding motifs of the PolyPatch model in the Γ*–η_c plane. We refer to this extension—distinctive in its ability to account for double-bonding motifs^20,21—as TPT1+. This theory can be solved to readily estimate the fractions of linkers that (1) form a loop, (2) remain dangling with a single end bonded to a colloid, (3) represent a free linker with no ends bonded, (4) make an effective bond between two colloids to extend the network, or (5) form a redundant bond between two colloids that already share an effective bond [Fig. 1(a)]. A full description of the theory, including expressions for quantities needed for computing the number of each linker motif per colloid, can be found in Refs. 20 and 21.

Wertheim’s perturbation theory decomposes the Helmholtz free energy density, a, into a sum of a bonding portion a_bond and a reference contribution a₀. The former accounts for directional attractions between particles with interaction sites (i.e., linkers or colloids), whereas the latter accounts for the free energy of the system in the absence of such attractions. The bonding free energy density is expressed in terms of a series of variables, $Xα(i)$⁠, each representing the fraction of component i [i.e., colloid (c) or linker (l)] not bonded at a set of sites α. Once the $Xα(i)$ is determined using TPT1+, one can compute the fraction of linkers in each of the five linker motifs.

In the parlance of TPT1+, if A₁₂ = {A₁, A₂} is the set of two (representative) neighboring bonding sites on a colloid and B₁₂ = {B₁, B₂} is the set of bonding sites on a linker, then the fraction of free linkers (i.e., with neither of two ends participating in a bond) is $χfree=XB12(l)$⁠. Likewise, the fraction of dangling linkers (i.e., with either one, but precisely one, of its ends bonded) is $χdangling=2(XB1(l)−XB12(l))$⁠. We denote the fraction of linkers with both ends bonded $χboth=1−χdangling−χfree=1−2XB1(l)+XB12(l)$⁠.

The fraction of linkers forming loops χ_loop and redundant bonds χ_redundant may be solved using relations formed in the derivation of TPT1+,²¹ 

where ν_c is the number of interaction site pairs on a colloid capable of participating in a loop or ring, ρ_i is the number density of component i, and Δ₂ and Δ₄ are bond volumes introduced to account for loop and ring structures.²⁰ The fraction of effective bonds χ_effective can be found with the relation,

These fractions can be transformed from a per linker basis to a per colloid basis by multiplying by Γ. One quantity of particular interest is ⟨n_eb⟩ = χ_effectiveΓ, the number of effective bonds per colloid, which is related to a colloid’s average coordination number ⟨CN₁⟩ = 2⟨n_eb⟩. CN₁ quantifies the number of colloids that a given colloid is directly bonded to via linkers.

By definition, a single effective bond joins each pair of linked colloids, connecting sites on their respective surfaces. Thus, although individual linkers can form loops on a colloid surface or form a redundant bond, TPT1+ treats the clusters of colloids linked by effective bonds as loop-less networks,²¹ consistent with the Flory–Stockmayer theory of gelation.^71,72 The presence of other linker bonding motifs (e.g., dangling, redundant-bonding, looping, and free linkers) only serves to reduce the probability p_b(βɛ_b, η_c, Γ*) = ⟨n_eb⟩/3 that a site participates in an effective bond. Flory–Stockmayer theory relates the value of this probability at the percolation line, $pbp$⁠, to the number of sites per colloid f via the expression $pbp=(f−1)−1$⁠. Thus, the TPT1+ and Flory–Stockmayer theory estimate for the percolation threshold of the six-patch PolyPatch model corresponds to the condition ⟨n_eb⟩ = 0.6, which can be used to locate the percolation boundary (e.g., in Γ* − η_c or βɛ_b − η_c planes).

As illustrated by the PolyPatch model's behavior, linker and colloid concentrations provide macroscopic handles for tuning network properties of linked-colloidal assemblies. A heat map of ⟨n_eb⟩ in the Γ*–η_c plane at βɛ_b = 18 predicted by TPT1+ shows a region of highly bonded network structures centered on the stoichiometric ratio Γ* = 1 for η_c of ∼0.10 and higher [Fig. 1(b)]. The correspondence between the colors inside of the open symbols and the background reflects the excellent agreement between the theoretical predictions and simulated results for ⟨n_eb⟩. At lower η_c and in the vicinity of Γ* = 1, effective bonding can be satisfied only by macroscopic phase separation, whereas at higher η_c, a single phase with sufficient bonding to percolate is found, consistent with the expectation of equilibrium gelation. Crossing contours of ⟨n_eb⟩ at constant η_c corresponds to re-entrant network formation that is approximately symmetric in ln Γ*. Lower values of ⟨n_eb⟩, and consequently the colloidal coordination number, that attain for conditions away from the stoichiometric ratio highlight that colloidal valence becomes restricted at low Γ* due to linker starving (i.e., too few linkers per colloid), whereas it is limited at high Γ* due to a lack of unoccupied colloidal sites to form effective bonds (i.e., linker saturation).

To further probe the microstructural ramifications of changes to linker concentration, we computed the number of each of the five linker bonding states shown in Fig. 1(a) per colloid as a function of Γ* at η_c = 0.10 [larger squares in Fig. 1(b)]. The changes in the various bond motifs per colloid predicted by TPT1+ are in good quantitative agreement with those computed via molecular simulation [Fig. 2(a)]. In the linker starved region (Γ* < 1), there is no shortage of available colloidal interaction sites, and thus, adding linkers increases the number of connections between colloids and extends the linked network. On the other hand, in the linker-saturated region (Γ* > 1), adding linker chains decreases the number of linkers bridging colloids, displacing them to form dangling linkers.

In addition to these modifications to the linker bonding motifs, changing Γ* has a pronounced effect on the colloidal coordination number CN₁. Histograms showing the simulated fractions of colloids with different coordination numbers, χ(CN₁), are unimodal and shift non-monotonically with Γ*. This trend reflects re-entrant structural modifications as the system passes through the stoichiometric ratio, where the colloids have the highest average coordination number [Fig. 2(b)]. The peak positions of the distributions closely track their average values ⟨CN₁⟩ = 2⟨n_eb⟩, further highlighting that ⟨n_eb⟩ captures key aspects of network connectivity.

To ascertain how much the bonding state of a colloid influences that of its neighbors, we computed ⟨CN₂⟩, the average coordination number of colloids linked to a central colloid with CN₁ [Fig. 2(c)]. To a first approximation, ⟨CN₂⟩ ≈ ⟨CN₁⟩, again following the re-entrant trends of ⟨n_eb⟩ in Fig. 2(a) as Γ* is increased. For each Γ*, there is a weak, positive correlation between ⟨CN₂⟩ and CN₁; neighbors of a more highly coordinated particle are themselves bonded to slightly more neighbors than the neighbors of low-coordinated colloids. This positive correlation, reflecting network cooperativity, is strongest at the stoichiometric ratio and flattens for Γ* in linker-starved or linker-saturated regions.

To study how linker concentration influences colloidal relaxation processes in the PolyPatch model, we computed the mean squared displacement ⟨(Δr)²⟩ and the self-intermediate scattering function F_s(k*, t) of the colloidal particles, where k*σ = 1 corresponds to the approximate position of the colloidal structure factor’s first peak (Fig. 3). At βɛ_b = 22 and Γ* = 0.33, i.e., a strong-bonding but linker-starved condition, ⟨(Δr)²⟩ shows equilibrium fluid-like behavior, crossing over from short-time dynamics to long-time diffusive scaling without exhibiting a plateau at intermediate times. The corresponding intermediate scattering function displays a fast, single-step decay. As Γ* is increased toward the stoichiometric ratio, plateaus form in ⟨(Δr)²⟩ and F_s(k*, t) at intermediate times, whereas a second, viscous relaxation process becomes evident at long times. The height of the plateau of F_s(k*, t), commonly referred to as the nonergodicity parameter A_s also shows pronounced increases as Γ* approaches one and the linked-colloidal network becomes more extensive. These trends are reminiscent of those observed experimentally in emulsions with increasing telechelic polymer concentration²⁶ or in DNA nanostars dispersions with increasing linker concentration.¹⁷ However, as even more linkers are added, creating linker-saturated conditions, the viscoelastic characteristics of the single-particle dynamics attenuate and the relaxation processes again approach simple fluid-like behavior.

Linker concentration and strength of attraction between linker ends and colloidal surface sites play distinctive roles in tuning dynamics of colloidal networks (Fig. 4). At fixed βɛ_b, changing Γ* from 0.67 to 1 increases both A_s and the characteristic slow relaxation time τ_s by over an order of magnitude, dynamic effects that are reversed upon increasing Γ* further to 1.67 [Figs. 4(a) and 4(b)]. Underlying this re-entrant behavior is a monotonic change in A_s with ⟨n_eb⟩ along this path; higher values of A_s emerging in strongly interconnected networks. Varying βɛ_b at fixed Γ* reveals an Arrhenius relationship for the slow relaxation, τ_s = τ₀ exp(Cβɛ_b), with activation energy proportional to the bond strength, i.e., Cɛ_b. These temperature-dependent dynamics are consistent with those observed in equilibrium gel-forming fluids of particles with microscopically restricted valence, e.g., DNA nanostars and patchy colloids,^13,73–75 as well as in strong glass-forming liquids.⁷⁵ 

Similar to the nonergodicity parameter, the activation energy prefactor C shows re-entrant behavior with Γ*, assuming lower values (C ≈ 0.9) in both the linker- and linker-saturated states and higher values (C ≈ 1.2) at the stoichiometric ratio. These trends mirror those observed in patchy colloids and other model network-forming fluids.^74–79 Activation energies in the latter can vary from 0.5ɛ to 2ɛ, with the larger values corresponding to more geometrically constrained networks where particle motion requires correlated, simultaneous breaking of multiple bonds.^{74,75,80–82}

To gain further insights into how changes in Γ* lead to different dynamic behaviors, a self-intermediate scattering function analysis was also carried out at βɛ_b = 22 for individual subpopulations of colloids with CN₁ = i with i ∈ {2, 3, 4, 5, 6} at t = 0. (Colloidal particles with CN₁ = 1 do not significantly contribute to two-step relaxation.) For simplicity and because, as discussed below, bonds persist for long times under these conditions, the initial subpopulation labels are maintained throughout this calculation. In addition to increasing the average coordination number of colloids and their neighbors [Figs. 2(b) and 2(c)], tuning Γ* toward the stoichiometric ratio also produces coordination-number sensitive relaxation times [Fig. 4(c)]. Similar coordination number-dependent single-particle dynamics were reported for the simulated tetrahedral network forming patchy particles.⁷⁴ In both the PolyPatch model at Γ* = 1 and the tetrahedral patchy particle fluid,⁷⁴ bond breaking is essential for colloidal translation (Fig. S4a). However, in linked-colloidal systems, modifying Γ* to linker-starved or linker-saturated conditions with fewer effective bonds per colloid substantially weakens the link between dynamics and coordination number [Fig. 4(c)], and at Γ* = 1.67, there are no detectable differences in relaxation times of colloidal particles with different CN₁. The reduced connectivity of the colloidal network under these conditions renders it less constrained, and significant colloidal displacements can occur without even rupturing bonds with the neighbors (Fig. S4b).

A comparison of continuous bond lifetime τ_b (i.e., time for a colloid–colloid bond to break) and bond persistence time τ_c (i.e., time that a colloid–colloid bond lasts allowing intermittent breaks) for the PolyPatch model provides insight about the microscopic processes underlying differences in structural dynamics as Γ* varies. Both times are found to exhibit Arrhenius temperature dependencies for the full range of conditions examined in this study [Figs. 5(a) and 5(b)]. However, the former is found to depend only on βɛ_b and not on Γ*, which determines the connectivity of the network and the constraints it imposes on the colloids. By contrast, the bond persistence time τ_c shows a Γ*-dependent activation energy, which is larger at the stoichiometric ratio than under linker-starved or linker-saturated conditions. This makes physical sense; at a given βɛ_b, looser, less interconnected networks allow colloids to diffuse away from a bonding partner after a linker end detaches from a colloidal site, although more geometrically constrained networks tend to promote bond reformation between the same pair. A direct comparison of data in Figs. 4(b) and 5(b) shows τ_c ∝ τ_s [Fig. 5(c)], highlighting that Γ*-dependent persistence of a colloid–colloid bond (allowing for intermittent bond breaks) controls the slow relaxation of the self-intermediate scattering function in these viscoelastic networks.

By studying a coarse-grained model, we explored how linked colloidal networks can exhibit dynamic hallmarks^12–14 of equilibrium-gel forming colloids with microscopically restricted valence. These dynamic properties include a characteristic two-step relaxation of the self-intermediate scattering function with a nonergodicity parameter that increases as gelation is approached within a single fluid phase (i.e., avoiding phase separation by spinodal decomposition) and a structural relaxation time with an Arrhenius temperature dependency. Unlike in single-component patchy particle networks, linker concentration provides a macroscopic handle in linked-colloidal gels for tuning network connectivity and viscoelastic dynamics. Re-entrant properties are observed as linker concentration is tuned away from the stoichiometric ratio where the number of linker ends matches the number of colloidal sites. Decreasing linker concentration reduces the connectivity of the network by starving the system of linkers that connect sites on neighboring colloids. Increasing linker concentration has a similar effect by displacing effective bonds connecting two colloids with dangling linkers. Through a bond lifetime analysis, we demonstrated that the colloid–colloid bond persistence time (with intermittent bond breaking) controls the characteristic slow relaxation of the self-intermediate scattering function.

Although qualitatively similar, the re-entrant network formation of the PolyPatch model is distinctive compared with the re-entrance exhibited by systems where the linking species nonspecifically binds to the colloid surfaces, bridging the particles.^37,38 In the latter, “stoichiometry” is satisfied when there are enough linkers present to saturate the surfaces of the colloids. Here, the number of functional sites on the colloid surface can be deliberately and discretely modified—either synthetically or by adding a capping species that binds to colloidal sites, displacing connecting linkers—to adjust the stoichiometric ratio where re-entrance occurs. Furthermore, the site–site separation on the colloid surface (relative to the linker’s end-to-end distance) can be tailored to influence the competition between linker bridging and loop formation, providing additional control over network formation.²⁰ Finally, in experiments, functionalizing the surface of colloids with chemically specific sites allows for customizing interactions to be complementary with specific linkers, e.g., to program assembly or co-assembly of gels.³³ 

In DNA nanostar gels, a capping strategy has been used to produce temperature-re-entrant networks.⁸³ Synthetic methods for creating dynamically activated (e.g., light-triggered) caps may further allow in situ control of network connectivity and properties in linked colloidal gels. Given the accuracy of TPT1+ for predicting the bonding motifs of the PolyPatch model, the theory should prove to be a powerful tool for predicting how the presence of caps influences the underlying equilibrium properties of such networks.

Following theoretical progress in patchy particle systems,^10,24,66,84 it would be interesting to use molecular simulation to systematically test TPT1+ and Flory–Stockmayer predictions for cluster-size distributions, and Flory theory estimates for the equilibrium gel transition across the broad Γ*–βɛ_b–η_c parameter space. It would also be helpful to extend the PolyPatch model to more accurately treat the dynamic covalent bonding strategies that are being implemented in experimental linked-colloidal networks.^30,32–35 Specifically, incorporating the thermodynamics and kinetics of the reactions that govern bonding between linker ends and colloidal surface sites (e.g., ligands) would allow the model to provide insights into how macroscopic rheological properties could be programmed using different dynamic covalent chemistries.^85,86

The supplementary material contains additional simulation data for analysis of the self-intermediate scattering functions, continuous bond lifetime distributions, bond correlation functions, and colloidal bond states during displacement.

This research was primarily supported by the National Science Foundation through the Center for Dynamics and Control of Materials: an NSF MRSEC under Cooperative Agreement No. DMR-1720595, with additional support from the Welch Foundation (Grant Nos. F-1696 and F-1848). We acknowledge the Texas Advanced Computing Center (TACC) at The University of Texas at Austin for providing HPC resources.

The authors have no conflicts to disclose.

Taejin Kwon: Conceptualization (equal); Data curation (equal); Formal analysis (equal); Investigation (equal); Methodology (equal); Validation (equal); Writing – original draft (equal); Writing – review & editing (equal). Tanner A. Wilcoxson: Conceptualization (equal); Data curation (equal); Formal analysis (equal); Investigation (equal); Methodology (equal); Validation (equal); Writing – original draft (equal); Writing – review & editing (equal). Delia J. Milliron: Conceptualization (equal); Funding acquisition (equal); Supervision (equal); Writing – original draft (equal); Writing – review & editing (equal). Thomas M. Truskett: Conceptualization (lead); Funding acquisition (lead); Project administration (lead); Supervision (lead); Writing – original draft (equal); Writing – review & editing (equal).

The data that support the findings of this study are available from the corresponding author upon reasonable request.