Grafted nanoparticles as soft patchy colloids: Self-assembly versus phase separation

Controlling the dispersion of nanoparticles in polymer melts is a challenging task, but one which promises advances in fields ranging from biomedical technology to nanoscale electronics.¹ The chemical dissimilarity between the nanoparticles and the polymer medium makes it difficult to disperse the nanoparticles homogeneously. This challenge is often overcome by covalently grafting polymeric components to the surface of the nanoparticles.^1,2 Although these polymer-grafted nanoparticles (GNPs) can be homogeneously dispersed in certain cases, the complex interplay between the thermal interactions of the cores and the entropy of the polymer coronas can also lead to unexpected self-assembly into anisotropic morphologies such as strings and sheets.^1,3,4 Because of these competing forces, GNPs are sometimes considered “nanoparticle amphiphiles” since their self-assembly can be interpreted as microphase separation.⁵ Temperature allows one to control the effective magnitude of the core-core attractive interactions, whereas the grafting density and length of the polymers in the corona control the entropically driven repulsion.

Recent theoretical work has shown that in the limit of strong core-core attractions, anisotropic morphologies of GNPs can be qualitatively predicted by considering only the packing of the nanoparticle cores and the subsequent entropic deformation penalty of their coronas.⁴ Early simulations qualitatively supported this theory, but their quantitative predictions differed by an order of magnitude from it. As this theory neglects other entropic contributions such as the translational entropy of the GNPs, this discrepancy is not surprising. More recent coarse-grained simulations⁶ treating the grafted polymers as Gaussian springs qualitatively confirm these morphological predictions but neither do they describe the corona in detail nor do they consider the full phase behavior of these systems in the thermodynamic limit. To our knowledge, no prior study has yet been undertaken to systematically study the phase behavior of GNP systems with an explicit corona. This is a significant shortcoming in light of another paradigm with which to consider GNPs, namely, as “patchy particles.”

So-called “patchy particles” are colloids with an inhomogeneous surface chemistry that provides anisotropic thermal interactions.^7–9 They have emerged as a natural extension to simpler two-sided “Janus” particles, serve as a useful model for various physical systems such as clays,^10,11 and are an attractive platform with which to produce different crystal lattices for catalytic and optical applications.^12,13 Sparsely grafted nanoparticles between the “mushroom” and “semi-dilute polymer brush” (SDPB) limits have clear analogies to such particles. In the mushroom limit, grafted chains are short enough, and their grafting density low enough, such that spheres swept out by their radii of gyration on the surface of a GNP do not significantly overlap.¹⁴ Whereas at higher grafting densities and chain lengths in the SDPB limit, chains overlap sufficiently to entirely cover the surface of the GNP core. However, between these limits, the chains will instead shield only certain regions of the surface, while isolating patches of unshielded surface area from one another. This means that numerous distinct portions of a GNP’s surface remain exposed, and are available to interact with other cores, while other portions are effectively shielded. Thus, isotropically decorated GNP surfaces could appear to have multiple isotropically placed surface patches. The coronas provide soft, deformable, repulsive patches while the unshielded surfaces attract one another. This immediately raises the question as to where on the spectrum between amphiphiles and currently accepted models of patchy particles the GNPs lie (cf. Fig. 1). Such an inquiry is made all the more intriguing as the thermodynamic properties at either ends are substantially different.

Amphiphilic self-assembly is a morphological transition which exhibits no discontinuities in any derivatives of the free energy. These transitions may sometimes appear first-order in simulations of a finite size, but in the thermodynamic limit, they are, in fact, continuous.⁵ On the contrary, patchy particles have been shown to undergo a true first-order phase separation between a low- and high-density fluid phase, separated by a critical point, when cooled sufficiently.¹⁵ When these patches are attractive, the equilibrium volume fraction of the high-density fluid decreases as the number of patches decreases, collapsing the two-phase region resulting in an “empty liquid” phase at densities above coexistence. Empty liquids are disordered phases with very long relaxation times, typically on the order of experimental timescales, which appear as a connected network dispersed in large regions of pure solvent.¹⁰ Despite their long relaxation times, they are ergodic on sufficiently long timescales and are thermodynamically stable, unlike non-equilibrium phases such as glasses or gels. Therefore, these empty liquids form equilibrium gels when cooled rather than phase separating. Clays including laponite have recently been shown to produce these empty liquid phases.^10,11 Empty liquids tend to have both structure and a relatively low density, which suggests the intriguing possibility that perhaps the anisotropic phases, most notable the string-like morphologies observed in GNP systems at low volume fractions,³ may in fact be low-density, yet GNP-enriched, empty liquids.

