Bottom-Up Colloidal Crystal Assembly with a Twist.

Globally ordered colloidal crystal lattices have broad utility in a wide range of optical and catalytic devices, for example, as photonic band gap materials. However, the self-assembly of stereospecific structures is often confounded by polymorphism. Small free-energy differences often characterize ensembles of different structures, making it difficult to produce a single morphology at will. Current techniques to handle this problem adopt one of two approaches: that of the "top-down" or "bottom-up" methodology, whereby structures are engineered starting from the largest or smallest relevant length scales, respectively. However, recently, a third approach for directing high fidelity assembly of colloidal crystals has been suggested which relies on the introduction of polymer cosolutes into the crystal phase [Mahynski, N.; Panagiotopoulos, A. Z.; Meng, D.; Kumar, S. K. Nat. Commun. 2014, 5, 4472]. By tuning the polymer's morphology to interact uniquely with the void symmetry of a single desired crystal, the entropy loss associated with polymer confinement has been shown to strongly bias the formation of that phase. However, previously, this approach has only been demonstrated in the limiting case of close-packed crystals. Here, we show how this approach may be generalized and extended to complex open crystals, illustrating the utility of this "structure-directing agent" paradigm in engineering the nanoscale structure of ordered colloidal materials. The high degree of transferability of this paradigm's basic principles between relatively simple crystals and more complex ones suggests that this represents a valuable addition to presently known self-assembly techniques.

The ability to intelligently direct the self-assembly of ordered materials at and below the micrometer length scale is the cornerstone of numerous modern technologies. It is a necessity in designing zeolite frameworks in the chemical industry,[1−4] biomimetic scaffolding used to grow artificial organs with stem cells,[5] and crystal lattices for optical applications.[6,7] Crystallization of colloidal suspensions often provides a direct route to some of these technologies; however, precise control during this phase transformation is usually difficult to achieve. The reason is that many crystals naturally form defects and exhibit polymorphism, where very closely related structures are separated by only a small free-energy difference. Two broad paradigms currently exist concerning control over colloidal crystal habit and polymorphism. In the “top-down” approach, a planar template or some other external guiding force is used to bias the formation of a specific crystal.[4,8,9] In such a case, the information about the final state of the system is “programmed” into the morphology of the template. Unfortunately, this is far from a panacea as this template must be made uniquely commensurate with a face of the desired crystal present in only that form and is therefore highly system-dependent.

Alternatively, this information can be encoded in the colloids themselves in the so-called “bottom-up” approach, where anisotropy in the shape, symmetry, and surface interactions of the colloids can be tuned to target a desired structure.[10−12] While this approach is attractive, its kinetics are subject to the smoothness of the free-energy landscape connecting competing forms and can be sensitive to small changes in the surface anisotropy.[10,13−17] Therefore, experimental realizations are subject to the fidelity of the synthesis protocol used to manufacture the constituent colloids.

Following queues from biological systems, one can imagine another alternative where this information is instead encoded into entirely separate molecular “chaperones” that guide the assembly process. Indeed, the utility of small cosolute molecules, or “structure-directing agents” (SDA), for the directed synthesis of zeolites and other microporous structures is already well-appreciated.[2,3,18,19] For instance, the introduction of organic polycations during the growth of siliceous ZSM-5 catalysts has been shown to control the morphology of the resulting crystallite.[2] This is believed to be a function of the way in which the nitrogen center on the amine positions itself at junctures between pores, while the aliphatic arms reside in the pores surrounding the juncture. By and large, the utility of SDAs has been linked to the chemical (thermal) interaction between the SDA and the porous material.[2,20−22] To our knowledge, the assembly of colloidal crystals via SDAs has not been realized or even appreciated.

Recent simulations have revealed how polymers can be used as SDAs to entropically bias the formation of specific hard-sphere polymorphs during crystallization, such as the face-centered cubic (fcc) or hexagonal close-packed (hcp) structures.[23−25] This mechanism exploits the difference in the interstitial void patterns between different crystals to amplify the free-energy difference between them. By connecting an otherwise pure colloid system to an athermal polymer reservoir at sufficiently high osmotic pressure, the polymer both induces crystallization of the colloids in a close-packed structure via depletion effects[26−28] and partitions itself into the resulting crystal (cf. top row of Figure ). Since the entropy cost of polymer confinement is strongly dependent on how complementary the polymer morphology and internal void structure of a crystal are, the polymer’s average shape offers a simple design parameter which can be used to tune the relative stability of one crystal over other competing forms. Thus, the efficacy of this approach hinges on the ability of the polymer to invade the interstices of the crystal phase without destabilizing it. In this approach to colloidal crystal assembly, the polymer cosolute, whose architecture and conformation can be finely controlled, effectively behaves as a SDA. However, its control is entirely entropically driven rather than enthalpically.

