Modeling colloidal interactions that predict equilibrium and non-equilibrium states

Colloidal suspensions, gels, and glasses form the majority of biological fluids andmany pharmaceutical preparations, and are prevalent in personal-care, agricultural, and industrial-coating materials. Even the simplest colloidal system—hard spheres suspended in a solvent—can display a wide range of behavior from Newtonian liquid-like flow to viscoelasticity to hard solid-like response, imparted by a hierarchy of microstructural relaxation length and time scales. Competition between Brownian relaxation and imposed flow can produce pronounced shear thinning, while hydrodynamic forces can induce shear thickening,^1,2 phenomena central to the flow of foodstuffs, such as mayonnaise,^3–5 or industrial fluids and polymeric solutions.⁶ Increased concentration of suspended particles can amplify such non-Newtonian behavior, while changes in particle-scale interaction can dramatically alter hydrodynamic shear thickening⁷ or induce discontinuous shear thickening.^8,9 Changes in concentration and particle interactions can also trigger thermodynamic phase transformation, such as phase separation or crystallization.¹⁰ A key point is that phase transformation can be manipulated to produce desirable microstructures that in turn impart desirable bulk properties. Foodstuffs provide many examples, including milk, a colloidal suspension that comprises fat droplets suspended in water to form an oil-in-water emulsion. Processing techniques, such as heating, acidification, and enzyme treatment, control the size and interactions between the droplets, leading to the formation of cream yogurt, and imbue unique textures to cheese.^11–13 Industrial coatings are another example of how colloidal microstructure is engineered during phase transformation, e.g., to cure a solid state for toughness or develop color and electronic properties. For example, phases of colloidal laponite,^14,15 pigment,^16,17 or metal oxide particles^18,19 can be fine-tuned to create colloidal films that impart corrosion resistant glazes, color, and magnetic properties for storage. Techniques, such as attraction-driven self-assembly or evaporation-driven convective assembly that sculpt crystalline, glassy, or bi-continuous network structures, have been developed for manufacturing paint,^16,17 ceramics,^20,21 magnetic storage,^18,19 and biomedical devices.^22,23

Relating interparticle potentials to equilibrium phase envelopes in colloids is routine for a wide range of systems.^24–26 The resulting crystalline or liquid states have properties that are straightforward to connect to underlying microstructure and tunable interparticle potential (IPP). However, there are exceptions for which conventional methods for relating IPP to phase behavior break down, for example, colloidal gels and glasses, which can form by the so-called arrested phase transitions.^27–29 Even the “phase” lines for gels and glasses cannot be predicted using conventional equilibrium approaches. Such colloidal states hold enormous promise for the development of novel materials, but techniques for fine tuning structure and rheology of novel materials are typically empirical.

The goal of the present study is to develop methods and tools for connecting material processing conditions, such as temperature to microstructure, that will enable detailed control over material properties and will reveal processing strategies that produce novel materials. Our method aims to apply to complex fluids in general, but to develop fundamental control principles, we develop it based on the workhorse model systems of spherical colloids suspended in a quiescent Newtonian liquid. In this case, two thermodynamic variables can control phase behavior: interparticle potential and colloid concentration. Although equilibrium phase transitions of hard-sphere colloids are routinely predicted by molecular theories, we are interested in colloidal gels and glasses, which are out-of-equilibrium states. Gels and glasses freeze in microstructures and exhibit rich properties that cannot be predicted by existing theories. The demand for design control of gel- and glass-based materials has opened a new frontier in the non-equilibrium colloids-processing space of materials science and engineering.

Here, we develop a predictive toolset that incorporates the impact of kinetics on the interparticle potential that predict non-equilibrium phases. Specifically, we identify and interrogate a competition between three kinetic processes as a material is “quenched” from weak attractions to strong attractions: the rate of the quench, the rate of Brownian diffusion, and the rate of bond formation. The first step in interrogating this competition is to develop an analytical form of the interparticle potential for an experimental model system; we select a colloidal system with IPPs that are thermo-sensitive because these hold great promise for future application to quench-and-anneal strategies for sculpting microstructure. We will demonstrate the basic mechanistic principles that govern structural connections to IPP, and we will identify how broadly these principles apply to other colloidal systems.

The first goal is a method that gives an analytical form of the interparticle potential between colloids. This analytical relation can be obtained directly by interrogating forces between particle pairs or indirectly by measuring particle distribution and inferring the underlying interparticle potential (IPP) that sets their spatial distribution. Examples of direct-measurement methods probe forces required to separate colloidal surfaces such as surface force apparatus (SFA) and³⁰ surface force balance (SFB) techniques,^31,32 atomic force microscopy (AFM),^33,34 and laser tweezer microscopy.^35–37 However, these methods are intractable when colloids are nanometers in size or have slippery surfaces.^33,34 Alternatively, scattering techniques^38,39 and optical microscopy^40–42 yield measurements of spatial distribution of colloids, from which the IPP can be inferred using liquid-state theory, and provide robust data for even nanometer-size particles and droplets.

Here, we develop a toolkit for extracting the thermodynamic variables that define an analytical IPP. For the widest applicability to many types of colloidal systems, we build it based on using measurements of particle configuration, rather than direct force measurements. We will expand the conventional approach, which is illustrated in Fig. 1: light, neutrons, or x rays are scattered through a sample and analyzed to produce the static structure factor, S_expt(q) [Fig. 1(b)], which is the wave-space microstructure. In parallel, an initial form of the IPP is proposed [see Fig. 1(c-i), 1(c-ii), or 1(c-iii)], which is inserted into a liquid-state theory model and solved iteratively to produce a predicted static structure factor, S_theo(q). This S_theo(q) is compared to the experimental S_expt(q); any mismatch indicates that the IPP needs to be refined. The refined IPP is again inserted into the theory model, and the process is repeated until S_theo(q) converges to S_expt(q). This procedure is illustrated as a flow chart in Fig. 2. While this approach accurately predicts equilibrium phase lines, it fails to predict gel or glass lines. As will be discussed in Sec. III A, this failure arises from two requirements of the conventional approach: diluteness and equilibrium.⁴³ We remark that closure approximations that account for many-body effects⁴⁴ augment the dilute structural model, permitting the conventional method to accurately predict dense equilibrium liquid states.^45,46 This approach is a powerful tool for integrating simulations with experiments and yields remarkable universal behavior across a wide range of physico-chemical origins of the IPP because the IPP is approximated by the second virial coefficient B₂. The extended law of corresponding states (ELCS) explains this success.⁴⁵ However, this approach breaks down in non-equilibrium conditions, such as gelation, or when the IPP is long-ranged.

There are two problems with relying on the ELCS to justify the conventional approach: it does not work with long-range potentials, and guarantees nothing regarding non-equilibrium conditions—it is an equilibrium-based law. For the former, the ELCS fails to accurately predict phase lines for intermediate or long-ranged potentials^47–49 because the detailed variation of potential with colloidal separation cannot be captured by B₂. For the latter issue, the liquid-state integral-equation theory⁴⁴ at the heart of the conventional approach is valid only at equilibrium. A quintessential example is the gelation of octadecyl-coated silica particles initially studied by Grant and Russel⁵⁰ and Verduin and Dhont²⁷ and later by Eberle et al.,^28,29 where small-angle scattering was used to extract the static structure factor and thence refine a Baxter potential. Under dilute equilibrium conditions, the theory accurately predicted measured structure. However, this approach significantly over-predicted the attraction strength necessary to trigger gelation. These results show that the IPP that triggers equilibrium phase transitions changes qualitatively during non-equilibrium conditions and that the second virial coefficient does not contain sufficient information to capture the how an IPP changes with crowding, temperature, pH, or other conditions. We believe this gap underlies the failure of conventional methods in predicting non-equilibrium phase behavior.

Nonetheless, attempts have been made to apply the conventional approach for predicting gelation. Lu et al. concluded that any potential profile (Asakura-Oosawa, Lennard-Jones, square-well) is equally suitable to describe the structural measurements needed to deduce a value of the second virial coefficient, ultimately suggesting that gelation of hard colloids with short-ranged attractions can be predicted via equilibrium theory with little sensitivity to the details of the theoretical form of the attraction profile.⁵¹ An equally important outcome of their study was the assertion that gelation occurs as a result of interrupted phase separation, a groundbreaking advance in the understanding of routes to gelation that suggests that the gel line and binodal are the same. However, the gel line is not always coincident with the binodal, as reported by Verhaegh et al.;⁵² Buzzaccaro, Rusconi, and Piazza;⁵³ and Zaccarelli et al.⁵⁴ In fact, studies by Gibaud et al.⁵⁵ and Foffi et al.⁵⁶ suggest that arrest lines depend sensitively on the detailed interaction potential even when equilibrium phase boundaries can be robustly predicted using only $B2*$ without regard to the detailed potential. This sensitive dependence is consistent with the much earlier findings of mode coupling theory, which showed that colloidal glass transition lines for various potential profiles with the same value of $B2*$ do not collapse upon rescaling with $B2*$⁠.⁵⁷ Overall, these foundational works indicate that gelation is a non-equilibrium phase transition that results from arrested phase separation and suggest that the location of the gelation line relative to the location of the binodal in a colloidal phase diagram is sensitive to the details of interaction potential profile. (For a more detailed discussion of these works, see  Appendix A.) Moreover, the current consensus paradigm concerning arrested phase separated is that that the gel line and binodal are distinct from one another. However, theoretical and modeling methods fail to predict the location of the gel line—and are still anchored in equilibrium metrics. We believe that there are two roots of the failure: first, the assumption that the interaction potential is invariant with colloid concentration and, second, utilizing strictly equilibrium metrics (i.e., B₂) to construct the gel line.

