Fluids with competing interactions. II. Validating a free energy model for equilibrium cluster size Free

Over the past century, colloidal aggregation has been observed and described in a wide range of contexts via progressively more powerful experimental techniques, phenomenological frameworks, and quantitative models.^1–6 Spanning processes including droplet nucleation and growth, gel and glass formation, various self-assembly processes, etc., an overarching goal has been to use statistical mechanical or molecular thermodynamic approaches adopted from atomic systems and, if necessary, empirical rules to relate the strength and lengthscale of particle interactions to resulting equilibrium (or non-equilibrium) structures and the thermodynamics (and kinetics) of their formation. These types of relations, especially when based on physically intuitive thermodynamic arguments, are not only of fundamental importance but also highlight pathways for engineering materials at the nano- to microscopic level.

In this article, we focus on fluids where interactions between primary particles (monomers) are characterized by attractions acting at small lengthscales close to contact that compete with repulsions acting at larger lengthscales. This class of interactions can drive the reversible formation of equilibrium cluster phases composed of self-terminating aggregates (droplets) of monomers. Such cluster phases have been the focus of much recent work, ranging from theoretical and computational studies of idealized colloidal or nanoparticle suspensions^7–19 to experimental demonstrations for both archetypal colloidal particles^20–23 and heterogeneous monomers with anisotropic interactions like proteins.^24–31 Despite the range of materials and lengthscales, the broad underlying physics appear universal: induce (or allow) depletion (dispersion) attractions between monomers to drive aggregation while simultaneously controlling electrostatic repulsions between the ionic double-layers of monomers such that they collectively build up to attenuate growth.

While this basic paradigm of frustrating interactions is well-accepted, it is not yet established how to best describe observed cluster phases in terms of their thermodynamics and phenomenology. For example, is it possible to develop a simple and physically motivated free-energy model that can generate accurate predictions of characteristic terminal cluster size N^∗ based on experimentally tunable variables governing monomer interactions? We address this question here by directly comparing free energy-based predictions against computer simulations for one of the most approachable and idealized cluster-forming models: the short-range attractive, long-range repulsive (SALR) pair potential.⁸ Once the behavior for this simple system that coarse-grains over many microscopic details of the short-range interactions, electrostatic double-layers, solvent, etc., is better understood, the goal is then to expand the framework to include more complex free energy contributions relevant for specific realizable colloidal suspensions.

First, we review the canonical a priori free-energy treatment for clustering colloidal suspensions due to Groenewold and Kegel,^7,22,32 where we compare its predictions for the cluster size N^∗ against a computational survey of phases comprising compact spherical aggregates. We take great care to clarify how this elegant and frequently cited model adapts the classical nucleation theory (CNT) of non-terminating droplets (or crystals)^33–35 for the SALR systems of interest by treating the latter as purely attractive reference fluids superimposed with perturbative effects due to charges on the monomers and in the suspending solvent. However, while frequently cited (and adapted for related systems, e.g., driven colloids³⁶), this model has not been systematically scrutinized against a large “test set” of cluster phases generated by gradually varying relevant independent variables, e.g., monomer surface charge Z. By conducting tests that align with the phenomenological assumptions underlying the model (e.g., apolar solvents, low cluster density), we readily find that the analytical predictive formula derived from their model exhibits a spurious scaling for the range of stable cluster sizes observable in systems governed by SALR pair potentials.

With this knowledge in-hand, we describe and validate an alternative free energy model that quantitatively predicts the characteristic cluster size N^∗ for approximately 100 different simulated SALR systems, which comprise compact spherical aggregates in the size range 6 ≤ N^∗ ≤ 60 for wide ranges in monomer packing fraction ϕ, attraction strength βε , monomer surface Z, and solvent screening length κ⁻¹/d (notably, even finite values outside the apolar limit). In essence, we find that a framework built on classical nucleation theory can indeed describe the thermodynamics of frustrated, finite-sized clusters provided one introduces a size-dependent enthalpic penalty of interface formation that accounts for the missing coordination bonds of “surface” particles in clusters. In justifying this approach, we also examine how the number of intracluster short-range bonds scales with size; interestingly, we find a superlinear crossover at our cluster sizes that bridges the previously established scaling regimes for very small sizes^37,38 (e.g., N^∗ ≤ 9) and larger, bulk-like droplets. Surprisingly, we also demonstrate that intercluster effects need not be considered to obtain correct predictions even for rather non-dilute conditions.

