Colloidal particle aggregation: mechanism of assembly studied via constructal theory modeling

The assembly of colloidal particles into ordered structures is of great importance to a variety of nanoscale applications where the precise control and placement of particles is essential. A fundamental understanding of this assembly mechanism is necessary to not only predict, but also to tune the desired properties of a given system. Here, we use constructal theory to develop a theoretical model to explain this mechanism with respect to van der Waals and double layer interactions. Preliminary results show that the particle aggregation behavior depends on the initial lattice configuration and solvent properties. Ultimately, our model provides the first constructal framework for predicting the self-assembly of particles and could be expanded upon to fit a range of colloidal systems.


Introduction
Constructal theory has been used to describe a number of naturally evolving processes/phenomena that include, but are not limited to, turbulent flow, heat and mass transfer, dendritic formation, and biological growth [1][2][3][4][5][6][7][8][9][10]. According to this law: "For a flow system to persist in time (to live), it must evolve in such a way that it provides easier access to the imposed (global) currents that flow through it" [1]. The main distinction between the classical laws of thermodynamics and constructal law is the former only describes the final equilibrium state of the system while the latter describes how and through what pathway(s) the system naturally progresses to achieve that final equilibrium state.
In concurrence with such constructal analyses, Bejan and Wagstaff showed that the natural coalescence of masses in space, with respect to attractive gravitational forces, will lead to a hierarchical distribution [11]. The internal tension of the system, which drives this natural coalescence, is decreased faster and more efficiently when the masses coalesce nonuniformly toward a state of equilibrium [10]. Similarly, Reis et al. used these principles to analyze how the electrostatic interactions between dust particles are minimized during particle agglomeration in air-cleaning devices [12]. It was concluded that the maximum collection rate was achieved when the aerosol particles first formed into clusters and then dendrites in order to balance their electrostatic forces. The aggregation of colloidal particles (small particles or droplets typically between 10 nm and 10 µm in diameter that can be suspended in a fluid via Brownian motion [13]) is, therefore, expected to adhere to a similar natural progression such that the total energy of the system is minimized. When analyzing colloidal systems, however, both attractive and repulsive forces must be taken into consideration.
The pairwise interactions between these particles are often described by traditional colloidal Derjaguin-Landau-Verwey-Overbeek (DLVO) theory [14]. Here, electrolyte ions in solution rearrange to form electric double layers near the charged surface of a particle. Although overlapping double layers result in repulsion between two particles, this force is constantly opposed by the attractive van der Waals force. The balance between these interparticle forces gives the total DLVO force and highly depends on system parameters, such as the electrolyte concentration and fluid dielectric constant, as well as the particle size, spacing, and surface charge. As a result, these system parameters can be tuned in order to control the balance between the DLVO forces.
Tunability, or the ability to direct the movement and assembly of particles into higher-order structures, is important for a wide variety of applications, including plasmonic/photonic devices, metamaterials, biomaterials and coatings, energetic materials, and high-strength ceramics [15][16][17][18][19][20][21][22][23][24][25]. Most applications, however, require precise manipulation of the particle placement, spacing, and arrangement in order to achieve desired macroscale properties. Thus, a theoretical understanding of the particle assembly mechanism is necessary to predict the nano/ microscale structures responsible for the desired or observable properties.
To our knowledge, constructal theory has never been applied/ implemented to describe the assembly of colloidal particles. In this article, we use said theory to predict the effect of varying system parameters on the magnitude of the total DLVO force for two scenarios, that is, uniform and non-uniform aggregation. The constructal approach presented here provides a fundamental framework for understanding and tuning the assembly behavior of colloidal particles, which could have implications across a broad range of fabrication techniques.