It is therefore important to consider carefully the behavior of GNP systems in the thermodynamic limit with a detailed description of the polymer corona in order to resolve the dichotomy presented in Fig. 1. In this work, we evaluate the thermodynamic behavior of model GNPs to establish if these nanoparticles can be engineered to behave as patchy particles, and thus provide a route to empty liquid phases. To this end, we present a lattice model of GNPs which reproduces the predictions of the previously developed mean-field treatment⁴ of GNP systems with a much higher degree of quantitative accuracy than previous models. We then employ extensive computer simulations to analyze the model’s true first-order phase transitions and micellar self-assembly behavior for a wide variety of GNPs, carefully investigating the dependence of the apparent location of phase boundaries on simulation system size, in order to properly distinguish between the two.¹⁶ 

The rest of this paper is structured as follows. In Sec. II, we describe our simulations and data analysis methodology. In Sec. III, we discuss our results comparing and contrasting the different types of thermodynamic behavior for a wide range of GNPs and offer conclusions about the validity of the patchy particle versus amphiphile paradigm in Sec. IV.

We performed grand canonical Monte Carlo (GCMC) simulations on a fine cubic lattice to obtain the thermodynamic properties of model polymer-grafted nanoparticles. Each lattice site is a cube with edge length σ, and the simulation box consists of L × L × L lattice sites. The lattice coordination number is z = 26. That is, a site’s neighbors consist of its 6 nearest neighbors separated by a distance σ, twelve next-nearest neighbors each at $2 σ$⁠, and 8 oblique neighbors each $3 σ$ away. This model has been shown to reproduce experimental thermodynamics of athermal colloid-polymer mixtures nearly quantitatively.^17–20 In these simulations, the solvent is implicitly represented by an unoccupied monomeric lattice site. A GNP consists of a core of sites which exclude all other entities, onto which f arms are isotropically grafted. Each of the polymer arms is composed of a flexible chain of M monomers, each comprising a single lattice site, connected according to the lattice’s coordination. When a site is occupied by a monomer, it also excludes all other entities from existing in that site. Thus, the monomers (or Kuhn segments) are fully excluded from one another, mimicking the good-solvent limit. All lattice sites within σ_c/2 of the center of the GNP core are considered part of the core. We fix the core diameter, σ_c = 5σ, for all simulations since at a two-body level this asymmetry has already been explored in some depth.²¹ Nanoparticle cores encompass 81 lattice sites, so we define the volume fraction of GNPs in a system as ϕ_GNP ≡ (81 + fM) N/L³, where N is the number of GNPs present. The location and relative orientation of each grafting point are fixed by a template at the beginning of a simulation regardless of the arm’s length, as depicted in Fig. 2.

Lattice site exclusion implies that all interactions in the system are purely athermal (repulsive) with the exception of the core-core attraction. To mimic the dispersion forces between the cores, we supplement the hard-sphere repulsion imparted by this exclusion with an attractive tail of the Yukawa form

We fix ϵ/k_B at unity so that the temperature reported throughout this text can be taken as the reduced temperature, T^∗ = k_BT/ϵ = 1/β. The range of interaction was fixed at κ = 0.4, r_cut = 2.5σ_c. The simulation box size was maintained well above 2r_cut and above three times the total size of an infinitely dilute GNP (⁠$σ c +4〈 R g,arm 2 〉 1 / 2$⁠).

Monte Carlo moves consist of insertion and deletion attempts for the entire GNP, displacement of the entire GNP by a single lattice site, displacement of a monomer on a polymer arm subject to connectivity constraints, regrowth of some fraction of a single arm, regrowth of some outer fraction of the entire corona, and rotation of the GNP. To accelerate the simulation, single polymer arms are always grown and moved via Rosenbluth sampling,²² whereas the rest of the moves are accepted with the standard Metropolis criterion. Rotation moves attempt to rotate the entire GNP by either ±90° or 180° in one of the three principle cartesian planes commensurate with the cubic symmetry of the underlying lattice. All moves which modify the corona, including insertion and deletion steps, also employ Rosenbluth weights to accelerate the sampling of important regions of phase space. The corona’s weight is assigned to be the product of the individual weights of each arm as described in more detail in Ref. 19,

