# Predicting the outcome of the growth of binary solids far from equilibrium

###### Abstract

The growth of multicomponent structures in simulations and experiments often results in kinetically trapped, nonequilibrium objects. In such cases we have no general theoretical framework for predicting the outcome of the growth process. Here we use computer simulations to study the growth of two-component structures within a simple lattice model. We show that kinetic trapping happens for many choices of growth rate and inter-component interaction energies, and that qualitatively distinct kinds of kinetic trapping are found in different regions of parameter space. In a region in which the low-energy structure is an ‘antiferromagnet’ or ‘checkerboard’, we show that the grown nonequilibrium structure displays a component-type stoichiometry that is different to the equilibrium one but is insensitive to growth rate and solution conditions. This robust nonequilibrium stoichiometry can be predicted via a mapping to the jammed random tiling of dimers studied by Flory, a finding that suggests a way of making defined nonequilibrium structures in experiment.

## I Introduction

Molecular self-assembly is the spontaneous organization of components, which move around under e.g. Brownian motion but are otherwise left undisturbed, into ordered structures. Self-assembly holds considerable promise for materials science Whitesides and Grzybowski (2002); Glotzer and Solomon (2007); Frenkel (2015); Cademartiri and Bishop (2015). The goal of molecular self-assembly in the laboratory is often to make an equilibrium structure, and the laws of statistical mechanics indeed dictate that components undergoing Brownian motion will eventually build themselves into the structure of least free energy. In practical terms, however, slow dynamical processes can prevent equilibration from happening on the timescale available to the process in question Whitelam and Jack (2015). In such circumstances the processes of nucleation and growth lead instead to the formation of kinetically trapped, nonequilibrium structures. Multicomponent systems, i.e. systems composed of more than one type of component, are particularly susceptible to kinetic trapping because the slow rearrangement of component types within a solid structure can prevent them from achieving their equilibrium arrangement as the solid structure grows. Frequently, the outcome of the nucleation and growth of multiple component types is an ordered crystal structure within which component types are arranged in a nonequilibrium way Sanz et al. (2007); Peters (2009); Kim et al. (2009); Scarlett et al. (2010); Whitelam et al. (2014). Such structures have potentially useful properties. However, predicting their component-type arrangements is not possible in general, because we cannot predict the outcome of self-assembly when that outcome is not the equilibrium structure.

Here we use simulation and analytic theory to study the component-type arrangements formed during the growth of two-component structures within a simple lattice model. In accord with several experimental results, growth results in the formation of nonequilibrium structures for a large range of growth rates and inter-component interaction energies. In some regions of parameter space the properties of nonequilibrium structures vary continuously with growth rate, while in other regimes of parameter space these properties are insensitive to growth rate. In a region of parameter space in which the low-energy structure is a binary ‘checkerboard’ we show that the grown nonequilibrium structure displays a component-type stoichiometry that is insensitive both to growth rate and to the abundance of component types in solution. We show that this robust nonequilibrium stoichiometry can be predicted via a mapping to the jammed random tiling of dimers studied by Flory. These findings suggest a route to the rational design of defined nonequilibrium structures in experiment.

## Ii Model and Simulation Methods

Kinetic trapping of component types within growing multicomponent structures has a simple physical origin – the slow dynamics of particles within a solid – and so can be reproduced by simple physical models that account for this slow dynamics Sanz et al. (2007); Peters (2009); Kim et al. (2009); Scarlett et al. (2010). Here we consider a lattice model of growth similar to the models used in Refs. Whitelam et al. (2014); Hedges et al. (2014); Sue et al. (2015). We focus on growth in a 2D system, but we will also present results for higher dimensions. As sketched in Fig. 1(a), lattice sites can be unoccupied (white) or occupied by a particle of one of two types (red or blue). Red and blue particles (or components) experience color-dependent nearest-neighbor interactions (see Appendix A) of energy , , and , in units of (which we shall set equal to unity). On a fully-occupied lattice (one without white sites) these interactions are equivalent to the Ising model with magnetic field and coupling constant Sue et al. (2015). White, blue, and red sites also receive energetic penalties , and , respectively. Here sets the relative abundance of colored and white sites in notional ‘solution’ (i.e. in the absence of energetic interactions), and is the notional solution fraction of colored blocks that are blue. We evolved this model using a grand-canonical Monte Carlo procedure that respects detailed balance, and that resolves the stochastic binding and unbinding of red and blue components. Unbinding dynamics is naturally slow when components possess many colored neighbors; we also imposed a kinetic constraint that prevents any change of state of a lattice site that possesses only colored neighbors. This constraint, which preserves detailed balance, is intended to model the fact that relaxation dynamics within solid structures is slow. In what follows we shall describe growth simulations done in the presence and absence of the kinetic constraint. The latter type of simulation represents a convenient way to assess the outcome of growth on timescales longer than we could otherwise access. In most simulations described below we used a 2D square lattice of lattice sites whose long edges were periodic and whose short edges were not. We began simulations in the presence of a ‘seed’ at the left-hand short edge of the simulation box, with the rest of the box left white, so that we could study growth without waiting for nucleation to happen. By varying we could change the rate of growth of the colored assembly. In what follows we refer to ‘growth’ simulations in which the simulation was stopped after 90% of the box become occupied by colored components, and ‘maturation’ simulations in which structures grown in this manner were allowed to evolve for an additional Monte Carlo cycles.

