## Abstract

Nucleation and growth of crystalline phases play an important role in a variety of physical phenomena, ranging from freezing of liquids to assembly of colloidal particles. Understanding these processes in the context of colloidal crystallization is of great importance for predicting and controlling the structures produced. In many systems, crystallites that nucleate have structures differing from those expected from bulk equilibrium thermodynamic considerations, and this is often attributed to kinetic effects. In this work, we consider the self-assembly of a binary mixture of colloids in two dimensions, which exhibits a structural transformation from a non–close-packed to a close-packed lattice during crystal growth. We show that this transformation is thermodynamically driven, resulting from size dependence of the relative free energy between the two structures. We demonstrate that structural selection can be entirely thermodynamic, in contrast to previously considered effects involving growth kinetics or interaction with the surrounding fluid phase.

## INTRODUCTION

As is well known because of Ostwald (*1*), the solid phase that, under given conditions, is most thermodynamically stable is typically not the one which first nucleates from a fluid. This behavior is common to many systems, and the explanation which is usually presented for this behavior invokes classical nucleation theory (CNT) (*2*). If one structure has a lower free energy barrier Δ*F** to formation than another, then its nucleation rate (*3*)

Although useful in many cases, the CNT approach is not always sufficient to explain certain observed crystallization phenomena (*4*, *5*). One issue not directly addressed by the classical explanation is the mechanism by which a metastable solid phase forming initially is converted to the final stable phase. An interesting manner in which this can occur is diffusionless transformation, as observed experimentally and in simulation in a three-dimensional (3D) binary colloidal system (*6*, *7*), as well as in simulation of a similar two-dimensional (2D) system (*8*). In these examples, the transformation in 3D is from a non–close-packed body-centered cubic (bcc) CsCl structure to a close-packed face-centered cubic (fcc) CuAu structure. Similarly, the transformation seen in the 2D system is between analogous square (fourfold coordinated) and hexagonal (sixfold coordinated) structures, respectively. Transitions from lower to higher coordinated lattices have been studied previously (*5*, *9*–*11*) but have been attributed to the solid phase matching its local structural order to that of a parent liquid phase. Symmetry considerations based on Landau theory suggest that bcc solids should be favored to nucleate from liquids because of configurational entropy favoring bcc arrangements close to the solid-liquid coexistence curve (*12*).

We show for the aforementioned 2D system (*8*), however, that such a transformation occurs in solid crystallites in a dilute fluid phase and even in isolation and is driven by the size dependence of the free energies of these finite-sized crystallites themselves. Although some previous works (*13*–*15*) use size in the discussion of transformations or other polymorphic behavior, a surrounding dense liquid phase and, in most cases, kinetic effects, are present and are deemed to play dominant roles in the observed behavior. Size-dependent thermodynamically driven selection of metastable polymorphs in inorganic solid phases has been considered (*16*). However, to our knowledge, this behavior also coupled with a diffusionless pathway between phases permitting transformation during crystallization as observed here has not been clearly illustrated elsewhere in the literature. The model system used in this work is a binary system with tunable short-ranged interactions, such as might be realized through DNA functionalization of micrometer-sized colloids (*6*, *8*). An understanding of the behavior in this system is therefore applicable to similar physical systems and to a more general understanding of the variety of nonclassical pathways for crystallization (*17*).

## RESULTS

