Low-density “equilibrium” gels that consist of a percolated, kinetically arrested network of colloidal particles and are resilient to aging can be fabricated by restricting the number of effective bonds that form between the colloids. Valence-restricted patchy particles have long served as one archetypal example of such materials, but equilibrium gels can also be realized through a synthetically simpler and scalable strategy that introduces a secondary linker, such as a small ditopic molecule, to mediate the bonds between the colloids. Here, we consider the case where the ditopic linker molecules are low-molecular-weight polymers and demonstrate using a model colloid–polymer mixture how macroscopic properties such as the phase behavior as well as the microstructure of the gel can be designed through the polymer molecular weight and concentration. The low-density window for equilibrium gel formation is favorably expanded using longer linkers while necessarily increasing the spacing between all colloids. However, we show that blends of linkers with different sizes enable wider variation in microstructure for a given target phase behavior. Our computational study suggests a robust and tunable strategy for the experimental realization of equilibrium colloidal gels.

## I. INTRODUCTION

Gels consisting of percolated, porous, and kinetically arrested networks of colloidal particles connected by chemical or physical bonds^{1,2} possess a remarkable combination of material properties: even at extremely low particle densities, they can be self-supporting and exhibit solidlike mechanical response.^{1–3} Low density colloidal dispersions are less prone to aggregation^{4} and more likely to be able to kinetically access deliberately assembled gel structures. In addition, many unique characteristics of nanoparticle dispersions are strongly perturbed in dense arrangements,^{5} motivating the development of strategies to control and limit their density even while assembling them. For example, coupling interactions between plasmonically active metal nanoparticles cause spectral shifts and broadening. Controlling the extent of such interactions in gels could yield a wide range of structure-dependent optical responses, e.g., a low-density metal oxide nanocrystal gel was recently shown to retain localized surface plasmon resonance characteristic of isolated nanocrystals but with a spectral shift triggered by the assembly process.^{6}

Most routes to gelation are nonequilibrium in nature, e.g., the spinodal decomposition of attractive colloids is kinetically arrested.^{7–9} Such gels tend to have uncontrolled local density variations and age as their microstructures slowly relax toward equilibrium phase separation. In contrast, “equilibrium” gels^{2,10–16} that are resilient to aging can be formed when the arrested structure is characteristic of a true equilibrium state (i.e., a snapshot of a single phase). Gelation occurs when the lifetime of the effective bonds between colloids in this state significantly exceeds the observation time.^{3} By restricting the number of these bonds, the spinodal boundary for spontaneous phase separation can be suppressed, opening a window of fluid states that can gel at small colloid volume fractions.^{2,13,17–19} This concept was originally demonstrated in Monte Carlo simulations for a model where the number of bonds between colloids was artificially constrained^{13} and for an anisotropic (“patchy”) colloid model.^{17–19} These computational studies motivated the discovery of several materials that can form equilibrium gels, including an anisotropic colloidal clay,^{20} DNA nanostars with complementary sticky ends,^{21–23} and a dipeptide with several aromatic rings that allowed for specific *π*–*π* stacking interactions.^{24}

Nevertheless, the preparation of equilibrium gels from patchy nanoparticles is challenged by difficulties in fabricating the latter with sufficient precision and in large enough quantities.^{25} This challenge motivated some of us to study an alternative route to equilibrium gelation that controls the extent of bonding between the (primary) colloids using a (secondary) linking agent.^{26} Many examples of such linking systems have been described, including colloid–polymer mixtures,^{27–31} mixtures of oppositely charged colloids,^{32,33} ion-linked functionalized nanocrystals,^{34–37} DNA-linked nanoparticles,^{38–40} and DNA nanostar mixtures.^{41–43} The degree of bonding between the colloids can be restricted by appropriately tuning the amount of linker added. Such macroscopic control is more readily achieved in experiments than careful microscopic control of the number of patches on a nanoparticle and so is a promising strategy for fabricating low-density, equilibrium colloidal gels.

We previously demonstrated using computer simulations that the phase behavior of a mixture of equally sized spherical colloids and linkers could be straightforwardly controlled through the linker concentration.^{26} However, at the time, we did not investigate how the design of the linker itself impacted gelation. Many linking agents, such as ditopic small molecules^{44} or DNA linkers,^{38,43} have synthetically tunable properties like their molecular size or flexibility that can be varied independently of the primary particle properties. The strength and range of the effective bonds formed by the linker between colloids should depend on these characteristics. We anticipate that the phase behavior, microstructure, and emergent macroscopic properties of a colloid–linker gel should be highly sensitive to the choice of the linker, yielding a rich design space for tailoring the gel properties.

