Modelling the role of interaction heterogeneity in the gelation of micron-scale colloidal systems

Abstract

Micron-scale colloids functionalized with supramolecular moieties represent a versatile platform for self-assembly. Single-stranded DNA functionalization, in particular, has been shown to be a highly tunable approach for inducing short-ranged attractive interparticle interactions that can be used to drive self-assembly into a wide range of crystalline structures. Recently, it has been noted that variability in the extent of surface functionalization across a population of colloids results in ‘interaction heterogeneity’, or IH, in which the binding strength between a pair of colloids varies according to the density of DNA strands on their surfaces. We have shown in previous work that IH strongly impacts colloidal crystal nucleation and growth but its impact on gelation further away from equilibrium conditions remains underexplored. In this study, we employ molecular simulations to systematically investigate the role of IH in colloidal gelation driven by thermal quenches. We consider four types of IH distributions: monodisperse (no IH), Gaussian, bidisperse, and uniform distributions, and analyze their effects on gel structure and gelation dynamics. Our results show that while IH minimally impacts macroscopic gel structure, it profoundly alters the local gel environment, as revealed by coordination number (CN) distributions. Principal component analysis of CN moments highlights distinct structural trends arising from the presence of IH, underscoring the sensitivity of local gel structure to IH. We also show that IH leads to a sequential aggregation of strong and weak binders, where strong binders first form a gel ‘backbone’ and weak binders subsequently decorate it. These findings highlight IH as a key parameter for modulating gel microstructure without significantly perturbing macroscopic organization.

---

1. Introduction

Micron-scale, short-ranged attractive colloids represent a broad class of self-assembling particle systems and are important constituents in various natural and engineered products, such as food,^1,2 cosmetics,³ and protein therapeutics.⁴ However, naturally derived colloidal formulations typically consist of particles with significant variability in size, shape, and poorly defined interactions, limiting their suitability for fundamental studies. To address these limitations, engineered model systems based on synthetic spherical colloids with well-characterized interparticle interactions have been developed. These systems have enabled detailed studies of diverse self-assembly phenomena, such as nucleation and growth processes,^5,6 crystallization,^7–12 gelation,^13–16 and even glass formation.^16,17 A widely-studied engineered example in this class is the DNA-functionalized, microscale colloidal particle.^5,11,12,18 Here, plastic spheres with diameters typically ranging from ∼200 nm to about a micron, are grafted with polymer brushes terminated with single-stranded DNA oligomers. The DNA oligomers are sequenced so that strands on neighbouring particles in solution are driven to reversibly hybridize at a certain temperature, producing an effective interparticle interaction that qualitatively resembles a simple pairwise interatomic potential, e.g., a Lennard-Jones or Morse function. A key advantage of the DNA-functionalized particle system is the high degree of specificity of DNA hybridization, enabling complex multicomponent, self-assembly processes with precisely tuneable interparticle attraction.^19–24

A large number of experimental and theoretical studies have provided substantial mechanistic insight into DNA-mediated colloidal assembly, with much of the emphasis being placed on the formation of complex, ordered assemblies.^{6,11,12,18,25–37} Yet some uncertainties remain, largely due to inherent ‘real world’ complexities. Among these is the role of hydrodynamic interactions between micron-scale particles during self-assembly, which can impact both crystal and amorphous assembly outcomes.^38–44 More recently, we have uncovered another such factor that can significantly impact the self-assembly behaviour of DNA-functionalized particles, namely interaction heterogeneity (IH). Unlike atoms, a population of colloids may exhibit variability in particle attributes, such as size or shape, resulting in heterogeneity in the interparticle interactions. While size, and to a lesser extent, shape polydispersity have been addressed in prior studies,^45–50 the DNA-functionalized colloidal system may also exhibit interaction heterogeneity as a consequence of the DNA grafting process itself. Specifically, it is known that different grafting methods can lead to different amounts of variability in the DNA coverage across a population of particles. Because the strength of the attraction between a pair of particles is directly governed by the number of DNA duplex bridges formed, this variability results in a variability, or heterogeneity, in the interparticle interactions. This type of population heterogeneity is fundamentally distinct from single-particle heterogeneity, e.g., due to patchiness in the grafting density over a single particle's surface. The latter type of heterogeneity has been addressed previously in the self-assembly literature, most notably in the context of mimicking directional bonding.^51,52 While we cannot exclude the possibility of natural patchiness in experimentally realized DNA-grafted particles, here we assume that each particle is uniformly coated with DNA strands.

Our previous studies demonstrated the potential of IH as a tuneable control parameter for colloidal self-assembly.^53,54 Using simulations we showed that unimodal IH modified crystallization dynamics by spreading out the nucleation process as temperature was lowered and attractive forces between particles increased. This ‘spreading out’ effect results from a stratification and sorting of the particle population according to binding strength and was found to enhance crystallization robustness by suppressing gelation under strongly driven conditions. Moreover, the effect of IH was found to be modulated by the width of its distribution, i.e., more heterogeneity led to increased crystallization robustness.⁵³ This is a useful result given that crystallization of micron-scale DNA-functionalized particles is often restricted to a very narrow temperature range due to the high temperature sensitivity of DNA duplex formation.^11,12,19 We have also shown that IH is important even under strictly equilibrium conditions. Specifically, unimodal IH shifts the fluid-side of the fluid-crystal coexistence boundary to lower volume fractions, thereby expanding the phase coexistence region in which crystallization occurs. These results were broadly consistent across different IH distributions, suggesting that intentionally engineered IH distributions could be used to optimize crystallization.

In the present study, we analyse the impact of IH on gelation of DNA-functionalized micron-scale particles. While colloidal gelation has received far less attention than crystallization, colloidal gels also offer the potential for creating assemblies with useful properties, such as structural coloration.^55–57

Gelation is notoriously hard to control because it is generally operational under highly non-equilibrium conditions. The mechanism of gelation in short-ranged attractive colloidal systems is multifactorial in nature and has been discussed in detail in a recent review by Royall et al.⁵⁸ In a pioneering study, Lu et al.¹⁴ demonstrated experimentally that gelation in these systems is often initiated with a phase separation in the form of liquid–gas spinodal decomposition. The liquid regions of high particle density subsequently become kinetically arrested, forming a gel phase. Nonetheless, certain systems, such as the patchy colloid system studied by Bianchi et al.,⁵⁹ have been shown to lead to gelation in the absence of spinodal decomposition. Moreover, other driving forces such as gravitational settling can also provide the necessary densification for gelation to proceed.⁶⁰ Recent work has also emphasized the emergence of percolation, mechanical rigidity, and aging across the gel phase as necessary elements of gelation.^61–63 Depending on the depth of the quench used to drive gelation, the arrested gel may be locally crystalline or fully amorphous, but precise control of gel structure is otherwise very hard to achieve. Moreover, the potential role of IH in gelation is also far from obvious because the highly non-equilibrium conditions and rapid kinetic arrest may prevent the particle sorting that lies at the heart of its ability to impact self-assembly.