In the present study, we propose a new method for determining the interaction potential that accurately predicts both the equilibrium phase lines and the arrest line for a model colloidal system: a nanoemulsion of oil-in-water droplets with adhering polymer chains for which phase instability is triggered via temperature changes that induce attractive forces between the non-deformable droplets. We select this system because it requires the use of a non-trivial analytical form of IPP, which demonstrates the applicability of our approach to a wide range of colloidal systems, and because our ultimate vision is to execute complex quench-and-anneal trajectories throughout the phase instability region to sculpt new materials, temperature is an easily accessible control variable to execute such trajectories. Finally, thermo-sensitivity permits gelation in a rheometer; rheology of a “pristine” gel avoids the issues of history required by shear rejuvenation.

The remainder of this paper is organized as follows: In Sec. II, we present the in vitro (experiments) and in silico (simulation) model systems and methods utilized in this study. In Sec. III A, we present the framework of our approach, outlining the iterative connection between experiments and simulations. The approach is demonstrated in Sec. III B by constructing an interparticle potential model of a thermoresponsive nanoemulsion system, in which gelation can be triggered from thermal quenches. We show via a comparison of gel lines constructed in vitro and in silico that a single IPP that accurately predicts both equilibrium and non-equilibrium behavior of a laboratory sample can be constructed. The study is concluded with a discussion of findings and prospects.

The experimental colloidal system is comprised of an oil-in-water nanoemulsion, in which polydimethylsiloxane (PDMS, viscosity = 5 cSt) nanodroplets [with an average radius of 50 nm and polydispersity index (PDI) of 0.28 determined by dynamic light scattering] stabilized by sodium dodecyl sulfate (SDS) surfactant are dispersed in an aqueous continuous phase containing telechelic oligomer, poly(ethylene glycol) diacrylate (PEGDA; M_w = 700 g mol⁻¹), and deionized water. In extensive previous works, it has been shown that this system exhibits temperature-responsive interdroplet colloidal attractions that produce colloidal gels at sufficient volume fraction of droplets, ϕ.^58–61 Figure 3 schematically illustrates the current consensus regarding the molecular mechanism of thermoresponsive colloidal attractions in this experimental model system. Specifically, hydrophobic effects associated with the acrylic ends of the PEGDA oligomers result in preferential bridging of the polymer between the nanoemulsion droplets at elevated laboratory temperatures, which gives rise to net attractions between droplets at moderate separations. In addition to the attractive interactions that arise from polymer bridging of the nanoemulsion droplets, short range repulsions from compression of bridging chains and longer range repulsive electrostatic interactions from the ionic SDS surfactant at the droplet interface are also present. These effects combine to produce a composite interaction—with a potential profile, which likely resembles those illustrated in the sketch in Fig. 3(c). Variation in the relative magnitudes and ranges of each individual contribution permits a wide range of potential profile shapes, as illustrated by the three curves in the figure. This complex combination of interactions were originally characterized in Helgeson et al.⁵⁸ through analysis of temperature-dependent small-angle neutron scattering (SANS) data collected on a dilute (ϕ = 0.01) sample. The SANS data were fit with a square-well potential model to determine how the strength of attractions, quantified by the square-well potential well-depth, changed with laboratory temperature. For more detailed SANS sample preparation and SANS data collection methods, see Helgeson et al.⁵⁸ 

Because the square-well model represents a gross oversimplification of the nanoemulsion interactions illustrated in Fig. 3(c), the SANS data from Helgeson et al.⁵⁸ were refit with a different interaction potential model, the Two-Yukawa (2Y) model. We note that several assumptions were made in the selection of the 2Y model to represent the coarse-grained description of the effective interdroplet interactions. The first assumption is that, regardless of the chosen shape of the interaction potential, this description of the system treats the continuous phase, including the polymer, as an effective medium, such that the sole contribution of the polymer is in determining the effective interaction potential between droplets. In reality, the polymer density need not be uniform throughout the continuous phase such that a full description would require modeling the polymer explicitly. However, in a recent work,⁶² it was found experimentally that the local polymer concentration is approximately insensitive to the droplet density, supporting the approximation of an effective medium to which the polymer implicitly contributes. The second assumption is that, although the quantitative magnitudes of the features in the interaction potential will depend on temperature, the qualitative shape (i.e., the number and relative order of local maxima and minima) does not. The final assumption is that the thermodynamics of the experimental system with the “true” interaction potential sketched in Fig. 3(c) can be faithfully approximated by the 2Y potential. Although the sketch in Fig. 3(c) is clearly better described by a potential with three or more contributions, the decision to limit the number of contributing terms to two is made in order to limit the number of adjustable parameters in the model interaction potential, thereby avoiding problems with uniqueness and convergence when optimizing the model parameters and also avoiding sharp features in the best-fit potentials that are challenging to replicate in numerical simulation. Under these assumptions, the 2Y model was used to describe the shape of the effective interdroplet interaction potential at different temperatures and was applied to refit the previously reported SANS measurements. Details for how the SANS data were refit with the 2Y potential can be found in Sec. II B. Nanoemulsion sample preparation, droplet size characterization, and rheological analysis for the nonequilibrium comparison testing are described in Sec. II A 2.

Polydimethylsiloxane (PDMS, viscosity = 5 cSt) is used for the dispersed droplet phase and PEGDA (M_w = 700 g mol⁻¹), sodium dodecyl sulfate (SDS) surfactant, and deionized (18.3 MΩ) water make up the continuous phase. Nanoemulsions of various dispersed phase volume fractions, ϕ = 0.2, 0.3, and 0.4, were made via homogenization or ultrasonication. The continuous phase contains a polymer volume fraction, P = 0.33, of PEGDA. The total SDS surfactant concentration was changed with ϕ in order to hold fixed the final “free” surfactant concentration in the bulk continuous phase at 20 mM for all ϕ. For ϕ = 0.2, 0.3, and 0.4, this corresponds to total SDS concentrations of 130, 180, and 240 mM SDS, respectively, which were determined according to the work of Pagenkopp and Mason.⁶³ In brief, for small ϕ, a lower total surfactant concentration is needed to obtain the final free surfactant concentration of 20 mM because less surfactant goes toward the stabilizing droplet interface. All chemicals mentioned earlier were purchased from Sigma-Aldrich and used without further purification.

The nanoemulsions were prepared according to previously published procedures.⁵⁸ A crude emulsion was prepared by mixing the continuous phase [PEGDA, SDS, de-ionized (DI) water, photoinitiator] with a stir plate set at 700 rpm followed by dropwise addition of the dispersed phase (PDMS). High pressure homogenization of emulsions with ϕ = 0.2 and 0.3 was performed using an Avestin Emulsiflex-C5 homogenizer operating at 15 kpsi. Between each pass, the samples were cooled in ice and 12–16 passes were needed until the desired droplet size distribution was achieved.

The high volume fraction emulsion with ϕ = 0.4 was processed using an ultrasonicator instead of the homogenizer because higher ϕ samples gel during processing in the homogenizer, which seizes the homogenization valve. The sample was stirred continuously during ultrasonication processing using a stir plate and cooled using an ice bath. The droplet sizes were periodically checked during processing until the desired droplet size was achieved. After synthesis, the nanoemulsions were stored at 8 °C to slow droplet growth. Subsequent characterization of the nanoemulsions shows that, under the appropriately controlled conditions reported, the two methods of nanoemulsification produce indistinguishable droplet size distributions.

Droplet sizes were measured upon ultrasonication and homogenization using dynamic light scattering (DLS). DLS was performed using a Brookhaven Instruments BI-200SM goniometer operating system equipped with a 500 mW dye-pumped solid state laser operating at a wavelength of 532 nm. Before measurements, 100 µl of mother nanoemulsion were diluted in 3 ml of DI water. Measurements were conducted at 20 °C at a scattering angle of 90°. The measured average diameter of the nanoemulsion droplets was 2a = 50 ± 3 nm with a polydispersity index (PDI) of 0.28 ± 0.02, where a is the average radius of the droplets.

Gelation of the nanoemulsions upon changes in temperature was monitored using small amplitude oscillatory shear (SAOS) rheological characterization performed on a TA Instruments ARG2 rheometer with a 60 mm, 2° upper-cone geometry and a Peltier bottom stage. Temperature ramp experiments were performed at various ramp rates spanning 0.05–15 °C/min, from an initial temperature of 25 °C and ramped until gelation was observed (i.e., a dominant elastic response, G′ > G″). SAOS experiments were carried out with an applied oscillation frequency of ω = 10 rad s⁻¹ and an oscillation strain of γ = 0.1%.