Molecular dynamics (MD) simulations of self-assembly similar to those presented by Song *et al*. (*8*), in which this specific 2D multiflavored system was first presented, were performed; Fig. 1 shows representative snapshots from such a simulation (see also movie S1). Under the conditions chosen, the like (A-A and B-B) attractive interactions have 15% of the strength of the unlike (A-B) interactions. These attractive strengths are denoted as *E*_{AA}, *E*_{BB}, and *E*_{AB}, respectively, and for convenience, we choose *E*_{AB} as a reference energy scale and work with the relative like attractive strengths *E*_{AA}/*E*_{AB} and *E*_{BB}/*E*_{AB}. Starting from a disordered mixture of particles (Fig. 1A, I), self-assembly proceeds to the compositionally ordered hexagonal structure (Fig. 1A, IV). However, instead of assembling directly into this structure, the first crystallites that appear have a square structure (Fig. 1A, II). On the basis of lattice energies alone, the hexagonal lattice should be favored under these conditions. Even at finite temperature (*k*_{B}*T*/*E*_{AB} = 0.1), where entropic effects could alter the relative stabilities of the two structures, hexagonal crystals are still eventually produced. The progression of assembly from a dilute fluid phase through this square intermediate raises the question of whether this behavior is similar to previously observed transformation effects, in the sense that it is related to interactions between the solid and the surrounding fluid phases, or whether it is driven by different underlying causes.

A few important observations can be made about the nature of the self-assembly process. First, because of the very low free energy barriers to nucleation and critical nucleus sizes (see fig. S1), aggregation of particles into clusters occurs rapidly, quickly reducing the density of the already dilute vapor phase (ρσ^{2} = 0.1, packing fraction ϕ = 0.083). Beyond this point, the predominant method for crystal growth is the merging together of smaller clusters. In addition, transformations of the type depicted in Fig. 1A were only observed to occur to square crystallites above a certain size. Furthermore, these crystallites moving freely through the simulation box occasionally collided and merged, as shown in Fig. 1B; these events were followed rapidly by transformation of the resulting larger cluster to a hexagonal structure. Figure 1C shows another example of the observed transformation from a square to a hexagonal arrangement. These events suggest that the transformations are size dependent, pointing toward a possible structural and therefore thermodynamic origin. The low-density conditions and lack of order in the fluid phase furthermore suggest that stabilization of the square structure is related to the solid crystallites themselves and is not due to their interaction with the surrounding fluid.

### Size dependence and reversibility of transformations