Numerous prior computational studies have investigated the role of quenching rate and/or particle volume fraction on gel structure^64–66 and mechanical response.^63,67,68 Here, we investigate how gelation structure is specifically impacted by a combination of the IH distribution and the applied thermal quenching profile. We systematically explore the influence of different types of IH distributions on gel structure, including monodisperse (no IH), Gaussian, uniform, and bidisperse distributions. The results show that although IH minimally impacts macroscopic gel structure (e.g., as determined by the static structure factor), the local gel environment, as measured by the distribution of particle coordination number, is found to be sensitive to IH. Most notably, the gel configurations obtained with IH are found to be otherwise inaccessible. Collectively, our findings confirm the relevance of IH in the setting of gelation and provide a practical potential pathway for realizing designer colloidal gel structures.

2. Methods

2.1. Simulation methodology

We employ standard techniques for simulating gelation in attractive particle systems whereby the attraction between particles is gradually increased in time to drive gelation.^63,64,69 All gelation simulations were performed using Langevin dynamics implemented in the LAMMPS⁷⁰ software package. Each simulation consisted of 4000 particles in a cubic domain with periodic boundary conditions and evolved in the NVT ensemble using a Langevin thermostat at a solvent viscosity of 1 cP. Selected cases were repeated with larger systems containing 16 000 particles to assess any finite size effects.

The pairwise interaction potential used to evaluate particle forces is based on a validated coarse-grained model for DNA-functionalized colloids developed by Rogers and Crocker.¹⁹ The interaction between DNA-functionalized colloids is derived from the summation of attractive and repulsive contributions. The attractive potential is derived from the bridging dynamics of DNA strands with complementary end sequences on the interacting colloids, while the hybridization entropy and enthalpy are derived from a nearest-neighbour model.⁷¹ The repulsive interaction is calculated by the entropy loss due to DNA strand confinement from adjacent colloids.

All parameters for the interparticle potential were fixed at the same values used in our previous study of IH-mediated crystallization.⁵⁴ Briefly, all DNA strands are assumed to be 40 nm in length, explicitly represented by 8 segments each with a Kuhn length of 5 nm. The enthalpy and entropy of hybridization between individual strands were fixed at −47.8 kcal mol⁻¹ and −0.1395 kcal mol⁻¹ K⁻¹, respectively. A reference potential, U_ref(r) was defined based on a colloid diameter of 1 µm and a DNA strand areal density of 1592 strands µm⁻², or 5000 strands per particle. In addition to these parameters, the DNA potentials are sensitively dependent on temperature, which was used to fix the depth of the attractive well, . A series of potentials with well depths ranging from U = 0 to U = −15k_BT at intervals of 0.1k_BT were precomputed and used to implement thermal quenches; see Fig. S1. Note that these well depths correspond to a small, implied temperature range of about 15 K centered around 300 K. Consequently, the actual simulation temperature employed in the NVT Langevin dynamics simulations was held fixed at 305 K for all simulations since particle diffusivity changes little over this temperature range.

Interaction heterogeneity was modelled by assigning a binding modulator, b_i, to each particle i to scale the interaction potential, so that the interaction potential between any two particles, i and j, is given by

U_ij(r) = b_ib_jU_ref(r) (1)

where U_ref(r) represents the reference DNA potential energy function, which was taken as the mean of the IH distribution. In other words, a given IH distribution is encoded entirely by the distribution of the binding modulator, for a given U_ref(r). Four different types of IH distributions were investigated, including monodisperse (no IH), Gaussian, bidisperse, and uniform distributions (see Fig. 1). To ensure consistency across different types of IH, the IH distribution was normalized such that the average squared binding modulator was unity, i.e. b² = 1. This normalization implies that the average interaction strength of the system corresponds to the average binding strength of two randomly selected particles. The IH distributions were specified as follows: the Gaussian distribution was fixed with unit mean and a standard deviation, σ_IH = 0.15. The bidisperse distribution contains two species, strong and weak, with distinct binding modulators, b_s and b_w, respectively. The ratio of binding modulators between the two species was used to represent bidisperse distribution, i.e. b_s/w = b_s/b_w. In this study, we fixed b_s/w = 1.50 for all simulations. Finally, the uniform distribution was defined with a unit mean and range [1 − Δb, 1 + Δb], where Δb assumed values of 0.1, 0.3, or 0.5.

Fig. 1 Examples of different IH distributions. (a) Monodisperse i.e. no IH. (b) Gaussian distribution with unit mean and standard deviation σ_IH. (c) Bidisperse distribution. Ratio of strong and weak binding modulator, b_s/w, is independently fixed with composition. (d) Uniform distribution with half range Δb.

Linear quench protocols, characterized by a constant rate of increase of the well-depth of the reference potential, dU/dt, were applied for all finite quench rate simulations. Due to technical constraints in LAMMPS, a stepwise approximation of the linear quench protocol was implemented with discrete steps in interaction strength of magnitude ΔU = 0.1k_BT. This interval size was determined to be functionally equivalent to a continuous quench over the quenching rates used in this study. Each quench simulation was first equilibrated with a purely repulsive potential, i.e., U = 0, followed by incremental increases of the well depth as described above, each maintained for a duration specified by the imposed quenching rate. A simulation was terminated when the system configuration became obviously kinetically arrested. Step quench simulations, equivalent to an infinitely fast quench, were implemented by a reference attractive interaction strength fixed at U = −15k_BT and held for 50 seconds following equilibration. To suppress crystallization in all simulations, particle size polydispersity was introduced in the form of a Gaussian distribution with unit mean and standard deviation of 0.05, independent of the IH distribution. The interparticle pair potential was accordingly shifted based on the sum of particle radii to reflect the effect of size polydispersity. Note that polydispersity in size would, in practice, result in an additional source of IH because of corresponding changes in particle curvature. We do not consider this source of IH here and instead focus solely on “direct” IH due to variations in grafting density across the particle population. In principle, a larger database of precomputed interaction potentials could be generated that accounts for both sources of IH but this would not qualitatively alter the conclusions of the present study. A summary of all quench simulations included in this study are listed in Table 1 and Table 2.

Table 1 Simulations with no IH and Gaussian distributions included in this study. dU/dt is in unit of k_BT s⁻¹. Circles: 4000 particles. Triangles: 16 000 particles

	ϕ = 0.20		ϕ = 0.30
	No IH	Gaussian, σIH = 0.15	No IH	Gaussian, σIH = 0.15