To systematically test the performance of free energy models for predicting equilibrium cluster formation, it is invaluable to be able to (1) rapidly generate aggregate configurations that can be analyzed in-depth and (2) unambiguously identify relevant free energy contributions. Thus, we consider one of the simplest and most frequently used models known to form self-limiting aggregates: the short-range attractive (SA), long-range repulsive (LR) pair potential.⁸ The combined SALR potential can be expressed as

where β = (k_BT)⁻¹ (k_B is the Boltzmann’s constant and T is temperature); x = r/d is the non-dimensionalized interparticle separation; and d is the characteristic particle diameter. We include the subscripts i and j to account for multiple particle types because we follow previous protocols^17,39 and examine size-polydisperse three-component mixtures that approximate colloids with 10% size polydispersity. (This favors the formation of amorphous fluid clusters over crystalline dynamically arrested clusters.¹⁷) In this context, x_i,j ≡ x − (1/2)(i + j)(Δ_d/d), where i (or j) = − 1, 0, 1 corresponds to small, medium, and large particles, respectively, and Δ_d/d is a perturbation to the particle diameter. Specifically, we study mixtures comprised of 20% small, 60% medium (characteristic size d), and 20% large particles with Δ_d = 0.158d.

The short-range attractions are expressed via a generalized (100-50) Lennard-Jones (LJ) model,

where βε is the reference monomer-monomer attraction strength and Δ_ε = 0.25k_BT is an energetic perturbation to promote mixing of the polydisperse particles. Given its simplicity, the contribution of Eq. (2) (similar to the contact attractions in the free energy model of Groenewold and Kegel⁷) does not specify the microscopic or chemical details; i.e., whether the attractions arise from depletion or other short-range interactions. Generally, the range of the attraction well is approximately 0.10d.

Long-range repulsions are calculated on the basis of the repulsive portion of the Derjaguin-Landau-Verwey-Overbeek (DLVO) potential,^2,3 which approximately captures the interactions of electrostatic double-layers formed around each monomer due to (homogeneously distributed) surface charge Z. This is expressed⁶ as

with

where βA_MAX is the maximum electrostatic barrier between particles at contact, κ⁻¹/d is the Debye-Hückel screening length, Z is the total surface charge per monomer, and λ_B/d is the Bjerrum length of the solvent. Crucially, this formulation neglects any long-range multi-body interactions^40,41 and any charge renormalization due to counterion condensation^42–45 or close monomer association.^18,29 As our goal here is to test how even the simplest clustering systems might be described from a free energy perspective, we reserve incorporation of these phenomena for future studies.

In using this model, we set the average monomer packing fraction ϕ = (π/6) ρd³ (where ρd³ is number density), charge Z, and screening length κ⁻¹/d and then independently tune the attraction strength βε to drive aggregation as if varying the amount of non-interacting depletant. In terms of the experimental control that one has over repulsive contributions to interactions in such systems, this picture is somewhat idealized: to wit, tunable repulsion-controlling parameters are more realistically (though still ignoring some possible interdependence) charge Z, solvent relative permittivity ϵ_R, and solvent ionic strength I. This is because even approximately, the screening length $κ−1/d=ϵ0ϵRkBT/(2d2NAe2I)$ and λ_B/d = e²/(4dπϵ₀ϵ_Rk_BT), where ϵ₀ is the vacuum permittivity, N_A is the Avogadro’s number, and e is the elementary charge. For simplicity, however, we universally fix the relative Bjerrum length λ_B/d, which means electrostatic effects are set via combinations of Z and κ⁻¹/d. With this experimental picture in mind, we also note that the repulsive strength in Eq. (4) can equivalently be written as $βAMAX=πdϵ0ϵRΨ02/(kBT)$⁠, where Ψ₀ is the surface potential on the monomer (often assumed to approximately equal the ζ-potential measured via electrophoresis).