We model the experimental system in silico by performing a simulated continuous ramp in laboratory temperature resulting in corresponding changes in the interaction potential. The computational model system is a colloidal suspension of 740 000 neutrally buoyant, hard Brownian spheres of average radius a dispersed in an implicit Newtonian solvent (representing the effective medium including the polymer and surfactant in the continuous phase) with viscosity η and density ρ. The Reynolds number, ρUa/η, where U is a characteristic particle velocity, characterizes the relative strength of inertial forces relative to viscous forces. The small particle size a ensures that the Reynolds number is small (Re ≪ 1) and fluid motion is thus governed by the Stokes equations. The colloids are size-polydisperse with 28% variance that matches the experimental model system. For dilute equilibrium structural characterization to match the experimental conditions in which neutron scattering experiments were conducted, we choose particle fraction ϕ = 4πa³n/3 = 0.01. For quench simulations that result in colloidal gelation, ϕ = 0.20, 0.30, and 0.40 were selected to match parallel experiments. The simulation cell is periodically replicated to model an infinite medium.

We select one of the simplest non-trivial models of attractions and repulsions of arbitrary strength and range is the Two-Yukawa model. Two variables K and Z set the strength and range of attractions, respectively. When a system exhibits both attractive and repulsion interactions, each variable has an attractive and repulsive form: K₁, K₂, Z₁, and Z₂. This workhorse analytical theory then gives the IPP as a function of colloid–colloid separation distance r and the interaction variables: V(r, K, Z).

For the thermoreversible system we study here, the temperature alters the attraction strength and range V(r; T) [Fig. 4(b)] arising from the hard-sphere repulsion and polymer-induced attraction [Fig. 3(c)]. The interaction between particles i and j is thus given by

where k_B is Boltzmann’s constant, T is absolute temperature, r_ij is the center-to-center separation distance between particles i and j, and a_i is the size of particle i (Fig. 4). The interaction variables K₁(T), K₂(T), Z₁(T), and Z₂(T) depend on the laboratory temperature, T, as described in Helgeson et al.⁵⁸ The selection of the variables K₁, K₂, Z₁, and Z₂ provides extremely broad utility of the Two-Yukawa potential. It is straightforward to identify values of these variables based on the physico-chemical origins of the attractive and repulsive forces as done widely in the literature to generate equilibrium phase lines.^64–67 

One can select a variety of values of the four variables and obtain identical values of $B2*$ and thus equivalent equilibrium phase lines. The fact that the selection made for K₁, K₂, Z₁, and Z₂ produces distinct interaction potential shapes is of little importance for predicting equilibrium phase envelopes (for short-range attractions). In fact, these combinations can be grouped into values of the ratios K₁/K₂ and Z₁/Z₂ to give a family of interaction profiles as sketched in Fig. 4(b). The blue, red, and green curves illustrate the three types of Two-Yukawa potentials. This versatility of the Two-Yukawa potential allows it to model a wide variety of systems.

However, as with any other potential, the value of $B2*$ associated with the IPP is insufficient to predict non-equilibrium behavior, such as gelation and gel lines. As a result, the details of the IPP will predict distinct gel lines. In the case of the Two-Yukawa potential, one will start with the three possible IPP profiles sketched in Fig. 4(b), but then a new method is required to tie one or more of those IPPs to match the gel line found in experiments to those predicted by the IPP, which will typically require dynamic simulation. In the present study, we will determine values of the variables K₁(T), K₂(T), Z₁(T), and Z₂(T) using scattering data selected temperatures (T = 10, 25, 40 and 60 °C) and the iterative refinement process as described in Sec. II B 3. Values at intermediate temperatures were computed via linear interpolation.

The motion of colloidal particles in the Stokes flow regime is governed by three forces: Brownian F^B, hydrodynamic F^H, and interparticle F^P. The Langevin equation describes the time-evolution of particle motion

where m is the mass or moment of inertia tensor, U is the particle velocity, and t is time.

The Brownian force acting on any particle i arises from thermal fluctuations of the solvent and obeys Gaussian statistics,

where the overbar denotes a time-averaged quantity over time scales longer than the solvent time scale, I is the identity tensor, and δ(t) is the Dirac delta function.

For a freely draining system, the hydrodynamic force acting on particle i by the surrounding solvent is described by Stokes’ drag law

where U_i is the velocity of particle i and ⟨u⟩ is the bulk flow velocity of the suspension. The angle brackets ⟨·⟩ indicate an average over fluid plus particles.

Finally, the interparticle force $FiP$ modeled here is derivable from a potential V_ij(r; T) between two particles i and j and is computed as a summation of pairwise interparticle interactions as

where $r̂ij$ is the unit vector along the line of centers of particles i and j. It is this potential V_ij(r; T) that we seek to model in the present work.

We use the Large-scale Atomic/Molecular Massively Parallel Simulator (LAMMPS) molecular dynamics package⁶⁸ with the implementation of an implicit solvent to construct the computational model system and conduct dynamic simulations. The simulation cell is first populated with particles distributed randomly throughout the domain. The simulation cell is initialized with only Brownian forces and hard-sphere repulsion to equilibrate the particle distribution. Following this, the attractive potential between particles is turned on, and the simulation commences. This in silico system evolves over time as particles diffuse due to Brownian motion, where an embedded Langevin “thermostat” enforces Brownian statistics, and particles undergo deterministic displacements due to attractive forces between them dictated by the potential V_ij(r; T) [Eq. (5)] as set by the Two-Yukawa potential [Eq. (1)]. In LAMMPS, the motion of each particle is computed by integrating the Langevin equation [Eq. (2)] forward over time using the velocity Verlet algorithm.⁶⁹ 

The system is brought to and held at simulated ambient laboratory temperature, V(r; T = 25 °C)/k_BT; because this system is in the liquid phase, particles still diffuse freely and equilibrate. Next, quenches are performed by implementing changes in the parameters in the Two-Yukawa potentials, as follows: An in silico temperature ramp is induced by adjusting the parameters of the interparticle potential K₁(T), K₂(T), Z₁(T), and Z₂(T) over time, which alters the strength of the interparticle potential and mimics experimental temperature ramps. The in silico temperature ramp induces a self-assembly into a scaffold-like bicontinuous network of particles. Figure 5 shows rendered images of a simulation cell during the ramp with magnified snapshots that highlight the strand structure.

To detect whether the material undergoes gelation during the ramp, a small amplitude oscillatory shear (SAOS) is imposed, with displacement γ(t) = γ₀ sin(ωt); this oscillation continues throughout the temperature ramp and is analogous to the procedure performed in experiments. To probe the linear-response regime, we impose strain amplitude γ₀ = 0.01 and oscillation frequency ω = 1/(a²/D₀), which produce a small dimensionless displacement, $Pe/(ωa2/D0)=γ̇a2/D0$⁠, where Pe is the Péclet number, that give the linear response. The resultant shear stress ⟨σ(t)⟩ is monitored continuously throughout the quench and comprises contributions from both the solvent and from the particle phase, ⟨Σ⟩,

where x is the distance from particle i to particle j and F^P = −∇V(r; T) [cf. Eq. (5)]. The angle brackets indicate an average over the entire particle phase, which deforms and relaxes over measurable time scales producing both a viscous response, G′, and an elastic response, G″. The viscoelastic moduli are averaged over several oscillation periods as described in Landrum, Russel, and Zia.⁷⁰ 

Particle positions and motion are tracked throughout each simulation to interrogate structure and dynamics. In addition, we calculate the static structure factor, S(q), the wave-space measurement of structure obtained in experiments,

where x_j is the position of particle j and q is the wave vector normalized by the average particle size a. We radially average the static structure factor to compute S(q), where q = |q|. Additional details about this calculation in our simulations are described in Zia, Landrum, and Russel.⁷¹ The next step is development of the interparticle potential that will change with temperature to mimic the experimental system.

The naïve approach for representing particle interactions of the experimental system in our computational model is to simply extract the reduced second virial coefficient, $B2*$⁠, by inferring the potential from structural configurations measured via SANS.⁵⁸ According to the extended law of corresponding states (ELCS), any short-ranged analytical potential adequately predicts the equilibrium phase envelope.⁴⁵ The premise of the present work is that this model-agnostic predictive power does not extend to predicting the conditions of out-of-equilibrium arrest, but it still provides the starting point for the formulation of valid candidate potentials. We will, using our new method, down-select to that which simultaneously predicts both the equilibrium structure (under conditions it can be sampled) and the arrest line (determined by the state points producing a transition to the arrested state). Thus, the first step is to construct this family of candidate potentials. Here we summarize our process for constructing interparticle potentials from scattering measurements (see the work of Helgeson et al.⁵⁸ for SANS sample preparation and data collection methods). The scattering intensity I(q) obtained from SANS encodes both intraparticle morphology (shape) and interparticle configuration: the shape of individual particles is quantified in the form factor intensity, P(q), and the spatial distribution or configuration of particles in the static structure factor S(q). Here, it is assumed that the particles are sufficiently dilute so that contributions to intra-particle and inter-particle interference are separable, such that I(q) = P(q)S(q). The form factor intensity P(q) is computed as part of the intensity calculation using size-dependent dynamic correlation data obtained separately via dynamic light scattering, also under semi-dilute conditions.⁵⁸ The static structure factor S(q) is then obtained directly once the scattering intensity I(q) is measured and P(q) factored out. Here, the wave vector magnitude q is a scalar, indicating a radial average of both the scattering intensity and static structure factor.