Model Details
First, we consider the pairwise DLVO interactions between two isolated spheres suspended in an electrolyte solution. The electrolyte ions arrange near the surface of the particles to form an electric double layer, thus screening the surface charge. The characteristic length or "thickness" of this double layer (which is a function of the ion concentration, I) is known as the Debye length, λ D , while the surface charge is parameterized by the zeta potential, ζ p . This repulsive double layer force is given by [14]: (1) where d is the particle separation distance, a is the particle radius, T is the temperature of the system, e is the elementary charge, k B is Boltzmann's constant, and ε 0 and ε c are vacuum permittivity and fluid dielectric constant, respectively. Conversely, the attractive van der Waals force is given by [14]: (2) where the characteristic energy scale is set by the Hamaker constant, A. It is noted that Equation 1 and Equation 2 assume spherical particles of equal radius and a sufficiently small separation distance (a ≫ d). A simplified schematic of this pairwise particle interaction is shown in Figure 1a. The resulting total DLVO force is then written as the sum of the double layer repulsion and the van der Waals attraction: When F DLVO < 0 or F DLVO > 0, negatively charged particles will experience a net attractive or repulsive force, respectively ( Figure 1b). Therefore, when F DLVO = 0, the particles experience zero net force. These fixed point locations are observed in Figure 1b where F DLVO intersects the equilibrium line. The location of the stable fixed point, which varies based on the particle and solvent properties, determines the preferred interparticle separation distance. Ultimately, this DLVO framework is used to describe the structural evolution of colloidal particle assemblies.

Aggregation in one dimension
Similar to the analysis by Bejan and Wagstaff [11], we begin by examining the one-dimensional (1D) assembly of particles. Initially, the particles are set at an evenly spaced distance (d) along an infinite line ( Figure 2a). In this initial state each particle feels a net DLVO force, F 11 , in either direction due to the forces from all other neighboring particles along the line. Therefore, based on Equation 2, the total sum of the van der Waals contributions is: Likewise, from Equation 1, the total double layer contribution is: As the system begins to evolve from this initial state, the particles can assemble either via uniform or non-uniform aggregation. With respect to uniform aggregation of identical particles, two neighboring particles uniformly come together to form doublets with a new spacing of 2d ( Figure 2b). Each particle separately feels the forces from the individual particles of the neighboring doublet, which results in an increase by a factor of four in the total force, where: (8) and (9) It is noted that the van der Waals force for uniform aggregation is equal to that of the initial particle configuration (i.e., F 11;vdW = F 12;vdW ). Therefore, the magnitude of F 12 compared to F 11 depends only on the double layer term. Alternatively, with respect to a non-uniform aggregation of identical particles, two particles combine with a single neighbor (Figure 2c), which results in a cluster of three particles with one lone particle remaining. Again, the spacing between each singlet/triplet pair is 2d. The total force, F 13 , on a single particle is, therefore, the sum of the contributions from alternating clusters and singlets, where: (b) The resulting uniform aggregation from (a) leads to a new square lattice of particle doublets with a uniform spacing of . (c) Particles in a uniform square lattice undergo non-uniform aggregation when four neighboring particles aggregate about a center particle (green arrows). (d) The resulting non-uniform aggregation from (c) leads to a new lattice of five-particle clusters surrounded by singlets. A single particle at the origin feels the net force from each infinite string of clusters and singlets. (10) and (11) Compared to the initial particle configuration, it can be deduced from Equation 10 and Equation 11 that the van der Waals effect decreases for non-uniform aggregation (F 13;vdW /F 11;vdW = 0.625), while the double layer effect exhibits a new exponential dependence from the formation of the particle clusters. Equations 5-11 can now be expanded in order to describe particle assembly in two dimensions.