In this article, we interrogate the phase behavior and microstructure of model equilibrium gels assembled from a mixture of spherical colloids and flexible-chain molecular (polymeric) linkers. Using a combination of thermodynamic perturbation theory (TPT) and computer simulations (Sec. II), we demonstrate how the linker length and concentration can be used to systematically modify the spinodal boundary for the mixture, which determines the colloid volume fractions at which equilibrium gelation can occur (Sec. III A). The simulations reveal that many linkers bind both ends to the same particle or form double bonds between particles, especially at a low colloid density, and this looping can inhibit percolation of the gel network. We also propose and demonstrate a strategy for simultaneously tuning the phase behavior and microstructure of the gel using blends of linkers with different lengths (Sec. III B). Macroscopic properties of the gel, like its equilibrium modulus, are shown to have a concomitant dependence on the distribution of linker lengths (Sec. III C).

## II. MODEL AND METHODS

To form an equilibrium gel, it is necessary to create a homogeneous, percolated network of linked colloids that subsequently undergoes an effective kinetic arrest.^{2,3} The region of state points where these conditions can be satisfied depends sensitively on the interactions between components in the mixture. The mixture may separate into two or more coexisting thermodynamic phases at low colloid concentrations or enter a single-phase, arrested state at high concentrations.^{1,2} Effective kinetic arrest is usually only achieved when thermal energy is sufficiently low (relative to bond strength) that the average lifetime of a bond in the percolated network becomes comparable to the experimental observation time; this point is typically gauged relative to the highest temperature for which phase separation occurs. Thus, it is important to determine the thermodynamic phase boundaries governing these transitions and to assess how they depend on the tunable linker properties.

To this end, we studied a simple model (Sec. II A) for a colloid–linker mixture that captures these salient physical behaviors using a combination of molecular simulations (Sec. II B) and thermodynamic perturbation theory (Sec. II C). The theoretical calculations predict equilibrium phase behavior under a reasonable set of assumptions, and these predictions can be validated by the more computationally demanding simulations. The simulations additionally provide detailed microstructural information and incorporate practical complications neglected by the assumptions of the theory.

### A. Model

The linkers were modeled as short linear chains of *M* segments with diameter *σ* and mass *m*, which we refer to as polymers. The segment interactions were purely repulsive to mimic good-solvent conditions for the polymer.^{45} The colloidal particles were represented by larger spheres of diameter *σ*_{c} = 5 *σ* and mass *m*_{c} = 125 *m*. The interactions between colloids and with the polymer segments were also purely repulsive, corresponding to good stabilization of the colloids in the (implicit) solvent and to neglecting polymer adsorption onto the colloid surface.

Polymer-mediated linking between colloids was incorporated by decorating the colloid surface with *n*_{c} = 6 attractive patches with nominal diameter *σ* and mass *m* in an octahedral arrangement (Fig. 1). The end segments of the polymer chains had a short-ranged attraction to these patches to mimic bond formation. In prior work that investigated the effect of limiting the number of bonds per particle using a single-component square-well fluid, it was shown that the phase behavior of particles that can form at least six bonds is similar to when the number of bonds per particle is not constrained.^{13} Additional details about the patch model and interaction are given in Sec. II B.

The colloid–polymer mixture is characterized by the temperature *T*, the volume *V*, and the number of colloids *N*_{c} and polymers *N*_{p}. The total number density of the mixture is *ρ* = (*N*_{c} + *N*_{p})/*V*, and the number density for each species is *ρ*_{i} = *x*_{i}*ρ*, where *x*_{i} = *N*_{i}/(*N*_{c} + *N*_{p}) is the mole fraction for component *i*. It is convenient to describe the composition using the colloid volume fraction $\eta c=\rho c\pi \sigma c3/6$ and the number ratio of polymers to colloids, Γ = *N*_{p}/*N*_{c}. Figure 1 shows an example snapshot of the full model, including the colloids with their patches and the polymer linkers, for *η*_{c} = 0.10 and Γ = 1.5.

### B. Simulations

We performed molecular dynamics simulations of the colloid–polymer mixtures at constant temperature *T* using lammps (16 March 2018).^{47} The interparticle interactions were modeled using a core-shifted Weeks–Chandler–Andersen potential,^{48}

where *r* is the distance between the centers of two particles *i* and *j*, *ε* = *k*_{B}*T* sets the energy scale of the repulsion with *k*_{B} being Boltzmann’s constant, and *δ* = (*σ*_{i} + *σ*_{j})/2 − *σ* shifts the divergence of the potential to account for the particle diameters, *σ*_{i} and *σ*_{j}. This potential was truncated at *r*^{*} = 2^{1/6}*σ* + *δ* to make the interactions purely repulsive.

Bonds between polymer segments were modeled using finitely extensible nonlinear elastic springs,^{49}

with the standard Kremer–Grest parameters for the spring constant (*k*_{p} = 30 *k*_{B}*T*/*σ*^{2}) and maximum bond extension (*r*_{0} = 1.5 *σ*).^{50} These parameters correspond to a polymer in a good solvent and ensure that there is no artificial chain crossing; the average bond length was 0.97 *σ*.

