1. Introduction
The formation dynamics and structures of colloidal clusters have been extensively studied because of their important role in a variety of applications, such as sensors, plasmonics and photonic crystals [
1,
2,
3]. Furthermore, colloidal clusters that contain a small number of constituent colloids can serve as building blocks for hierarchically organized superstructures [
4,
5]. Manoharan et al. [
6] developed an assembling route based on the use of emulsion droplets as a three-dimensional (3D) template to obtain identical colloidal clusters with high yields. In this approach, colloidal microspheres absorbed at the interface of emulsion droplets reduce the interfacial free energy. During the subsequent evaporation of the oil droplets, colloidal particles are forced to pack into clusters by capillary forces and strongly bound together by van der Waals attractions.
Following the pioneering work of Manoharan and coworkers [
6], numerous studies based on emulsion droplet evaporation have been performed. Some of these focused on different types of materials/solvents [
7,
8,
9,
10], others on submicron-sized particles [
11,
12,
13] and others still on clusters with large sizes [
14,
15]. Interestingly, a common feature of the small colloidal clusters is that their geometric structures are in many cases unique for a given number of constituent colloids,
. The uniqueness of clusters for each
suggests that they are optimal packings and appear to be controlled by minimization of the second moment of the mass distribution in the cluster [
16,
17]. Manoharan et al. [
18,
19] examined the structure and free energy of colloidal clusters of hard spheres with a short-range attraction, albeit without geometrical confinement. The authors found that except in trivial cases with
, the clusters possess multiple configurations whose occupation probabilities are determined by geometrical rules. In addition, rotational entropy disfavors high-symmetry clusters (octahedron occurs less frequently than poly-tetrahedron for a cluster of six particles).
Despite the strong interest in 3D cluster structures, little work has been done to investigate the emulsion-assisted formation of two-dimensional (2D) colloidal clusters. Perry and co-workers [
20] studied, both experimentally and theoretically, transitions between ground states and excited states of 2D clusters of spherical colloids bound by depletion interactions. The ground state clusters adopt three distinct configurations that are dictated primarily by symmetry (permutational entropy). Because the number of ground and excited states increases rapidly with
, the authors considered only the case of clusters made of six particles. The experimental studies of Hurd and Schaefer [
21] and Oneda [
22] on colloidal suspensions, confined to two dimensions at the air–water interface, showed, depending on the strength of attractive forces, either an irreversible or reversible clustering process. However, the authors did not identify the observed cluster structures. Iwamatsu [
23] analyzed the lowest energy structure of 2D colloidal clusters as a function of
using a realistic model potential, which consists of a short-ranged van der Waals attractive interaction and a long range dipolar electrostatic repulsion. Although all of the predicted cluster configurations match exactly those of Lennard–Jones 2D clusters, they are not always consistent with experimental structures [
24,
25]; for example, for
the predicted cluster adopts a regular hexagonal structure, whereas the experimentally observed cluster forms a circular shell structure. Note that the particles self-assemble spontaneously into such clusters, which are just intermediate states during the formation/dynamics of larger assemblies, because there is no restriction on the size of a cluster driven by attractive forces.
Polyiamonds are plane figures formed by joining congruent equilateral triangles edge to edge. Polyiamonds are considered distinct (free polyiamonds) if they are different in shape—more specifically, they are invariant under reflection, translation and rotation. The number of free polyiamonds without holes as a function of
, where
is the number of equilateral triangles, is given by the Online Encyclopedia of Integer Sequences (OEIS): A070765 [
26]. Polyiamonds, one problem in combinatorial geometry, have applications in statistical mechanics of macromolecules [
27] and percolation [
28]. However, to our best knowledge, there is no report about the occupation probability of these polyiamonds of a given
, and a connection between polyiamonds and colloidal clusters.
In this work, similar in spirit to a model originally developed for 3D systems [
29,
30], we formulate a model of colloid-droplet mixtures in 2D. The structures and occurrence probabilities of small clusters of colloidal particles assembled in the spherical confinement of emulsion droplets are analyzed. We find that each cluster of
particles possesses a unique configuration for
, but multiple, degenerate configurations for
. We suggest that the known ground state clusters for systems of particles interacting through the weak depletion interaction [
20] are also relevant to colloidal clusters assisted by evaporating emulsion droplets. In particular, our cluster structures exactly correspond to those of polyiamonds.
This paper is organized as follows. In
Section 2 we introduce the model and simulation methodology a binary colloid-droplet mixture. Our results are presented in
Section 3. Conclusions are given in
Section 4. A detailed description of polyiamonds is given in
Appendix A.
2. Model and Simulation Method
We consider a binary mixture of
colloidal particles of diameter
and
droplets of diameter
in the plane. Similar to three dimensional (3D) systems that have been studied [
30], the colloids interact via a long-range screened Coulomb repulsion (Yukawa potential) and short range attraction, which can be expressed as
where
is the depth of the square-well potential and
with
the width of the attractive well. Two any colloids form a bond when their center-to-center separation is smaller than or equal to
. We chose the simulation parameters to have a sufficiently large attractive component of the potential function, so that physical bonds between colloids, once formed, would be irreversible, to mimic short-range van der Waals. The parameter
measures the strength of the Yukawa repulsion and the inverse shielding length
controls the range of the Yukawa interaction. Models of this type, as illustrated in
Figure 1a, have proven useful in studies of the formation of colloidal shells [
31] and the clustering of colloidal particles [
29,
30,
32,
33].
Emulsion droplets stabilized by colloidal particles, so-called Pickering emulsions, are considered to be highly stable, due to a significant reduction in the surface free energy [
34]. The loss of the interfacial energy, dependent on the particle size, surface tension and equilibrium three-phase contact angle, is typically of the order of hundreds or millions of
[
35]. To model this effect of the interfacial energy, Schwarz et al. [
29] first introduced a colloid-droplet energy function which is related to the surface tension
by
, where
S is the droplet surface that is covered by the colloid (details in [
29]). In an analogous fashion, for two-dimensional systems, we consider a line tension instead of the surface tension and a contact line instead of the contact surface. Thus, the colloid-droplet energy can be expressed for
, as
and
,
where
is the line tension and
l is the length of the contact line, given by
The interfacial energy for various droplet-colloid size ratios is pictured in
Figure 1b, from which it can be seen that the potential shape is very similar to that in the three-dimensional case [
29]. Moreover, for simplicity, we assume that the coalescence of the droplets is negligible. The droplet–droplet interaction has a hard-core repulsion of effective diameter
to avoid one colloid being shared by two droplets.
We performed kinetic Monte Carlo (MC) simulations of colloid-droplet mixtures in the canonical ensemble. The trial state was generated by giving the particles small maximum displacements,
for the colloidal particles and
for the droplets. This ensured that kinetic Monte Carlo simulations approximated Brownian dynamics simulations [
36]. The evaporation dynamics were introduced by forcing the droplets diameter
to shrink at a fixed rate. As a result, the resulting transient structures during evaporation were out-of-equilibrium. The evaporation rate was chosen so that all the droplets vanished after one-half of the simulation time (
MC sweeps). This left another
MC sweeps to investigate the (quasi-)stability of the clusters against thermal fluctuations. During the span of the MC simulations, a cluster, a series of colloids connected to each other by a path of bonds, can be formed spontaneously or formed via the evaporation of the emulsion droplets. Therefore, the collective modes of motion of particles in the cluster, i.e., collective translational and rotational cluster moves [
37,
38], are taken into account. In a MC cluster move (and its reverse), the clusters are rotated around a random rotation axis with a maximum angle
and translated with a maximum linear displacement
with
number of constituent colloids in the cluster and
the effective diameter of the cluster taken to be
. The choice of translational and rotational motions that satisfies the conditions for the translational diffusion constant (as discussed later) is to approximate a hydrodynamic damping of a circular cluster. Moreover, these cluster moves are proposed and accepted according to the principle of detailed balance; i.e., any cluster–cluster aggregation or cluster–particle aggregation is rejected [
29].
The simulations were performed in a cubic box with periodic boundary conditions for 500 colloidal particles with packing fraction
, and at a fixed droplet packing fraction
for a number of droplets of 6–12. For a given set of parameters, the statistical data of the sampled quantities were analyzed by running 30 independent simulations. To study the role of the interaction range on the cluster size distribution and cluster structure, we chose the width of the attractive well
in the range 0.09–0.21. Additionally, it is known that the cluster size distribution is driven by various factors, including (i) the droplet (colloid) packing fraction, (ii) the initial droplet size, (iii) the interaction parameter and (iv) the rotational and translational dynamics of clusters. Therefore, for a large parameter space, we restricted ourselves to simulation parameters (
,
,
and the interaction parameters given in the caption of
Figure 1) that provide a distribution of small-size clusters (
).
3. Result and Discussion
Here we roughly estimate the physical time per MC sweep via the diffusion coefficient because this quantity provides a link between the mean square displacement of the particles and the simulation time. The self-diffusion coefficient of clusters
is defined by the Einstein relationship [
39]:
where
n is the number of MC sweeps and
is the physical time per MC sweep. The mean square displacement of the clusters after
n sweeps,
, is given by [
30]
where
is the number of clusters containing
colloids, and
is the displacement of a cluster from its center-of-mass after
n sweeps.
For two-dimensional liquids, unlike three-dimensional liquids, the validity of the Stokes–Einstein relation that relates the diffusion coefficient of a Brownian particle and the fluid shear viscosity
over a wide range of density or temperature is still unclear [
40]. Hence, instead of using the two-dimensional Stokes–Einstein relation [
40],
, which is only valid at high densities [
41], we calculate the self-diffusion coefficient from the Enskog equation in the first Sonine approximation; that is [
42,
43],
Here
is the number density;
is the colloid mass; and
is the equilibrium pair correlation function at contact. The relevant mesoscopic time scale is the Brownian time
, which is the time required for an isolated cluster to diffuse over a distance equal to its diameter, where the diameter of the circular cluster
. An estimate of
in terms of
is, therefore, given by
with
At room temperature, for the clusters containing six colloidal particles and typical values of
, we found that
s and
s; those values are much smaller than the realistic long-time dynamics in experiments (usually the order of minutes or hours). Despite the difference in timescale between the simulation and experiments, we do not expect the timescale to affect the final cluster structures [
30].
Figure 2 shows representative snapshots of the simulation trajectory for the binary mixture of colloids and droplets. The initial configuration consists of non-overlapping of free colloids (purple) and droplets (pink) (
Figure 2a). After
MC sweeps, accompanied by the slow evaporation of the droplet phase, a statistical number of particles is strongly confined in the droplet surface (
Figure 2b). Since the droplet diameter
changes continuously during the simulation, the system is not at equilibrium. After this out-of-equilibrium assembly, a process towards equilibrium starts, which leads to the rearrangement of the particles to a cluster. In reality, stable clusters can be achieved only when the system is close to its equilibrium structure, a structure that forms on a very long timescale. The presence of the droplet-induced clusters as a final phase with respect to thermal fluctuations is shown in
Figure 2c at the end of the simulation (after
MC sweeps). As a test, we performed an additional run with
MC sweeps and observed no change in the number of resulting clusters and none in their structures (
Figure 2d).
Another feature is that the minimum separation between colloids at the initial stage of the simulations is larger than the bond length. In order to form a bond, the colloids must be provided with sufficient thermal energy to cross the repulsive energy barrier () of the colloid–colloid interaction. This probability is very small, and consequently, in the restricted simulation timescale, the colloidal clusters can be considered to be (quasi-)stable. However, for longer times, further clustering might occur.
The geometric structure of clusters as a function of the number of constituent particles in the plane has been studied theoretically. The optimal packings are obtained using the Lennard–Jones potential [
23], minimal second moment of the mass distribution (
, where
is the position of the particle
i and
is the position of the cluster center-of-mass) [
44,
45] and depletion interaction [
20].
The geometry of two-dimensional Lennard–Jones clusters where the particles interact via the Lennard-Jones potential was predicted based on genetic algorithm [
23]. The first few Lennard–Jones clusters are moniamond for
, diamond for
, triamond for
, tetriamond (chevron) for
, pentiamond (
), and hexiamond (hexagon) for
. Notably, all of predicted cluster structures, which can be regarded as a spreading of triangular networks, seem to be the same as those of
-minimal clusters for a given value of
[
23,
44]. The minimum colloid–colloid pair interaction, which is the sum of van der Waals attraction and dipolar repulsion, also produces configurations identical to that of both the Lennard–Jones cluster and
-minimal cluster [
23]. However, for hard-discs interacting via a short-ranged depletion attraction, the geometric structure of clusters, for example, for
, adopts three different ground state configurations; namely, parallelogram, chevron and triangle [
20]. We note that three such configurations are known as the tetriamond (4-iamond). A unique set of
-iamond sequences for
up to 5 are illustrated in
Figure 3; each case is indexed by a pair
, where
is the number of equilateral triangles and
m is the index of the configuration in the
-iamond sequence. Hence, for a given
, the maximum value of
m corresponds to the number of different configurations in each
n-iamond. As shown in
Figure 3, there is only one configuration for the moniamond (
), diamond (
) and triamond (
), while the tetriamond (
), pentiamond (
) and hexiamond (
) possess three, four and twelve possible configurations, respectively. The number of distinct polyiamonds of size
increases rapidly as
increases. Elements of the sequence for
and 8 are given in
Appendix A. It should be noted that the polyiamonds mentioned above are considered distinct if they have different number of equilateral triangles or different shapes. However, if account is to be taken of the number of vertices or point group rather than the number of equilateral triangles, then a few configurations in the sequence, e.g., the configurations (6,4), (6,6) and (6,9), are equivalent to each other. In the present work, we classify the colloidal clusters obtained via emulsion droplet evaporation by referring to the models mentioned above.
Figure 4 shows all clusters observed in simulations (
,
). For very small clusters, i.e.,
= 2–5, the structures are unique depending on
. These structures are identical to those of Lennard–Jones clusters,
-minimal clusters and clusters of attractive hard discs. This indicates that the nature of the colloid–colloid potential has little effect on the optimal structure of small clusters. However, for
, in addition to a familiar configuration of a "chevron" shape (labeled as (4,2) in
Figure 3), we obtain two other isomers: one is a (4,1) parallelogram and the other a (4,3) triangle. Four distinct isomeric structures are observed together for
, as pictured in
Figure 4. Surprisingly, they are exactly the same as the configuration of pentiamonds (
Figure 3). Of these isomeric structures, the (5,3) configuration, an equivalent version to the (6,11) configuration with a sixfold axis of symmetry, corresponds to the theoretically predicted Lennard–Jones cluster and
-minimal cluster. In the eighth clusters (
), we find seven isomeric structures compared to nine possible structures of hexiamonds (only nine different configurations are counted because the (6,11) and (5,3) configuration; the (6,4), (6,6) and (6,9) are equivalent to each other under point symmetry). Two missing structures are (6,3) and (6,4), in which case (6,4) has same as the structure as the Lennard–Jones cluster. For higher order clusters (
), all cluster structures obtained from computer simulations are also members of the set of polyiamonds (see
Appendix A), although the number of cluster structures is much less than that enumerated in the sequence of polyiamonds. For higher order clusters (
), all cluster structures obtained from computer simulations are also members of the set of polyiamonds (see
Appendix A), although the number of cluster structures is much less than that enumerated in the sequence of polyiamonds. This is because the number of higher order clusters, or equivalently the number of different structure types possibly formed, directly relates to cluster size distribution, which in turn is driven by various factors (i) a droplet (colloid) packing fraction
, (ii) relative size ratio
.
As an illustration of the role of colloid packing fraction, consider the cluster size distribution shown in
Figure 5, in which
is the number of clusters with
colloids—normalized by dividing by the number of single particles
. From
Figure 5a–d, the yield for high order clusters increases with increasing colloid packing fraction, signaled by a shift in cluster size distribution toward larger clusters. A similar trend can be observed in the cluster size distribution (not shown) when the interaction parameter
decreases from 1.18 to 1.09. This result can be explained by a larger probability of capturing colloids by droplets as the colloid concentration increases. Likewise, the bonding between the constituent colloids of clusters becomes more stable against thermal fluctuations as
decreases (see illustration in
Figure 1a), leading to an increase in the yield of large clusters.
A previous study [
20] enumerated a wide variety of cluster structures of six colloidal particles with short-ranged depletion interactions. By using complicated free energy calculations, the authors determined the occupation probabilities of ground-states and excited states, and the structural rearrangements of clusters. The cluster formation of planar triangular particles is expected to be influenced by a face-to-face contact alignment, which is quite different from a point-to-point contact alignment for spherical or disc-like particles. However, in our theoretical model, we explore the occupation probability of the clusters simply based on the number of ways to add a single particle to a given
-colloid cluster and with constraint that a point-to-face contact alignment is obeyed. We made this assumption based on the fact that the growth of cluster as a spreading of triangular networks was observed by the computer simulation [
23] and experiments [
24,
25]. As an illustration,
Figure 6 shows the possible configurations obtained for the three configurations (4,1), (4,2) and (4,3) by addition of one particle. In the case of (4,1), there are six ways to add one particle (labeled 1–6), but only three of them ((5,1), (5,2) and (5,4)) are inequivalent under rotational symmetries. Similarly, insertion of one particle in (4,2) can have one of three possible inequivalent configurations: (5,2), (5,3) and (5,4), whereas insertion of one particle in (4,3) produces a unique configuration; i.e., (5,2). From the number of occurrences, we calculate theoretically the probability of each inequivalent configuration (
p (theo.)), and the results are given in
Table 1. For comparison, we show the corresponding probabilities,
p (simul.), which are directly determined from the number of times a configuration is observed through MC simulations in the final stage.
Table 1 shows that the theoretical probabilities are in good agreement with the simulation observations for configurations from (1,1) to (4,4), or, alternatively, clusters containing three to seven particles. This is intuitively reasonable, since such cluster configurations are known as ground states whose occupation probabilities are determined primarily by the rotational entropy (permutational entropy) [
20]. At
, the overall agreement between the theoretical and simulated probability is generally good. However, in contrast to a high occupation probability of the (6,4) configuration predicted by the theory (
), we have not observed the isomeric cluster of (6,4) from the present simulation. This may be explained as follows. Theoretically, there exist three types of equivalent configurations mentioned above, i.e., (6,6), (6,9) and (7,18), that contribute to the total occupation probability of the configuration (6,4). Thus, the predicted probability of (6,4) is high. However, our theoretical calculation is based on the assumption that each
-disc cluster containing particles is formed by addition of single particle at each triangular lattice site to clusters containing
particles. Hence, the probability of finding certain cluster shape depends only on the number of ways to generate a packing. Within this model, we neglect the influence of long-ranged repulsive interaction and the presence of the droplet phase. In simulations and experiments, the repulsive interactions hinder spontaneous clustering and the droplets confine a certain number of particles at their surface. Despite a lack of consistency, our theoretical model is able to predict the simulation results quite accurately. However, it should be noted that this theory fails to provide the occurrence probability of clusters with one particle located at the center of the cluster—the internal particle, surrounded by the other particles. This is because in both experiments [
6] and our simulations, while the colloids adsorbed at the droplet interface, there does not appear to be a route to an internal colloid upon collapse. A constraint subject to this route may hamper the rearrangement of particles. For
, the absence of many isomeric clusters in the simulation compared to the theoretical calculation comes from the fact that most of them have relatively low probabilities of being found, and therefore can easily be missed in the simulation data.
4. Conclusions
We have considered a two-dimensional suspension in which the spherical colloids and droplets are confined at the interface, making the system quasi-two-dimensional (q2D). In such a q2D system, the particle translational motion is restricted to a 2D surface and rotational motion is restricted about a line normal to the plane of the interface. We note the difference between the q2D that we studied and the q2D confinement of other studies. For example, in the work by Manahoran et al. [
20], particles were confined between two coverslips with a certain spacing to generate 2D templates. Another q2D monodisperse colloid suspension was performed by construction the chamber with a thickness only slightly greater than one colloid particle diameter [
47], thereby restricting the centers of the colloid particles to lie close to a plane, but allowing small amplitude deviations from that plane.
By combining the Einstein relationship and Sonine approximation, we obtaiedn a rough estimate of the simulation time scale (of order 1 s), which is typically three orders of magnitude smaller than the experimental time scale.
On the time scale over which the simulation is performed, we observed stable clusters that ranged from moniamonds to 8-iamonds. From a comparison of the six-particle ground state clusters in [
20] to three possible configurations of tetriamonds, we suggest that all higher order clusters of polyiamonds (
) may be expected to be ground states. In particular, using a simple geometrical rule, on the basis of permutational degeneracy, we have calculated the occupation probabilities of polyiamonds and found good agreement with those obtained from simulations. This is intuitively reasonable, since the occupation probability of ground state clusters is largely controlled by symmetry [
19,
20].
Unlike 3D systems, in which emulsion droplets play the role of a driving force to generate the uniqueness in cluster structures, in 2D systems droplets have little effect on stable cluster structures except that they still significantly reduce the probability of clusters that contain an internal particle. We interpret a coexistence of different structures in 2D clusters as a result of a smaller free-energy difference between those structures compared to that in 3D clusters.
Our simulation also shows that the model potentials predict the growth of clusters of a span of triangular lattices, broadly consistent with that obtained for the global minimization of the Lennard–Jones potential, the second moment of the mass distribution and an attractive van der Waals interaction. Furthermore, we found no change in the cluster structures for a variety wide of the tunable parameters of the model. These results suggest that the structures of 2D colloidal clusters are insensitive to the nature of the potential.
It is worth noting that the simulations described here could lead to the formation of uniform assembly of colloidal particles characterized by a range of well-defined sizes, shapes and structures. The availability of well-controlled assembly of colloidal particles will provide an opportunity to experimentally probe the hydrodynamic, aerodynamic and optical properties. For example, the 2D structure of polystyrene microspheres could be put to use for microlenses with controllable focal lengths [
48]. In addition, some of the colloidal structures, e.g. the dimers, trimers and tetramers, can be further explored as a class of new building blocks for self-assembly to generate mesostructured systems that may find uses in practical applications, such as photonics, electronics, plasmonics and condensed matter physics [
3,
48].