Aggregation in two dimensions
Next, we consider two variations of an infinite, two-dimensional (2D) array of particles, that is, a uniform square lattice and a hexagonal (triangular) lattice. First, we focus on the uniform square lattice (Figure 3a). The net DLVO force, F 2 , on a single particle positioned at the origin (O) results from the sum of each infinite string of particles acting in the positive x-direction. Symmetry about the x-axis causes the horizontal components of the force to be additive while the upward (y > 0) and downward (y < 0) vertical components of the force will cancel each other out. For example, the x-component of the van der Waals force along the 45° line in Figure 3a is with a total contribution of 2F 21x;vdW . Thus, F 2 is given by the sum of all van der Waals and double layer contributions in the horizontal direction (x > 0): The first summation on the right-hand side of Equation 13 and Equation 14 accounts for the force contributions for F 20 < θ < F 21 , while the second summation in each equation accounts for the force contributions for θ > F 21 , where θ is the angle of the force vector.
From this basic square lattice arrangement, the particles can again assemble by either uniform or non-uniform aggregation. When neighboring particles combine via uniform aggregation (blue arrows in Figure 3a), the result is a new square lattice of particle doublets with a uniform spacing of ( Figure 3b) and an increase by a factor of four in the total force (F = 4F 2 ( )). For clarification, the pathway outlined by the blue arrows in Figure 3a is just one example of several possible uniform assembly pathways. Contrarily, when four neighboring particles combine with a center particle via non-uniform aggregation (green arrows in Figure 3c), a new lattice of five-particle clusters surrounded by singlets is formed (Figure 3d). Here, a single particle at the origin feels a net force, F 3 , from each infinite string of clusters and singlets, given by: where B = [(4n) 2 + 4], and C = [3 n + (−1) n ]. Similar to the result for F 2 , the first and second summations on the right-hand side of Equation 16 and Equation 17 account for the force contributions from F 30 < θ < F 31 , while the third and fourth summations in each equation account for the force contributions from θ > F 31 .
The approach outlined above can also be used for the analysis of the hexagonal lattice structure (Figure 4a) where the net DLVO force, F 4 , is again given by the sum of all van der Waals and double layer contributions in the horizontal direction (x > 0): The resulting uniform aggregation from (a) leads to a new hexagonal lattice of seven-particle clusters with a uniform spacing of . (c) Particles in a hexagonal lattice undergo non-uniform aggregation when non-neighboring particle aggregates form about a center particle (green arrows). (d) The resulting non-uniform aggregation from (c) leads to a new lattice of seven-particle clusters surrounded by single particles. A single particle at the origin feels the net force from each infinite string of clusters and singlets.
The first summation on the right-hand side of Equation 19 and Equation 20 accounts for the outer force contributions where θ < F 43 and θ > F 42 , while the second summation in each equation accounts for the more central force contributions where F 43 < θ < F 42 .
Particles in this hexagonal lattice configuration undergo uniform assembly when neighboring particles aggregate about a shared, central particle (blue arrows in Figure 4a). A new uniform hexagonal lattice is formed with clusters of seven particles separated by a triangular spacing of (Figure 4b), while the total force of the system increases by a factor of 49 (F = 49F 4 ( )). Alternatively, non-uniform aggregation from the hexagonal lattice occurs when non-neighboring particle aggregates form (green arrows in Figure 4c) leading to a new lattice of seven-particle clusters surrounded by singlets (Figure 4d). The net force on a single particle at the origin, F 5 , is given by: Equations 12-23 can now be used to predict the effects of varying system parameters on the particle assembly behavior. It should be noted that Figures 2-4 are merely snapshots of a transient system -not the final state of the system -, which will continue to evolve and form new aggregates until equilibrium is reached. The approaches outlined above are recursive and can continue to be carried out as the particle spacing and aggregate size changes.