−0.005	O	O	O	O
−0.001	O	O	O	O
−0.020	O, ▲	O, ▲	O	O
−0.050	O	O	O	O
−0.200	O	O	O	O
−0.500	O	O	O	O
−2.000	O, ▲	O, ▲	O	O
−5.000	O	O	O	O
Step	O, ▲

---

Table 2 Simulations with bidisperse and uniform distributions included in this study. dU/dt is in unit of k_BT s⁻¹. Circles: 4000 particles. Triangles: 16 000 particles

	ϕ = 0.20
	Bidisperse, bs/w = 1.50			Uniform
	p(bs) = 0.20	p(bs) = 0.50	p(bs) = 0.80	Δb = 0.10	Δb = 0.30	Δb = 0.50
−0.020	O	O, ▲	O	O	O	O
−0.200	O	O	O	O	O	O
−2.000	O	O, ▲	O	O	O	O
Step	O	O	O			O

---

2.2. Structural analysis methods

We employed several different order parameters and functions to assess the impact of quench rate and IH on gel outcomes. Below, we briefly summarize the various measures; many of these are standardized and implemented in LAMMPS and other software packages.

2.2.1. Coordination number (CN) distribution. Particle coordination number distributions were computed from particle configurations at each time point. Neighbours were identified for each particle as those within a cutoff distance defined as the sum of the particle radii plus a shell thickness of 0.04 µm. The CN distribution represents the fraction of particles corresponding to each CN. Due to size polydispersity, CN as high as 13 may occasionally occur using a fixed cutoff radius.

2.2.2. Radial distribution function. The radial distribution function, g(r), was calculated using the Python package Freud.^72,73 To obtain g(r), distances between particle pairs were histogrammed and normalized by the average number density, further adjusted by the volume of spherical shells at each radial bin, i.e.,where ρ is the number density, N is total number of particles, |r_i − r_j| is the distance between particles i and j, δ is the Dirac delta function, and 〈·〉 denotes ensemble average. Here, g(r) was computed over a distance range from 0 to 6 µm, using a bin width of 0.1 µm.

2.2.3. Static structure factor. The static structure factor, S(k), was calculated using the Python package Freud.^72–74 S(k) quantified spatial correlations among particles in reciprocal (Fourier) space, typically obtained by evaluating the ensemble average of the Fourier transform of particle density fluctuations. S(k) is defined by
where k is the wavevector, (r_i − r_j) is the displacement vector between particles i and j, and 〈·〉 denotes ensemble average. Isotropic averaging was performed over all wavevectors sharing the same magnitude, i.e., k = |k|. In this study, S(k) was calculated at discrete wavevector magnitudes defined as k = 2πn/L, where L is the edge length of the cubic simulation box, and integer values n ranging from 1 to 40.

2.2.4. Fractal dimension. The fractal dimension of gels were calculated using the box-counting method.^75,76 In this approach, the simulation box is divided into smaller cubes of varying side length, s. For each cube size, the number of cubes containing portions of the gel structure, denoted as N_c(s), is counted and the fractal dimension, d_f, is given byIn practice, d_f is determined by fitting the linear region of a plot of log N_c(s) versus log(1/s), thereby characterizing the scaling behaviour and self-similarity of the gel structure across different length scales. Cube lengths within the range of 1.2 to 2 µm were specifically used for this linear fit.

2.2.5. Topological cluster classification (TCC). Topological cluster classification (TCC) was performed using the software package developed by Royall and coworkers.^77–79 Briefly, TCC is used to quantify local structure in the gel by identifying bonded particle networks and classifying their topology into a library of well-defined motifs. Here, the bond network was constructed using a cutoff radius of 1.19 µm. The TCC algorithm first identifies all 3-, 4-, and 5-membered shortest-path rings; these are subsequently used as the basis for identifying more complex topological motifs and larger clusters. To avoid double-counting, the fraction of particles associated with each cluster size accounts only for the largest identified topological cluster to which a particle belongs.

2.2.6. Principal component analysis (PCA). PCA was applied to the histogrammed properties of the final gel configurations using the Scikit-learn package in Python.⁸⁰ PCA is a statistical technique that reduces the dimensionality of multivariate data by transforming the original variables into a new set of orthogonal axes, or principal components (PCs), which are ranked according to the amount of variance they capture from the input data. The PCA procedure used in this study included several steps. First, the raw data were organized into a matrix X with dimensions (n × p), where n is the number of IH cases, and p is the number of histogrammed bins of a specific property. To ensure all features contribute equally, each column of X was normalized to zero mean and unit standard deviation, i.e.,where X_ij is the original data for the ith IH cases and jth bin, µ_j and σ_j are the mean and standard deviation of the jth column, and X̃_ij is the normalized X_ij. Next, the covariance matrix C was computed asyielding eigenvalues λ_k and eigenvectors u_k, satisfying the relation

Cu_k = λ_ku_k, (7)

where k = 1,2,3, …, p. The eigenvectors u_k represent principal directions in feature space, while the eigenvalues λ_k measure the variance along each corresponding principal component. To quantify the relative importance of each principal component, the eigenvalues were normalized by their sum, providing the fraction of total variance explained (FVE), i.e.,

3. Results and discussion

3.1. Influence of quench rate on gel structure without IH

The influence of quench rate on monodisperse (no IH) systems was investigated first to provide a baseline. Shown in Fig. 2(a) are the kinetically arrested configurations at varying quench rates at a fixed volume fraction, ϕ = 0.20. Visually, the structures transition from open, dendritic networks at high quench rates to more compact clusters as the quench rate decreases. The impact of quench rate on gel coordination number (CN) distribution is shown in Fig. 2(b). As expected, the distributions shift towards higher average CN at lower quench rates because slower quenching provides sufficient time for particles to settle into energetically favourable configurations before becoming kinetically arrested. While all the distributions remain unimodal across the cooling rates considered here, closer inspection reveals a transition in the coordination number distributions at a cooling rate of approximately of dU/dt = −0.20 k_BT s⁻¹ (green curve in Fig. 2(b)). Faster quench rates generally produce distributions peaked at a coordination number of 6, while slower quenches result in a peak located between 7 and 8.

Fig. 2 (a) Final gel configurations obtained at different quench rates (in units of k_BT s⁻¹) without IH at ϕ = 0.20. (b) CN distributions at different quench rates (unit: k_BT s⁻¹) with no IH at ϕ = 0.20. (c) Sorting length and the corresponding multiples of mean free path (see text) as a function of quench rate. Orange dash line: L/L_MFP = 1. Green dash line: L/L_MFP = 3.

To understand the transition in peak coordination number and obtain a quantitative measure of “fast” versus “slow” thermal quenching, we evaluated the quench rate in the context of particle sorting during gelation. The sorting length, L, was determined through a scaling analysis that compares the particle diffusion timescale, τ_D = L²/D, with the gelation timescale,