After an equilibration period (of at least 1 × 10⁹ steps) to allow properties to converge, a histogram of total energy and number of GNPs, h(U, N), was collected during a production phase lasting at least as long. Standard reweighting techniques can then be applied to obtain the critical points, binodals, and pressure isotherms.²³ After patching together overlapping histograms from different simulations in a self-consistent manner,²⁴ the underlying entropy surface can be extracted. This technique allows us to perform simulations at high T^∗, where ergodic sampling is easily achieved, to obtain information about the system at lower T^∗ where direct simulations generally suffer from poor sampling. Critical points were identified by optimizing the mixed-field order parameter, m = N − sU, such that the probability distribution, P(m), at critical conditions was fit to the universal three dimensional Ising distribution. The temperature, chemical potential, and field mixing parameter, s, were each adjusted during the multidimensional optimization to yield the best possible agreement. Subsequent simulations were then performed at the identified critical point to verify the results of reweighting. Errors at the critical point were obtained by taking the standard deviation of three separate simulations with different random number sequences. Tabulated values and details are relegated to the supplementary material.²⁵ Binodals were then computed using the equal peak weighting criterion to determine coexistence conditions at temperatures below the critical point. Pressure on the lattice can be computed from the partition function as follows. At each state point (β, μ), the values in the histogram provide the unnormalized probability of a given state according to

The entropy, $S= k B ln Ω ( U , N , V )$⁠, is only known up to some run-specific constant, c,

Similarly, the partition function, which is obtained by the sum,

is only known up to some unknown multiplicative constant, C, related to the reference state of the system. Once the entropy surface Ω(U, N, V) has been obtained by self-consistently patching multiple histograms together to cover as much of phase space as possible, this uncertainty can be calculated by making use of the ideal gas reference state at low density, ρ = N/V. Since the grand canonical partition function is related to the pressure of a system by

and in the ideal gas limit Pβ → ρ, we have the condition that

When this is satisfied up to approximately 1% error or less, C′ is known and Eq. (6) is completely defined. Pressure isotherms are a key component of our analysis of these systems because they reflect the system’s finite-size dependence away from any critical point. We emphasize that our ability to precisely determine whether or not a system undergoes a true first-order phase transition, based on finite-size scaling, is the key feature distinguishing our work from previous investigations of similar systems.

We began by computing the phase behavior for the reference case of f = M = 0. These bare attractive colloids phase separate below a critical temperature of $T c ∗ ≈1.08$ into an isotropic GNP-lean phase (gas-like) and an isotropic GNP-rich (liquid-like) phase as depicted in Figs. 3 and 5. We then grafted the surface with athermal polymer arms to shield the cores from one another. Figure 3 illustrates all combinations of f and M we explored. When the surface of the GNP is entirely shielded by its corona, which corresponds to the transition from the mushroom to the semi-dilute polymer brush limit, the GNP should not appear patchy. This limit is found by estimating the projected surface coverage from a single grafted chain by using its ungrafted, infinitely-dilute radius of gyration, $〈 R g 2 〉 1 / 2 =0.508 M 0 . 588$⁠, as previously determined for this lattice.^17,20 The crossover to the SDPB limit is indicated by the magenta curve on Fig. 3 using

Below this limit, we expect to see signatures of patchy particle behavior, if it is to be found anywhere. Overall, we observed three basic types of thermodynamic behavior upon cooling each homogeneous system from a reference temperature of T^∗ = 1.08: first-order fluid-fluid phase separation, continuous self-assembly into planar sheets, and continuous self-assembly into string-like morphologies. We compared the types of thermodynamic behavior observed in our simulations to the predictions of the mean-field theory in Ref. 4 which has previously been shown to rationalize this behavior based on the geometry of different GNP morphologies. As we illustrate in Fig. 3, the present simulation model produces results which agree with the theory quite reasonably. Since GNPs at moderate f and M were found to undergo either phase separation or self-assembly over a range 1/10 ≲ T^∗ ≲ 1/3 (cf. supplementary material²⁵ and Fig. 5), the background of Fig. 3 was generated from theory at T^∗ = 1/7. Discrepancies between the theory and simulations come partly from the fact that the theory predicts density-independent behavior at a fixed T^∗, whereas the simulations are classified based on the nature of their observed behavior upon cooling from a homogeneous fluid state; morphological and phase transitions occur at different temperatures for different combinations of f and M.