To provide better evidence for a thermodynamic driving force behind the transformations, we now consider MD simulations starting not from dilute fluid but from predefined lattices of differing structures and sizes. We selected ranges of values of interaction strengths *E*_{AA}/*E*_{AB} and *E*_{BB}/*E*_{AB}, numbers of particles *N*, and temperatures *k*_{B}*T*/*E*_{AB}. For each set of parameters, 10 replicate simulations were performed from each starting structure to demonstrate the reproducibility of the transformations. Figure 2A shows the predominant structures as a function of *E*_{AA}/*E*_{AB} = *E*_{BB}/*E*_{AB} and *N* resulting from simulations starting with crystallites in each structure. It can be seen (Fig. 2A, left) that square crystallites transform to hexagonal structures for sufficiently high like attraction strengths. This is an expected result, given the presence of two additional like particle contacts present for each particle in the hexagonal lattice, which are absent in the square lattice. However, note that as the size *N* becomes smaller, square crystallites remain stable at larger interaction strengths; in other words, conditions that may favor the hexagonal structure for large crystallites favor square arrangements for small crystallites. This is consistent with observations from self-assembly simulations.

Furthermore, the right panel of Fig. 2A, showing results for crystallites starting in hexagonal arrangements, demonstrates that these transformations do not only occur in the square to hexagonal direction but also are reversible in the sense that sufficiently small hexagonal crystallites will spontaneously transform to square structures. This can also be seen in individual trajectories. Figure 2B shows a square crystallite with *N* = 100 transforming to a hexagonal crystallite; a smaller hexagonal crystallite with *N* = 64 under otherwise identical conditions transforms in the opposite direction, as shown in Fig. 2C. Last, sufficiently small crystallites near the square-hexagonal boundary will oscillate continuously between the two structures, as shown in Fig. 2D, clearly demonstrating the reversibility of these transformations and confirming a thermodynamic origin for the effect. For visualizations of transformation trajectories, see movies S2 to S4 corresponding to Fig. 2 (B to D), respectively. Although there is some difference between the results starting from square and from hexagonal structures, this becomes substantial only for larger *N* and can ultimately be attributed to the finite lengths of the MD simulations (see fig. S2).

Transformations can be observed for temperatures and ranges of interaction strengths other than those presented in Fig. 2 alone. Size and temperature dependence of the favored structures, as well as the effect of asymmetric attractive interactions (*E*_{AA}/*E*_{AB} ≠ *E*_{BB}/*E*_{AB}), are reported in fig. S3. It can be seen that the square lattice is stabilized for larger values of the like interactions as the temperature increases. This is likely due to a greater vibrational contribution to its entropy (*18*) compared to the hexagonal lattice: As shown by its radial distribution function (RDF) in fig. S4B and by movies S1 to S4, the square lattice is substantially less rigid. As interactions become progressively more asymmetric, it is also possible to stabilize a “rhombic” phase, as shown in fig. S4A, where particles of both types have the four unlike contacts common among all three structures but one type has an additional like contact corresponding to the more self-attracting component.

Note that the shapes chosen for these simulations, represented at the top of Fig. 2A, were selected to be physically consistent with respect to transformation (see fig. S4). In addition, square clusters with surfaces oriented mostly along (10) rather than (11) planes predominate during assembly, and (10) surfaces have a lower surface energy in the *T* → 0 limit, motivating the use of the selected shapes, which present these particular surfaces. Additional tests, as shown in fig. S5, show that the transformation effect is still present and qualitatively similar when more complicated shapes are chosen; however, the square structure can also be stabilized over the hexagonal structure by increasing the crystallites’ perimeters while holding *N* fixed. This suggests that the thermodynamic effect at play involves a difference in line tensions between the surfaces of crystallites with different structures, which drives the transformations. Hence, *N* is not the only order parameter which governs the transformation behavior, but it appears to characterize the effect well in assembling compact crystallites in addition to the isolated systems used here.

Keep in mind that although these crystallites are simulated in isolation while a dilute fluid is present during MD simulations of self-assembly, the two kinds of conditions should be comparable. The fluid quickly becomes less dense than the starting ρσ^{2} = 0.1 as crystals assemble. For comparison, approximate densities of the square and hexagonal lattices are ρσ^{2} = 0.946 and 1.093, respectively. Note also that although, at higher temperatures, particles occasionally come detached from the surfaces of the isolated crystallites, these events are rare and do not substantially affect crystallite geometries. Last, during assembly, shapes vary because of the stochastic nature of crystallization, continuous surface rearrangements, and merging of clusters. However, crystallites still increase in *N* and decrease in perimeter-to-area ratio as they assemble, regardless of their precise shapes.

A natural question that arises from these results is whether these reversible transformations appear only because of the specific chosen pairwise interactions or whether they are more generally realizable in a physical system. Figure 3 (A and B) shows the effect occurring in a different simulated system of multiflavored DNA-functionalized particles (DFPs) with explicitly modeled DNA strands, while Fig. 3C demonstrates the dependence on the average number of nearest neighbors based on trajectories starting from both square and hexagonal arrangements. Given that the effect appears in both the model system with short-range interactions and this DFP system with much longer-range interactions (and explicit strand binding and unbinding kinetics), it is reasonable to expect that it should be present for a variety of interaction ranges and in many physically realizable systems with similar interaction schemes.

### Free energies of finite-sized crystallites

Although the MD simulation results provide strong evidence supporting the thermodynamic origin for size-dependent stabilization of different structures, the finite lengths of the MD simulations may still have an effect on the apparent coexistence conditions for finite-sized crystallites. Free energy calculations can provide a more quantitative look at the transformation phenomenon and can also completely separate any such dynamical issues from the underlying thermodynamic driving forces. These calculations were performed through a modification of the Einstein molecule method (*19*) for nonperiodic crystals (see Materials and Methods and the Supplementary Materials for more details). In this way, transformations occurring in crystallites with specified shapes and exposed faces can be studied, allowing us to find free energy differences associated with specific transformation events. For comparison, we can first look at lattice energies for finite clusters. Figure 4A shows a representative set for nonperiodic square and hexagonal crystals, along with their respective periodic (*N* → ∞) limits. On the basis of MD data, the crossing point for transformation between the two structures is expected to appear near *N* = 144. However, the data of Fig. 4A suggest that the hexagonal structure should be favored not only in the periodic limit but also for nonperiodic crystals of all *N* in the considered range. Free energies for this range are presented in Fig. 4B. Here, a crossover at *N* ≈ 100 in the free energy curves is seen, indicating that the square lattice clusters smaller than this size should be favored over the hexagonal clusters and the reverse should be true for larger *N*. Although the agreement between MD results in Fig. 4C and expected equilibrium probabilities in Fig. 4D as calculated by *P _{j}* = exp(−

*A*/