Figure 1

Schematic of the assembly paradigm in which a suspension of otherwise repulsive colloids is connected to a polymer reservoir where the polymer’s morphology is selected such that, when partitioned into the colloidal crystal at sufficient osmotic pressure, its morphology is commensurate with the local void symmetry found exclusively in one polymorph or another. (Top) This assembly route is purely based on the entropy loss of the polymer when adsorbed in a given crystal and has been previously shown to strongly bias the formation of one of the competing cubic and hexagonal polymorphic forms of close-packed crystals. (Bottom) By introducing thermal interactions between the colloids, it is plausible that a combination of colloidal energy and polymer entropy in the resulting crystal phase could be used to produce much more complex crystal morphologies in an intelligent fashion. In this work, we consider tetrastack crystals which are akin to close-packed ones in the sense that they have an identical topological distribution of octahedral voids (OV) and tetrahedral voids (TV) between different polymorphs.

This technique seems thermodynamically versatile as it depends only on the conformation of the polymer, and its kinetics are controlled simply by the polymers’ bulk concentration (degree of partitioning into the crystal phase), offering what would appear to be complete control over the crystallization and annealing process. However, the generality of this SDA approach hinges on its ability to be deployed to bias the formation of complex, arbitrary crystal structures into which there is no guarantee the polymer will partition without destabilizing the structure. Previous work has only demonstrated control over the relative stability of close-packed crystals.[23,25] In this work, we show that this approach can, indeed, be employed to stabilize non-close-packed structures, thereby suggesting that this method could provide the foundations of a different paradigm for the assembly of colloidal crystals. Due to the entropic nature of polymer-induced depletion interactions in these mixtures, in previous work, the polymer SDA was responsible for inducing both crystallization (which exclusively results in close-packed morphologies) as well as polymorph selection. Here, we decouple the driving forces for each phenomenon, combining an entropic polymorphic selection mechanism with an energetically driven crystallization one. The latter is in the spirit of the second paradigm previously discussed but is often incapable, on its own, of producing colloidal crystals with sufficient long-range order to be technologically relevant. This hybrid approach amounts to a linear combination of the two factors contributing to the free energy of the state of the system, energy and entropy, where their relative contributions can be modulated via temperature (cf. Figure ). This work demonstrates the fruitfulness of this combined approach in producing open colloidal crystals of a single polymorphic form. We illustrate how observations from close-packed systems can be easily extended to understand the assembly of more complex crystals, opening alternative avenues for engineering assembly pathways and experimental protocols to produce large-scale colloidal crystals relevant for a variety of applications.

Results and Discussion

Open colloidal crystals have lattice constants comparable to the wavelength of light, making them good candidate materials for photonic applications.[29−31] In this work, we consider the crystallization of triblock Janus colloids, that is, colloids decorated with two circular polar patches. Experimental work has demonstrated that the synthesis of these building blocks can be precise enough to make the self-assembly of crystals possible.[12] However, this system is a prototypical instance of a crystallizing colloidal system whose resulting morphology is difficult to control via patch chemistry alone. Indeed, self-assembly is strongly affected by defects arising from polymorphism between cubic and hexagonal forms and has degenerate noncrystalline structures (clathrates).[16] Romano and Sciortino[16] have shown that specific polymorphs can be stabilized by finely tuning the shape of the polar patches, but patterning the surface of the particles with the required degree of precision is currently beyond experimental capabilities. Therefore, this system is an attractive platform on which to rigorously test the merit of polymers as an experimentally realizable SDA paradigm. In the pure colloidal system, each colloid has two spherically symmetric patches placed on the poles. If the size of the patches is chosen in such a way that each particle can have a maximum of six bonded partners, at low temperature, these colloids crystallize into tetrastacks, depicted in Figure . These crystals exhibit a cubic and hexagonal polymorphism analogous to that of hard spheres that form close-packed structures. The interstices in these close-packed crystals can be described as being either octahedral voids (OV) or tetrahedral voids (TV), which can be visualized by connecting the centers of mass of neighboring colloids. The edge length of these platonic shapes is the diameter of a single colloid. In the fcc crystal, all OVs are entirely surrounded by TVs and vice versa, whereas in the hcp, the OVs share faces and stack in columns (cf. Figure ). In fact, the voids of the tetrastack polymorphs have a topology identical to that of their close-packed counterparts; however, the edge length of these polyhedra is now twice as long. Because the OVs enclose roughly 6 times the volume as the TVs, in the hexagonal polymorph, an adsorbed linear polymer of sufficient size will spread itself between pairs of cavernous OVs. By contrast, in the cubic polymorph, the same polymer will experience more confinement as it must spread between an octahedron and a tetrahedron. In close-packed crystals, this confinement can strongly bias the formation of the hexagonal phase by amplifying the free energy difference between the two polymorphs by 3 orders of magnitude over a system of pure colloids.[23,32,33]

