Binary nanoparticle superlattices of soft-particle systems

Significance The phase diagram of a system of two species of particles with different diameters interacting with a soft (1/r12) potential is provided. The results provide a general framework to predict the crystalline phases observed in nanoparticle superlattices. The approach is particularly relevant for nanoparticles with hydrocarbon capping ligands where theoretical descriptions have been based on hard-sphere systems and is general enough to consider supercrystals obtained by DNA programmable self-assembly, thus providing a conceptual unification for general nanoparticle superlattice systems. The solid-phase diagram of binary systems consisting of particles of diameter σA=σ and σB=γσ (γ≤1) interacting with an inverse p = 12 power law is investigated as a paradigm of a soft potential. In addition to the diameter ratio γ that characterizes hard-sphere models, the phase diagram is a function of an additional parameter that controls the relative interaction strength between the different particle types. Phase diagrams are determined from extremes of thermodynamic functions by considering 15 candidate lattices. In general, it is shown that the phase diagram of a soft repulsive potential leads to the morphological diversity observed in experiments with binary nanoparticles, thus providing a general framework to understand their phase diagrams. Particular emphasis is given to the two most successful crystallization strategies so far: evaporation of solvent from nanoparticles with grafted hydrocarbon ligands and DNA programmable self-assembly.

The solid-phase diagram of binary systems consisting of particles of diameter σ A = σ and σ B = γσ (γ ≤ 1) interacting with an inverse p = 12 power law is investigated as a paradigm of a soft potential. In addition to the diameter ratio γ that characterizes hard-sphere models, the phase diagram is a function of an additional parameter that controls the relative interaction strength between the different particle types. Phase diagrams are determined from extremes of thermodynamic functions by considering 15 candidate lattices. In general, it is shown that the phase diagram of a soft repulsive potential leads to the morphological diversity observed in experiments with binary nanoparticles, thus providing a general framework to understand their phase diagrams. Particular emphasis is given to the two most successful crystallization strategies so far: evaporation of solvent from nanoparticles with grafted hydrocarbon ligands and DNA programmable self-assembly.
phase separation | superlattices | crystalline phases | stoichiometry A rrangement of nanoparticles (NPs) into structures with long-range order encompasses a fundamental new type of materials with potential revolutionary applications in optics, photonics, catalysis, or novel fuel energy sources, just to name a few. Over the recent years there has been a spectacular success in the assembly of nanoparticle superlattice (NPS), with the two most successful strategies consisting of evaporation of a solvent from NPs with grafted hydrocarbon chains (1)(2)(3)(4) [solvent evaporation (SE) systems] or the programmable self-assembly of DNA grafted NPs (5)(6)(7) in water (DNA systems).
Although there are different models available to investigate DNA systems (8)(9)(10)(11)(12)(13), studies of SE systems have been almost (except, for example, in ref. 14) exclusively based on the hardsphere (HS) model, following the pioneering work of Murray and Sanders on micrometer-sized colloidal systems (15). However, HS models do not provide a satisfactory description of experiments, as clear from the fact that crystals isomorphic to MgZn 2 , CaCu 5 (1), body centered cube AB6 (bcc-AB 6 ) (4) [also known as Cs 6 C 60 (7)], or quasi-crystals (16), just to name a few, have not been reported as equilibrium phases for HS (17,18). It has also been observed that different binary systems with the same NP hydrodynamic radius (with different hydrocarbon chain length, for example) do not exhibit the same equilibrium phase (19), clearly pointing to a phase diagram that depends on other parameters besides the ratio of the two NPs diameters, which completely determines the phase diagram of the HS system (15).
The interaction between two NPs is far more complex than a HS because the polymer shell (consisting of grafted hydrocarbons or DNA) is flexible. In the limit where the grafted polymers are infinitely long (f-star limit) such potential is known analytically (20) and does reveal a very soft tail. In SE evaporated systems, however, the grafted hydrocarbon chains contain between 9 and 18 hydrocarbons (3), which are too short to be described by the f-star limit. Motivated by the partial success of HS models, namely the imperfect but clear correlation between experimental equilibrium structures and those with high packing fraction (PF) and the need, for the reasons exposed (see also ref. 19), to consider a soft potential, we examine particles of different diameters interacting with an inverse power law (IPL): where σ jh = 1=2ðσ j + σ h Þ and« jh is a dimensionless number. For p = ∞, the HS model is formally recovered, but at any finite p the model is continuous (in r) and thus provides a generic example of a soft interaction; Fig. 1. In this paper, only the p = 12 case will be considered; this is a common approximation to describe a soft core, as done, for example, in Lennard-Jones systems. The phase diagram for other sufficiently short-range interactions follows similar qualitative trends, but with longer range potentials (p = 6) significant differences, which will be discussed elsewhere, are found. IPL potentials with p in the range 3-12 have been used in models of soft colloids before, such as, for example, in characterizing the phase diagram of microgel particles (see review in ref. 21), thus showing that IPL potentials provide a generic description of soft particles, although, admittedly, more refined potentials are needed for more rigorous quantitative studies (22).