## Iii Growth simulations

Growth carried out using different choices of the inter-component energetic parameters and ^{1}^{1}1The Ising parameters and are a convenient way to describe the three blue-red interactions , , and , but they provide only a partial description of the model’s parameter space, which includes chemical potential terms and interactions with white sites. For instance, one can add constant terms to the parameters that leave and unchanged but change the energetics of the model., shown in Fig. 1(b), is similar in the following respects (see Fig. 2). At vanishing rates of growth a structure resembling the equilibrium one is generated; at very large rates of growth a ‘solid solution’ is obtained, i.e. red and blue components are arranged randomly on the lattice in proportion to their solution proportions; and at intermediate rates of growth one obtains nonequilibrium structures that differ from both of these limiting cases. These nontrivial nonequilibrium structures can be different in different parameter regimes. For instance, using the ‘ferromagnetic’ energy scale hierarchy shown by the symbol on Fig. 1(b), nonequilibrium structures include ‘critical’ arrangements in which red and blue component domains of a broad size distribution are present Whitelam et al. (2014) (this behavior may be related to that seen in certain irreversible cellular automata Ausloos et al. (1993); Candia and Albano (2008)). At the parameter combination labeled , nonequilibrium structures consist of large domains of the blue component within which a small red impurity fraction is found (see Fig. S2). This impurity fraction is only weakly sensitive to growth rate over some range of growth rates, a result that reproduces the qualitative outcome of growth in experiments and off-lattice simulations done by other authors Kim et al. (2009). The reproduction of these results by the present model suggests that it captures key physical aspects of real growth processes. Finally, at the parameter combination labeled on Fig. 1(b), growth results in a nonequilibrium component-type stoichiometry identical (on an nbo lattice) to that seen in a certain metal-organic framework; this stoichiometry is insensitive to growth rate and component solution stoichiometry over some range of those parameters Sue et al. (2015). As we shall show, robust nonequilibrium stoichiometry is also seen in other parameter regimes left of the dotted line in Fig. 1. Here we aim to provide a partial physical understanding of this behavior.

We start by showing in Fig. 2 the outcome of growth simulations done in the aforementioned regime of parameter space, where the red-blue energetic interaction is lower in energy than both of the like-color interactions. Here the low-energy structure, and the thermodynamically stable structure for the parameter values we shall consider, is an alternating red-blue ‘antiferromagnet’ or ‘checkerboard’. In the top panel of Fig. 2 we show results for the parameter combination (meaning that red-red and blue-blue interactions are of equal strength), while in the bottom panel we consider an asymmetric energy hierarchy for which . At low rates of growth, in both cases, the structure generated dynamically, upon 95% filling of the simulation box (panel (a)), is close in nature to the equilibrium structure, and so possesses a blue fraction (fraction of colored components that are blue) . For large rates of growth the structures obtained are kinetically trapped arrangements of components whose blue fraction is related to that of the notional solution (we considered three different solution stochiometries, shown as black, red, and blue lines). At intermediate rates of growth the blue fraction of the grown structure interpolates smoothly between these limiting cases. Configuration snapshots are shown in Fig. S3. However, when allowed to further evolve or ‘mature’ for Monte Carlo sweeps (panel (b)) (in the absence of the kinetic constraint so as to effectively allow access to longer timescales), structures generated at intermediate growth rates did not evolve to equilibrium, but instead became kinetically trapped in configurations whose blue fractions display plateaux as a function of growth rate. That is, the stoichiometry of those nonequilibrium structure is insensitive to growth rate. Furthermore, in the case of the asymmetric energetic hierarchy (bottom panel) this stoichiometry was also insensitive to solution stoichiometry. Near-equilibrium and far-from-equilibrium regimes are separated by a regime of large fluctuations of color and energy (panels (c) and (d)), suggesting the existence of a nonequilibrium phase transition similar to that seen in the ‘ferromagnetic’ regime of parameter space Whitelam et al. (2014).