Near and above the SDPB limit, we generally find that these GNPs self-assemble in a thermodynamically continuous fashion to produce either string- or sheet-like morphologies. This self-assembly is reflected by the characteristic shape of the pressure isotherms which rise quickly at volume fractions where these morphologies begin to form and plateau shortly thereafter. Figure 4 illustrates these isotherms for progressively lower temperatures obtained for a GNP with f = 4, M = 15 which forms lamellar sheets. From simple geometry, we can calculate the volume fraction where n sheets form in a cubic simulation box,

In general, the effective surface area of each sheet will be slightly higher than what is expected in an idealized geometry due to fluctuations. Usually, agreement between this calculation and simulation results was improved by increasing the effective surface area, and therefore volume fraction, by approximately 5% but never more than 10%. The three vertical lines in Fig. 4 correspond to n = 1, 2, and 3, respectively. At higher temperatures, translational entropy inhibits low density systems from forming low n structures. However, at higher density, the entropy difference between an amorphous morphology and an ordered sheet-like structure is reduced, allowing higher n structures to form. As T^∗ is reduced, the effects of translational entropy are as well stabilizing progressively lower n structures. Of course, as the box size increases, the volume fractions at which these lamellae form will decrease (cf. Eq. (9)). Similar geometric considerations do not apply for systems which form strings, but their assembly is still visible in the pressure isotherms which show no signatures of first-order phase separation (cf. supplementary material²⁵).

However, well below the SDPB limit, poorly shielded GNP cores undergo first-order fluid-fluid phase separation. This regime must be carefully examined for signatures of novel behavior. Figure 5 illustrates representative binodals we obtained for select (f, M) indicated as undergoing phase separation in Fig. 3. A table with numerical values for the critical points of all phase-separating systems, as well as their respective binodals, can be found in the supplementary material.²⁵ For f ≳ M, we exclusively observed first-order phase separation into isotropic fluid phases for GNPs that approach and surpass the SDPB limit starting from low f (cf. Fig. 5(a)). This is easily rationalized as the corona in this portion of (f,  M) space is composed of many short chains which cannot be deformed significantly from their grafting points. Thus, GNP-GNP interactions are expected to be soft but radially isotropic. As the corona becomes more dense, this phase separation is delayed to progressively lower temperatures but the basic shape of the binodal curve does not change. However, we expect GNPs to mimic patchy particles best when the corona is below the SDPB limit and composed of sparse, highly deformable chains as they can more easily stretch away from the surface leading to poorer shielding.⁶ 

Although we considered numerous values of f (cf. supplementary material²⁵), we focus here on the case where f = 4. For an isotropically decorated GNP, f = 4 corresponds to a tetrahedrally coordinated patchy particle. Although the polymer grafts repel other GNP surfaces, the self-dual nature of a tetrahedron implies that the remaining exposed (attractive) surface patches are also arranged in a tetrahedral fashion. Within the patchy particle paradigm, increasing M at a given f should eventually create grafted “mushrooms” that overlap enough to isolate regions of the exposed surface into independent patches. Indeed, Fig. 5(b) illustrates that in addition to lowering $T c ∗$⁠, increasing M at f = 4 seems to cause the high-density portion of the binodal to collapse, apparently commensurate with patchy particle behavior. However, this conclusion is incorrect for the following reasons. As the temperature decreases, it is clear from Fig. 6(a) that the pressure appears to increase quickly at volume fractions just above the computed GNP-rich branch of the binodal around ϕ_GNP ≈ 0.30. Moreover, the crossing of these isotherms illustrates the anomalous sign of the system’s isobaric expansion coefficient,

However, as the box size is increased, this anomaly begins to disappear as shown in Fig. 6(b). Furthermore, at no condition does the order parameter distribution for this system match the distribution for the three dimensional Ising model at its critical point. The most optimized agreement we obtained for each simulation box size is presented in Fig. 6(c), but it is clear that deviations become more pronounced as the box becomes larger. Even at L = 30, for which the “binodal” was computed in Fig. 5(b), the Ising distribution cannot be perfectly recovered implying this is not a true critical point. In supercooled water, a first-order liquid-liquid phase transition, metastable with respect to crystallization, is often linked to anomalies including negative expansion coefficients.²⁶ However, in this case, the anomaly emerges as a result of competition between phase separation and self-assembly. The order parameter distribution cannot be mapped onto the Ising model since the high-density phase for these GNPs is no longer an isotropic fluid, but is instead a self-assembled phase of planar sheets.