Thermodynamic Parameters
The potential Eq. 1 is parameterized by with A being the particle with the largest diameter and B the smallest. These parameters are the ratio of the two diameters and the strength of the A-B and B-B interactions in units of the A-A interaction (« AA ≡ 1). An important advantage of power-law potentials is that all excess quantities do not depend on (T,V) or (T,P) but rather on a single parameter, either ξ p or its conjugate

Significance
The phase diagram of a system of two species of particles with different diameters interacting with a soft (1=r 12 ) potential is provided. The results provide a general framework to predict the crystalline phases observed in nanoparticle superlattices. The approach is particularly relevant for nanoparticles with hydrocarbon capping ligands where theoretical descriptions have been based on hard-sphere systems and is general enough to consider supercrystals obtained by DNA programmable selfassembly, thus providing a conceptual unification for general nanoparticle superlattice systems.
Author contributions: A.T. designed research, performed research, analyzed data, and wrote the paper.
The author declares no conflict of interest.
This article is a PNAS Direct Submission.
The dimensionless variables areT = k B T=« AA ,ρ = Nσ 3 A =V. The parameter N B is the number of particles of type B and N = N A + N B the total number of particles.
It should be noted, however, that, without loss of generality it can be set that« BB ≡ 1 [by defining γ′ = γð« BB Þ 1=p and« AB ′ = « AB ðð1 + γÞ=1 + γ′Þ p ]. So, in what follows, we will assume« BB ≡ 1. In summary, the phase diagrams will be presented in terms of the following four parameters: [4] All binary lattices considered (Table 1) were chosen following two criteria: high PF (see Discussion below and plot of the respective PF in SI Appendix, Fig. S1) and being experimentally relevant. All lattices are parameterized in terms of an equivalent binary HS model d 1 ≡ d AA and d 2 ≡ γ L d AA . We note that the diameters ðd 1 , d 2 Þ parameterizing the lattice do not need to be the same as ðσ A , σ B Þ in Eq. 2, but test calculations showed that allowing (γ L , d 1 ) to differ from (γ, σ A Þ did not lead to any significant new prediction. All calculations presented are for d 1 = σ A and γ L = γ. The fact that the underlying HS radii parameterizing the lattice match so well with the two radii defined by the potential Eq. 2, as well as the strong correlation between stability of the phase and HS PF, shows that the parameter γ has a clear physical interpretation as the ratio of the two soft-particle radii.

Results
The general phase diagram is a function of four coordinates, defined in Eq. 4: two that parameterize the interaction (Eq. 2) and two purely thermodynamic variables (Eq. 3). The excess free energy per particle in units of k B T for a given crystal phase in terms of the variable ξ p (described in Eq. 3) can be written as (25) The first term on the right-hand side is the free energy at zero temperature (lattice sums), the second terms are the harmonic dynamical lattice theory (DLT) corrections, and the last term encompasses the anharmonic contributions. How to obtain the explicit values of að p, γÞ and bð p, γÞ (dependence onê AB is omitted) from Eq. 11 has been extensively discussed in ref. 25. Concrete values for these coefficients are given in SI Appendix , Tables S2  and S3). The symbol x refers to the reference state (as described in SI Appendix, section SII, Eq. S14). The anharmonic function f anh ðxÞ is finite at x = 0 and can be computed from a thermodynamic integration (25,26). This calculation is beyond the scope of this paper and f anh ðxÞ ≡ 0 is assumed. As discussed in Materials and Methods as well as in ref. 25, the error introduced by this approximation is less than 0.1% right at the melting temperature, becoming virtually negligible deep in the solid phase, which is the main region of interest in this paper. The equation of state relates λ p and ξ p The chemical potential is given from Eqs. 5 and 6 as [7] In the limit of low temperatures and/or high density (ξ p → ∞), the chemical potential has the analytical form Phase diagrams for large ξ p , as determined from this equation, provide exact limiting results that guide the phase diagrams of the more general cases. The calculation of a phase diagram at a given value for the potential parameters Eq. 2 proceeds by first establishing the stable phases, determined by the condition that all eigenvalues of the dynamical matrix (DM) except for the three acoustic modes, are positive. Then, the aðp, γÞ and bðp, γÞ coefficients are calculated for all these stable phases (examples are provided in SI Appendix , Tables S2 and S3), and with it, the chemical potential Eq. 7. The stability range for each phase (for the phase diagrams reported in this study) is provided in SI Appendix, Table S1 for all cases considered; a given lattice is not always stable, but if it is, it  is always around a finite interval around its highest PF (η) and becomes unstable around η ∼ 0.5.
The next step consists of establishing whether the stable pure A k B q phase, which exists only at stoichiometric value x B = q=ðk + qÞ, does not phase separate into two coexisting A r B s and A t B u phases. The chemical potential for these coexisting phases follows from basic thermodynamics (for completeness, a derivation is provided in SI Appendix, section S1): [9] In terms of the chemical potential μ k,q ðλ p Þ of the pure A k B q phase, such phase is unstable toward phase separation if μ k,q ðλ p Þ > μðλ p , x B = q=ðk + qÞjr, s, t, uÞ.
Two illustrative examples are provided in Fig. 2. The first consists ofê AB = 0.2 at γ = 0.5275: out of the 15 possible phases, only 5 are stable: fcc/hcp, NaCl, AlB 2 , and NaZn 13 . The second consists ofê AB = 0.8 at γ = 0.558, and here six stable phases: fcc/hcp, CsCl, AlB 2 , CaCu 5 , and NaZn 13 . Both pure A and B phases are fcc. The resulting phase diagram is shown in Fig. 3: In the first case all stable phases are equilibrium at least within a limited pressure range. In the second case the CsCl is not an equilibrium phase and CaCu 5 is only stable over a very limited pressure range. In both cases, the phase diagram is different from the one obtained by direct application of the large ξ p formula, Eq. 8.