Each colloid and its patches were modeled as a rigid body. The patches were placed in an octahedral arrangement on the colloid surface at a radial distance *r*^{*} from the center of the colloid so that a segment bonding with the patch was not repelled by the colloid (*δ* = 2 *σ*). The end segments of the polymer chains were attracted to these patches through a short-ranged potential,

where *ε*_{b} sets the bond strength. The length scale of attraction for *u*_{b} was chosen so that two polymer segments could not bond to the same patch. We additionally allowed the patches to interact with each other through Eq. (1) to ensure that multiple patches could not coordinate with the same end segment. This computationally convenient scheme mimics the formation of an exclusive covalent bond between the patch and a segment, especially as *ε*_{b} becomes large compared to the thermal energy *k*_{B}*T*, which is the limit relevant to typical experimentally realizable linker chemistries.^{43,44,51}

We prepared configurations containing *N*_{c} = 1000 colloids and *N*_{p} = 1500 polymers (Γ = 1.5) in a cubic, periodic simulation box of edge length *L* at colloid volume fractions ranging from *η*_{c} = 0.01 to 0.15. We initially equilibrated each mixture for 0.5 × 10^{4} *τ* without any patch bonding (*βε*_{b} = 0), where $\tau =\beta m\sigma 2$ is the unit of time and *β* = 1/*k*_{B}*T*. The integration time step was 0.001 *τ*, and the temperature was controlled by a Langevin thermostat with friction coefficient 0.1 *m*/*τ* applied to each particle. We then followed a slow annealing protocol to switch on the patch attractions in an attempt to produce equilibrium gel structures. For each volume fraction, the bond energy was first linearly ramped to *βε*_{b} = 10 over 10^{4} *τ*, followed by a 10^{4} *τ* equilibration period. The bond energy was then increased by *k*_{B}*T* over a 0.5 × 10^{4} *τ* period, and the mixture was allowed to re-equilibrate for 1.5 × 10^{4} *τ*. The final configuration from this period was saved, and the incrementing process was repeated until *βε*_{b} = 20. We took the configurations generated for each (*η*_{c}, *ε*_{b}) pair and used them to sample the structure of the mixture during production simulations of 10^{4} *τ*, saving configurations every 10 *τ* for subsequent analysis.

### C. Theory

Wertheim’s thermodynamic perturbation theory (TPT)^{52–54} is unrivaled in its ability to predict the equilibrium phase behavior of strong, short-ranged, highly directional, bond-forming interactions in particle-based systems. Its predictions are often in close qualitative and quantitative agreement with the exact simulations of patchy particle models.^{10,12,14–17,26} Formally, TPT is derived from a graph-theoretic expansion and a partial resummation centered around a few key physical assumptions:^{52–55} (1) each patch on a particle can only bond with one other patch, (2) bond formation is statistically uncorrelated between patches on a particle, and (3) bonding occurs in a treelike structure (i.e., closed loops are not permitted). As a result of these assumptions, TPT has a relatively simple closed-form expression for the free energy, rendering thermodynamic stability and equilibrium phase coexistence calculations computationally inexpensive.

Within TPT, the Helmholtz free energy *a* is decomposed into a reference contribution *a*_{0}, devoid of patchy interactions, and a bonding contribution *a*_{b} that adds in the effect of patches,

Here, the free energies are made intensive by the total number of colloids and polymers, *N* = *N*_{c} + *N*_{p}. We approximate the reference free energy *a*_{0} as that of a hard-chain mixture assembled from a dissociated hard-sphere fluid.^{53} The linear chains of species *i* are made up of *M*_{i} hard spheres with diameter *σ*_{i}. The free energy, itself derived using TPT, is

where *a*_{hs} is the free energy of an equivalent hard-sphere mixture obtained by fully dissociating all chains, and $gii(\sigma i+)$ is the corresponding contact value for the radial distribution function of species *i* in the hard-sphere fluid. Practically, we employ Boublík’s well-known and accurate equation of state^{56} for *a*_{hs} and $gii(\sigma i+)$.

The bonding contribution to the free energy is

where *n*_{i} is the number of patches and *X*_{i} is the fraction of patches that are not bonded for species *i*. The values of *X*_{i} are determined by a coupled set of chemical equilibrium equations,

with Δ_{ij} being an effective bond partition function. For our model, where only polymers bond to colloids,