As illustrated in Fig. 1, we conduct a wide survey of Z-κ⁻¹/d combinations designed to span the weakest repulsions that produce self-limiting aggregates (i.e., near the boundary of macrophase separation) to repulsions with strengths up to A_MAX ≈ 2.0k_BT. Here, note that to examine this range of repulsion strengths referenced against any plausible relative Bjerrum length λ_B/d (e.g., λ_B/d = 0.014, corresponding to d = 50 nm monomers suspended in room temperature water with λ_B = 0.7 nm), one must consider monomers with very low effective charge density. Throughout the publication, we reference Z-values based on the choice of λ_B/d = 0.014, though choosing a different reference λ_B/d simply renormalizes the range of Z under consideration, with an example of this rescaling given in Fig. 1. All of the parameter combinations (ϕ, Z, κ⁻¹/d) we examine are listed in Table I by their respective critical attraction strengths (discussed below). Finally, note that throughout the remainder of the publication, we notate $βui,jSALR(xi,j)$ as βu(r) for aesthetic simplicity unless otherwise indicated.

We generate the configurations of cluster phases via three-dimensional (3D) MD simulations of the ternary SALR mixtures described above, where we generate trajectories using LAMMPS.⁴⁶ We perform simulations in the NVT ensemble with periodic boundary conditions using an integration time step of $dt=0.001d2m/(kBT)$ (taking the mass m = 1) and fix temperature via a Nosé-Hoover thermostat with time-constant τ = 2000dt. As outlined in Table I, we consider many combinations of charge Z and screening length κ⁻¹/d at four different packing fractions: ϕ = 0.015, 0.030, 0.060, and 0.120 (where we simulate N_box = 1920, 2960, 6800, and 6800 particles, respectively). Beginning with randomized initial configurations, we equilibrate systems at ϕ = 0.015, 0.030, 0.060, and 0.120 for 3 × 10⁷, 1 × 10⁷, 3 × 10⁶, and 2 × 10⁶ steps, respectively, and confirm that they are equilibrated on the basis of energy convergence and visualization, where the latter shows that the systems are ergodic (aggregates undergo frequent intra- and intercluster rearrangements and exchanges). We cut off the pair potential for a given Z and κ⁻¹/d such that the interaction strength at distance $xi,jc$ (note the explicit use of the mixture notation) is $βui,j(xi,jc)≤2e−3$ and the force is simultaneously $−d[βui,j(xi,jc)]/dxi,j≤1e−3$⁠.

To make our analysis tractable, we focus on states corresponding with the onset of clustering (i.e., at the cluster transition locus), where the phases are composed of fluid aggregates with characteristic size N^∗, but the systems have not yet begun to form percolated phases or become dynamically arrested. To characterize the size of equilibrium aggregates, we calculate cluster-size distributions (CSDs), which quantify the probability p(N) of observing clusters comprising N particles. Here, we follow the established convention^8,14,15,17 of considering two monomers as part of the same cluster if they are directly bonded to one another (i.e., within the range of the attractive well) or each directly bonded to a shared neighbor (i.e., are connected via some percolating pathway).

In turn, to locate the cluster transition locus, we make sweeps in attraction strength βε (at increments of Δε = 0.05k_BT) and identify states at the onset of clustering based on the following criteria: (1) the p(N) distribution exhibits a visibly apparent local maximum (mode) at some 1 < N^∗ ≪ N_box, where the corresponding local minimum between N = 1 and N^∗ is notated as N_min; and (2) that 80% of the particles in the system participate in aggregates of size N ≥ N_min, i.e., $0.80=∑n=NminNboxp(N)$⁠, where p(N) is appropriately normalized. Taken together, these conditions correspond to the emergence of meaningful bimodality (coexistence) in p(N) between N = 1 and the cluster mode N^∗. In this way, we obtain the characteristic cluster size N^∗ associated with a particular combination of ϕ, Z, and κ⁻¹/d and the corresponding critical attraction strength βε^∗. All of the parameter combinations we consider in our analysis are listed by their respective βε^∗ values in Table I.

Before discussing free energy models for characteristic cluster size N^∗, we begin by briefly describing the cluster morphologies under examination: for the approximately 100 different combinations of packing fraction ϕ, surface charge Z, and screening length κ⁻¹/d that we consider (listed in Table I), we observe phases at the corresponding critical attraction strengths βε^∗ that comprise compact spherical clusters with characteristic sizes in the range 6 ≤ N^∗ ≤ 60, as plotted in Fig. 2. In terms of cluster shape, we find that by measuring the radius of gyration R_G/d and plotting it versus cluster size N^∗, our results obey the relation