where D is defined as the particle diffusivity at a given volume fraction, and ΔU is the interaction strength range over which gelation occurs. Equating the two timescales gives an approximate sorting length ofAssuming particle binding becomes essentially irreversible for |U| > 5k_BT and that particle sorting does not become significant until |U| ∼ 3k_BT, suggests that the relevant interval for particle sorting during a quench is about ΔU ∼ 2k_BT in width; see Fig. S2. Given a particle diffusivity at ϕ = 0.20 of approximately 0.30 µm² s⁻¹, the sorting length is estimated to be L ≈ 5.5 µm at a quench rate of dU/dt = −0.20 k_BT s⁻¹, or about five particle diameters. The impact of the sorting length was assessed by comparing it to the particle mean free path, L_MFP, which is defined as⁸¹where v_p is the single particle volume and R_p is the particle diameter. Note that the expression in eqn (11) is derived by assuming particles are spatially well mixed and interact purely entropically, i.e., via hard-sphere interactions. In the quench simulations, L_MFP therefore represents the average distance that a particle travels before colliding with other particles. For ϕ = 0.20 and R_p = 1 µm, L_MFP is estimated to be 1.85 µm. Consequently, the ratio between sorting length L and L_MFP represents the average number of collision events (through which sorting occurs) experienced by each particle during gelation. On this basis, we identify L ∼ 3 L_MFP as the threshold for a “slow quench”, corresponding to a quench rate of −0.02k_BT s⁻¹, (Fig. 2(c)). This estimate is fully consistent with the observed transition in the coordination number peak observed in Fig. 2(b).

3.2. Influence of IH on macroscopic gel structure

We begin our analysis of IH impact on gelation by considering various types of IH distributions at a fixed quench rate of −0.020k_BT s⁻¹, which as described in the prior section, corresponds to a sorting radius of about five particle diameters. Shown in Fig. 3 are the final, kinetically arrested configurations obtained with different IH distributions. The upper row displays the gel structures with uniform particle colouring, revealing no apparent differences in gel structure produced by the various IH types. The lower row shows particles coloured by their binding modulator values: red represents strong binders (b > 1), blue represents weak binders (b < 1), and white corresponds to b = 1. Compared to the monodisperse case (no IH), configurations incorporating IH exhibit localized regions enriched with strong binders (indicated by green arrows). These local aggregations suggest that the presence of IH may influence the gelation process by promoting spatial heterogeneity in particle binding, potentially altering the resulting gel structure in subtle ways that are not immediately apparent by visual inspection.

Fig. 3 Final colloidal gel configurations obtained at dU/dt = −0.020k_BT s⁻¹ and ϕ = 0.20 with various IH distributions. The upper and lower rows show identical configurations with different colour schemes. (Upper row: particles shown single colour. Lower row: particles coloured by binding modulators. Green arrows: local regions enriched with strong-binding species).

The potential impact of IH on macroscopic gel structure was probed further using the radial distribution function, g(r), and structure factor, S(k), as shown in Fig. 4 for the same quench rate of dU/dt = −0.020k_BT s⁻¹. Neither g(r) nor S(k) show significant impact arising from IH. However, gels without IH exhibit subtly elevated S(k) values within the intermediate range (k ≈ 2 − 4), suggesting slightly higher structural correlations at length scales corresponding to roughly 2–3 particle diameters. In contrast, at even lower k values (k < 1), gels with bidisperse and uniform IH distributions show slightly higher S(k) values compared to gels without IH, reflecting enhanced aggregation over larger length scales due to IH. Overall, the differences in S(k) between systems with and without IH are subtle, suggesting a lack of significant large-scale structural distinctions with IH.

Fig. 4 Structural analysis of colloidal gels formed with various IH distributions using (a) radial distribution function, g(r), and (b) structure factor, S(k) at dU/dt = −0.02k_BT s⁻¹ and ϕ = 0.20. (c) Fractal dimension, d_f, as a function of quench rate at ϕ = 0.20. Squares correspond to data from larger systems (16 000 particles). Data from step quenches are plotted at a quench rate value of dU/dt = −10.0k_BT s⁻¹ and denoted by the arrow.

Another metric commonly used to characterize gel structures in the literature is the fractal dimension, d_f. The fractal dimension characterizes the spatial complexity and self-similar arrangement of colloidal gel networks. Higher d_f values indicate denser and more compact structures, while lower d_f values correspond to more open, porous morphologies. Typically, diffusion-limited cluster aggregation (DLCA), governed by rapid and irreversible particle binding limited only by diffusion, yields lower fractal dimensions in the range d_f ≈ 1.7 − 1.9. In contrast, reaction-limited cluster aggregation (RLCA), characterized by slower binding kinetics that allow extensive structural rearrangement, produces gels with higher fractal dimensions around d_f ≈ 2.0 − 2.2.^82–88 Shown in Fig. 4(c) are the fractal dimensions d_f as a function of quench rate for various IH distributions. The value of d_f increases from approximately 1.8 to 2.3 as the quench rate becomes slower, indicating a transition of the gel structure from DLCA to RLCA. At lower quench rates, the fractal dimensions across different IH conditions exhibit minor differences. While at intermediate quench rates gels formed in systems with IH appear to lead to higher d_f, results from larger system simulations containing 16 000 particles (denoted by square symbols) suggest that this is likely due to statistical noise or systematic finite size effects.

3.3. Influence of IH on local gel structure

Next, the effect of IH on the local gel structure as represented by the coordination number (CN) distribution was investigated. In Fig. 5(a), the CN distribution is shown for a step quench, where the impact of IH is expected to be negligible and indeed no effect is observed. Interestingly, however, for a slow quench rate of dU/dt = −0.020k_BT s⁻¹, the influence of IH on CN distribution is quite pronounced; see Fig. 5(b). In the absence of IH, the CN distribution is a unimodal peak, a feature that persists across all quenching rates as shown previously in Fig. 2(b). The imposition of Gaussian IH broadens and flattens the distribution, forming a plateau between CN = 4 and CN = 8. Even more notable is the effect of bidisperse IH here fixed at p(b_s) = 0.50, representing equal fractions of weak and strong binders, which generates a distinct bimodal CN distribution. Finally, uniform IH produces a CN distribution also exhibiting two peaks but somewhat broader than those resulting from bidisperse IH. All three IH distributions lead to higher fractions of particles at the extremes of the CN distribution compared to the monodisperse reference case. Notably, IH minimally impacts the average CN of the final structures across different quench rates, as shown in Fig. 6. Select CN distributions were also examined for finite size effects using additional simulations with 16 000 particles; as shown in Fig. S3, the results are insensitive to size.