_{j}*k*

_{B}

*T*)/∑

*exp(−*

_{i}*A*/

_{i}*k*

_{B}

*T*) is not perfect, the thermodynamic origin of this transformation effect is once again confirmed. Any observed differences may arise from the longer time necessary to observe transformation in MD as

*N*increases (see fig. S2) and from shape changes resulting from rearrangements of crystallites due to the sliding of layers, which occurs in normal MD simulations but may be hindered during thermodynamic integration.

The free energy data also confirm other trends observed from the MD simulations. Figure 5A shows the square-hexagonal phase boundaries for various values of *E*_{AA}/*E*_{AB}, *E*_{BB}/*E*_{AB}, and *N* at *k*_{B}*T*/*E*_{AB} = 0.08, while Fig. 5B demonstrates the temperature dependence of these boundaries. The trends here are in agreement with those described previously. The movement of the contours toward higher like attraction strengths with the increasing temperature, as well as the increased spacing between them, suggests that entropic effects are at the root of the stabilization of the square lattice for small sizes. Furthermore, the increased spacing between the contours as *N* becomes smaller points toward a surface-related phenomenon. To better understand the underlying cause, we decompose free energy into its component parts in two ways. First, by calculating free energies at many temperatures and extracting their temperature derivatives, we can obtain the internal energy as *U* = *A* − *T*(*∂A*/*∂T*)* _{V}*. In addition, a bulk contribution to some extensive quantity ξ should scale with

*N*, and a surface contribution should scale with

*N*against

*N*≥ 100 for extrapolation because of nonlinearity at lower

*N*and perform calculations with

*N*≥ 36 as numerical issues can result from the instability of very small crystallites, especially at higher temperatures.

Figure 5C demonstrates this decomposition for the static lattice energy *U*_{latt} under the same conditions shown in Fig. 4. In this form, it is clear that the hexagonal lattice has a lower bulk energy and that the surface energies are nearly identical; thus, no crossover is observed in the range of sizes shown. Turning to the results from the free energy calculations, Fig. 5 (D and E) shows the internal and free energies at *k*_{B}*T*/*E*_{AB} = 0.06. The surface contribution to the internal energy appears greater for the hexagonal lattice, but it alone is insufficient to cause a crossover. The free energy, including the −*TS* term, shows the expected crossover. In this particular case, the lower surface entropy of the square lattice is seen to yield the crossover by overcoming the free energy penalty in the periodic limit, although the thermal contribution to the surface internal energy (coming from the heat capacity) also plays a role. Examples under two additional representative conditions are provided in fig. S6. Overall, the effect can be seen to come from a combination of the contributions of heat capacity and entropy to the overall free energy.

## DISCUSSION

In summary, we have demonstrated that the transformation effect appearing in a 2D colloidal system during crystallization is a consequence of the size-dependent thermodynamics of the involved solid phases. Although initially developed to study a specific system of micrometer-sized DFPs, the primary model used here is based on simple pairwise interactions. Meanwhile, the explicit DFP model also discussed accounts for individual interactions between grafted DNA strands, but both types of systems show a similar transformation effect. The observed phenomenon may therefore be applicable to a variety of systems and may help to explain transformation and polymorph selection mechanisms during their nucleation and growth processes. For example, in 3D, the micrometer-sized particle system studied here produces CsCl and CuAu structures (*20*) analogous to the square and hexagonal structures considered in this work. Experimental and simulation evidence in a similar system (*6*) provides evidence of transformations, indicating that they occur in a diffusionless fashion similar to the transformations considered here. Although there are ultimately more pathways for transformation to occur in 3D due to the existence of multiple close-packed arrangements, a similar mechanism, which is responsible for their appearance in 2D, may drive them in 3D by stabilizing sufficiently small crystallites of the non–close-packed CsCl structure.

This particular transformation effect, as it appears during self-assembly, is remarkable in the sense that it is of a different nature from previously observed structural transformations in other systems. We find that, under conditions where barriers to nucleation are very low and crystal growth occurs through attachment of clusters, this effect still appears and typically manifests itself in crystallites whose sizes are much greater than critical. In this case, instead of a by-product of nucleation kinetics, it can be considered as a direct consequence of the relative thermodynamic stabilities of the lattices when accounting for the presence of their surfaces. Furthermore, it is present and can be explained in complete isolation from a surrounding fluid phase, contrary to other observations of similar behaviors. The initial appearance of the square structure in place of the hexagonal structure, which is stable in the bulk, is consistent with Ostwald’s step rule, but the transformation effect itself does not necessarily correspond with and cannot be fully explained by the associated CNT interpretation. In the end, these results highlight the importance of considering thermodynamic size and surface effects in understanding self-assembly processes.

## MATERIALS AND METHODS

MD simulations were performed using the Large-scale Atomic/Molecular Massively Parallel Simulator (LAMMPS) package (*21*). Interactions between “A-” and “B”-type colloids were modeled through a Fermi-Jagla pairwise potential (*22*). Details are available in Supplementary Text (see the “Pairwise interactions” section). All simulations were performed at constant NVT with a Langevin thermostat. Assembly simulations were performed with at least 1000 particles in a 1:1 ratio of A to B particles at a number density ρσ^{2} = 0.1. Those starting from isolated crystallites used clusters of sizes 36, 64, 100, 144, 196, 256, and 324 in square, hexagonal, and rhombic configurations. Details on equilibration and other simulation parameters are available in Supplementary Text (see the “MD simulation details” section).

To determine whether crystallites remained in or transformed between square (fourfold coordinated), hexagonal (sixfold coordinated), or rhombic (four- and fivefold coordinated) structures, a method based on RDFs was used. This is a simplified and slightly modified version of a method proposed by Oganov and Valle (*23*). The two given structures α and β and their RDFs *g*_{α}(*r*) and *g*_{β}(*r*) can be used to calculate a similarity parameter defined as follows

*S* is equal to 0 and 1 for completely dissimilar and identical structures, respectively. Pairwise RDFs were not used, and the parameter cannot distinguish between structures with identical underlying lattices but different compositional ordering. However, it was sufficient for the structures of interest here. RDFs were sampled from trajectories to give time series of 400 independent *g*(*r*) samples for each simulation. The interactions (*E*_{AA}/*E*_{AB} and *E*_{BB}/*E*_{AB}) for the reference RDFs were set to prevent transformations. The integrals were approximated by sums with *g*(*r*) binned by Δ*r* = 0.01 σ from *r* = 0 to *r _{c}* = 5 σ and smoothed with a Gaussian filter having an SD σ

*= 5 Δ*

_{g}*r*. Sample reference RDFs and details on the reference interactions are available in fig. S4.

An extension to Frenkel-Ladd integration (*24*, *25*), specifically to the Einstein molecule method (*19*), was devised to perform free energy calculations on nonperiodic crystals having exposed surfaces. As these crystals do not interact with themselves through periodic boundaries, constraints must be applied to remove both translational and rotational degrees of freedom. This was most straightforwardly done by not only fixing entirely one particle as in the standard method but also constraining additional particles in particular directions (*26*). These calculations were carried out using LAMMPS, similar to the method described by Aragones *et al*. (*27*), but with the additional rotational constraints. The same systems and conditions were studied as in the MD simulations, except for temperature, which spanned the range of *k*_{B}*T*/*E*_{AB} = 0.02 to 0.11 at intervals of 0.01. Gauss-Legendre quadrature with 32 sample points was used for the thermodynamic integration. Sizes *N* ≥ 36 were used as smaller systems tended to become unstable near one end point of the thermodynamic integration, leading to numerical issues. More details are available in Supplementary Text (see the “Free energies of nonperiodic crystals” section).

The simulations of DFPs with results presented in Fig. 3 used a model described by Ding and Mittal (*28*) in which single-stranded DNA is represented explicitly using a coarse-grained DNA model. Sixty single-stranded DNA strands, T_{6}C_{4} and T_{6}G_{4}, were grafted onto each particle using multiflavoring (*8*) with a blending ratio of γ = 0.05 to generate attractive like and unlike particle interactions. Simulations were performed with the HOOMD-blue package (*29*, *30*) at constant NVT using a Langevin thermostat (*k*_{B}*T*/ϵ = 0.37). Particles were constrained to a quasi-2D region using parallel Weeks-Chandler-Andersen (*31*) repulsive walls separated by 40 σ.

## SUPPLEMENTARY MATERIALS

Supplementary material for this article is available at http://advances.sciencemag.org/cgi/content/full/5/9/eaaw5912/DC1

Supplementary Text

Fig. S1. Temperature dependence of nucleation barriers and critical sizes at *E*_{AA}/*E*_{AB} = *E*_{BB}/*E*_{AB} = 0.15.

Fig. S2. Mean lengths of time that the initial structures in MD simulations of individual crystallites persist before either transformation occurs or the simulation ends.

Fig. S3. Predominant structures in transformation MD simulations for asymmetric interactions.

Fig. S4. Shapes and RDFs for crystallites of different structures.

Fig. S5. Predominant structures in transformation MD simulations for various crystallite shapes and aspect ratios.

Fig. S6. Comparison of components of free energy for additional cases of interactions and temperatures.

Fig. S7. Interaction potentials for different DFP models.

Fig. S8. Interaction dependence of the system crystallization temperature.

Fig. S9. Illustration of the thermodynamic pathway used for the Einstein molecule method applied to nonperiodic systems.

Movie S1. Animation showing a trajectory of self-assembly from which the snapshots presented in Fig. 1 are taken.

Movie S2. Animation corresponding to the transformation trajectory shown in Fig. 2B.

Movie S3. Similarly, animation corresponding to Fig. 2C.

Movie S4. Similarly, animation corresponding to Fig. 2D.

Reference (*32*)

This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial license, which permits use, distribution, and reproduction in any medium, so long as the resultant use is **not** for commercial advantage and provided the original work is properly cited.

## REFERENCES AND NOTES

**Acknowledgments:**

**Funding:**This work was supported by the U.S. Department of Energy, Office of Basic Energy Science, Division of Material Sciences and Engineering under award DE-SC0013979. This research used resources of the National Energy Research Scientific Computing Center, a DOE Office of Science user facility supported under contract no. DE-AC02-05CH11231. Use of the high-performance computing capabilities of the XSEDE (Extreme Science and Engineering Discovery Environment), which is supported by the National Science Foundation (project no. TG-MCB120014), is also gratefully acknowledged.

**Author contributions:**E.P. and J.M. designed the research and wrote the paper. E.P., R.M., H.Z., M.S., and Y.D. performed the research and analyzed the results.

**Competing interests:**The authors declare that they have no competing interests.

**Data materials and availability:**All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials. Additional data related to this paper may be requested from the authors.

- Copyright © 2019 The Authors, some rights reserved; exclusive licensee American Association for the Advancement of Science. No claim to original U.S. Government Works. Distributed under a Creative Commons Attribution NonCommercial License 4.0 (CC BY-NC).