Stability Predicted from Monte Carlo Simulations

Polymer SDAs bias the formation of a specific crystal polymorph by having a lower specific free energy, or chemical potential, when confined in that crystal over competing forms.[23] Thus, we begin by considering the total difference in the free energy, Δμ/kBT = (μcub – μhex)/kBT, of a linear, fully flexible homopolymer confined in each of the two tetrastack polymorphs that form spontaneously at low temperature. The total difference in the free energy of a single confined polymer between the two crystals plateaus to roughly 0.1 kBT, in favor of the hexagonal polymorph, but decreases by up to an order of magnitude as the colloid to monomer size ratio increases (cf. Supporting Information Figure S1). In close-packed crystals, this difference is much larger, Δμ/kBT > 2.[23] While it is plausible that such a difference could yield a detectable shift in the relative amount of the hexagonal polymorph produced in an experimental realization, we would like to improve upon this. In order to amplify Δμ/kBT, the polymer must be able to better detect differences between the local interstitial environments in the crystals. We achieved this by making the polymer more compact, that is, by converting the linear polymers into star polymers. Indeed, Figure illustrates that as we increase the functionality, or number of arms, on the polymer from f = 2 (linear homopolymer) to f = 8, the difference in free energy of an adsorbed polymer between the polymorphs is amplified by more than an order of magnitude. This maximum difference is achieved near q = 2Rg/σc ≈ 1.3, where Rg is the star polymer’s radius of gyration and σc is the colloid diameter. Thus, q quantifies the relative size of the polymer and the colloids. In close-packed crystals with confined linear homopolymers, this maximum has been found at q ≈ 0.71.[23] Again, this difference is due to the distribution of the polymer between neighboring voids, which in the hexagonal polymorph corresponds to a pair of large neighboring OVs, while in the cubic polymorph, it is always a set of one large OV and one small TV (cf. Supporting Information Figure S2). Compared to close-packed crystals, tetrastacks have an identical topological distribution of voids but the edge length of each void, and therefore distance between the centers of each, is twice as large. Consequently, in tetrastacks we a priori expect this maximum difference to occur at q ≈ 2 × 0.71, which is indeed observed.

Figure 2

Relative stability of tetrastack polymorphs in the presence of star polymers as predicted by Monte Carlo simulation. (a) Depiction of the cubic tetrastack crystal color-coded to illustrate the “ABC” stacking pattern. Each layer comprises tetrahedral sets of colloids packed in a hexagonal arrangement stacked out of the plane of the image in an alternating fashion, which is successively repeated to produce the macroscopic crystal. (b) Hexagonal polymorph is shown color-coded to illustrate the “ABAB” stacking pattern out of the plane. The visible gaps in the crystal correspond to octahedral voids which are stacked out of the plane shown. (c) Difference in the chemical potential of a given polymer confined in the two polymorphs, where Δμ = μcub – μhex, thus a positive Δμ corresponds to a lower free energy (higher stability) of the polymer when confined in the hexagonal polymorph. These simulations were performed with full monomeric resolution so that each point corresponds to a star, with a total of f arms, where the length of each arm, Marm, is increased by one monomer. (d) Δμ/kBT for coarse-grained high functionality star polymers as the overall size of the star, q, increases. Both the local maxima and minima in q are observed for all stars, as in the case of fully detailed simulations of low functionality stars in close-packed polymorphs.[25]