where α(ϕ) is a ϕ-dependent prefactor of magnitude approximately 1/2 (hereafter notated α). Together with the fractal dimension d_f = 3, this signifies that the aggregates are compact objects, and visual inspection of the MD trajectories confirms that the clusters are indeed highly packed amorphous droplets that are spherical on average and undergo frequent intracluster rearrangement and intercluster exchange (seen previously in Refs. 17 and 39). As shown in the inset of Fig. 2, the clusters do become slightly less packed with increasing ϕ, which is attributable to an increasing frequency of intercluster exchange. (These transfer events tend to instantaneously but, on average, isotropically distort the clusters, effectively expanding them.) We discuss trends in the cluster size and shape from a different perspective (and in more detail) in Paper I.⁶³ 

In terms of cluster number size N^∗, there are two important observations from Fig. 2: (1) characteristic cluster size depends only weakly on the packing fraction for the range of 0.015 ≤ ϕ ≤ 0.120; and (2) the morphologies associated with unscreened electrostatic repulsions (i.e., κ⁻¹/d → ∞) are effectively generated when the screening length approaches κ⁻¹/d ≈ 4.0. As shown by considering Figs. 2(a) and 2(b) simultaneously, increasing the packing fraction ϕ (given fixed Z and κ⁻¹/d) does not systematically shift N^∗ but does slightly inflate the cluster radius R_G/d. (We do note that the CSD peaks at N^∗ also become wider with increasing ϕ due to more frequent intercluster contacts.) The second point is apparent based on Fig. 2(a), which demonstrates that for the larger screening lengths κ⁻¹/d tested, cluster sizes N^∗ at fixed ϕ and Z have already nearly reached asymptotic values, i.e.,

This ability to access the Coulombic limit at finite κ⁻¹/d is important for Secs. III B and III C.

We now begin our discussion of the canonical framework for cluster formation due to Groenewold and Kegel⁷ (with subsequent follow-ups^22,32), with an emphasis on making clear important concepts and assumptions underpinning the model. The model aims to predict characteristic cluster size $N∞∗$ for large and perfectly monodisperse aggregates governed by short-range attractions (SA) and long-range (LR) unscreened Coulombic interactions between monomers (the subscript alludes to the κ⁻¹/d → ∞ limit). This prediction necessarily begins with an expression for the extensive free energy βΔF of cluster formation as a function of N (agnostic to $N∞∗$⁠),

where βF_N and βF₁ are the free energies of the N-sized clusters and monomers, respectively.

The free energy change βΔF is broken into reference and perturbative contributions: the reference portion is taken to be the free energy of aggregate formation for a SA (i.e., purely attractive) fluid, which can be described via the classical nucleation theory (CNT) for large droplets (or crystals).^33–35 Meanwhile, the perturbations are any contributions to the free energy due to the electrostatic effects. This is simply expressed as

where we detail these (reference) attractive and (perturbative) repulsive free energy differentials in order below.

The CNT-based free energy contributions of the reference SA system comprise two terms, which capture competing effects that scale with aggregate volume and surface area, respectively. The first term accounts for the transfer of monomers from the low-density dispersed phase to the dense (bulk) fluid or crystal phase corresponding to the cluster interior. This transfer is characterized by a favorable change in chemical potential per particle with the magnitude $βΔμ0SA$⁠. The second term is an enthalpic penalty^47,48 characterized by surface tension βγ^SAd², which accounts for the relative number of “missing” intracluster coordination bonds z_c,m of the particles at the droplet surface relative to, e.g., the bulk-like coordination number z_c,0 of the cluster interior.⁴⁹ These contributions can be written as

where, reflecting our observed morphologies, we incorporate the expression for cluster surface area assuming spherical droplets with radius R_c/d. Going forward, this radius is considered interchangeable with the radius of gyration within some O(1) prefactor, i.e., R_c ≈ R_G.

In turn, the perturbative electrostatic contributions are treated as arising from unscreened repulsions acting between all intracluster pairs of particles (i.e., N(N − 1)/2 ≈ N²/2 interactions), which can be written as

where 〈βu^LR〉 ≈ Z²(λ_B/d)/(R_c/d) is the Coulombic limit (κ⁻¹/d → ∞) of the DLVO-type potential of Eqs. (3) and (4) evaluated at r = R_c/d, which assumes that the characteristic (average) intracluster pair distance is simply the cluster radius.⁵⁰ The form of Eq. (10) implies that the repulsive free-energy contribution of each monomer in the dispersed phase is truly negligible compared to the intracluster contribution, which is consistent with the choice of Groenewold and Kegel to ignore intercluster interactions, i.e., consider the limit of very low ϕ. Note that Groenewold and Kegel also originally include a term (see Eq. 18 in Ref. 7) that roughly accounts for counterion condensation,^42–45 which could occur for strong bare surface charges. However, we neglect this contribution because their approximation naturally drops out of the subsequent analysis, and the coarse-grained SALR potential considered here only captures a constant net-effective charge.