We next produce interparticle potentials V(r) that are consistent with the static structure factor S(q) obtained via scattering measurements. The connection between V(r) and real-space structure g(r), which is the Fourier transform of S(q), is a liquid-state integral-equation,⁴⁴ namely, the Ornstein–Zernike (O–Z) equation. The O–Z equation posits that pair correlations g(r) are separable: they can be split into direct correlations between pairs, c(r), and indirect correlations, where the correlation between particle pairs is influenced by the background ensemble configuration as follows:

where r_ij is the center-to-center separation distance between particles i and j. Because Eq. (8) contains g(r) and c(r) that are unknown, a closure relation is needed to relate the configuration integral [last term in Eq. (8)] to these two unknown functions. For example, the Percus–Yevick approximation is often used as a closure for many-body correlations; it approximates all interactions using a Boltzmann-like distribution, with direct and indirect potentials,

where w(r) is the potential of mean force that accounts for the free energy contribution arising from the surrounding particles. Substituting Eq. (9) into Eq. (8) gives the well-known Percus–Yevick equation,

In theory, Eq. (10) provides a direct connection between V(r) and g(r), and therefore S(q), meaning, one should be able to measure S(q), compute its Fourier transform g(r), and then solve Eq. (10) for V(r). Yet, analytically solving the integral expression in Eq. (10) to obtain V(r) from g(r) [or S(q)] is only tractable for monodisperse and purely repulsive hard spheres (or a simple potential, e.g., the square-well potential). A standard workaround is to propose an “initial” potential of known analytical form (ideally one that has features that represent physical sources of repulsion and attraction in the experimental system, e.g., the Two-Yukawa potential), whence it is tractable to numerically compute g(r) from Eq. (10). When a numerical solution to Eq. (10) is not tractable for the selected potential, other closure approximations may be necessary to further approximate Eq. (10). However, the selection of closure relation—e.g., Percus–Yevick (PY) approximation or hypernetted-chain (HNC)—is just as important as the O–Z equation because each involves approximations that ultimately propagate to the detailed shape of the interparticle potential. For example, a commonly used truncation of the PY closure, the mean-spherical approximation (MSA), underestimates the attractive strength $B2*$ required to induce the formation of large clusters or induce phase separation. To avoid this pitfall, we utilize the full Percus–Yevick equation, for which we developed a new numerical solution algorithm, and we combine it with the O–Z equation to predict a structure. This structure is compared with that obtained in experiment; any difference between the two is used to refine the proposed V(r), and this iterative process is continued until Eq. (10) produces a structure S(q) that matches experimental results. The details of how this process is carried out in the present study for the Two-Yukawa potential are given in  Appendix B.

This process leads to the construction of a family of several interparticle potential profiles [constrained to have matching values of $B2*(T)$] via use of different initial guesses that, as we will show, produce indistinguishably accurate predictions of the experimentally observed S(q). In previous applications of the conventional approach, this degeneracy is ignored for sufficiently short-ranged potentials by invoking the ELCS. While prior studies have attempted to utilize S(q) obtain at non-dilute ϕ, these efforts faced significant challenges in the collection and analysis of scattering data due to multiple scattering and emergence of arrested states (see  Appendix C for a more detailed discussion on these attempts to resolve the degeneracy).^28,29

These shortcomings of the conventional approach motivate the need for our new approach to selecting model potentials among candidate shapes. We propose to down select a single potential from multiple degenerate candidates that faithfully reproduces both equilibrium and non-equilibrium states. We utilize this novel approach, which builds upon the conventional approach as described in the sections to follow to directly incorporate non-equilibrium property determination (in this case, time-dependent SAOS rheological response) into the down-selection process.

When colloidal interactions are long-ranged or undergoing non-equilibrium transitions, the ELCS does not apply, which is borne out in the fact that a given value of $B2*$ alone cannot predict gel lines, for example. We described in Sec. II B that in such conditions, one must select some analytical form of a potential that is responsive to the actual colloidal system and then begin the iterative process. This new process is illustrated schematically in Fig. 6. Here, we present the results of this approach, starting with selection of the Two-Yukawa potential.

In both the new and the conventional process, the first goal is to obtain a value of $B2*$ to constrain the initial guess for the variables that set the proposed interaction potential. This goal is achieved by the first two steps in Fig. 6. First, the experimental system is interrogated to measure scattering intensities [using small-angle light, neutron or x-ray scattering (SALS/SANS/SAXS), described in Sec. II] that are then used to compute the static structure factor S(q) at equilibrium, where q is the wave vector. In the present study, we used neutron scattering (cf. Sec. III B) for the nano-emulsion system under dilute or semi-dilute conditions to isolate the contributions of pair-level potential interactions (these measurements can be repeated for several volume fractions to determine if they depend on volume fraction even in dilute conditions). The interaction potential may be sensitive to, e.g., temperature or depletion or concentration, for example; this dependence must be encoded into the potential, which can be done by systematically altering the control variables to obtain a set of scattering data. In the present study, we collected scattering data for a range of values of laboratory temperature. Following the procedure, we propose a simple potential (here, a square-well potential), insert it into the Ornstein–Zernicke equation as described (here) to predict structure; this structure is compared to scattering data and then iteratively refines the V_SW(r) until the predicted structure matches the experiments. We then insert this V_SW(r) into the definition of the reduced second virial coefficient to obtain $B2*$⁠,

We now have the $B2*$ needed to constrain the selected potential [here, the Two-Yukawa potential, Eq. (1)], i.e., we used the square-well potential as a bridge to develop an initial guess of the more sophisticated potential. We selected the Two-Yukawa potential [Eq. (1)] for our nano-emulsion system because it readily models temperature-dependent attractions and repulsions. This approach is general; one can apply the same procedure for obtaining an initial guess for a Lennard-Jones or DLVO (Deraguin-Landau-Verwey-Overbeek) potential. As detailed in Sec. II B, a single value of $B2*$ provided three equally valid forms of the Two-Yukawa potential as sketched in Fig. 4(b); any of the three will accurately predict the binodal, but not gel lines.

The goal of the new method is to extract which of the equilibrium-equivalent potentials also predicts gelation. To make progress, we recognized that the degeneracy problem involves the kinetic behavior of the IPP, which is most observable during non-equilibrium transformations, such as gelation or vitrification. Dynamic simulations are free of assumptions of diluteness or equilibrium constraints; instead, they enforce conservation of momentum to capture colloidal-scale diffusive dynamics and detailed particle interactions directly. As a result, simulations directly model the competition between interparticle attractions and Brownian motion that underlies colloidal gelation or arrested phase separation. The three candidate Two-Yukawa potentials each exhibit distinct changes in attraction and repulsion during a quench.

Thus, the key aspect of our method is the use of dynamic simulations to implement the set of potentials that satisfy equilibrium conditions and also induce gelation and vitrification; we implemented these new steps in the method (steps 3–6 in Fig. 6) as follows. We constructed an algorithm for translating laboratory temperature to the model parameters K₁(T), K₂(T), Z₁(T), and Z₂(T), and the analytical expression used to calculate the interaction between all pairs encodes the impact of relative size, producing a computational model that represents the diffusive and interparticle colloidal-scale dynamics of the experimental system, and how it varies with laboratory temperature, size distribution, and volume fraction. This step of connecting the in silico model to the laboratory conditions is essential and generalizable to other model systems where one might instead translate pH or salt content, for example, to the in silico potential. One must then validate the model in dilute conditions to assure recovery of the equilibrium behavior. Our results for this step of the process are discussed in Sec. III B.

The next step is to subject the suspension—both experimental and computational—to conditions of high concentration and/or stronger attractions, limits in which the ELCS does not apply, and defines the matching condition. In the present system, we aim to trigger gelation over a range of volume fractions and identify which of the candidate potentials modeled in silico accurately predicts it. Our thermo-sensitive system will produce a gelation line unique to each potential. Each gel line is a set of state points (in $B2*$ and ϕ) that separates equilibrium fluids from gels. The candidate interparticle potential that produces the best agreement between the in silico gelation line and the experimentally constructed gelation line will be selected as the correct potential, i.e., the potential that passes the non-equilibrium test as well as the equilibrium test.

In this section, we demonstrate how the details of an interparticle potential exert a discernible influence on gelation and therefore demonstrate how the proposed algorithm enables development of an interparticle potential model that accurately predicts equilibrium and non-equilibrium phase transformation. The method we present is agnostic to the origin of colloidal scale interactions. We follow the steps in the process in Fig. 6.

We first performed step 1(a) in the process (Fig. 6), working with two qualitative features of a Two-Yukawa potential: attraction and repulsion. Next, we perform step 1(b): measure the static structure factor as a function of wave vector at semi-dilute equilibrium conditions in laboratory samples of our model system. To do so, we conducted Small Angle Neutron Scattering (SANS) experiments (described in Sec. II A) on samples of colloids volume fraction ϕ = 0.01, measuring neutron scattering intensities at laboratory temperatures T = 10, 25, 40, and 60 °C to capture the temperature dependence of the polymeric bridging interactions the underlie attractive and repulsive colloidal interactions.⁵⁸ The radially averaged static structure factor S(q) for each temperature is plotted in Fig. 7(a).