Structures generated dynamically in the presence of certain energetic interactions therefore display a stoichiometry that is different to the 1:1 equilibrium one, but that is robust with respect to changes of growth rate and solution stoichiometry over a considerable range of those parameters. We call these robust nonequilibrium stochiometries ‘magic numbers’. The existence of magic numbers has potential application for materials science, because it suggests that one can grow two-component solids, out of equilibrium, in a predictable manner. Magic number materials may have already been synthesized. One particular two-component metal-organic framework (MOF), called MOF-2000, displays a stoichiometry that is robust to solution stoichiometry over a considerable range Sue et al. (2015). The numerical value of this stoichiometry can be reproduced by the growth of magic-number structures using the present model on a 3D framework whose topology is appropriate to the crystal structure of MOF-2000 Sue et al. (2015).

In physical terms, nonequilibrium structures emerge because microscopic contacts that are not the equilibrium or ‘native’ red-blue one (such as red-red contacts) can appear stochastically during growth and can become trapped by the arrival of additional material. In some regions of parameter space the properties of the resulting kinetically trapped structures vary continuously with growth rate and solution stoichiometry; in the magic number regime they do not. As shown in Fig. 3, magic number configurations, reached upon ‘maturation’ of a structure after its initial growth, are long-lived: evolution to the equilibrium checkerboard structure does not happen on the timescale of simulation (Fig. 3d; Fig. S5). Magic number structures can also be generated by zero-temperature, single-spin-flip Monte Carlo sampling of fully-occupied lattices; that is, magic number structures are also inherent structures of the lattice Hamiltonian. These inherent structures are accessible from a wide range of initial conditions (they have a large basin of attraction), and the numerical values of the associated magic numbers are dependent upon lattice connectivity and dimensionality (Fig. S6).

## Iv Mapping to jammed dimer systems

The interaction energies used to obtain the magic numbers seen in Fig. 2 satisfy the hierarchy . In words, the blue-red contact is the ‘native’ or equilibrium one; red-red contacts are higher in energy but can occur during growth; and blue-blue contacts are so unfavorable that they cannot form at reasonable rates of growth. In one dimension this energetic hierarchy results in the growth of red-blue arrangements, such as those shown in Fig. 4(a), that map to a tiling of dimers with voids. ‘Dimers’ are red-blue pairs, and ‘voids’ are red particles. The particle-to-dimer mapping produces one of two equivalent void arrangements, as shown. The long-time outcome of our growth process then becomes equivalent to that of random sequential absorption (RSA) of dimers on a one-dimensional lattice. This problem was studied by Flory Flory (1939), who computed the dimer filling fraction to be . The mean blue fraction of our equivalent red-blue structure is then half of this value, i.e. . This value is indeed the mean value of the stoichiometry of inherent structures of our lattice model (recall that in this regime of parameter space the ‘matured’ growth configurations are also inherent structures of the model with no white sites), which we shown in Fig. 4(b). Thus, the nonequilibrium stoichiometry resulting from growth can be predicted via a mapping to a jammed system of dimers.

The equivalence between our growth process and dimer deposition does not hold in dimensions greater than one. Nonetheless, we can use the Flory result to estimate numerically the magic number ratio seen in our growth simulations in 2D and 3D. Consider a periodic hypercubic lattice that possesses lattice sites or nodes in each dimension, and so has nodes in total. Each node may be occupied by one red or one blue particle. Let be the number of void sites that exist on a connected row of nodes in any given dimension, and assume that Flory (1939). Assume further that each dimension is independent, so that each void site connects to a continuous chain of void sites that extends independently into each of the remaining dimensions; see Fig. 4(c,d). Thus, each void region contains in total voids. Summed over all independent dimensions there therefore exist voids in total, meaning that the void density is . We therefore predict the nonequilibrium ‘magic number’ stoichiometry of our red-blue structure, grown in dimensions, to be

(1) |