Given these expressions for the free energy contributions, one can proceed to the crux of the analysis: identifying the characteristic cluster size $N∞∗$ at which the driving force to associate per monomer is at its largest magnitude (or energetic minimum), i.e., $βΔf(N∗)≡minN[βΔf(N)]$ where βΔf(N) ≡ βΔF(N)/N. Of course, here one requires a βΔf(N) function where the sole dependent variable is N. By combining Eqs. (9) and (10) with the known relation between cluster radius and number size R_G/d = αN^1/3 for compact spherical aggregates, one can readily write

and evaluate its derivative to find the global minimum

which, dropping prefactors, gives the scaling relation

This states that the cluster size is simply governed by the strength of the surface energy relative to the characteristic strength of electrostatic repulsion.

To write Eq. (13) completely in terms of experimentally tunable parameters, one then approximates^22,47,48 the surface tension of the SA reference fluid βγ^SAd² as scaling like the attraction strength βε multiplied by the aforementioned number of missing bonds per surface particle z_c,m (divided by a “surface area” per monomer A_m), i.e.,

Because z_c,m is considered constant with respect to N for large, low-curvature droplets, combining Eqs. (13) and (14) leads to the master a priori scaling relation,

Reintroducing prefactors, Eq. (15) is written as $N∞∗=αν0βε/[Z2(λB/d)]$⁠, where α remains from the repulsive term in Eq. (11) and ν₀ is a prefactor that is the product of z_c,m and some conversion factor to arrive at a surface energy per area.

Given our wide survey of compact spherical cluster morphologies, we can perform the first systematic test of the master scaling law given by Eq. (15) for SALR pair potentials by plotting measured cluster sizes N^∗ for systems with sufficiently large screening lengths κ⁻¹/d at various ϕ, Z, and βε. Specifically, in Fig. 3, we plot N^∗ values observed at critical attraction strengths βε^∗ and screening length κ⁻¹/d = 4.0, where the latter corresponds to effectively unscreened systems (see Section III A) as assumed in writing Eq. (15). Here, we note that we use the version of Eq. (15) that incorporates prefactors α and ν₀, which shift predicted sizes approximately in line with the measured N^∗ values (of course, including or excluding these prefactors does not affect scaling itself).

In Fig. 3, we do indeed observe a master ϕ-independent relation between the N^∗ values measured in simulations and the relative strength of attractions and repulsions between monomers, i.e., the ratio βε/[Z²(λ_B/d)]; however, the observed scaling does not reflect the exponent of 1 that is expected based on the free energy model underlying Eq. (15). Instead, we clearly observe the empirical relation,

for various cluster sizes and packing fractions. This immediately begs the questions: what alternative (and, ideally, comparatively simple) free energy model for SALR systems results in this softer master scaling? Furthermore, can this alternative model also readily predict N^∗ for finite screening lengths κ⁻¹/d?

To ascertain what new model can capture the empirically observed scaling in Fig. 3 (and be extended for generic κ⁻¹/d), we first ought to identify which of the current free energy terms in Eqs. (9) and (10) correctly (or incorrectly) describe the energetics of cluster formation in the MD simulations. Given its simplicity, the most straightforward candidate to consider is the repulsive free energy contribution of Eq. (10), which we can test against MD configurations by adding up the total repulsive energies (between all intracluster pairs of monomers) of simulated clusters as a function of characteristic size N^∗.

As shown in Fig. 4, we observe that the repulsive free energy contribution of Eq. (10) quantitatively describes MD results in the unscreened limit and, with a simple extension, also works for finite screening lengths κ⁻¹/d; in other words, the current perturbative free energy term capturing electrostatics is self-consistent and should be retained. In Fig. 4(a), we see that βΔF^LR measured in simulations, when normalized by the maximum repulsion barrier βA_MAX = Z²(λ_B/d) (corresponding to the κ⁻¹/d → ∞ limit of Eq. (4)), scales as N^5/3. Of course, this N^5/3 scaling is expected given N² intracluster pair interactions occurring on the lengthscale of the cluster radius, which scales as N^1/3 (see Eq. (10)). Meanwhile, Fig. 4(b) demonstrates that the same scaling holds for finite κ⁻¹/d away from the Coulombic limit provided one appeals to the more generalized form of Eq. (4) for the maximum repulsive barrier energy, i.e., βA_MAX = Z²(λ_B/d)/[1.0 + 0.5/(κ⁻¹/d)]².