We then completed step 2 in the algorithm: we computed the temperature-dependent reduced second virial coefficient $B2*(T)$ for each temperature using Eq. (11) and plotted it in Fig. 7(b).

We now commence step 3, the first of the new algorithm steps: we use the second virial coefficient to constrain the four variables in the Two-Yukawa potential. This produces a family of Two-Yukawa potentials consistent with B₂: depending on the selection of the K₁, K₂, Z₁, and Z₂—or rather the ratios K₁/K₂ and Z₁/Z₂ as discussed above—the potential can have any of three distinct shapes as sketched in Fig. 4(b): (1) short-range repulsion plus longer short-range attraction (blue curve); (2) hard-sphere repulsion plus short-range attraction (red curve); or (3) short-range attraction and a longer short-ranged repulsion (green curve). (A similar process is applied for other types of systems, for example, gelation via depletion interactions, where one would systematically alter the colloid-to-depletant size ratio Δ instead of K and Z, and depletant concentration instead of temperature.) Starting with the SR-LSRA profile [blue curve, Fig. 4(b)], we produced an initial estimate for the potential V(r; T = 10 °C) for the range and strength variables K₁, K₂, Z₁, and Z₂ using a seed estimate based on the $B2*$ just calculated. We then iteratively refined this potential until convergence—i.e., the predicted S_theo(q) from Ornstein–Zernike equation matched S_expt(q) from experimental SANS data, as shown in Fig. 7(c). Next, to construct the potential for higher temperatures (T = 25, 40, and 60 °C), we used the potential at T = 10 °C as an initial guess for the iterative refinement process to ensure that the shape remains constant as attraction strength changes with temperatures. The resulting final potential at each temperature is shown for the SR-LSRA potential in Fig. 7(d) (see Sec. II B 3 and  Appendix B for additional details). We repeat step 3 for the other two potentials, the Hard-Sphere Repulsion plus Short-Range Attraction (HS-SRA) [red curves in Fig. 4(b)] and Hard-Sphere Repulsion plus Short-Range Attraction and Longer Short-Range Repulsion (HS-SRA-LSRR) [green curves in Fig. 4(b)]; the structural comparison is provided in Fig. S2 of the supplementary material, and the final converged potentials are plotted in Figs. 7(e) and 7(f). (In another type of system, e.g., depletant gelation, one would induce deeper and deeper quenches by systematically increasing depletant concentration.) This validation of the equilibrium condition recovers behavior predicted by the extended law of corresponding states and justifies our use of the $B2*$⁠.

For step 4, we next implement these potentials into our computational model and perform a detailed in silico validation at equilibrium by comparing simulation-generated structures to experimentally measured structure.

We first simulated the same semi-dilute conditions and temperatures of the SANS experiments and measured S_sim(q). We then compare the in silico structure to S_theo(q) and S_expt(q). The data are presented in Fig. 8 and show excellent agreement between all three, validating our computational model. Furthermore, these results confirm that the three potentials constructed from equilibrium methods accurately predict equilibrium behavior. Finally, the results in Fig. 8 reaffirm the ambiguity of equilibrium methods: three distinct shapes of the interparticle potential produce identical equilibrium structure over a range of laboratory temperatures. Distinguishing the three potential shapes—i.e., deciding which accurately predicts both equilibrium behavior and gelation—requires one of two approaches. One could increase the volume fraction of the sample but, as discussed in Sec. I, multiple scattering undermines the accuracy of structural characterization. The second approach, taken here, is to connect each potential to its unique non-equilibrium state.

We now prepare to commence step 5 of the method. In our thermo-reversible system, we trigger non-equilibrium conditions by subjecting the laboratory sample and computational sample to temperature ramps. (Similar ramps should be employed in other systems, e.g., adding depletant or salt, for example.) Kinetics are expected to play an important role in how the system responds to increased attractions, for example, the rate at which the temperature is increased may change the temperature at which gelation commences, thus altering the gel line. We emphasize that since kinetics will play a role, one must first establish and validate a method for interrogating the rheology and structure during the quench-ramping procedure.

To mimic the laboratory temperature ramp, we developed an in silico algorithm that ramps up the interactive potential following the same ramp rate. We implemented the ramp [Fig. 9(a)] and interrogated the material throughout the temperature ramp. To extract the elastic and viscous moduli G′(ω) and G″(ω), we imposed small-amplitude oscillatory shear (SAOS) with strain amplitude γ₀ = 0.01 and fixed frequency ω = 1/(a²/D₀). This moderate frequency was selected because it will give viscous response in a suspension but elastic response in a gel. Simultaneously, we monitored the static structure factor S(qa) over a range of wave vectors qa at key instants in time during the quench. Finally, we monitored the contact number N_c, which is the number of neighboring particles within center-to-center distance 2 ≤ r/a ≤ 2.2. We evaluated each of these metrics to identify the gel temperature. We compare the simulation results to experiments at the same initial temperature T = 25 °C and colloid volume fraction ϕ = 0.30 at a ramp rate of 0.05 °C/(a²/D₀). Here, a is the average colloid radius in the computational model, and D₀ is the diffusivity of a particle of radius a alone in solvent.

Results are plotted in Fig. 9 for the HS-SRA-LSRR potential (Figs. S3 and S4 for other potentials in the supplementary material). Simulation images of the structure are shown for five instants in time (and instantaneous temperatures) during the ramp [Fig. 9(b)]. The images show the formation of particle strands and solvent-filled pores that coarsen during the ramp, both hallmarks of gelation and aging.^59,71 The storage and loss moduli are plotted in Fig. 9(c). We find that the sample is predominantly viscous (G″ > G′) at the beginning of the temperature ramp, matching expectations for the dispersed structure in the image in Fig. 9(b). As temperature ramps up, the material becomes dominantly elastic, commencing with a crossover of the linear viscoelastic moduli G′ and G″ at Δt = 420a²/D (T = 46 °C). This emergence of elasticity matches the apparent emergence of bi-continuous structure in the corresponding simulation image. As the final temperature (65 °C) is approached, both moduli approach a plateau; structural coarsening [Fig. 9(b)] also appears to slow down.

These rheological signatures of gelation were confirmed by structural measurements in Figs. 9(d) and 9(e). The static structure factor S(qa) is plotted in Fig. 9(d) for x ≤ qa ≤ y at several times during the quench. The structure curve is nearly flat at early times (corresponding to 38 °C) except a weak peak at large qa, consistent with of a fluid of transient clusters. A small-qa peak first emerges at Δt = 420a²/D (corresponding to 46 °C) at the cross-over temperature from the modulus plot. The location of the peak corresponds to long-range structure—a hallmark of gel morphology—that emerges at precisely the same time as the crossover of G′ and G″ and when the images reveal an apparent bi-continuous network. As the temperature ramps up over time and particles become more attractive, the low-qa peak shifts to lower-qa and grows, in agreement with the structural images, signaling gelation with concomitant phase separation.^59,71

We next investigate the particle-scale structure as it evolves during the quench,⁷¹ as shown in Fig. 9(e), a plot of the distribution of contact numbers P(N_c) at the same instants in time during the quench at which rheology and mesoscale structure were presented in Figs. 9(c) and 9(d). At the beginning of the ramp (Δt = 0), a peak in P(N_c) at N_c = 2 indicates very little bonded structure. A rightward shift and broadening of the peak during the quench shows that clusters have formed by 270a²/D and that even larger structures subsequently form. While the rightward shift indicates a growing volume-to-surface-area ratio, this shift slows down, indicating the arrest of structural evolution.⁷¹ We remark that the presence of N_c > 12 arises from the 28% polydispersity of the sample. These results showing emergence of arrested phase separation at approximately T = 46 °C are consistent with our prior experiments elsewhere.^58,59 This completes validation of the computational model.

We now commence step 5 in the method. Our validated non-equilibrium test is gelation. Cognizant that kinetics are an important part of the picture, we repeated the temperature ramp just discussed at several temperature ramp-rates [Fig. 10]. As expected, changing the quench rate changes the crossover time [Fig. 10]. In addition to different potentials producing different gel lines, we have different ramp rates of a given potential producing different gel lines. Ideally, we would like to find a gel line that is independent of the ramp rate. We approached this condition by quenching the material at a systematically slower ramp rate.