which depends on both the radial distribution function *g*_{cp} between the colloid (c) and polymer segment (p) and the bonding Mayer f-function^{52,53,57}*f*_{b} averaged over all angular orientations of the colloid and polymer segment at a center-to-center distance *r* = |**r**|. As a further simplification, we assume that the range where *f*_{b} is nonzero is short and take $gcp(r)\u2248gcp(\sigma cp+)$ for a reference hard-sphere fluid at contact distance *σ*_{cp} = (*σ*_{c} + *σ*)/2. Details regarding the evaluation of Eq. (8) are provided in the Appendix.

To divide the parameter space into homogeneous (single-phase fluid) and phase-separated regimes, we computed the spinodal and binodal boundaries for the mixture. The spinodal demarcates the conditions under which the single-phase fluid becomes unstable toward spontaneous decomposition into multiple phases; true equilibrium phase separation at the binodal always pre-empts the spinodal, and the single-phase fluid is only metastable in the region between the binodal and spinodal.^{58,59} At fixed *N*, *V*, and *T*, the spinodal is determined by the locus of points for which the determinant of the stability matrix **H** is zero, with *H*_{ij} = *∂*^{2}*a*/*∂x*_{i}*∂x*_{j} being the second derivative of the Helmholtz free energy with respect to composition.^{58,59}

The binodal is a complex, multidimensional construct containing multiphase information. In this work, we computed it in the plane of constant composition by postulating the existence of a two-phase system and minimizing its total free energy subject to particle and volume conservation between phases. (This procedure is equivalent to solving for temperature, chemical potential, and pressure equality between the two phases.) The binodal is defined by the conditions where an infinitesimal amount of an incipient phase emerges in coexistence with the homogeneous phase, also known as a cloud curve.^{12,26} Further inside this multiphase region, the incipient phase becomes finite and demixing occurs, leading to two phases with different compositions than that of the corresponding homogeneous state. Additional details of this procedure are outlined in the Appendix of Ref. 26.

## III. RESULTS AND DISCUSSION

### A. Phase behavior and linker length

We first investigated the phase behavior of mixtures of colloids with linkers of a single length. Unlike our prior work with spherical linkers,^{26} the polymer linkers are flexible and possess internal degrees of freedom, leading to softer effective interactions between linked colloids. Because the linkers in the prior work occupied the same excluded volume as the primary colloids, the polymer linkers considered here also occupy a lower volume fraction. Hence, we hypothesized that the colloid–polymer mixtures might exhibit different thresholds for network percolation and phase separation than the sphere mixtures and that the phase behavior may depend on the polymeric linker length.

The formation of a percolated, self-supporting network is the first requirement for gelation. Percolation only occurs once a sufficient amount of linker has been added at cold enough temperature (high enough bond energy) and can be predicted using TPT as the conditions under which there is unbounded growth of the network of effective bonds between colloids. In the limit β*ε*_{b} → ∞, relevant to strong bonds between the linkers and colloids, the percolation threshold depends simply on the patch stoichiometry and is estimated to be Γ = 0.6 for our model. In previous simulations, though, we found that percolation did not occur until at least Γ = 0.9 for the spherical linkers.^{26}

Based on this, we initially performed simulations at Γ = 1.0 but were unable to form a percolated network even at colloid volume fractions predicted to be outside the phase separated region. In the simulations, percolation was defined as the formation of a network spanning the periodic simulation box in at least one dimension. We speculated that this discrepancy was caused by the linkers forming “loops” in the colloid network, which are neglected by TPT. Redundant linker bonds create cycles in the network graph having length equal to the number of linkers in the loop and lead to a smaller number of effective connections between colloids than expected. When we increased the linker concentration to Γ = 1.5, we were able to obtain a percolated network.

We can characterize the number of linkers that do not form new connections in the colloid network directly from these simulations. We considered the two shortest looping motifs, which we expected to occur in largest number and so to have the greatest impact on the phase behavior of the system: a cycle of length one, where a linker bonds back onto the same colloid, and a cycle of length two, where two linkers form a double bond between two colloids. We counted the fraction of linkers *f*_{p} that formed redundant bonds from these motifs for linker lengths *M* = 2 and 8 at Γ = 1.5 as a function of *ε*_{b} (Fig. 2). Only the *M* = 8 linkers were long enough to bond two patches on the same colloid, but both linkers were able to form double bonds. We considered both low (*η*_{c} = 0.01) and high (*η*_{c} = 0.15) colloid densities, which TPT predicted to be phase separated and homogeneous, respectively.

For low *ε*_{b}, there was a negligible number of redundant bonds in the network, although the total number of bonds was also lower. As *ε*_{b} increased, *f*_{p} slowly increased until a sharp rise occurred near *βε*_{b} ≈ 15. This rise in *f*_{p} happened at lower *ε*_{b} for the *M* = 8 linkers, where most of these redundant bonds initially came from the linkers bonded to the same colloid. The total number of redundant bonds at *βε*_{b} = 20 was much larger in the low density mixture (20%–40% of linkers) than in the high density mixture (10% of linkers) and was also larger for the longer linkers at *η*_{c} = 0.01. However, because *f*_{p} was similar for both mixtures when *βε*_{b} was smaller and also at *η*_{c} = 0.15, this increased tendency to form redundant bonds at larger *ε*_{b} and *M* for lower *η*_{c} may be due to the onset of phase separation. The redundant bonds likely inhibited percolation compared to TPT predictions and are a real effect that must be considered in an experimental realization of linker gels.

We then determined approximate phase diagrams of the colloid–polymer mixtures from the simulations in order to compare with TPT’s predictions. It is challenging to rigorously compute these diagrams from the free energies of competing phases, and so we instead looked for signs of phase separation emerging within our finite simulation box and finite observation time. We computed the partial static structure factor for the colloids,^{57}

where **q** = 2*π***n**/*L* is the wavevector, **n** is a vector of integers, and **r**_{j} is the position of the *j*th colloid. We determined *S*_{cc} for the 22 smallest nonzero wavevector magnitudes *q* = |**q**|, averaging over the wavevectors corresponding to each *q*. We then extrapolated the structure factor to zero-wavevector, *S*_{cc}(0), by fitting *S*_{cc} to a Lorentzian form^{60}