Given that intracluster repulsions scale as expected (with N^5/3), the simplest extensive free energy expression (resembling that of Groenewold and Kegel) that readily leads to the empirically observed scaling in Fig. 3 is one where the surface-energy penalty, rather than scaling as N^2/3, instead effectively scales with a lesser exponent,

Here, ν₁ is some (as yet undetermined) dimensionless prefactor distinct from the ν₀ above. In turn, it is easily shown that solving Eq. (17) for $βΔf(N∗)≡minN[βΔf(N)]$ results in the generalized scaling N^∗ ∝ {βε/[βA_MAX]}^3/4 or, in the unscreened limit, $N∞∗∝{βε/[Z2(λB/d)]}3/4$⁠.

During the remainder of this section, our ultimate goal is to demonstrate that this reduced exponent for the surface energy term naturally emerges for our clustered systems because the effective energy penalty is dependent on cluster size N^∗ in the range 6 ≤ N^∗ ≤ 60. Conceptually, this size-dependence for the surface energy echoes the long-established notion that the generalized surface-tension of a liquid droplet with high curvature γ(R) will depart from the reference surface tension γ^∞ of a planar liquid-vapor interface (or very large droplet with low curvature). Indeed, starting with the pioneering work by Tolman,⁵¹ a vast number of studies have been dedicated to measuring first- and/or second-order corrections for γ(R)/γ^∞ (the classic first order correction depends on the “Tolman length”) to better model, e.g., homogeneous nucleation, but this topic continues to be an active and challenging area of research even for model systems like the LJ fluid.^52–58 Compared to these studies, which are especially difficult given their general focus on critically unstable droplet formation (usually droplets with radius R ≈ 4d at the smallest), the following analysis is notable because we consider stable droplets with effective surface tensions dominated by short-range attractive bonds (much shorter than, e.g., the LJ attraction range) and radii of less than three particle diameters.

Specifically, to capture this size-dependent surface energy, one ought to account for an N-dependent number of missing co-ordination bonds z_c,m(N) for the surface particles relative to the reference bulk (interior) co-ordination number z_c,0. The surface energy penalty in Eq. (17) can then be written as

with the (to be demonstrated) scaling

where we still assume that the number of surface particles at least roughly scales as N^2/3, i.e., proportional to the squared cluster radius (R_G/d)² = α²N^2/3, though making a formal distinction between interior and surface particles is difficult for small N (as discussed later). To demonstrate that the scaling in Eq. (19) is reasonable, we show in Figs. 5 and 6 that this size-dependence for z_c,m(N) = z_c,0 − z_c(N) originates based on the coordination number of (surface) particles z_c(N) measured from MD configurations, which we calculate from the extensive number of intracluster bonds n_B(N). Given our measurement of n_B(N) is at the root of much of this analysis, we consider its behavior first and proceed backwards to the scaling of Eq. (19).

Looking towards estimating z_c,m(N), consider in Fig. 5(a) the extensive number of intracluster bonds n_B(N) measured from MD simulations, where we observe a superlinear growth rate over the range of cluster sizes that we generate. Interestingly, this superlinear behavior contrasts with known small- and large-cluster limits, which are linear in N. Here, n_B(N) is nicely captured at each packing fraction for 6 ≤ N^∗ ≤ 60 by the empirical expression,

where k is a ϕ-specific O(1) prefactor⁵⁹ and we include a division by 2 for aesthetic purposes. This superlinear regime contrasts with the small cluster regime (3 ≤ N ≤ 9), where it is known^37,38 that colloidal clusters dominated by SA bonds maximize their extensive bonding number according to the expression n_B(N) = 3N − 6. Likewise, in the limit of large droplets, the number of bonds must scale increasingly like in the corresponding bulk fluid, i.e., n_B(N) → (z_bulk/2) N, where z_bulk is the co-ordination number of the reference fluid (or crystal) phase.