Starting with simulations at a colloid volume fraction ϕ = 0.30 and initial temperature T = 30 °C, we subjected the system to the four ramp rates shown in the bottom row of Fig. 10, looking for a crossover temperature as indicated by G′ > G″ using the same SAOS protocols. G′ and G″ are plotted as a function of ramp time in Figs. 10(a) and 10(b) for simulations and experiments, respectively. The temperature ramp was continued, but the crossover is marked with a vertical dashed line from the modulus to the temperature ramp, denoting where gelation occurred in the ramp. In simulations, this process was repeated for all three IPP profiles; Fig. 10 shows the data for the HS-SRA-LSRR potential, and Fig. S5 of the supplementary material shows the results from other potentials. In both experiments and simulations, slower quenches lead to slower gelation, but the crossover occurs at lower temperatures (shallower B₂ quench) as illustrated by the heavy curve in the T vs Δt plots. That is, we find that lowering the quench rate permits gelation at lower temperatures (shallower B₂ quenches). This trend can be understood by recalling the microscopic forces that contribute to particle dynamics during phase transitions: attractive forces and crowding, which hinder particle rearrangements, and diffusion, which promotes particle migration. Attractive forces and crowding promote condensation into particle-rich structures, but diffusion is required for such rearrangement to occur. However, attractions hinder this diffusive rearrangement, and this underlies the arrest of phase separation.⁷¹ A slow quench gives the system more time in the phase separation process where particle mobility is strong, allowing the formation and coarsening of the structure into large, dense domains. Steric crowding is then able to contribute and combined, and these produce an elastic network. In contrast, a fast quench plunges the system quickly to a regime where particles bond together with little rearrangement, allowing little coarsening of the newly formed bi-continuous structure. The resulting tenuous network must rely almost entirely on the bond strength to acquire elastic strength. Overall, systems subjected to fast ramps form thin strands of particles due to rapid arrest of phase separation and thus require stronger attractions (higher temperature, larger $B2*$⁠) to produce an elastically dominant mechanical response. In contrast, systems evolving during a slow ramp have time to undergo some phase separation before bonds become strong enough to arrest the phase separation.

Getting back to our goal of finding a ramp-rate independent gel line: connecting the gel point at each ramp rate (heavy curve in the T vs Δt plots) reveals a plateau for slower and slower quenches, as indicated by the dashed projection line. This projected plateau suggests an asymptotic limit—a quench above which (temperature below which) gelation will never occur; that is, there exists a minimum temperature of gelation. We suggest that this minimum gelation temperature sets the fundamental emergence of gelation in the system, that there is a “fundamental” gel line that is a unique identifier of an interparticle potential that robustly predicts gelation, and that this fundamental gel point is that quench temperature below which (or quench depth B₂ above which) gelation would never occur, i.e., that phase separation will not arrest.

We repeated step 5 in silico for several volume fractions, in pursuit of a gel line, for each of the three candidate potentials. In parallel, the same temperature quenches were applied in experiments to reveal the in vitro gelation behavior. Based on the findings discussed in the previous paragraph (Fig. 10), each quench was repeated for several quench rates. The quenching and interrogation protocols remained identical to that discussed earlier. The resulting data are plotted for ϕ = 0.20, 0.30, and 0.40 in Figs. 11(a), Figs. 11(b), and 11(c), respectively. There are two key features of the three plots. First, gelation temperature varies with the quench rate; we will need a rationale for picking of the “true” gel point from the data points. Second, the curves within each plot are separated from one another, showing that the variation with the quench rate differs from one potential to another, confirming the suspicion that despite their common value of $B2*$⁠, the three potentials predict different gelation behavior.

For the first point, the quench-rate dependence of the gel point: there are at least two distinct regimes. For fast quenches, variation is approximately linear, in contrast to gel temperature that is very sensitive to changes in the quench rate at slower quenches. The sensitivity to quench rate is most pronounced when the attraction range is short and the volume fraction is low [Fig. 11(a)], and, conversely, gelation temperature is less sensitive to the quench rate when attractions are longer-ranged and volume fraction is high [Fig. 11(c)]. Mechanistically, reduced particle mobility reflects an “effective temperature” that, when decreased, induces gelation because bond formation and crowding arrest dynamical rearrangements required for phase separation to complete. Longer-range attractions arrest dynamics at lower volume fractions, where higher crowding sterically reduces mobility that can arrest dynamics with shorter attractions. Consequently, there is a crossover in the curves at the transition from the fast- to the slow-quench regime, most prominent in Fig. 11(a). However, as crowding grows [Figs. 11(b) and 11(c)], this effect becomes less pronounced because particles are brought closer together by high volume fraction. The slow-quench regime is the most important because there the gelation temperature is highly sensitive to changes in the quench rate.

The slowest quenches provide data that extrapolate (quadratically) to the infinitely slow quench rate axis, which suggests that there is a limit below which gelation will not occur—a fundamental gelation temperature independent of the quench rate. We propose that this minimum gelation temperature sets the fundamental non-equilibrium behavior of the system, and that this regime corresponds to a quasi-equilibrium quench that, in principle, should allow phase separation. Furthermore, each potential reaches the vertical axis at a different point, indicating that each potential has a unique gelation temperature. Now that we have identified a unique gelation temperature for each candidate potential at a range of volume fractions, we compare the experimental and simulation results to find the best match among the models that best describe the experimental system.

We now complete step 6: Plot the phase diagram of the computational model and the experimental system and compare the gel lines. For the in silico system, we extract the fundamental gelation temperature at each volume fraction studied and plot it as a function of volume fraction in Fig. 11(d). Below this gel line (the locus of fundamental gelation temperatures at various volume fractions), the system will equilibrate and never arrest, regardless of the quench rate. We sketch in an extension of the arrest line and an equilibrium phase envelope for each potential (the equilibrium fluid–fluid phase boundaries are presented in our companion article⁶²). The extension of arrest line delineates the expected kinetically arrested glassy behavior outside of the equilibrium phase envelope.^72,73 The sketched phase envelopes are shown as shaded regions to emphasize that these are estimates. Dashed lines suggest low-volume fraction percolation behavior and high-volume fraction vitrification lines (also explored in our companion article⁶²).

Each in silico potential gives a distinct gel line. Overall, a shorter-ranged potential (HS-SRA-LSRR) produces gelation that is more sensitive to changes in volume fraction than a longer-ranged potential (as evidenced by the more pronounced slope of the gel line). The shorter the attraction range, the higher the concentration must be in order to gather particles together into a fully connected network (hence the steeper slope at lower volume fraction). In contrast, longer bonds can span a greater space between particles, such that the bonded network can be less dense but fully connected (thus the flatter gel line across a wider range of volume fraction). Comparison of the three in silico gels with the laboratory sample [black in Fig. 11(d)] reveals that the HS-SRA-LRR potential is the best predictor of gelation. Indeed, the in silico predictions for the gel line for this potential [green filled triangles in Fig. 11(d)] is statistically nearly indistinguishable from that obtained in the experiments, without any additional adjustment in the model parameters.

Steps 1–6 identify the best qualitative match between the predicted and the experimental gel line. The selected potential can be further refined (step 7) by adjusting the range of attraction and repulsion (in our example, the variables K and Z) and repeating steps 5 and 6 to obtain exact quantitative match. If steps 1–6 produce predictions with poor qualitative match between simulation and experiment away from equilibrium, the recourse is to examine whether the originally selected potential in step 1(a) was the best choice (in our study, the Two-Yukawa potential). In the present case, because predictions from the selected potential [green filled triangles in Fig. 11(d)] lie nearly within the associated uncertainty in the gel point extrapolated from experiments (black filled diamonds), steps 1–6 are sufficient to demonstrate the success of the method.

In this study, we devised and validated an algorithm that extracts the detailed interparticle potential from an experimental system and replicates it in a computational model that successfully recovers the equilibrium phase envelope as well as the gel line. This result is novel because prior approaches to modeling interparticle potential could only predict the family of potentials that give the correct equilibrium phase envelope but could not predict the gel line. We showed that the potential must encode non-equilibrium phenomena that embody the competition between quenching kinetics and particle dynamics to faithfully match the experimental system.

Our framework can be applied to any system in which arrested states can be formed via temporal changes in interparticle interactions or colloid density, whether such changes are induced by temperature, pH, or concentrations of other components. The requirement is the capability to perform such temporal changes and detect the emergence of arrested states.

The results of the non-equilibrium test provide a deeper understanding of the arrest line in the gelling system. The locus of state points defining the fundamental transition to the gelled state is set by the interplay of quench kinetics with the steric and attractive hindrance (ϕ and 1/T) that slow particle dynamics sufficiently to interrupt phase separation. Although the second virial coefficient is essential to constructing the gel line, it predicts multiple possibilities. Moreover, even if the correct interparticle potential was known in analytical form, the rate of the quench also predicts multiple possible locations of the gel line. The latter behavior is somewhat similar to vitrification in molecular liquids, where a “fictive temperature” is often used to communize experiments but does not resolve the issue of kinetics-induced ambiguity. Our algorithm extracted mechanistic differences between the multiple potentials; in the case tested, changes in the gel line with variation in the ratios K₁/K₂ and Z₁/Z₂ pointed to the relative influence of steric and attractive hindrance on gelation temperature. From this, we found that the quench can further tip this balance.

We deduced that this path dependence can be removed by an infinitely slow quench that will predict the deepest quench possible that does not subsequently arrest. Using this part of our algorithm, we identified the time-independent limit. From this comparison, we identified a single potential that captures both the equilibrium and non-equilibrium behavior of experiments. We claim that this process reveals the fundamental arrest line for a given potential because, “above” this arrest line, phase separation always proceeds unobstructed.