with *ξ* being a correlation length. In the thermodynamic limit, *S*_{cc}(0) should diverge at the spinodal, but, in practice, it remains finite in the simulation and *S*_{cc}(0) > 10 is usually considered to be phase separated.^{13,26} We qualitatively confirmed that this criterion held for our model by visually inspecting configurations that it identified as phase separated, finding that they had clear density variations.

We compared the measured structure factor to the TPT predictions of phase coexistence for two chain lengths, *M* = 2 and *M* = 8, with Γ = 1.5 (Fig. 3). The theoretical spinodal boundary (solid line) indicates the state points for which the system will spontaneously demix, while the binodal boundary (dotted line) fully encompasses the region where phase separation can occur, including metastable regions. Phase separation should nearly always occur in the simulations for state points inside the spinodal, but may or may not be observed in the metastable region between the binodal and spinodal due to the finite simulation time.

Overall, the simulations and TPT were in good agreement, particularly with respect to *η*_{c}. Accurately predicting the phase boundary with respect to *η*_{c} is important for determining the colloid volume fractions that can form a low-density equilibrium gel. We additionally confirmed that the structures that were not phase separated formed percolated networks once *ε*_{b} was sufficiently large, typically around *βε*_{b} ≳ 16 for *η*_{c} = 0.15. There was a discrepancy of about 15% between the simulations and TPT in the minimum *ε*_{b} needed for phase separation to occur. We speculate that this is due to the presence of linker loops in the colloid network. The discrepancy is qualitatively consistent with the number of redundant bonds in the network near the onset of phase separation (Fig. 2), but it is not obvious how to renormalize the TPT to account for this or cycles of longer lengths. Ultimately, reliably predicting the minimum *ε*_{b} needed for phase separation is not as crucial for experimental realization of the gels as predicting the minimum *η*_{c} at which single phase gel assemblies may be realized. TPT can be reasonably used to qualitatively assess trends in the phase behavior as functions of linker length and mixture composition and to suggest experimental strategies to realize such colloidal gels at low density.

In our previous work,^{26} we showed how Γ controls the onset of percolation and also the range of colloid volume fraction over which phase separation is expected for colloids joined by spherical linkers. These mixtures exhibited re-entrant phase behavior as a function of Γ, as the spinodal initially broadened after percolation with a maximum width at the stoichiometric ratio of linkers with the 6-patch colloids (Γ = 3), before shrinking and eventually vanishing at higher Γ. In this work, we are primarily interested in the low-linker (Γ ≤ 3) regime, where the phase behavior is not as strongly influenced by the number of bonding sites on the colloids. Moreover, restricting the linker concentration may help circumvent linker solubility issues^{45} and avoids the introduction of additional interactions to the system (e.g., significant depletion attraction^{4,61} between colloids due to free linkers^{6}) that may affect the phase behavior. For polymeric linkers with *M* = 8, analogous to the spherical linkers, the spinodal broadens as Γ is increased from 0.9 to 3.0 (Fig. 4). However, the shift in the high density branch of the spinodal with respect to *η*_{c} is only a factor of 2.1. Moreover, Γ is additionally constrained by the need to form a percolated network and so can be decreased over only a limited range.