To quickly understand why n_B(N) growth should be superlinear over this size range, we show in Fig. 5(a) extensions of the small- and large-cluster linear regimes (to large and small N where they should, respectively, fail) to demonstrate that the function n_B(N) = (k/2) Nln(N) connects these otherwise disparate limits, while quantitatively overlapping with the upper reaches of the small cluster trend at N ≈ 10. To wit, notice that the characteristic slope of the small-N regime, (m = 3), differs meaningfully from the typical slope in the large-N regime of a very dense bulk fluid or crystal, which we estimate as m = z_bulk/2 = 6 with z_bulk = 12 because it is the sphere kissing number in three dimensions⁶⁰ (this is justified later). Thus, provided z_bulk is decidedly larger than 6, a superlinear regime allows for a smooth continuous growth in n_B(N) with respect to N.

This connectivity between very small and large cluster sizes is clearly echoed by the next necessary quantity we must calculate: the average co-ordination number z_c(N) = 2n_B(N)/N = kln(N), which we show in Fig. 5(b) for the same selected clustered states.⁶¹ Here, we plot z_c(N) values calculated from MD configurations, which begin to bridge the gap (up to the highest cluster sizes we observe) between the highly bond-restricted regime at small N and the bulk regime at large N , where the co-ordination number approaches z_c(N) → z_bulk. Notably, z_c(N) varies by approximately a factor of 2 over the size range of interest 6 ≤ N^∗ ≤ 60, which underlines that the conventional practice (for larger droplets) of assuming that surface effects are size-independent is problematic for these smaller aggregates.

With z_c(N) in hand, we can proceed to calculate the average number of missing bonds per particle z_c,m(N), which indeed collapses onto a master curve scaling as N^−1/3 (shown in Fig. 6) when the magnitude of the reference (fitting) co-ordination parameter z_c,0 is set—in line with measurements of cluster interiors—at values appropriate for highly packed bulk fluids. To do this, we use the expression

where the only as-yet undetermined value is z_c,0, which is the co-ordination number of the reference bulk SA fluid that represents the idealized cluster interior. For our immediate purposes, we treat this parameter as tunable and verify our choices as reasonable below. As shown in Fig. 6(a), our data approximately collapse onto a master curve with characteristic N^−1/3 dependence when z_c,0 = 12.0, 11.5, and 10.5 for ϕ = 0.015, 0.030, and 0.060, respectively. All of these values—especially for the lowest-density case—are reflective of bulk fluids dominated by short-range attractions, especially here given that energetic gains from bonding occur within attractive wells beyond surface contact that are approximately 0.1d in width.

In Fig. 6(b), we demonstrate that these z_c,0 values are appropriate based on direct measurements of the locally averaged co-ordination number z_c(r) as a function of radial position within clusters (relative to cluster center-of-mass). Here, we specifically show results from some of the largest clusters observed (50 < N^∗ < 60), which are most likely to possess bulk-like interiors as r → 0; indeed, it is evident that z_c(r → 0) ≈ 12 for these larger clusters, though the limiting value (as above) slightly decreases as ϕ increases, presumably due to the previously discussed trend in intracluster density. We also observe in Fig. 6(b) that the z_c(r → 0) limit is similar even for smaller clusters, e.g., N^∗ ≈ 20, where central particles can still be surrounded by a packed shell of intracluster neighbors.

Taken altogether, the results of Figs. 5 and 6 nicely justify the choice to quantify the surface energy penalty of cluster formation from the perspective of a size-dependence in the relative number of missing bonds z_c,m(N). Before moving on to consider the impact of the scaling relationships in Eqs. (18) and (19) for predicting cluster size, we pause to note that in the analysis above, we approximate z_c,m(N) for surface particles based on an average measurement of n_B(N) for all cluster constituents. We take this somewhat imprecise approach because it draws upon relatively unambiguous measurable quantities and bypasses the fraught process of definitively distinguishing between surface and interior particles (consider, e.g., Fig. 6(b)). Our approximation is sufficient for the proof-of-concept analysis here, but we imagine a more exacting analysis in this vein would be a worthwhile future endeavor.⁶² 

Based on the new scaling for the free energy surface penalty justified in Section III D and the generalized free energy term for repulsive contributions in Section III C, we can return to the extensive free energy model of Eq. (17) and readily derive a new master equation for predicting cluster size N^∗ based on experimentally tunable parameters,

where, as before, α is the known ϕ-dependent prefactor relating the cluster radius and number size and ν₂ is a constant similar to those above that scales the surface energy penalty, which we treat as an empirical tuning parameter. Eq. (22) is the central result of this publication.