Overall, the new algorithm provides a platform for constructing models that can be rapidly tuned to faithfully represent an experimental colloidal system undergoing equilibrium and non-equilibrium phase transformations. We expect that this framework will also be useful for predicting percolation, vitrification, and other self-assembled states. The nearly hard sphere parameters selected in the present work do not place a restriction on the regime of validity of our method; conversely, our method will be particularly useful for handling soft repulsions and deformability. Size polydispersity is also an important factor to be considered because many real colloidal systems present great variation from one type to the next or even one batch to the next in technologically relevant materials. In the present work, both the in situ and in silico particles had 28% size polydispersity. There is no limitation on size polydispersity inherent in the model or the method. Alterations in polydispersity do shift the average diffusivity of particles in the system, which we expect will produce quantitative shifts of the fundamental arrest line, with the volume fraction adjusted accordingly. These effects will emerge naturally from the model. Our method can be applied to colloidal systems where interactions are patchy,^74,75 long-ranged,⁷⁶ or between irregularly shaped particles—including globular proteins, depletion systems, electrostatically interacting systems, and would be especially powerful for industrial suspensions where little is known about particle interactions. For systems with extremely complex interactions, our method can be applied with a more expressive potential, such as a Three-Yukawa potential. Clearly, the use of a more complex interaction potential model will increase difficulty of optimizing the potential. If such added complexity is warranted, one might elect to construct a potential de novo using molecular-dynamics simulations.

In silico quench studies in this work were conducted in the freely draining limit. Hydrodynamic coupling between colloidal particles is expected to reduce particle mobility, with an effective influence that is O(k_BT), which is weaker than the influence of attractive forces of order a few to several k_BT during a quench. Although hydrodynamic interactions do not influence equilibrium properties, several viewpoints on the role of hydrodynamic interactions on the location of the non-equilibrium arrest line have been suggested by prior literature.^77–82 These effects produce quantitative shifts in linear-response rheology and structure⁸³ that are not fundamentally different from changes in attraction strength and range. We expect the impact of hydrodynamics to emerge directly in the asymptotic limit of the gel temperature as the quench rate is decreased, which is where the quantitative position of the arrest line is set.

The seven-step method can be carried out just one time in order to construct a model that can be used over and over for a given laboratory system. We described the computational process in detail, which was carried out very efficiently using the highly optimized and parallelized LAMMPS algorithm.⁶⁸ Now that the procedure is established, routine computational activities are sufficient to implement it. Implementations of the process for a totally new potential can be completed in less than a week, from start to finish: laboratory experiments, model set up, validation, and selection of the potential. Once completed, the computational model can be used for as many experiments or projects as desired with no further algorithmic intervention.

Future studies should widen the predictive power to percolation and vitrification. Within the phase-separated region, the new understanding of how quench rate competes with diffusion and attractions can be leveraged to reveal a “goldilocks zone” where structure is forming but malleable, ideal for intervention to sculpt the forming microstructure. This approach would be very useful in designing new soft materials with novel structural, mechanical, and optical properties.

See the supplementary material for the complete predicted scattering intensity profiles I_theo(q), figures demonstrating influence of square-well potential depth on extracted $B2*$⁠, and time evolution of structural and rheological features.

This research was supported by NSF (Grant No. 1729017) from the program for Designing Materials to Revolutionize and Engineer our Future (DMREF). B.K.R. acknowledges the support of an NSF Graduate Research Fellowship (Grant No. 1656518). R.N.Z. acknowledges computational support from the National Science Foundation Extreme Science and Engineering Discovery Environment (XSEDE) Research Award (Grant No. OCI-1053575) and its resources at the Texas Advanced Computing Center (TACC) at the University of Texas at Austin. We acknowledge the support of the National Institute of Standards and Technology, U.S. Department of Commerce, in providing the neutron research facilities used in this work. This work utilized facilities supported by the Center for High Resolution Neutron Scattering, a partnership between the National Institute of Standards and Technology and the National Science Foundation under Agreement No. DMR-2010792. We thank Dr. Yun Liu (NIST) for helpful discussions.

The authors have no conflicts to disclose.

All datasets underlying the conclusions of the paper are either published directly in the figures or are available to readers upon request to the corresponding author. Raw data were generated at the TACC Stampede large scale facility, as well as NIST. Derived data supporting the findings of this study are available from the corresponding author upon reasonable request.

In Sec. I, we discussed examples that aimed to apply the conventional approach for equilibrium gelation. These studies clearly show that the gel line is not coincident with the binodal. For example, in a study of globular protein solutions, Gibaud et al. interrogated three ideas: first, globular proteins behave as colloids in the context of phase separation and gelation; second, phase separation is a route to gelation;⁵⁵ and third, the binodal and gel line are coincident. The authors established that the proteins exerted short-ranged attractions and verified that the ELCS held for the equilibrium phase envelopes; they found that the proteins did indeed exhibit colloidal phase transitions. In addition, they found that gelation followed the onset of phase separation, but the resulting gel line did not overlay the binodal. Thus, Gibaud et al. provide strong evidence that suggests that gelation is a non-equilibrium phase transition by connecting it to the glass line outside the binodal. That study is particularly convincing owing to the fact that the gel line and binodal were both generated by direct measurements in experiments, although without regard to absolute values of the potential. A hint for the mechanistic origin of this non-equilibrium behavior can be gleaned from the fact that, while changes in the induced attractive force between proteins could be scaled out for the equilibrium phase transition, that approach failed for the gel line: the reduced second virial coefficient could collapse all equilibrium phase separation onto a single binodal, but the resulting gel lines were distinct and could not be collapsed onto a single gel line using the reduced second virial coefficient $B2*$ alone.⁵⁵ In fact, earlier work provides further insight into a possible mechanistic origin of the separation between binodal and arrest line: Foffi et al.⁵⁶ found that the arrest line is set by the intersection point of the glass line (outside the phase envelope) and the binodal; even earlier simulation studies⁸⁴ demonstrate that this intersection point depends sensitively on the details of the attraction potential. From this, one can conclude that arrest lines depend sensitively on the detailed interaction potential even when equilibrium phase boundaries can be robustly predicted using only $B2*$ without regard to the detailed potential. This sensitive dependence is consistent with the much earlier findings of mode coupling theory, which showed that colloidal glass transition lines for various potential profiles with the same value of $B2*$ do not collapse upon rescaling with $B2*$⁠.⁵⁷ Overall, these foundational works indicate that gelation is a non-equilibrium phase transition that results from arrested phase separation and suggest that the location of the gelation line relative to the location of the binodal in a colloidal phase diagram is sensitive to the details of interaction potential profile. Moreover, the current consensus paradigm concerning arrested phase separated is that that the gel line and binodal are distinct from one another. However, theoretical and modeling methods fail to predict the location of the gel line—and are still anchored in equilibrium metrics.

We believe that there are two roots of the failure: first, the assumption that the interaction potential is invariant with colloid concentration and, second, utilizing strictly equilibrium metrics (i.e., B₂) to construct the gel line. This can be explained as follows. The study by Lu et al.⁵¹ determined the interaction potential from a system with a given depletant concentration via an iterative refinement procedure: the structure required to compute the potential is first extracted experimentally by measuring cluster-size distribution at equilibrium using microscopy, restricted to 16% maximum volume fraction. Separately, a set of interparticle potentials V(r) was implemented in simulation to induce cluster formation in a colloidal suspension, also at low volume fraction. Each simulation potential was iteratively refined until the in silico cluster size distribution matched the experimentally measured cluster size distribution. Fortunately or unfortunately, all of the candidate potentials (Asakura-Oosawa, Lennard-Jones, and square-well) produced a matching structure. Moreover, the value of B₂ calculated directly from the in silico V(r) was identical at fixed depletant concentration. Thus, the authors were free to select the simplest potential. However, the fact that the calculations were restricted to equilibrium cluster formation at low volume fraction all but guarantees that the phase transition predicted will be the binodal rather than the gel line. The problem with this approach is the assumption that the potential V(r) is independent of volume fraction—that the potential predicted at lower volume fraction is sufficient to describe interactions at high concentration—a dubious supposition if one assumes that gelation follows phase separation into a dense phase. Validation of their approach should have included implementing the potential in silico to generate gelation and, from it, construct a gel line. We expect this would have yielded not only one gel line distinct from the binodal, but three distinct gel lines, and a conclusion that the interparticle potential changes with volume fraction and that the non-equilibrium nature of arrested phase separation requires non-equilibrium tests to identify the unique potential that accurately predicts both the binodal and the arrest line for a given colloidal system.

In the case of the Two-Yukawa potential, there are four independent parameters that are intractable to fix by iteration without an initial estimate for each, i.e., a means to constrain them. To address this added challenge, we devised an iterative algorithm for constructing an interparticle potential that consists of two steps: first, we extract a $B2*$ of interparticle interactions using a square-well potential; second, we utilize this $B2*$ as a constraint to construct a Two-Yukawa potential.