In contrast, the spinodal boundary is far more sensitive to the length of the linkers. Figure 4(b) shows that the high-density branch of the spinodal can be shifted by roughly a factor of 5.6 with respect to *η*_{c} as *M* is increased from *M* = 2 to *M* = 16. This trend is also born out in the simulations (Fig. 3). In patchy particle systems, an analogous compression of the spinodal is achieved by limiting the number of patches on the particles, but here we are able to keep the primary particles the same and simply adjust the linker length. These results demonstrate that the linker length can be used to control and expand the low-density window for producing equilibrium gels.

### B. Microstructure and linker mixtures

Although linker length provides effective control over the phase behavior of the system, it is expected to simultaneously affect the microstructure of the gel. Longer polymers occupy and exclude a larger volume than shorter polymers, which should increase the average distance between the colloids they connect; a similar effect has been observed in DNA-linked colloidal crystals.^{62,63} This increased interparticle spacing may have consequences for the macroscopic properties of the gel. For example, near-field coupling between plasmonic nanoparticles rapidly decays with increased surface separation.^{64} Hence, it may be necessary to simultaneously consider the phase behavior and microstructure of the gel when designing the linker.

Given that the average interparticle spacing for a given linker is likely intimately related to its length, we propose one strategy to partially decouple the phase behavior and microstructure by using mixtures of linkers having different lengths. The linker mixture is expected to lead to different phase behavior and microstructure than its pure components, and these should be sensitive to both the linker lengths and concentrations. To demonstrate this idea, we considered an equimolar mixture of *M* = 2 and *M* = 8 linkers. We computed the spinodal boundary for this system using TPT and also ran equivalent simulations [Fig. 5(a)]. The spinodal of the mixture lies in between that of the systems with only *M* = 2 or *M* = 8 linkers (Fig. 3), and the TPT prediction was again in reasonable agreement with the simulations. We further compared the phase behavior of this mixture to an effective single-linker system having the same number average molecular weight, *M* = 5 [Fig. 5(b)]. Interestingly, the spinodal predicted by TPT is nearly the same for the mixture and the effective linker, and the phase behavior measured in the simulations shows similar boundaries in *η*_{c}, with some discrepancies again arising in *ε*_{b}.

Despite this similarity in phase behavior, the simulations revealed that the linker mixture had a significantly different microstructure than the single linker with the same average length. We computed the radial distribution function *g*_{cc} between the colloids for the three single-linker systems studied (*M* = 2, 5, and 8) and also the mixture at different *η*_{c} (Fig. 6). As the linker length increased from *M* = 2 to *M* = 8, the first peak in *g*_{cc} was pushed outward to larger *r*, from near contact at *r* ≈ 5 *σ* for *M* = 2 to *r* ≈ 7 *σ* for *M* = 8. This increase in the colloid spacing with the linker length was expected because the preferred size of a linear chain increases with *M* (polymer radius *R* ∼ *M*^{0.588} in a good solvent), and there is a corresponding free energy penalty to confine the linkers between two closely spaced colloid surfaces.^{45}

The microstructure of the mixture was dissimilar to the equivalent single-linker (*M* = 5), instead showing signatures of particle spacings typical of both the *M* = 2 and *M* = 8 component linkers and a more uniform long-ranged structure on average. The split peaks close to contact, which also arise for the *M* = 2 linkers, are a key indicator of the underlying local heterogeneity of the microstructure in the mixture. These peaks are due to colloids in configurations rotated around the short *M* = 2 (dimer) bond, from contact to separation by the bond length, and are absent for the longer linkers.

We measured *g*_{cc} at both *η*_{c} = 0.10 [Fig. 6(a)] and *η*_{c} = 0.15 [Fig. 6(b)], obtaining the same qualitative trends for the first peaks in *g*_{cc} at both volume fractions. *g*_{cc} was largely unchanged by *η*_{c} for the *M* = 8 linkers, which is consistent with a sufficiently low-density fluid interacting through a potential of mean force determined by the linker. As the linker length decreased, the first peaks in *g*_{cc} became more pronounced at the lower *η*_{c}. This additional structuring is likely due to the approach to the spinodal boundary for the *M* = 5 linkers and the mixture (Fig. 5), where increased long-ranged structuring was also detected by *S*_{cc}(0). In fact, at *η*_{c} = 0.10, the *M* = 2 linker has crossed into the spinodal region [Fig. 3(a)].

These results demonstrate how the linker composition can be used to design the phase behavior and microstructure of colloidal gels at a target volume fraction. The *M* = 2 linkers promote a microstructure with close contacts between particles, but the colloid–linker mixture phase separates at *η*_{c} = 0.10. The *M* = 5 and *M* = 8 linkers can form homogeneous equilibrium gels at *η*_{c} = 0.10, but the first neighbors of the colloids are farther apart. By blending the *M* = 2 and *M* = 8 linkers, an equilibrium gel can be produced at *η*_{c} = 0.10 with a microstructure that shows signatures of both linkers. In fact, *g*_{cc} for the mixture was effectively given by a composition-weighted average of *g*_{cc} for the pure single linker fluids (dotted line in Fig. 6).