Fig. 5 CN distributions with various IH distributions at ϕ = 0.20 under (a) step quench and (b) slow quench at dU/dt = −0.020k_BT s⁻¹ and ϕ = 0.20.

Fig. 6 Average CN as a function of quench rate for different IH conditions at ϕ = 0.20. Data at dU/dt = −10k_BT s⁻¹ correspond to step quench results.

The local configurational impact of IH was investigated further with TCC analysis. The results of the TCC analysis are shown in Fig. 7 for the same selected cases shown previously in Fig. 5. Note that the data in Fig. 7 has been grouped so that only cluster size, and not cluster geometry, is shown; the raw data which resolves distinct cluster geometries is shown in Fig. S4. Also note that each particle is only assigned to the largest cluster that it belongs to. Generally, IH is seen to induce a shift to larger TCC cluster sizes (∼10–13) relative to the no IH case, which is relatively enriched in intermediate sizes (∼5–8). This broad trend suggests that IH promotes further aggregation and consolidation of intermediate motifs into larger, more compact structures. Among the different types of IH, the uniform case (U) appears to be most effective at converting intermediate clusters to larger motifs (∼10–11), although interestingly, the Gaussian case (G) results in the highest number of size 13 clusters. This is notable because size 13 clusters may represent highly coordinated, low energy local packings. Overall, we find that interaction heterogeneity shifts the gel structure from intermediate-sized motifs toward larger, more compact topological clusters, with uniform distributions promoting general densification and Gaussian heterogeneity selectively stabilizing highly coordinated (size 13) structures. This suggests that IH enables continued structural evolution within the arrested gel.

Fig. 7 Grouped cluster size distribution (see text) by topological cluster classification (TCC) with various IH distributions at ϕ = 0.20 under (a) step quench and (b) slow quench at dU/dt = −0.020k_BT s⁻¹ and ϕ = 0.20. For particles identified in multiple topological clusters, only the largest cluster size is counted.

PCA was performed to further analyse and classify gelation outcomes based on CN distributions. We applied PCA on the first through tenth moments of the CN distributions, according to the methodology outlined in Section 2.2. As shown in Fig. 8, the fraction of total variance explained (FVE) captured by the first two PCs is over 99%, indicating that the first two PCs successfully capture the characteristics of CN distributions. A two-dimensional projection of the gelation outcome data in the space defined by the first two principal components is shown in Fig. 9. The effects of quench rate and IH in the PCA space are both evident and distinguishable. The impact of quench rate is most easily discerned by considering the homogeneous cases with no IH, as represented by the filled symbols and the black dashed line. As the quench rate is decreased, the data is observed to shift in a straight line in the PC1 + PC2 direction.

Fig. 8 Fraction of variance explained (FVE) by each principal component.

Fig. 9 Two-dimensional projection resulting from PCA applied to statistical moments of the CN distributions. Symbols indicate distinct IH conditions and particle volume fractions (ϕ), while colors represent varying quench rates. Solid and dash lines are trendlines with the effect of quench rate and IH, respectively, as guide to the eye.

Note that the homogeneous data includes contributions from gelation runs at two volume fractions, ϕ = 0.2 (circles) and ϕ = 0.3 (squares). While there is a modest shift in the PCA space due to changes in the volume fraction, the overall trends are consistent. While not a primary focus of the present analysis, the effect of volume fraction on IH-mediated gelation is expected to be multifactorial. First and most obviously, particle diffusivity in the fluid generally decreases with increased volume fraction, leading to a decrease in the sorting distance available to the gel phase. On the other hand, increasing the volume fraction also implies that additional particles would be available for sorting for a given distance. Finally, the overall time available for sorting during the gelation process is itself a function of overall volume fraction.

The impact of IH is also well delineated in the (PC1, PC2) space (open symbols represent cases with IH). At a given quench rate, IH drives the data along the PC1–PC2 direction, perpendicularly to the quench rate trendline as shown by the colored solid lines in Fig. 9. The excursions in (PC1, PC2) space due to IH are consistent across quench rates as evidenced by the fact that the linear guides are roughly parallel. However, as expected, the impact of IH is limited at high quench rates due to short sorting length scales. Notably, all the locations corresponding to various types and extents of IH are inaccessible in homogeneous systems, confirming the notion that IH may provide a tuneable parameter for creating new gel structures. Finally, we note that PCA of selected Steinhardt structural order parameter distributions fail to capture the effect of IH and/or quench rate; see Fig. S5.

3.4. Influence of IH on gelation dynamics

Finally, we study the impact of IH on gelation dynamics using CN as a metric to quantify the evolution of gel structure as a function of time. Shown in Fig. 10 is the evolution of the average CN as a function of interaction strength at a fixed quench rate of dU/dt = −0.020 k_BT s⁻¹ under different IH conditions. In general, IH is seen to broaden the gelation process—the average CN begins to increase at lower interaction strengths and reaches its final value at a higher interaction strength. Recall from Fig. 6 that the final value of the average CN is relatively unaffected by IH so the impact of IH is limited to the transient evolution of average CN. The observed effect can be broadly attributed to early aggregation driven by strong binders followed by late attachment of weak binders; together these extend the duration of the overall gelation process.

Fig. 10 Evolution of the average CN during gelation as a function of reference potential well-depth, U, at quench rate of dU/dt = −0.020k_BT s⁻¹and ϕ = 0.20 for various IH types.

Shown in Fig. 11 are example particle configurations during the gelation process under different IH conditions. Particles with CN ≥ 6 and their neighbours are identified as solid-like and coloured by their binding modulators. The remaining particles are identified as fluid-like and represented by smaller yellow spheres. Without IH, the gelation process occurs between interaction strengths U = −3.5k_BT and U = −6k_BT, resulting in complete particle aggregation by U = −6k_BT. Conversely, systems with IH exhibit a broader interaction strength range for gelation process, reflecting more gradual gelation behaviour, similar to our previous findings in the case of crystal nucleation.⁵³

Fig. 11 Configurations at selected U (in units of k_BT) during quenching at a rate of dU/dt = −0.020k_BT s⁻¹ under varying IH conditions. Potential nuclei (CN ≥ 6) and their bound particles are coloured with binding modulator b. Particles in No IH case are distinctly colored in purple. Fluid particles are depicted as smaller yellow spheres. (Parameters for different IH types: Gaussian: σ_IH = 0.15; Bidisperse: b_s/w = 1.50 and p(b_s) = 0.50; Uniform: Δb = 0.50).