Discussion
In this way, by performing the calculations outlined in the previous section, the phase diagram for all γ-values at 12 different values ofê AB are compiled. Four representative phase diagrams (the rest are found in SI Appendix, Figs. S2-S13) are shown in Fig. 4. The interval where each phase is not just stable, but an equilibrium phase is displayed, with its corresponding PF highlighted. At a given γ, many phases can exist, either because they occur with different stoichiometry (such as NaCl and AlB 2 in Fig. 3) or because there is a phase transition at a given λ p (such as the NaZn 13 phase). If this latter case occurs, some loss of information results from the plots, because a full phase diagram like the one in Fig. 3 cannot be fully reconstructed from those plots as the critical pressure λ p (for the corresponding γ) is not provided.
Rather interestingly, phase diagrams forê AB > 1.1 are trivial (consist of A + B). As an example in Fig. 4, the caseê AB = 1 has only MgCu 2 as the only nontrivial equilibrium phase. Nontrivial phases are found forê AB < 1 (Fig. 4). In almost all cases, the equilibrium phases exist around the maximum of the PF; the opposite is not true: there are many examples where a high PF is not enough to ensure that a phase will be in equilibrium. The ZnS phase, for example, is very stable and has high PF at γ ∼ 0.22, but was never found to be an equilibrium state. Fig. 5 summarizes all phase diagrams. In this plot, the regions inê AB where each phase is an equilibrium state at the corresponding γ-values are shown.