The magic number structures seen in the 2D growth processes whose results are reported in Fig. 2(b, bottom) and Fig. 3 have a stoichiometry (magic number blue fraction) of . The estimate of Eq. (1) is , which agrees closely with our inherent simulation results (Fig. 5) and with the plateaux seen in our growth simulations (Fig. 3 and Fig. 2(b) bottom). Thus the expression (1) can be used to predict the nature of a kinetically trapped structure generated in 2D by two-component growth. In 3D the estimate (1) predicts a magic number blue fraction of ; our inherent structure simulations done in 3D display a very similar magic number ratio of ; see Fig. 5. In dimension 4 and up the predictions of Eq. (1) depart from the results of our numerical simulations (see Fig. 5), signaling the breakdown of the approximations we used to derive the equation. But in dimensions relevant to laboratory self-assembly, for square and cubic crystal structures, we can rationalize the stoichiometry that results from two-component kinetic trapping by analogy to a jamming problem.

We have shown that of a lattice model of two-component growth displays a rich range of phenomenology, key aspects of which reproduce behavior seen in experiments Kim et al. (2009); Sue et al. (2015). Growth can result in near-equilibrium structures and far-from-equilibrium structures. In certain regimes of parameter space the component-type stoichiometry of these nonequilibrium structures is independent of growth rate and solution stoichiometry, and the numerical value of this stoichiometry can be predicted via a mapping to a jammed tiling of dimers. These observations suggest that one can grow, far from equilibrium, defined two-component structures in experiment.

## Appendix A Further details of simulation methods

Our lattice model has energy function

(2) |

The first sum runs over all distinct nearest-neighbor interactions, and the second sum runs over all sites. in Eq. (2) can be either w (white), b (blue) or r (red), depending on the color of node ; is the interaction energy between colors and ; and the chemical potential is , and for w, b and r, respectively. In the absence of pairwise energetic interactions (i.e. in notional ‘solution’), the likelihood that a given site will be white, blue or red is respectively .

Monte Carlo simulations were done as follows. We started with a simulation box that is 400 sites wide and 40 sites high, with the first six columns populated with the equilibrium checkerboard structure. We selected a node at random, and proposed a change of color of that node. If the chosen node was white, we attempted to make it colored; if the chosen node was colored, we attempted to make it white. If the chosen node was white, then we proposed to make it blue with probability ; otherwise, we proposed to make it red. No red-blue interchange was allowed, mimicking the idea that unbinding events are required in order to relax configurational degrees of freedom. To maintain detailed balance with respect to the stated energy function, the acceptance rates for these moves were as follows ( is the energy change resulting from the proposed move):

The notional solution abundances of red and blue are controlled by the chemical potential term that appears in Eq. (2) and therefore in the term . Our choice to insert blue particles with likelihood does not by itself result in a thermodynamic bias for one color over the other (because this bias in proposal rate is countered by the non-exponential factors in the acceptance rates). Instead, we bias insertions so that the dynamics of association is consistent with the thermodynamics of the model. For instance, if blue particles are more numerous in solution than red ones, we consider it to be physically appropriate to insert blue particles into the simulation box more frequently than red particles. Consider the limit of large positive : the ‘solid solution’ that results as the box fills irreversibly with colored particles will have a red:blue stoichiometry equal to that of the notional solution only if blue particles are inserted with likelihood . (As a technical note, the chemical potential term present in ends up simply canceling the non-exponential factors in the acceptance rates, but we have chosen to write acceptance rates as shown in order to make clear which pieces are imposed by thermodynamics, and which pieces we have chosen for dynamical reasons).

We also imposed a kinetic constraint that prevents any change of state of a lattice site that possesses only colored neighbors. This constraint, which respects detailed balance, is intended to model the fact that relaxation dynamics within solid structures is slow. In some simulations we omit the kinetic constraint in order to assess the outcome of slow internal evolution on timescales longer than we could otherwise access. In some regimes such constraint-free evolution leads rapidly to equilibrium, while in others it does not. For instance, for the inter-component interaction energies used to obtain Fig. S2, grown structures evolve quickly to equilibrium if the kinetic constraint is not used. The kinetic constraint is therefore needed in order to capture the physical character of growth seen in experiments. By contrast, for the interaction energies used to obtain Fig. 2, grown structures evolved even in the absence of the kinetic constraint fail to reach equilibrium, because of the deep kinetic traps associated with interaction energies large on the sale of . There we can omit the kinetic constraint in order to effectively simulate longer, and still obtain nontrivial results.

The parameter values obtained from Refs. Whitelam et al. (2014), Kim et al. (2009) and Sue et al. (2015) and marked on Fig. 1 are , (also see Table S1), and (,,).