Vertical solid black lines in Fig. 6 show the predictions of Eq. (9) for n = 1, which scale as ∼1/L. Clearly, the GNP-rich phase detected at L = 30 by histogram reweighting is actually the point at which the system has self-assembled into a single lamellar sheet. Without carefully considering the order parameter distribution and finite-size effects, it is clear how such a transition could easily be mistaken as being first-order, suggesting an erroneous analogy to patchy particle behavior. The system’s isobaric expansion coefficient anomaly, which also occurs at this point, is simply a reflection of the elasticity of the sheet as it buckles. Since a sheet is stabilized by energetic interactions, cooling makes it more pliable, whereas heating entropically favors dissolution of the sheet, making it more brittle. As the box size increases, the anomalous expansion behavior decays owing to the fact that longer wavelength fluctuations are permitted in the system allowing it to progressively relax elastic deflections.

As Fig. 3 illustrates, we investigated other combinations of (f,  M) in detail as well, however, in all cases, we found unambiguous first-order phase separation or continuous self-assembly. At no point did we obtain a trend where the modulation of either f or M led to a collapsing high-density volume fraction for the binodal as would be expected during the formation of empty liquids or equilibrium gels. Therefore, we conclude that GNPs comprise an entirely separate class of materials from patchy particles. The latter are prone to form novel empty liquid phases, whereas after extensive searching, we found no evidence that our model GNPs can exhibit the same phenomenon. This supports the classification of GNPs wholly as nanoparticle amphiphiles, and the conclusion that the anisotropic morphologies which are found in these systems are the result of continuous self-assembly processes.

This firmly places GNPs on the left side of Fig. 1, and reveals that GNPs near the SDPB limit are not equivalent to particles with multiple attractive surface patches. We attribute this difference to the softness and deformability of the repulsive corona, which allows the shape and orientation of the unshielded attractive surface “patches” to change. Currently, both simulation models of patchy particles⁹ and their experimental realizations¹⁰ fix the size, number, and shape of the patches since small changes to these parameters are known to affect their thermodynamic behavior.²⁷ Indeed, introducing flexibility and softness into the interparticle potential for certain classes of tetrahedrally coordinated colloidal particles has been shown to stabilize or destabilize a myriad of different first-order transitions.²⁸ 

More recently, a quantitative analogy between Janus particles and GNPs has been found, suggesting that the self-organization of GNPs can be understood entirely in terms of the self-assembly behavior of particles with a single attractive patch.²⁹ To the contrary, Fig. 1 posits that the addition of more arms to the surface of a GNP changes the number of attractive surface patches which remain unshielded, suggesting the analogy to particles with many surface patches (which can produce empty liquids). However, as f and M are systematically explored, we have shown that our model GNPs do not exhibit the novel phases these particles are associated with. This supports the remarkable position of Ref. 29 that, regardless of the corona, a GNP never behaves as if it has multiple surface patches. Instead, GNPs always remain effectively amphiphilic in nature.

In this work, we studied model GNPs using grand canonical Monte Carlo on a fine lattice to elucidate the thermodynamics of these systems when both the corona’s internal degrees of freedom and translational entropy are taken into account. Since GNPs have characteristics of both patchy particles and of amphiphiles, it is important to establish which model more accurately reflects the behavior of GNP systems. We found that, consistent with recent theoretical work,²⁹ GNPs exhibit behavior thermodynamically consistent with amphiphilic nanoparticles. That is, we observed thermodynamically discontinuous first-order phase separation into two isotropic fluid phases at certain conditions and thermodynamically continuous self-assembly at others. The fact that the morphologies of these self-assembled structures are fully consistent with previously established scaling theories⁴ supports the recently proposed analogy between GNPs and Janus particles.²⁹ That is, that GNPs behave as amphiphilic patchy particles with only a single patch. The patchy particle paradigm of GNPs would suggest that by tuning the number of core-shielding surface grafts on the GNP when below the SDPB limit, the number of exposed attractive surface patches could also be tuned. Our results firmly contradict such a position, revealing that GNPs do not belong to a class of multi-patch colloid models, which have the potential to produce such novel phases as equilibrium gels, but should instead be rigorously considered “nanoparticle amphiphiles.”

This publication is based on work supported by Grant No. CBET-1403049 from the U.S. National Science Foundation (NSF). Additional support has been provided by the Princeton Center for Complex Materials (PCCM), a U.S. National Science Foundation Materials Research Science and Engineering Center (Grant No. DMR-0819860).