Implications for NPS
In addition to the diameter ratio γ that fully characterizes the phase diagram of HS models, IPL potentials have one additional parameterê AB , which allows one to relate them to the experimental systems mentioned.
DNA Systems. NPS in DNA systems are driven by hybridizations between A and B particles. Such hybridizations can be described by the interaction energy between A-B particles« AB being much smaller than those for A-A (« AA = 1) and B-B (« BB = 1). In this way, association of A and B particles is energetically favored over AA or BB, thus mimicking the role of hybridizations. Therefore, DNA systems should be described by the region« AB 1 (see Fig. 4 for an example).
As is clear from Fig. 5, those systems display the CsCl, AlB 2 , bcc-AB 6 (known as Cs 6 C 60 in this context) found in experiments Fig. 2. Difference of chemical potential between a pure crystal phase and a mixture of A + B, as calculated from Eq. 9 for (γ = 0.5275,ê AB = 0.2) and (γ = 0.558,ê AB = 0.8). Fig. 3. Phase diagram as calculated from Eq. 9 for (γ = 0.5275,ê AB = 0.2) and (γ = 0.558,ê AB = 0.8). The shadowed area is in the liquid phase, as estimated from the Lindemann criterion. (7,27). There are four additional equilibrium phases in this region: NaCl, AuCu, AuCu 3 and NaZn 13 . Explicit calculations (7) (see also ref. 28) have ruled them out as they do not optimally hybridize the DNA shells, but those calculations do not consider the role of entropy and the possibility of coexistence. As the results presented here reveal, they could become equilibrium phases, at either small γ (NaCl), γ ∼ 1 (AuCu 3 ), or when there is an abundance of B particles, that is, x B ∼ 1 (NaZn 13 ). Another possibility to target these phases is to modify the DNA shell to include additional AB repulsion forces that compete against hybridizations (28)(29)(30), which is achieved by inclusion of neutral (those that do not form hybridizations) DNA strands into the NP shell (28,30). In those systems, AuCu phases have been found (28,30). SE Systems. Experimentally, the NP radius is defined as half the smallest NP-NP in a 2D hexagonal lattice of same NPs (2). The PF broadly used to interpret experimental results corresponds to the value γ obtained by the ratio of the smallest to the largest radii obtained this way. Here, I introduce the softness asymmetry (SA) parameter as s AB = effective thickness of the ligand shell for the small NP effective thickness of the ligand shell for the large NP . [10] The effective thickness is obtained by subtracting the NP core radii from the total radii. Data compiled from refs. 2, 4, 19 comprising 36 experiments are plotted in Fig. 6. The NP radius (defined from the 2D hexagonal separation) was taken from ref. 2, table 1 or by the Optimal Packing Model formula (19, 31), depending on whether the hydrocarbon capping ligands are unsaturated or saturated. The plots in Fig. 6 reveal a clear correlation both vertically (PF) and horizontally (SA). In fact, such correlation is so remarkable that rescaling the SA by a factor c = 0.6 allows a semiquantitative agreement between theory and experiments. There are some differences, namely: the Laves MgCu 2 is predicted to be the equilibrium phase over the competing MgZn 2 , which is the one experimentally reported. The calculations show a preponderance of CsCl over the AuCu competing phase, and finally, the CaB 6 are not found to be equilibrium phases (and are stable over a very limited PF).
The difference in free energy between MgZn 2 and MgCu 2 , as clear from the a, b coefficients provided in SI Appendix, Tables S2 and S3), is extremely small (but always favorable to the latter). The reasons for the preponderance of CsCl against AuCu seem less clear, although for the values of γ reported, both phases are structurally very close (they become identical for γ < 0.73). Furthermore, in the very narrow range where CaB 6 is stable, the free-energy difference with the bcc-AB 6 is very small.
A closer analysis of Fig. 6 shows that the results of ref.    6. Data compiled from the three references (including 36 different experiments and phases not analyzed in this study) with γ and s AB calculated as described in the text, superimposed to the theoretical result (Fig. 5). The coefficient is c = 0.6. In those cases where many phases exist at a point, the symbols are split horizontally for visualization. however, should a longer effective thickness for oleic acid be used, the results would be within consistency. In any case, given the simplicity of the model, the overall agreement with the theoretical phase diagram is remarkable.

Conclusions
The phase diagram for IPL potentials has been described and provides a paradigm for any general short-ranged soft potential; introduction of soft repulsive potential allows a characterization of a significant morphological diversity in self-assembled superlattices. In the region where phase separation into A + B can be overcome, the phase diagram is much richer than the one reported for HS (15,17,18). It remains an open question as to whether more general HS models (where, for example, the additivity of the radius is dropped) can account for these additional phases. More importantly, the study provides a semiquantitative description of the equilibrium crystals observed in NPS experimental systems by considering two parameters only: PF and the SA; Fig. 6. How robust predictions based on SA are is a subject for further exploration. This paper will be followed by further studies where more sophisticated potentials or free energies are considered. There have been many important developments relating the relation of potentials to given lattice structures (32)(33)(34) or characterization of optimal PFs (35-37) that will definitely synergize with this work.
It is my contention that the presented calculations capture the conceptual relevant aspects of the phase diagram with soft interactions, thus providing at least a qualitative description of mesoscale NPS. Furthermore, the methods presented are completely general and can accommodate any continuous potential. In addition, this paper highlights the clearly unappreciated value of DLT: Calculation of these phase diagrams by standard thermodynamic integrations would have required enormous computational resources.

Materials and Methods
DLT (38) provides an exact low-temperature expansion to investigate the thermodynamic properties of any phase with long-range order (23,24). Quantitative comparisons of excess free energies in soft potentials have demonstrated agreement within less than 0.1%, extending in some cases to the liquid-solid transition (25,26). Within DLT, the potential energy is expanded to quadratic order U = NU 0 + 1 2 X n X n′ u α,s ðnÞD α, sjβ,t ðn − n′Þu β,t ðn′Þ, [11] where U 0 is the potential energy at zero temperature (lattice sums), D is the DM, and u α,s ðnÞ are the displacement of the particle at the s-lattice basis at the unit cell n along the α = 1,2,3 direction. All thermodynamic quantities can be calculated from the DM and U 0 . Both the lattice sums and the DM were calculated using the HOODLT software, which is described elsewhere (25). The calculations were run in a CPU cluster using Message Passing Interface (MPI) as implemented in the Python package PyPar. A calculation of a phase diagram at all γ for given values of ð« AB Þ using 16 cores takes of the order of 6 h.