Since computationally it is intractably expensive to accurately sample all monomeric degrees of freedom of large star polymers when confined in a crystal, we now turn to a coarse-grained model, which treats the star polymers as a single soft sphere. Qualitatively identical results are obtained with this model, but the magnitude of the effect is relatively underestimated. Comparing the maximum in Δμ for f = 8 from the detailed simulations with that of f = 10 from the coarse-grained model, it is apparent that the coarse-grained model underpredicts the difference by roughly a factor of 2. This is not entirely unexpected. Even in close-packed crystals, simulations with explicit monomer detail have shown that a star’s corona deforms significantly and anisotropically under such confinement due to simultaneous interaction with multiple colloids.[25,34,35] The coarse-grained model struggles to quantitatively capture this effect because it was derived in the limit that the corona has a homogeneous angular monomer distribution. This causes the maximum in Δμ to shift to slightly higher q, an effect which increases with f. Coarse-grained stars do not fully account for volume conservation of the monomer corona when invaded by overlapping colloids; thus, when confined in a crystal, they must be made marginally larger to reach a size that fully detailed monomer simulations already account for. This effect is rather small, and the qualitative behavior of Δμ as a function of q is well captured, displaying a prominent local maximum where the hexagonal polymorph is strongly favored due to neighboring OVs, followed by a curious minimum at q ≈ 3. In fact, this minimum has been observed for star polymers confined in close-packed crystals, as well.[25] Extrapolating the minimum in Δμ at q = 3 reveals that Δμ becomes negative at f ≈ 8 (cf. Supporting Information). This is consistent with the results found for low functionality stars in close-packed crystals, where this change in the sign of Δμ corresponds to a translocation of the star polymer’s core from an OV to a TV.[25] As a star grows, eventually the entropic considerations of its corona outweigh that of its core, causing the core to shift from a large OV into a smaller TV if confined in the cubic polymorph. This is because all TVs in the cubic polymorph are surrounded by large OVs, meaning the corona gains more free volume when the core can position itself to access this pseudospherical symmetry (cf. Figure ). The hexagonal form does not have such a fortuitous arrangement of OVs and thus becomes less stable than the cubic form.

Thus, we find that observations from close-packed systems in the presence of linear homopolymer SDAs are reproduced in tetrastack systems with star polymers. The qualitative maximum in the bias for the hexagonal polymorph, Δμ, occurs when the SDA spreads between pairs of neighboring voids, due to the increased accessible volume in the hexagonal polymorph due to its arrangement of large OVs. This is followed by a minimum, which can become negative at very low f, due to the competition between the core and corona entropy of the star polymer. This competition is minimized by the way the cubic polymorph surrounds each TV completely by large OVs. The magnitude of these effects depends on the properties of the SDA; however, the location of these effects is dependent only on the size and separation of the voids, which scales linearly with the edge length of the void polyhedron. This complete transferability demonstrates that this SDA paradigm for polymorph selection is not a specific effect observed only in simple close-packed crystals but can be fruitfully applied to more general systems.

Molecular Dynamics Realizations

The Monte Carlo calculations performed in the previous section are based on single-particle insertions and hence do not take into account many-body effects that could play a substantial role in determining the relative stability of bulk phases. Going one step further, we studied the crystallization of triblock Janus particles in pure systems as well as in mixtures with star polymers using constant-temperature molecular dynamics simulations. We modeled the colloid–colloid interaction with a continuous version of the model employed by Romano and Sciortino[16] (see Figure a), the star–colloid interaction with the coarse-grained potentials developed in ref (36), and the star–star interaction with the coarse-grained potentials derived by Likos and co-workers.[37] We analyzed the resulting configurations by classifying colloids according to their local symmetry, as determined by Steinhardt–Nelson bond order parameters (cf. Methods).[38] In what follows, we will focus on stabilizing the hexagonal tetrastack polymorph which, as predicted by Monte Carlo results, can be strongly favored by mixing colloids with large star polymers.

Figure 3