The bidisperse IH case is particularly useful for additional analysis of the impact of IH on gelation. The bidisperse configuration at U = −3.5k_BT in Fig. 11 shows that gelled particles are almost entirely composed of strong binders despite the system containing equal proportions of strong and weak binders. Only later configurations (U = −6.0k_BT and later) include a significant number of weak binders in the gel phase. Although somewhat different in terms of the overall interaction matrix across particle types, this sequential gelation behaviour is qualitatively similar to that reported experimentally in ref. 15 A quantitative representation of the sequential gelation effect is shown in Fig. 12, which shows the initial increase in average CN for strong binders at approximately U ∼ −3.0k_BT (dotted line), followed by subsequent increase in average CN for the weak binders (dashed line) at U ∼ −5.0k_BT. While strong binders continue to increase their average CN during weak binder aggregation, the rate at which they do so is much lower than the initial burst. The difference in final average CN between strong and weak binders is significant—strong binders reach an average CN ∼ 9, while weak binders achieve a much lower average CN ∼ 4.5. The total system average CN (solid line) therefore exhibits a distinct two-step behaviour, plateauing at an intermediate average CN ∼ 6.5. For the same system parameters reported in the caption of Fig. 12, representative gel configurations at two points, U = −5k_BT and U = −12k_BT, during the quench are shown in Fig. 13. Here, the first column shows the overall process in which all particles are visible, while the second and third columns show only the strong (red) and weak (blue) binders, respectively. The entire gelation process is shown in an animation in Movie S1 (SI).

Fig. 12 Average CN as a function of average interaction strength for a quench rate of dU/dt = −0.020k_BT s⁻¹ and ϕ = 0.20 with bidisperse IH (b_s/w = 1.50 and p(b_s) = 0.50). Solid line – all particles, dash line – weak particles, and dotted line – strong particles.

Fig. 13 Configurations of gel structures for bidisperse IH with b_s/w = 1.50 and p(b_s) = 0.50, shown at (a)–(c) U = −5k_BT and (d)–(f) U = −12k_BT (final structure). Panels show (a, d) all particles, (b, e) strong-binding species only, and (c) and (f) weak-binding species only. (Red: strong-binding species. Blue: weak-binding species.) An animation of the complete gelation process is provided in Movie S1.

4. Conclusions

Unlike colloidal crystal nucleation and growth which tends to occur under weakly non-equilibrium conditions, gelation is most often a highly non-equilibrium process that is rather difficult to control. The most common control parameter is often the rate at which the driving force for aggregation is increased, usually imposed by temperature or concentration change, or by the modulation of an external field. While this can alter some broad structural features of the resulting gel, it is still rather limited. For example, faster thermal quenches typically reduce the gelation time and lead to less compact, more branched gel structures but the distribution of particle coordination remains relatively unchanged.

In this study, we identify and explore interaction heterogeneity (IH) as a new modality for controlling gelation in addition to thermal quench rates. We systematically investigated the influence of IH on colloidal gelation, focusing on the combined effects of varying quench rates and several types of IH distributions, including monodisperse (no IH), Gaussian, bidisperse, and uniform distributions.

At the overall system level, the gel structure was found to be only weakly influenced by IH. Analyses of the radial distribution function, g(r), and static structure factor, S(k), show subtle differences with the presence of IH, indicating minor changes in large-scale structural organization. Although fractal dimension measurements at intermediate quench rates suggest a potential influence of IH, this trend diminished in larger systems, indicating that the effect may not be robust.

In contrast, IH exhibits a profound impact on local gel structure. Broadly speaking, these effects originate from the tendency of IH to spread out the aggregation process in time, analogously to our previously reported observations of IH-mediated colloidal crystallization. More specifically for gelation, coordination number distributions are significantly altered by the specific type of IH under slow quench conditions. Most visibly, certain types of IH were found to generate bimodal coordination number distributions, which are otherwise inaccessible in monodisperse systems. Principal component analysis of coordination number distribution moments effectively distinguish structural outcomes due to quench rate variation and IH, underscoring the unique role IH plays in determining gelation pathways. Topological cluster classification (TCC) analysis further supports this picture, showing that IH enhances the formation of connected motifs in the gel phase. TCC also confirms that the specific type of IH is important in establishing the resulting particle network.

Finally, dynamical analysis of the impact of IH on gelation using the relatively simple bidisperse IH system shows that gelation is discretized into two sub-processes by IH—strong binders aggregate first to form a gel backbone, followed by decoration of the backbone with weak binders. Importantly, strong and weak binders exhibit markedly different coordination number environments, suggesting a pathway for assembling gels with tailored properties.

A potential application of using IH as a control parameter for colloidal gels is the fabrication of amorphous photonic crystals (APCs), which exhibit short-ranged order structure while maintaining macroscopic disorder. Unlike photonic crystals, APCs function as isotropic optical media, meaning that their photonic band gaps are independent of direction.^89,90 Our findings suggest that IH can be exploited to sequentially control the gelation process by the distinct behaviours among species with different binding modulators. A potential pathway to APC formation involves initiating local assembly through the early nucleation of strong binding species, resulting in short-ranged order domains. Subsequently, weak binding species can fill the voids between these small clusters, forming a disordered macroscopic gel structure. At least in concept, therefore, IH may offer a promising strategy to direct colloidal gelation toward structures characteristic of APCs.

Overall, our work demonstrates that IH modulates the sequence and kinetics of particle binding during gelation, enabling control over local structural features without significantly altering macroscopic organization. These insights suggest IH as a key design parameter for tuning the structure and functionality of colloidal gels, particularly for applications in isotropic photonic materials. Establishing detailed quantitative connections to experiment with particle-resolved techniques, such as 3D confocal microscopy, along with reliable measures of IH, will be required to confirm the mechanistic predictions of this work and enable systematic design of gelation protocols for desired optical responses.

Author contributions

Po-Ting Wu: Data curation (lead); formal analysis (equal); investigation (equal); methodology (equal); writing – original draft (equal); writing – review & editing (equal). Ying-Shuo Peng: formal analysis (equal); methodology (equal); writing – review & editing (equal). John C. Crocker: conceptualization (equal); formal analysis (equal); writing – review & editing (equal). Talid Sinno: Conceptualization (equal); Formal analysis (equal); Investigation (equal); methodology (equal); writing – original draft (equal); writing – review & editing (equal).

Conflicts of interest

The authors have no conflicts to disclose.

Data availability

Raw data for this article, consisting of particle coordinates and binding modulators for the various runs described in this study, are available at Scholarly Commons at https://repository.upenn.edu/handle/20.500.14332/13.

Supplementary information (SI) is available. See DOI: https://doi.org/10.1039/d6sm00117c.

Acknowledgements

This work was partially supported by the Department of Chemical and Biomolecular Engineering at the University of Pennsylvania.

References