Using a mixed linker composition not only provides facile control over gelation through the linker concentration but also offers flexibility in extending the low-density window for gelation and tuning the microstructure of the resulting gel; this paradigm can be applied to other systems with a broad range of linking chemistries. Moreover, the linker gel clearly offers greater tunability compared to a single-component gel, including those based on conventional patchy particles, because linker properties like concentration, length, and flexibility can be controlled independently of the primary colloid.

### C. Macroscopic properties

Having established how the linker can be used to design the phase behavior and microstructure of the gel, we finally considered potential impacts of the linker on macroscopic properties that depend on the underlying gel structure. As a rough proxy for the impact of the linker length on the optical properties of a gel composed of plasmonically active nanoparticles, we computed the average fraction of colloids in a given configuration having a certain number of neighbors with surface separation less than 1 *σ* (Fig. 7). This separation corresponds to 20% of *σ*_{c}, which is the distance at which significant coupling between plasmonic nanoparticles may be expected to influence the macroscopically observable optical properties of such a gel.^{64} Hence, the optical spectra of plasmonic nanoparticle gels are expected to depend strongly on the average number of closely contacted neighbors.

Consistent with Fig. 6, the shorter linkers created microstructures with more colloids having at least one close contact (88% for *M* = 2), while the longer linkers had far fewer of these (57% for *M* = 8), as shown in Fig. 7(a). The mixture had more colloids with close contacts (77%) than its effective single-length linker (68%), and its coordination was intermediate between its component linkers. The differences in coordination between the mixture and its component linkers are visually apparent in representative snapshots from the simulations [Figs. 7(b)–7(d)]. This example clearly highlights an advantage of the linker mixtures compared to the single-length linker: the microstructure can be reasonably controlled while simultaneously maintaining the desired phase behavior. However, further study using chemically detailed simulations and a quantitative electromagnetic model will be required to fully understand how the linker can be used to modulate, e.g., the absorption spectrum of a gel composed of plasmonically active nanoparticles.

We also expected that the linker length might influence the gel’s mechanical properties, which have important implications for how the assembled gel responds to applied stresses during use. The linkers create effective bonds between the colloids that can store energy. The gel can then be thought of as an elastic network^{65} analogous to classic theories for polymers, with the primary colloids serving as effective cross-link junctions. In the polymer models, the shear modulus is dominated by entropy and so is roughly *k*_{B}*T* per network strand per volume.^{45} It was recently shown, however, that the degree of bonding controls the rheology of DNA nanostar hydrogels^{66} through both the elasticity of the effective bonds in the network and the constraints imposed on the junctions themselves. Moreover, real chain molecules have additional interactions like excluded volume that prevent chain crossing, which can cause the modulus to increase.^{67} Hence, we characterized the linear rheology for our linker model directly from the simulations.

To simulate long-lived network connections relevant to experiments (*βε*_{b} → ∞) without sacrificing numerical accuracy, we permanently bonded the end segments of the polymers to the colloids. (In practice, this was done in lammps by removing and replacing the overlapped colloid patches in the rigid body by the end segment of the polymer.) We additionally switched off the attractions to the remaining patches from any free polymer ends to completely freeze the network configuration. We then measured the equilibrium stress tensor ** σ** over trajectories generated from 21 statistically independent network configurations for each linker length at

*η*

_{c}= 0.15. We computed the shear-stress relaxation function from the off-diagonal (

*i*≠

*j*) components of the stress,

for *t* < 100 *τ*. Because of the permanent network bonds,^{45} *G* did not decay to zero at long times but rather reached an equilibrium modulus *G*_{e}. We determined *G*_{e} by averaging *G* over times *t* ≥ 25 *τ*.

To compare the stress relaxation of the different network configurations and linker lengths, we defined Δ*G*(*t*) = *G*(*t*) − *G*_{e} (Fig. 8). We averaged Δ*G* and *G*_{e} across the independent configurations of a given linker length. There were clear differences in Δ*G* for the different linker lengths, with *M* = 2 showing a faster initial decay in the stress followed by a stronger anticorrelation than the longer linkers. The differences between the *M* = 5 and *M* = 8 linker were less pronounced. Interestingly, the mixture of *M* = 2 and *M* = 8 linkers showed a stress relaxation intermediate between the systems with one of those two linkers. This relaxation was initially consistent with the arithmetic mean of the two, although small differences are apparent for longer times; this deviation may be partially due to statistical uncertainty in the simulations. These results suggest that linker mixture composition can then be used to tune the macroscopic effective stress relaxation in the gel, which is connected to the frequency-dependent moduli.