Molecular dynamics simulation snapshots of triblock Janus particles. (a) Triblock Janus colloid modeled as a particle with two spherical polar patches (in blue). When two patches of different colloids are close to each other, the two colloids feel a net attraction whose magnitude depends on their relative distance and orientations (cf. Methods for the functional form of the potential). (b–d) Simulation snapshots where only crystalline colloids are shown (cf. Methods). Particles are color-coded according to the local environment: green (red) particles are in a cubic (hexagonal) local environment, while yellow particles are in a mixed environment. Patches are depicted in blue. (b) Snapshot of a pure colloidal system. The majority of the particles are in a crystallite formed by randomly stacking hexagonal and cubic layers. (c,d) Snapshot of simulations of a binary mixture of patchy particles and star polymers with f = 30 and q = 1.5. In (c), the relatively small star–star repulsion makes the system undergo a demixing transition, which pushes the colloids close to each other, forming a close-packed crystal. In (d), the same simulation was performed with charged star polymers, which experience a higher mutual repulsion that stabilizes open structures and biases the formation of the hexagonal tetrastack polymorph. As a consequence, nearly all crystalline particles are in a locally hexagonal environment.

Simulating the pure colloidal system at low temperatures and moderate densities always yielded a polymorphic crystal, that is, a crystal with stacking faults, as shown in Figure b. The snapshot exhibits a clear lack of long-range order. As suggested by Monte Carlo results, star polymers could be used to select the hexagonal polymorph. However, the choice of q and f is a delicate matter. On one hand, Monte Carlo results clearly indicate that the higher the f the better, whereas an optimal q value should lie in the range [1.4,1.8]. On the other hand, the biasing provided by the star polymers relies on them being able to explore the void distribution of the crystals, meaning that the colloid–star binary mixture should not be prone to demixing. As noted in ref (36), the tendency to demix can be estimated by performing a mapping of the system to an effective nonadditive hard-sphere binary mixture. Under this assumption, the stability of the star–colloid system against demixing is qualitatively controlled by the nonadditivity parameter Δ, defined aswhere σ is the effective hard-sphere diameter between particles of species i and j, which we obtained through the Barker–Henderson approximation.[39] Negative values of the nonadditivity parameter favor mixing between different species, while positive values of Δ may indicate a tendency to demix, provided that the densities of both species exceed certain thresholds that are, in turn, Δ-dependent.[40,41] We note that the results contained in ref (36) show that mixtures of colloids and stars with q ≈ 1.5 and f ≳ 20 always have Δ ≳ 0, indicating a possible instability with respect to demixing. However, this approach is known to be quantitatively inaccurate especially for relatively small f and cannot be relied upon to estimate the parameters for which demixing occurs.[35] A numerical check is thus required. We performed simulations of a mixture of colloids and star polymers with f = 30 and q = 1.5, for which Monte Carlo results predict a free-energy difference between the hexagonal and cubic tetrastacks, Δμ/kBT ≈ 3.5 per star polymer adsorbed. In principle, such a large value makes it possible to clearly detect the predicted entropic effect in simulations. Unfortunately, we observed that as soon as the density of the stars becomes appreciable, the system demixes and the colloids assemble into close-packed crystals, as shown in Figure c. The net effect of the star polymers is thus to suppress the formation of low-density phases characteristic of patchy particles. Increasing f only further increases Δ, while increasing q provides only a marginal reduction.[36] In all cases, the energetically degenerate clathrate structures were also suppressed in favor of one or more crystals.

Therefore, to observe this SDA-driven polymorphic biasing, one must have a system in which Δμ/kBT is of the desired sign and the nonadditivity parameter is sufficiently small to avoid demixing and subsequent collapse of the low-density crystal phase. Since σsc and σcc are fixed, the former by the polymorph free-energy difference one wants to enforce and the latter by the size ratio q, the only way to enhance the tendency to mix of stars and colloids is to increase the star–star interaction (cf. eq ). It is advantageous that the class of soft colloids is very rich, and it includes a variety of polymer-based colloids, offering flexibility in the control of their interactions. In particular, by employing polymers with ionizable units and changing the pH of the solution, one can control the charge on the polymer chains, creating thereby polyelectrolyte stars.[42] By charging the stars, while at the same time keeping the colloids neutral, it is possible to enhance the star–star repulsion without significantly affecting the other interactions. The interstar repulsion in the presence of charge also exhibits a Yukawa tail at non-overlapping distances, which is controlled by salinity and can be thus tuned to yield a Yukawa decay length the same as that of the neutral stars. The amplitude of the repulsion is controlled essentially by the overall valency Z of the star.[43] Accordingly, charging up the stars is equivalent to multiplying their mutual repulsion by a factor of α, which depends solely on the charge Z and the star polymer’s diameter, and can be used to tune Δ. In the q = 1.5, f = 30 case, the nonadditivity parameter, which is a monotonically decreasing function of α, takes the value Δ ≈ 0.22 in the neutral case and quickly drops to zero for α ≳ 7.8. As an example, we considered mixtures with α = 6.5, for which Δ ≈ 0.016 and which can be realized by employing stars that are polyelectrolytes containing a few hundred elementary charges, depending on the diameter (cf. Supporting Information). Under such conditions, we successfully observed the formation of hexagonal tetrastack crystals as predicted by Monte Carlo calculations. A typical snapshot of the final configuration of a crystallized sample is presented in Figure d, which shows a crystallite that does not contain any randomly stacked layers, with nearly all the particles in a locally hexagonal environment.