In the first step of our algorithm, we estimate the parameters of a square-well potential V_SW(r; T) that models the thermoresponsive bridging interactions. This V_SW(r; T) is inserted into Eq. (10) to output the predicted structure and its Fourier transform is computed to obtain at theoretical prediction of the structure factor, S_theo(q). The predicted structure S_theo(q) is compared to the experimentally obtained static structure factor S_expt(q) as follows. If |S_theo(q) − S_expt(q)| is greater than a prescribed threshold, we refine V_SW(r; T). If the error is smaller than the threshold, we “accept” the square-well parameters and compute the $B2*$ via

where a is the average size of a particle. This $B2*$ is utilized in the next step as an additional constraint for a Two-Yukawa potential. We note that the choice of the square-well model—made by convenience due to its minimal number of parameters and featureless attraction shape—is arbitrary, and we found in initial testing of the algorithm that this choice, including the square-well width, has negligible influence on the extracted values of $B2*$ for sufficiently dilute values of ϕ, consistent with the dilute, small-angle approximation of the virial expansion S(q → 0, ϕ → 0) ≈ 1 + B₂(T)ϕ + O(ϕ²) (see Fig. S1 of the supplementary material for the negligible influence of square-well width on $B2*$⁠).

In the second step of the algorithm, we estimate the parameters of a Two-Yukawa model V(r; T) that is consistent with the value of $B2*(T)$ from the square-well model. This estimated V(r; T) is inserted into the O–Z equation to predict S_theo(q), which is again compared to S_expt(q). When comparing the two structure factors, we account for instrument resolution smearing effects and variable uncertainty in S_expt(q) by “smearing” S_theo(q) and appropriately weighing the difference between S_theo(q) and S_expt(q) (see  Appendix E for details of this process). We find that if these effects are not accounted for, this comparison and refinement process may result in erroneous values of $B2*$⁠. The potential V(r; T) is again refined until suitable agreement is achieved between S_theo(q) and S_expt(q). For construction of a family of Two-Yukawa potentials, we select several different initial estimates of K₁(T), K₂(T), Z₁(T), and Z₂(T) and repeat this part to obtain Two-Yukawa potentials with different detailed profiles for different values of T. The two-step process above is necessary because the iterative refinement process becomes under-constrained for complex models, such as the Two-Yukawa potential, in which a multiplicity of potential shapes with its own $B2*$ are permissible from a set of potential parameters (K₁, K₂, Z₁, and Z₂).

In the present work, the attractive interactions between particles change with temperature. Temperature dependence is extracted and implemented in our Two-Yukawa potential model as follows. Beginning with a “Soft Repulsion plus Longer Short-Range Attraction” (SR-LSRA), we produced an initial estimate for the profile V(r; T = 10 °C) for the range and strength parameters K₁, K₂, Z₁, and Z₂ using a seed estimate based on the $B2*$ just calculated. We then iteratively refined this potential until convergence to the correct value of $B2*(T=10°C)$ was obtained and the predicted I(q) from Ornstein–Zernike equation matched I(q) from experimental SANS data. Next, to construct the potential for higher temperatures (where the attractions are stronger and, presumably, higher-order interactions can become more important), one path forward is to “bootstrap” from one temperature to the next using all four parameters from one temperature as seed estimates at the next temperature. However, doing so permits the attraction strength parameter K₁ to change non-monotonically with temperature—which is both unphysical and difficult to represent analytically or computationally. Instead, for higher temperatures, (T = 25, 40, and 60 °C), we held three of the parameters fixed: the attraction range, Z₁, repulsive strength, K₂, and repulsive range, Z₂. The interparticle potential was then constructed by iteratively refining the attraction strength K₁ until the correct value of $B2*(T)$ was obtained and S(q) from theory matched S(q) measured from experiments. The justification of this process is as follows. The attraction range Z₁ can be fixed because the range of polymer bridging attractions is set primarily by the length of the PEGDA oligomer and does not change significantly with temperature.⁵⁸ Next, the repulsion strength K₂ and range Z₂ are only weakly dependent on temperature, because they are governed by screened electrostatic repulsions, which are determined primarily by the free surfactant (SDS) concentration in the bulk phase,^85–87 which is held fixed. The resultant set of V(r) curves are shown for the SR-LSRA in Fig. 7(d). For the reader’s reference, the theoretical predictions for the scattering intensity I_theo(q) at each temperature are shown in the inset in Fig. 7(a) (solid lines) for the SR-LSRA potential, which show excellent agreement with the experimental data (open symbols). This process is repeated for a “Hard-Sphere Repulsion plus Short-Range Attraction” (HS-SRA) potential and “Hard-Sphere Repulsion plus Short-Range Attraction and Longer Short-Range Repulsion” (HS-SRA-LSRR) potentials to produce the sets of curves in Figs. 7(e) and 7(f).

The conventional method (Fig. 2) produces a degeneracy of interaction potentials V(r) accurately predicts equilibrium structures. In previous applications of the conventional approach, this degeneracy can be ignored for sufficiently short-ranged potentials by invoking the ELCS. Otherwise, it is typically resolved by two approaches. The first is to make empirical guesses regarding the qualitative shape of V(r), which rules out certain candidate families of potential shapes. The second is to extend the analysis, including fitting of S(q), to sufficiently large ϕ where the high-q features of S(q) are influenced by the detailed shape of V(r) so that the predicted S(q) begins to differ significantly between different candidate shapes with otherwise equal $B2*$⁠. Although such ϕ-dependent analysis has been applied in colloidal systems involving short-ranged attractions,^28,29 this approach is challenging due to the multiple scattering events and also because the appropriate choice and accuracy of a closure approximation to the O–Z equation becomes critical due to differences in the way the different closures approximate the many-body interaction term encoded in the configuration integral in Eq. (8).

There is another, more fundamental issue with the extension of the conventional approach to larger ϕ as typically applied, which is that the analysis above, starting with the O–Z equation and continuing through the various closure approximations for g(r), assumes that the sampled S(q) is the one corresponding to thermodynamic equilibrium. Although this will always be true for systems with sufficiently weak interactions under sufficiently dilute conditions, in many cases, the equilibrium analysis under the conventional approach is applied to conditions, both in terms of sufficiently large ϕ and sufficiently strong attractions in V(r), in which the structure measured in experiment corresponds to a state of the system that is arrested out-of-equilibrium. For example, in Eberle, Wagner, and Castaneda-Priego,²⁸ interaction potentials were extracted under conditions where the system formed colloidal gels that do not represent the equilibrium state of the system, thereby invalidating the analysis, resulting in a fictive V(r) that may accurately describe the measured non-equilibrium structure, but does not accurately represent the “true” underlying V(r) that drives the system from equilibrium to begin with. Thus, the approach produces a situation in which a unique interaction potential (or temperature-dependent family thereof) cannot simultaneously predict or reconcile the observed equilibrium properties and phase behavior of the system for some regions of state space with the appearance of non-equilibrium arrested states in others.

These shortcomings of the conventional approach motivates the need for our new approach to selecting model potentials among candidate shapes, which is not based on subjective decisions and, more importantly, can simultaneously reconcile the prediction of equilibrium behavior with the appearance of non-equilibrium arrested states across the phase diagram, by down-selecting a single potential from multiple degenerate candidates that faithfully reproduces both equilibrium and non-equilibrium states. We utilize this novel approach, which builds upon the conventional approach as described in the sections to follow to directly incorporate non-equilibrium property determination (in this case, time-dependent SAOS rheological response) into the down-selection process.

Table I shows the Two-Yukawa potential parameters constructed and utilized in this study. Parameter uncertainties were estimated by computing ranges that increase the difference in predicted and experimentally measured structure |S_theo(q) − S_expt(q)| by 5%.

Intensity measurements from small-angle neutron scattering (SANS) experiments are typically smeared due to the resolution of the instrument. The SANS instrument utilized in this work for collection of dilute structure is pinhole collimated and uses a detector with finite pixel size.³⁸ While the effect of resolution smearing is often slight especially for dilute samples without sharp structural features in scattering intensities, we account for pinhole-smearing effects as described as follows. Upon detection, the incident scattering intensity, I_unsmeared(q), is smeared as

where I_smeared(q′) is the smeared intensity and R(q′, q) is the resolution function.⁸⁸ The resolution function can be described as a Gaussian-like function to good approximation as

where σ_q is the standard deviation of the Gaussian approximation to the q-resolution. In theory, the integral in Eq. (E1) can be solved to correct the measured, smeared intensity (i.e., “desmear”). However, such calculation is typically challenging. Thus, when construction of an interparticle potential, smearing the structure predicted from theory is preferred rather than desmearing the measured intensity.⁸⁹ In this work, we smear the theoretically calculated structure predicted by Eq. (8) in the iterative process described in Sec. II B 3 as

This smeared and predicted structure S_{theo, smeared}(q) is then compared to the experimentally obtained structure for iterative refinement of V(r).

Radially averaged SANS intensity measurements typically produce variable uncertainty at different q-values because of differences in area of the detector that is binned to each q-value—error bars are smaller for small-q data and larger for high-q data (see error bar sizes of markers in Fig. 8). Error bars are computed in the iterative refinement process by normalizing the error in static structure factor,

where σ_q(q) is the standard deviation of the intensity for each radially averaged intensity data point. In our V(r) refinement process described in Sec. II B 3, we iteratively minimize the error in Eq. (E4) to accurately account for the influence of resolution smearing and data uncertainty.