The equilibrium modulus *G*_{e} also systematically increased with the linker length (inset of Fig. 8), although there was significant variability due to the different network configurations. Flory’s elastic theory for cross-linked polymer networks^{68,69} predicts a modulus of *ξ*(*k*_{B}*T*/*V*), where *ξ* = *B* − *J* − *C* is the cycle rank of a network having *B* bonds, *J* junctions, and *C* connected components. If all linkers were incorporated into our network uniformly and without any defects, *ξ* = (Γ − 1)*N*_{c} − 1, giving a modulus of 0.001*k*_{B}*T*/*σ*^{3} at *η*_{c} = 0.15. Clearly, *G*_{e} was much larger than this for most of the points, especially as the linker length increased. We confirmed that *ξ* did not deviate significantly (<5%) from its idealized value for any of the networks. We accordingly speculate that the large differences in *G*_{e} are due to approximations made by mapping the simulated linker gel onto the elastic network model, namely, the neglect of excluded volume interactions between monomers, which may cause the polymer chains to deviate from the theoretically assumed Gaussian statistics and prevent chain crossing, and from the large colloids that form the network junctions, which have both excluded volume and rotational degrees of freedom. Further work is required to quantitatively understand how and why *G*_{e} and other rheological properties depend on *M* for polymer-linked colloidal gels.

## IV. CONCLUSIONS

We have demonstrated how the phase behavior and microstructure of an equilibrium linker gel assembled from a model colloid–polymer mixture can be systematically tuned through the linker concentration and length using a combination of thermodynamic perturbation theory and computer simulations. The spinodal boundary for phase separation was dramatically compressed toward low colloid volume fractions upon increasing the linker length, predicting the potential for realizing equilibrium gels at low density. However, these low density gels consisted of colloidal particles with an increased spacing in the microstructure, which may have macroscopic consequences on the gel’s optical properties or rheology. To offer flexibility in designing gels that reach the targeted properties of interest, we proposed using blends of linkers of different lengths to allow independent control over both the phase behavior, which was consistent with an effective single-length linker, and the microstructure, which was consistent with an average of the constituent linkers. Together, these results constitute a robust and facile strategy for fabricating equilibrium gels with tunable densities and microstructures by straightforward macroscopic choices that are experimentally accessible, namely, the length of polymeric linkers and their concentrations relative to the colloids undergoing assembly.

Our simulations revealed that a substantial fraction of linkers did not form effective bonds in the colloidal gel network, either due to bonding both ends to the same colloid (looping) or by forming multiple bridging linkages between the same two neighboring particles, i.e., a double bond. Both conditions represent a “waste” of linkers that are no longer available to contribute to percolation of the colloidal gel network. We expect that the fraction of linkers forming such bonds would increase with the number of binding patches on the colloidal particles. Since these wasted linkers may inhibit percolation of the gel network, this effect should be taken into account in the linker design and formulation of experimental systems. For example, it may be possible to optimize the linker itself or its interactions with the colloids to limit the number of such bonds that form.

In this article, we restricted ourselves to consider only the effect of linker length on the phase behavior and microstructure, but other linker properties like internal flexibility or valency^{40,42} may similarly impact the gel. Moreover, we also neglected solvent-mediated interactions,^{70} including effective interactions between the linkers and the colloids and also hydrodynamic interactions between the various components. We plan to investigate these effects with a particular focus on how the linker influences macroscopic properties like the gel rheology in future work.

## ACKNOWLEDGMENTS

We thank Zachary Sherman for helpful comments on this manuscript. 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 and also by 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.

### APPENDIX: EVALUATING THE BOND OVERLAP FUNCTION

Without loss of generality, we place the colloid at the origin and assume that the polymer segment lies a distance *r* along the *z*-axis. The polymer patch is effectively at the center of the segment from Eq. (3), and so rotation of the segment can be neglected from the average in Eq. (8). The colloid patch is offset a distance *r*^{*} from its center, which is a point on the surface of a sphere with coordinates (*r*^{*}, *θ*, *ϕ*), where 0 ≤ *θ* < 2*π* is the azimuthal angle and 0 ≤ *ϕ* < *π* is the polar angle. The vector pointing to the colloid patch from the polymer segment is Δ**r** = [*r*^{*}cos *θ* sin *ϕ*, *r*^{*} sin *θ* sin *ϕ*, *r*^{*} cos *ϕ* − *r*]. The average of the Mayer f-function is then

where the integrand depends on both *θ* and *ϕ* through Δ**r**; this integral is easily computed by quadrature. However, for efficiency in the spinodal and binodal calculations that require many evaluations of this expression, we first evaluated Eq. (A1) at numerous values of *βε*_{b} between 0 and 20 and then fit the entire functional form to a basic neural network model that produced virtually exact results over the fit range.