In order to prove that the hexagonal polymorph is indeed the most stable one and not merely a metastable intermediate, we also performed annealing simulations. These simulations were started from configurations in which all the colloids were initially on a perfect cubic tetrastack lattice and the star polymers were randomly placed throughout the sample. The results, as well as representative snapshots, are shown in Figure . During the full course of the simulation, we observed complete annealing of the cubic tetrastack into the hexagonal one. Indeed, ⟨q4⟩ decreases monotonically to the value associated with the hexagonal polymorph, and its initial rapid descent demonstrates that star polymers make the cubic tetrastack very unstable with respect to the hexagonal one because no nucleation event is required to start the annealing process. As shown in the snapshots, the hexagonal polymorph grows directly on the surface of the cubic crystallite, which slowly but steadily acquires a hexagonal morphology. The final crystalline cluster, which has a long-range hexagonal symmetry, is stable and does not undergo any further evolution with time. Similar (partial or complete) annealing runs were obtained at different star polymer densities and with stars containing a different number of charges (shown in the Supporting Information), proving the robustness of this polymorph stabilization mechanism. As an additional test, we also performed annealing simulations in the presence of smaller stars (f = 30, q = 0.75) which, according to Monte Carlo results, should not affect the stability of either polymorph. The resulting ⟨q4⟩ time series, plotted as the green and red curves in Figure , are constant throughout the simulation, confirming that the hexagonal and cubic crystallites are not affected by the stars.

Figure 4

Evolution of patchy particle crystals in molecular dynamics simulations. The color coding is the same as that in Figure . Coarse-grained star polymers with f = 30, q = 1.5 and a star–star repulsion increased by a factor of 6.5 due to electrostatic repulsion can be employed to anneal a cubic tetrastack crystallite into a hexagonal one. The snapshots (each of which refers to the point in the plot indicated by the accompanying letter) demonstrate how the annealing process takes place: the simulation starts from an initial configuration where all colloids are arranged on a cubic tetrastack lattice and the star polymers (which are shown only in the top-left picture) are randomly distributed throughout the sample. As the simulation proceeds, particles continuously detach from and reattach to the crystallite, slowly annealing into a hexagonal crystal. The black line in the plot shows the time dependence of the average bond order parameter ⟨q4⟩, which starts from the cubic tetrastack value (green line, ⟨q4⟩c ≈ 0.513), decreases monotonically, and eventually plateaus to a value very close to the hexagonal tetrastack one (red line, ⟨q4⟩h ≈ 0.172). The small difference between the two values is due to a few (2–3) particles that are still in a mixed environment. The red (green) curve was obtained by simulating a hexagonal (cubic) crystallite with smaller stars (q = 0.75) which, as predicted by Monte Carlo calculations, do not affect the stability of either polymorph.

Conclusions