1. E. Dickinson, Soft Matter, 2006, 2, 642–652.
2. R. Mezzenga, P. Schurtenberger, A. Burbidge and M. Michel, Nat. Mater., 2005, 4, 729–740.
3. L. Xu, Y. Zhou and S. Amin, Polym. Colloids, 2019, 9, 399.
4. N. Devi, M. Sarmah, B. Khatun and T. K. Maji, Adv. Colloid Interface Sci., 2017, 239, 136–145.
5. A. Hensley, W. M. Jacobs and W. B. Rogers, Proc. Natl. Acad. Sci. U. S. A., 2021, 119(1), e2114050118.
6. W. B. Rogers, W. M. Shih and V. N. Manoharan, Nat. Rev. Mater., 2016, 1, 1–14.
7. A. van Blaaderen, R. Ruel and P. Wiltzius, Nature, 1997, 385, 321–324.
8. A. D. Dinsmore, J. C. Crocker and A. G. Yodh, Curr. Opin. Colloid Interface Sci., 1998, 3, 5–11.
9. U. Gasser, E. R. Weeks, A. Schofield, P. N. Pusey and D. A. Weitz, Science, 2001, 292, 258–262.
10. M. E. Leunissen, C. G. Christova, A.-P. Hynninen, C. P. Royall, A. I. Campbell, A. Imhof, M. Dijkstra, R. van Roij and A. van Blaaderen, Nature, 2005, 437, 235–240.
11. A. J. Kim, P. L. Biancaniello and J. C. Crocker, Langmuir, 2006, 22, 1991–2001.
12. Y. Wang, Y. Wang, X. Zheng, É. Ducrot, J. S. Yodh, M. Weck and D. J. Pine, Nat. Commun., 2015, 6, 7253.
13. E. Zaccarelli, J. Phys.: Condens. Matter, 2007, 19, 323101.
14. P. J. Lu, E. Zaccarelli, F. Ciulla, A. B. Schofield, F. Sciortino and D. A. Weitz, Nature, 2008, 453, 499–503.
15. L. Di Michele, F. Varrato, J. Kotar, S. H. Nathan, G. Foffi and E. Eiser, Nat. Commun., 2013, 4, 2007.
16. C. P. Royall, S. R. Williams and H. Tanaka, J. Chem. Phys., 2018, 148, 044501.
17. L. J. Kaufman and D. A. Weitz, J. Chem. Phys., 2006, 125, 074716.
18. P. L. Biancaniello, A. J. Kim and J. C. Crocker, Phys. Rev. Lett., 2005, 94, 058302.
19. W. B. Rogers and J. C. Crocker, Proc. Natl. Acad. Sci. U. S. A., 2011, 108, 15687–15692.
20. W. Benjamin Rogers, T. Sinno and J. C. Crocker, Soft Matter, 2013, 9, 6412–6417.
21. R. Dreyfus, M. E. Leunissen, R. Sha, A. V. Tkachenko, N. C. Seeman, D. J. Pine and P. M. Chaikin, Phys. Rev. Lett., 2009, 102, 048301.
22. R. Dreyfus, M. E. Leunissen, R. Sha, A. Tkachenko, N. C. Seeman, D. J. Pine and P. M. Chaikin, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2010, 81, 041404.
23. S. Angioletti-Uberti, B. M. Mognetti and D. Frenkel, Phys. Chem. Chem. Phys., 2016, 18, 6373–6393.
24. P. Varilly, S. Angioletti-Uberti, B. M. Mognetti and D. Frenkel, J. Chem. Phys., 2012, 137, 094108.
25. D. Nykypanchuk, M. M. Maye, D. van der Lelie and O. Gang, Nature, 2008, 451, 549–552.
26. S. Y. Park, A. K. R. Lytton-Jean, B. Lee, S. Weigand, G. C. Schatz and C. A. Mirkin, Nature, 2008, 451, 553–556.
27. A. J. Kim, R. Scarlett, P. L. Biancaniello, T. Sinno and J. C. Crocker, Nat. Mater., 2009, 8, 52–55.
28. D. Sun and O. Gang, J. Am. Chem. Soc., 2011, 133, 5252–5254.
29. R. J. Macfarlane, B. Lee, M. R. Jones, N. Harris, G. C. Schatz and C. A. Mirkin, Science, 2011, 334, 204–208.
30. M. T. Casey, R. T. Scarlett, W. B. Rogers, I. Jenkins, T. Sinno and J. C. Crocker, Nat. Commun., 2012, 3, 1209.
31. W. B. Rogers and V. N. Manoharan, Science, 2015, 347, 639–642.
32. F. Lu, K. G. Yager, Y. Zhang, H. Xin and O. Gang, Nat. Commun., 2015, 6, 6912.
33. W. Liu, M. Tagawa, H. L. Xin, T. Wang, H. Emamy, H. Li, K. G. Yager, F. W. Starr, A. V. Tkachenko and O. Gang, Science, 2016, 351, 582–586.
34. Y. Wang, I. C. Jenkins, J. T. McGinley, T. Sinno and J. C. Crocker, Nat. Commun., 2017, 8, 14173.
35. É. Ducrot, M. He, G.-R. Yi and D. J. Pine, Nat. Mater., 2017, 16, 652–657.
36. H. Lin, S. Lee, L. Sun, M. Spellings, M. Engel, S. C. Glotzer and C. A. Mirkin, Science, 2017, 355, 931–935.
37. M. Song, Y. Ding, H. Zerze, M. A. Snyder and J. Mittal, Langmuir, 2018, 34, 991–998.
38. J. Harting, M. Hecht, H. J. Herrmann and S. McNamara, Multifield Problems in Solid and Fluid Mechanics, Springer Berlin Heidelberg, 2006, pp. 113–143.
39. B. Dünweg and A. J. Ladd, Lattice Boltzmann simulations of soft matter systems, in Advanced computer simulation approaches for soft matter sciences III, Springer Berlin Heidelberg, Berlin, Heidelberg, 2009, pp. 89–166.
40. R. G. M. van der Sman, Soft Matter, 2009, 5, 4376–4387.
41. D. S. Bolintineanu, G. S. Grest, J. B. Lechman, F. Pierce, S. J. Plimpton and P. R. Schunk, Comput. Part. Mech., 2014, 1, 321–356.
42. S. Poblete, A. Wysocki, G. Gompper and R. G. Winkler, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2014, 90, 033314.
43. Y.-S. Peng, Y.-J. Sheng and H.-K. Tsao, Soft Matter, 2020, 16, 5054–5061.
44. Y.-S. Peng and T. Sinno, J. Chem. Phys., 2024, 160, 174121.
45. P. G. Bolhuis and D. A. Kofke, Phys. Rev. E: Stat. Phys., Plasmas, Fluids, Relat. Interdiscip. Top., 1996, 54, 634–643.
46. S. Henderson, T. C. Mortensen, S. Underwood and W. Megen, Phys. A, 1996, 233, 102–116.
47. S. Auer and D. Frenkel, Nature, 2001, 413, 711.
48. P. Sollich and N. B. Wilding, Phys. Rev. Lett., 2010, 104, 118302.
49. P. Sollich and N. B. Wilding, Soft Matter, 2011, 7, 4472–4484.
50. R. Rice, R. Roth and C. P. Royall, Soft Matter, 2012, 8, 1163–1167.
51. Y. Wang, Y. Wang, D. R. Breed, V. N. Manoharan, L. Feng, A. D. Hollingsworth, M. Weck and D. J. Pine, Nature, 2012, 491, 51–55.
52. F. Sciortino and E. Zaccarelli, Curr. Opin. Colloid Interface Sci., 2017, 30, 90–96.
53. I. C. Jenkins, J. C. Crocker and T. Sinno, Phys. Rev. Lett., 2017, 119, 178002.
54. P.-T. Wu, J. C. Crocker and T. Sinno, J. Chem. Phys., 2024, 161, 244902.
55. G. H. Lee, J. B. Kim, T. M. Choi, J. M. Lee and S.-H. Kim, Small, 2019, 15, e1804548.
56. Z. Cai, Z. Li, S. Ravaine, M. He, Y. Song, Y. Yin, H. Zheng, J. Teng and A. Zhang, Chem. Soc. Rev., 2021, 50, 5898–5951.
57. T. Liu, B. VanSaders, J. T. Keating, S. C. Glotzer and M. J. Solomon, Langmuir, 2021, 37, 13300–13308.
58. C. P. Royall, M. A. Faers, S. L. Fussell and J. E. Hallett, J. Phys.: Condens. Matter, 2021, 33, 453002.
59. E. Bianchi, J. Largo, P. Tartaglia, E. Zaccarelli and F. Sciortino, Phys. Rev. Lett., 2006, 97, 168301.
60. J. M. Kim, J. Fang, A. P. R. Eberle, R. Castañeda-Priego and N. J. Wagner, Phys. Rev. Lett., 2013, 110, 208302.
61. S. Zhang, L. Zhang, M. Bouzid, D. Z. Rocklin, E. Del Gado and X. Mao, Phys. Rev. Lett., 2019, 123, 058001.
62. F. Bonacci, X. Chateau, E. M. Furst, J. Fusier, J. Goyon and A. Lemaître, Nat. Mater., 2020, 19, 775–780.
63. M. Bantawa, B. Keshavarz, M. Geri, M. Bouzid, T. Divoux, G. H. McKinley and E. Del Gado, Nat. Phys., 2023, 19, 1178–1184.
64. C. P. Royall and A. Malins, Faraday Discuss., 2012, 158, 301–311 ; discussion 351-70.
65. B. K. Ryu, S. M. Fenton, T. T. D. Nguyen, M. E. Helgeson and R. N. Zia, J. Chem. Phys., 2022, 156, 224101.
66. A. D. Smith, G. J. Donley, E. Del Gado and V. M. Zavala, ACS Nano, 2024, 18, 28622–28635.
67. M. Bouzid and E. Del Gado, Langmuir, 2018, 34, 773–781.
68. H. Bhaumik, T. B. Liverpool, C. P. Royall and R. L. Jack, Phys. Rev. E., 2025, 111, 055412.
69. S. M. Fenton, P. Padmanabhan, B. K. Ryu, T. T. D. Nguyen, R. N. Zia and M. E. Helgeson, Proc. Natl. Acad. Sci. U. S. A., 2023, 120, e2215922120.
70. A. P. Thompson, H. M. Aktulga, R. Berger, D. S. Bolintineanu, W. M. Brown, P. S. Crozier, P. J. in ’t Veld, A. Kohlmeyer, S. G. Moore, T. D. Nguyen, R. Shan, M. J. Stevens, J. Tranchida, C. Trott and S. J. Plimpton, Comput. Phys. Commun., 2022, 271, 108171.
71. J. SantaLucia Jr, Proc. Natl. Acad. Sci. U. S. A., 1998, 95, 1460–1465.
72. V. Ramasubramani, B. D. Dice, E. S. Harper, M. P. Spellings, J. A. Anderson and S. C. Glotzer, Comput. Phys. Commun., 2020, 254, 107275.
73. D. Frenkel and B. Smit, Understanding Molecular Simulation: From Algorithms to Applications, Elsevier, 2023.
74. E. Fransson, M. Slabanja, P. Erhart and G. Wahnström, Adv. Theory Simul., 2021, 4, 2000240.
75. S. Buczkowski, S. Kyriacos, F. Nekka and L. Cartilier, Pattern Recognit., 1998, 31, 411–418.
76. S. Griffiths, F. Turci and C. P. Royall, J. Chem. Phys., 2017, 146, 014905.
77. S. R. Williams, arXiv, 2007, preprint, arXiv:0705.0203.
78. A. Malins, S. R. Williams, J. Eggers and C. P. Royall, J. Chem. Phys., 2013, 139, 234506.
79. C. P. Royall, J. Eggers, A. Furukawa and H. Tanaka, Phys. Rev. Lett., 2015, 114, 258302.
80. F. Pedregosa, G. Varoquaux, A. Gramfort, V. Michel, B. Thirion, O. Grisel and E. Duchesnay, J. Mach. Learn. Res., 2011, 12, 2825–2830.
81. J. W. Rohlf and P. J. Collings, Phys. Today, 1994, 47, 62–63.
82. D. A. Weitz, J. S. Huang, M. Y. Lin and J. Sung, Phys. Rev. Lett., 1985, 54, 1416–1419.
83. M. Mellema, J. H. J. van Opheusden and T. van Vliet, J. Chem. Phys., 1999, 111, 6129–6135.
84. W. R. Heinson, C. M. Sorensen and A. Chakrabarti, J. Colloid Interface Sci., 2012, 375, 65–69.
85. P. J. Lu and D. A. Weitz, Annu. Rev. Condens. Matter Phys., 2013, 4(1), 217–233.
86. S. Lazzari, L. Nicoud, B. Jaquet, M. Lattuada and M. Morbidelli, Adv. Colloid Interface Sci., 2016, 235, 1–13.
87. S. Jungblut, J.-O. Joswig and A. Eychmüller, Phys. Chem. Chem. Phys., 2019, 21, 5723–5729.
88. J. Opdam, M. Tateno and H. Tanaka, ACS Nano, 2025, 19(23), 21515–21524.
89. L. Shi, Y. Zhang, B. Dong, T. Zhan, X. Liu and J. Zi, Adv. Mater., 2013, 25, 5314–5320.
90. L. Bai, V. C. Mai, Y. Lim, S. Hou, H. Möhwald and H. Duan, Adv. Mater., 2018, 30(9), 1705667.

---