Results and Discussion
The constructal model outlined above provides a useful framework for predicting, tuning, and optimizing colloidal particle assembly. As a proof of concept, we start with a representative particle/solvent combination, where a 0 = 25 nm, ζ p = −25 mV, I = 0.011 M, and ε c = 80 (water), with an initial particle lattice spacing of d = 5 nm. The effect a given system parameter has on the overall mechanism of assembly with respect to lattice configuration can now be evaluated by simply adjusting the desired parameter while holding the others constant.
The predicted results of the square and hexagonal lattice arrangements are shown in Figure 5 and Figure 6, respectively. The total DLVO forces from both uniform (left) and nonuniform (right) aggregation are plotted as functions of particle spacing or particle radius and electrolyte concentration or fluid dielectric constant. The same general trends were observed for both lattice configurations and are describe herein. Figure 5a,b and Figure 6a,b show the results of varying particle spacing as a function of the electrolyte concentration. Here, the magnitude of the total DLVO force is represented by the color intensity while the color itself represents whether attractive or repulsive forces dominate. Red indicates areas where the repulsive double layer force dominates while blue indicates areas where the attractive van der Waals force dominates. When the double layer and van der Waals forces are equal, the total DLVO force is zero leading to a stable fixed point (cf. Figure 1b). The dashed lines in Figure 5 and Figure 6 denote the locations of these stable fixed points. As expected, when the particles are moved closer together the magnitude of the DLVO force increases. Likewise, as the electrolyte concentration is increased (for both uniform and non-uniform aggregation) the particles become less repulsive and transition to an attractive regime where aggregation is favorable. Figure 5c,d and Figure 6c,d show a similar result for the effect of particle spacing as a function of the fluid dielectric constant such that lower dielectric solvents lead to an increased attractive force and particle aggregation. The effect of the particle radius on the assembly behavior in each lattice arrangement, with respect to the solvent properties, is shown in Figure 5e-h and Figure 6e-h. As expected from Equation 1 and Equation 2, the size of the particle only influences the magnitude of the DLVO force, which scales linearly with particle radius, that is, larger particles are more attractive/ repulsive. Ultimately, only the solvent properties are responsible for determining whether or not aggregation will be favorable. Similar to varying the particle spacing, aggregation is most favorable at high salt concentrations (Figure 5e  where N A is Avogadro's number, and e is the elementary charge. A sufficiently thin double layer will allow the particles to come into closer contact where van der Waals attraction can occur. Consequently, increasing the electrolyte concentration or decreasing the fluid dielectric constant will decrease the Debye length and promote aggregation. The results in Figure 5 and Figure 6 reveal that the magnitude of the total DLVO force, regardless of nature (attractive or repulsive), is expected to be greater for uniform aggregation than for non-uniform aggregation. The nature of this force does, however, naturally influence particle behavior. For example, at a fluid dielectric constant value of ε c = 10, the total DLVO force is negative (Figure 5c,d,g,h and Figure 6c,d,g,h). Since the magnitude of this force is greater in the uniform scenario than in the non-uniform scenario, uniform aggregation is expected. For the square lattice, specifically, this increased magnitude is accompanied by a shift in the location of the stable fixed points, which suggests that uniform aggregation of particles is typically preferred with respect to van der Waals and double layer interactions. For comparison, a value of ε c = 80 yields a positive total DLVO force and separation is expected. Since the magnitude of this force is smaller for the non-uniform scenario, the system is "less repulsive" such that a non-uniform pathway could potentially be more conducive for aggregation.
Unlike the square lattice, where predicted particle movement (aggregation/separation) is the same regardless of the assembly mechanism (uniform/non-uniform), the hexagonal lattice shows regimes where the assembly mechanism seems to influence the particle behavior. For example, when the particles are in close proximity (d = 1 nm) at an electrolyte concentration of I = 0.08 M (Figure 6a,b) or a fluid dielectric constant value of ε c = 35 (Figure 6c,d), the uniform scenario leads to a net repulsive force. Conversely, for the non-uniform scenario these system parameters predict a net attractive force. Therefore, the solvent properties could be tuned by changing either the salt concentration or dielectric constant in order to induce a desired assembly mechanism. Moreover, the overall magnitudes of the DLVO forces in the hexagonal lattice ( Figure 6) are greater than the magnitudes of the forces in the square lattice ( Figure 5), that is, the driving force for aggregation (or separation) is greater for the hexagonal system. In the square lattice structure, a single particle is surrounded by four neighbors within one d-spacing, while a particle in the hexagonal lattice structure is surrounded by six neighbors within one d-spacing. Thus, the increase in driving force for the hexagonal system is likely due to this increase in the coordination number.
We should note that the results presented here assume the selfassembly of monodisperse spherical particles suspended in a monovalent salt solution arranged in a uniform starting configuration. Given the Brownian nature of colloidal particles, future simulations should account for the randomness of colloidal suspensions. These simulations should also include the effects from polydisperse particle mixtures or divalent salt solutions, and they should be able to handle the addition of an external stimulus that could facilitate and drive particle assembly.

Conclusion
We have developed the first model based on constructal theory to explain the self-assembly mechanisms of charged colloidal particles. This theoretical model, based on traditional DLVO interactions, can be broadly applied to many colloidal systems and is used here to predict the particle aggregation and equilibrium behavior as a function of varying system parameters. Our modeling results indicate that uniform aggregation of a square lattice is preferred over non-uniform aggregation with respect to the overall magnitude of the total DLVO force. When arranged in a hexagonal lattice, however, the preferred mechanism of assembly is determined by the solution properties and can be tuned by varying the salt concentration or fluid dielectric constant. Possible expansions to this theory include accounting for various particle mixtures, solvent properties, and particle configurations, as well as adding external driving forces (i.e., electric or magnetic fields) and, finally, a third dimension. Nonetheless, our model provides a fundamental framework for understanding and tuning the assembly behavior of colloidal particles and can have broad implications on fabrication techniques.