Inherent structures in used to make Fig. 4b and Fig. 5 were obtained using zero-temperature single-spin-flip moves Sue et al. (2015) (also see Fig. S6) starting from initial conditions in which all 2000 sites of the periodic system are randomly colored red or blue, with equal likelihood (Fig. S6). This procedure was carried out until no more spin flips occurred. At least 100 independent inherent structures were obtained for each datapoint. As shown in Fig. S6d, the average value of the resulting is insensitive to system size.

###### Acknowledgements.

This work was done at the Molecular Foundry at Lawrence Berkeley National Laboratory, supported by the Office of Science, Office of Basic Energy Sciences, of the U.S. Department of Energy under Contract No. DE-AC02-05CH11231.## References

- Whitesides and Grzybowski (2002) George M Whitesides and Bartosz Grzybowski, “Self-assembly at all scales,” Science 295, 2418–2421 (2002).
- Glotzer and Solomon (2007) Sharon C Glotzer and Michael J Solomon, “Anisotropy of building blocks and their assembly into complex structures,” Nature Materials 6, 557–562 (2007).
- Frenkel (2015) Daan Frenkel, “Order through entropy,” Nature Materials 14, 9–12 (2015).
- Cademartiri and Bishop (2015) Ludovico Cademartiri and Kyle JM Bishop, “Programmable self-assembly,” Nature Materials 14, 2–9 (2015).
- Whitelam and Jack (2015) Stephen Whitelam and Robert L. Jack, “The statistical mechanics of dynamic pathways to self-assembly,” Annual Review of Physical Chemistry 66, 143–163 (2015).
- Sanz et al. (2007) Eduardo Sanz, Chantal Valeriani, Daan Frenkel, and Marjolein Dijkstra, “Evidence for out-of-equilibrium crystal nucleation in suspensions of oppositely charged colloids,” Physical Review Letters 99, 055501 (2007).
- Peters (2009) Baron Peters, “Competing nucleation pathways in a mixture of oppositely charged colloids: Out-of-equilibrium nucleation revisited,” The Journal of Chemical Physics 131, 244103 (2009).
- Kim et al. (2009) Anthony J Kim, Raynaldo Scarlett, Paul L Biancaniello, Talid Sinno, and John C Crocker, “Probing interfacial equilibration in microsphere crystals formed by dna-directed assembly,” Nature materials 8, 52–55 (2009).
- Scarlett et al. (2010) Raynaldo T Scarlett, John C Crocker, and Talid Sinno, “Computational analysis of binary segregation during colloidal crystallization with dna-mediated interactions,” The Journal of Chemical Physics 132, 234705 (2010).
- Whitelam et al. (2014) Stephen Whitelam, Lester O Hedges, and Jeremy D Schmit, “Self-assembly at a nonequilibrium critical point,” Physical Review Letters 112, 155504 (2014).
- Hedges et al. (2014) Lester O Hedges, Ranjan V Mannige, and Stephen Whitelam, “Growth of equilibrium structures built from a large number of distinct component types,” Soft matter 10, 6404–6416 (2014).
- Sue et al. (2015) Andrew C.-H. Sue, Ranjan V. Mannige, Hexiang Deng, Dennis Cao, Cheng Wang, Felipe Gándara, Fraser Stoddart, Stephen Whitelam, and Omar M. Yaghi, “Two-component metal-organic framework displaying compositional robustness to solution constitution.” PNAS 112, 5591–5596 (2015).
- (13) The Ising parameters and are a convenient way to describe the three blue-red interactions , , and , but they provide only a partial description of the model’s parameter space, which includes chemical potential terms and interactions with white sites. For instance, one can add constant terms to the parameters that leave and unchanged but change the energetics of the model.
- Ausloos et al. (1993) Marcel Ausloos, Nicolas Vandewalle, and Rudi Cloots, “Magnetic eden model,” EPL (EuroPhysics Letters) 24, 629 (1993).
- Candia and Albano (2008) Julián Candia and Ezequiel V Albano, “The magnetic eden model,” International Journal of Modern Physics C 19, 1617–1634 (2008).
- Flory (1939) Paul J Flory, “Intramolecular reaction between neighboring substituents of vinyl polymers,” Journal of the American Chemical Society 61, 1518–1521 (1939).

Supplementary Materials

Parameters from Ref. Kim et al. (2009) | Energies from Ref. Kim et al. (2009) | Our parameter equivalents | ||||||
---|---|---|---|---|---|---|---|---|

hierarchy # | = | = | = | |||||