As demonstrated in Fig. 7, Eq. (22) successfully predicts characteristic cluster sizes for the vast majority of our ≈100 cases over various Z-κ⁻¹/d combinations and near order-of-magnitude ranges in both size 6 ≤ N^∗ ≤ 60 and bulk monomer packing fraction 0.015 ≤ ϕ ≤ 0.120. This wide applicability is notable as the underlying free energy framework remains very simple: to wit, intercluster effects can evidently be neglected even as conditions become less dilute (e.g., ϕ ≈ 0.120), though the current model cannot predict more subtle trends known for the SALR model³⁹ like the growing polydispersity of aggregates with increasing ϕ. Meanwhile, the biggest deviations between measured and predicted N^∗ (larger than 20%) occur for states combining large charge (e.g., Z = 15.0) and small screening length (e.g., κ⁻¹/d = 0.70), which results in rather non-idealized repulsions that are both strong relative to k_BT and far from the Coulombic limit. Finally, note that the value for prefactor ν₂ that shifts the (already collapsed) predictions into the correct range is $ν2≈2π$⁠, which we expect should apply rather generally for compact colloidal clusters as it simply converts between measurements of the cluster surface-size based on population and radius.

We have validated a new and readily applied formula (Eq. (22)) that can predict the characteristic cluster size N^∗ for idealized SALR suspensions as a function of the variables controlling monomer-monomer interactions (including attraction strength βε, surface charge Z, and screening length κ⁻¹/d). Eq. (22) and its underlying free energy model represent a semi-empirical adaptation and extension of the canonical free energy model due to Groenewold and Kegel,^7,22,32 where we find that the latter exhibits a spurious scaling of N^∗ away from the large-droplet limit with respect to the ratio of attractive and repulsive interaction strengths driving aggregation. We subsequently find that Eq. (22) performs excellently based on direct comparisons of predicted cluster sizes and measurements of N^∗ from MD simulations of approximately 100 different systems for very wide ranges in ϕ, Z, and κ⁻¹/d, where we examine states at the onset of clustering that exhibit compact spherical aggregates in the size range 6 ≤ N^∗ ≤ 60.

The predictive quality of Eq. (22) demonstrates that a simple free energy model, which treats SALR systems as reference SA fluids (via classical nucleation theory) with additive repulsive perturbations due to electrostatic effects, can be applied down to extremely small cluster sizes (N^∗ < 10) provided one properly corrects for surface effects at small N^∗. Conceptually, this is in the spirit of long-standing investigations regarding size-dependent surface tensions in small droplets,^51–58 and practically, we find that one can treat the energy penalty of “interface” formation as a function of an N-dependent number of missing co-ordination bonds z_c,m(N) for surface particles (referenced against the co-ordination number in the bulk fluid). Here, this picture is validated in part by configurational analysis of the number of extensive intracluster bonds n_B(N), which exhibits a superlinear scaling regime over the size range 6 ≤ N^∗ ≤ 60. Meanwhile, based on the form of the free energy model, we confirm that intercluster effects can be neglected even for rather non-dilute conditions (e.g., ϕ = 0.120), which are reflected by our observation that cluster size N^∗ exhibits little variability with respect to ϕ given otherwise fixed conditions.

We look forward to testing the predictive capability of Eq. (22) for real colloidal suspensions that exhibit equilibrium cluster phases, which could help bolster whether SALR pair potentials are a sufficient (if idealized) description of experimental systems. For instance, there has been recent discussion^16,18,29 in the literature as to whether accounting for charge renormalization during aggregation is necessary for describing cluster behavior; likewise, Groenewold and Kegel initially postulated that non-trivial charge effects^42–45 could affect the free energy picture in certain limits. Of course, these effects are not captured by the canonical SALR pair potential examined here, but we now possess a free energy model known to describe this simpler system. Thus, ascertaining whether cluster size N^∗ scales in experiments similarly to the empirical scaling of Eq. (22) would help clarify the degree to which phenomenology and interpretive guidelines derived from the pairwise model are appropriate for real systems. Similarly, it is fascinating to consider how accounting for size-dependent surface effects, here so crucial for producing quantitative predictions, might change for less compact (e.g., elongated) aggregates than those considered here.

This work was partially supported by the National Science Foundation (Grant No. 1247945) and the Welch Foundation (No. F-1696). We acknowledge the Texas Advanced Computing Center (TACC) at The University of Texas at Austin for providing HPC resources.