We have demonstrated the validity of an alternative approach for the directed self-assembly of high-fidelity complex colloidal crystals that exploits polymers as structure-directing agents to thermodynamically assign the most stable crystal structure from a suite of potential competitors. This can be achieved purely by engineering the polymer’s interactions with each other and with the colloids. Monte Carlo simulations revealed that correlations between different voids due to different symmetries between crystal structures are sufficient to amplify weak inherent free-energy differences between polymorphs by several orders of magnitude if the polymer is chosen wisely. We revealed how qualitative inferences from Monte Carlo simulations on simple, close-packed crystals are transferable to more complex ones and suggest how to make such choices. Molecular dynamics simulations fully support Monte Carlo predictions, validating the use of polymers as SDAs to entropically stabilize complex, open colloidal crystals on the basis of their internal void structure. The only caveat is that one must take care to avoid destabilizing open crystals with respect to their close-packed counterparts because, as the polymer density is increased, depletion effects tend to stabilize more compact structures. To this end, we illustrated that the most stable crystal structure can be tuned by two factors: the difference in polymer free energy between different crystal structure environments, Δμ/kBT, and the nonadditivity parameter, Δ. The former reflects the relative amount of free volume a polymer can access inside each crystal, while the latter reflects the capacity of a polymer to invade the interstices of the crystal phase in the first place. For close-packed crystals, the latter could be effectively controlled by polymer density (osmotic pressure); however, additional polymer–polymer interaction tuning is necessary to stabilize low-density crystals over their close-packed analogues. Both Δμ/kBT and Δ can effectively be tuned by changing the polymer size and charge, for a given colloidal system, neither of which affect the intrinsic colloid–colloid interaction. This opens additional pathways to engineering colloidal crystals as this SDA-based mechanism can be deployed, not only independently but also in parallel to other paradigmatic approaches which rely on energetic interactions between the colloids themselves or between colloids and some external medium. Although this approach does not solve all colloidal crystal design problems, it serves to illustrate an alternative with which to supplement existing approaches in cases where other strategies might produce unsatisfactory results. We emphasize that the success of the design logic itself for this set of crystal polymorphs is as important as the results themselves, and we expect the heuristics outlined here will serve as a guide to experimental realizations of this phenomena in the near future.

Methods

Monomer-Resolved Monte Carlo Simulations

Our simulated star polymers were composed of f linear chains of equal-sized monomer beads which at one end were bonded to a common central monomer. Each monomer bead had a diameter of σm = 1.0, and each arm contained Marm beads, each of which was bonded in sequence by the finitely extensible nonlinear elastic potential:[44]

We employed the well-known Kremer–Grest model, wherein k = 30ϵ, r0 = 1.5 were used to minimize bond crossings. A linear homopolymer is simply a star polymer with a functionality of f = 2. All beads also interacted with each other, and with colloids, through a shifted Lennard-Jones potential, which is cut and shifted at the minimum. For equal-sized beads this is known as the Weeks–Chandler–Andersen potential (WCA):[45]

This potential has been shifted by a factor, Δ = (σ + σ)/2 – σ, which accounts for the difference in diameters between dissimilar species such that the slope of the potential remains identical for all pairs as they begin to overlap regardless of their individual size. For the detailed MC simulations, colloids had a diameter of σc = 6.0σm and were treated as purely repulsive spheres without any patches since these patches do not interact with monomer beads.

For these simulations T* = kBT/ϵ = 1, however, since the interactions are purely and steeply repulsive, the difference in the polymer’s chemical potential between the two polymorphs was found to be invariant when normalized by kBT, hence we always express the difference as Δμ/kBT to be consistent with the temperature used in the molecular dynamics portion of this work, where kBT is reduced to crystallize the colloids when attractive patches are present. Colloids were held fixed for the duration of these simulations and were initialized on either the cubic or hexagonal tetrastack lattice. The former was generated in a periodic box with dimensions (L,L,L) = , which represents one complete “ABC” stacking pattern in the z direction containing 48 tetrastacks or 192 total colloids. The latter was generated in a periodic box with dimensions (L,L,L) = , which represents a pair of “AB” stacks in the z direction containing 64 tetrastacks or 256 total colloids. The additional spacing of 0.12 ≈ 21/6 – 1 between nearest neighbors was added because the colloids interacted via the repulsive WCA potential rather than as perfectly hard spheres. This additional gap results in crystals with zero internal energy and configurational pressure, as in previous work.[23−25]

To measure the chemical potential, or free energy, of a single star polymer confined in a given crystal, we employed incremental-growth Monte Carlo simulations[46] described in more detail elsewhere.[25] We briefly summarize the technique here. Polymers were grown incrementally by attempting to insert an additional “ghost” monomer onto a fully inserted portion of the star at the end of each arm, from arm 1 to f. The fully inserted portion was relaxed in the canonical ensemble between insertion attempts, and relaxation moves included local displacements of monomers, displacements of the entire polymer’s center of mass, and regrowth of all or parts of a randomly chosen individual arm via Rosenbluth sampling.[47] These moves typically occurred with a 8:1:1 ratio. The ensemble-averaged Boltzmann factor to insert an additional ghost monomer on the end of an arm was then used to calculate the incremental excess chemical potential:

Monomer insertions were generally attempted after several thousand relaxation moves. Between 1 × 103 and 3 × 107 relaxation moves were performed for each monomer fully inserted in the system before the polymer was appended with an additional monomer. All simulations were repeated between 10 and 40 times to obtain reliable statistics. At the end of the relaxation period, an additional bead was formally inserted on the end of the arm being grown, and the process was repeated for the next monomer on the next arm. The star polymer’s core was similarly added via a series of test insertions. The total chemical potential of the star polymer is given bywhere μcoreex and μex are the incremental chemical potentials of the core and of the ith monomer of the jth arm, respectively.

Coarse-Grained Monte Carlo Simulations

A Widom insertion scheme is used to obtain the chemical potential of a single coarse-grained star inside these polymorphs using previously validated colloid–polymer potentials. A star is modeled as a single soft sphere where the form of the interparticle potential is dependent on the relative size of the colloids and polymers, q = 2Rg/σc.[36,48] The size of a star under dilute conditions is estimated from f and Marm using previously validated scaling on a fine lattice.[34] The effective force between a star and a colloid is determined by considering the osmotic pressure exerted by the star on the colloid’s surface to obtain a mean force, F(r), at discrete interparticle separations, r, which is accurately splined and integrated from infinitely far away to r. Again, these are purely repulsive interactions. These MC simulations were performed at kBT = 0.155, where crystallization of the patchy colloid counterparts in the MD simulations was observed. This coarse-grained potential allows us to consider stars at much larger q than in the detailed monomer MC simulations but effectively capture identical trends.

Molecular Dynamics

Triblock patchy particles were modeled as nearly hard-sphere particles decorated with two polar patches shaped as truncated spherical cones. The core–core repulsion was given by[49]The patch–patch attraction was provided by a continuous version of the Kern–Frenkel model.[50] It is a square-well-like attractive potential modulated by an orientation-dependent function. Namely, the interaction between patch i on particle α, identified by the versor α̂, and patch j on particle β, identified by the versor β̂, was given bywhere r⃗ = r⃗α – r⃗β, r = |r⃗|, r̂ = r⃗/r, and Ω is a steep modulating function that takes into account the orientation of a patch with respect to the versor connecting the particles centers. It has the form of a generalized Gaussian function:Star–star interactions were modeled through the coarse-grained potentials developed in ref (37). Using the original potential resulted in a collapse of the open crystal into a close-packed one. In order to overcome this limitation, we enhanced the penetration of the crystallite by the star polymers by rescaling the potential by a factor α = 6.5. We note that this can be realized by employing charged star polymers. Depending on the size of the star, this enhanced repulsion requires a few hundred elementary charges per star or, equivalently, tens of charges per arm (cf. Supporting Information). A possible realization of this system could be provided by triblock patchy particles fabricated by glancing angle deposition, which can be used to produce patchy particles with diameters of down to 10 nm.[51] For example, σcc = 65 nm would require star polymers with σss ≈ 100 nm and roughly 10 elementary charges per arm. These numbers are perfectly compatible with experimental realizations of star polymers. We used the same star–colloid interactions as in the Monte Carlo simulations.

We ran constant-temperature molecular dynamics simulations on GPUs with the simulation package oxDNA.[52] All simulations were performed at a reduced temperature T* = kBT/ϵ = 0.155, enforced through a modified Andersen-like thermostat,[49] for at least 109 time steps of length .

During the analysis, we classified crystalline colloidal particles according to their local environment by using the well-known Steinhardt–Nelson bond order parameters.[38] First, following Romano and Sciortino,[16] we labeled as crystalline all those colloids having exactly six bonded neighbors (two particles are considered bonded if their pair interaction energy is less than −0.1 ϵ). Then we described the local structure around each crystalline colloid i with the quantitywhere N(i) is the number of i’s neighbors, r̂ is the versor connecting colloid i and its jth neighbor, and Y is the lth degree spherical harmonic of order m. Each colloid was then assigned a locally averaged bond order parameter of order l given bywherewhere N̅(i) is the number of colloids in the neighborhood of i, including itself. Both q̅4 and q̅6 were able to distinguish between hexagonal and cubic local environments, yielding equivalent results. In order to follow the annealing to the hexagonal polymorph, we used the globally averaged bond order parameter ⟨q4⟩, defined aswith Nc being the total number of crystalline colloids.