M Calandra

C Attaccalite

G Profeta

F Mauri

à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Comment on "Electronic Structure of Superconducting KC 8 and Nonsuperconducting LiC 6 Graphite Intercalation Compounds: Evidence for a Graphene-Sheet-Driven Superconducting State from electron-phonon interaction." (Li stacking is AA) [2] and cuts over the kz = 0 (b) and kz = π/c (c) planes within the local density approximation (LDA). In (b) empty blue dots are ARPES data from [1]. (d) DFT Fermi surface (red) of a Li-intercalated graphene bilayer (C6LiC6) compared with ARPES data (blue dots) from ref. [3]. (e) DFT and GW bulk LiC6 electronic-structure. The Fermi level is at zero. IL labels the metal-interlayer state. In all plots we use the LiC6 Brillouin zone where the graphene Dirac point is folded at Γ.

1 Å -1 (b) 1 Å -1 (c) 1Å -1 (d) Γ M K Γ A L H A
In [1] it is claimed that the charge transfer from the metal to the graphene layers is larger in KC 8 than in LiC 6 . This relies on a two-dimensional analysis of photoemission (ARPES) data that, a priori and without experimental verification, totally neglects the electronic band dispersion perpendicular to the graphene planes. Pan et al. could have verified this assumption by performing ARPES measurements as a function of the photon energy, as it was done for the out-of-plane dispersion in graphite [4]. Here we demonstrate that, in LiC 6 , the small interlayer distance results in a very high banddispersion along the k z direction, invalidating the analysis and the conclusions of [1].

The density-functional-theory (DFT) band-structure in Fig. 1 (e) is highly dispersive along k z (ΓA direction) by more than 1.2 eV. This results in the threedimensional Fermi surface shown in Fig. 1 (a),(b) and(c). Note, in particular, that the Fermi surface cuts in the k z = 0 and k z = π/c planes are completely different, the k z = 0 being closer to the experimental data in [1].

Correlation effects slightly increase the k z dispersion in bulk LiC 6 as shown by our GW [5] calculation. Both GW and DFT find an empty interlayer metal state in LiC 6 [6,7] (complete charge-transfer from Li to graphene) and a partially occupied one in KC 8 (incomplete chargetransfer from K) [7,8]. Thus the calculated chargetransfer is larger in LiC 6 than in KC 8 , in disagreement with Ref. [1].

To validate the accuracy of DFT on Li intercalated graphite systems, we calculate [9] the Fermi surface of a Li-intercalated graphene bilayer, obtaining results in excellent agreement with ARPES data [3] (Fig. 1 (d)). The experimental Fermi surface splitting of the bilayer measures the hopping between graphene planes separated by Li atoms and, indirectly, the k z dispersion in the bulk. Indeed, the hopping between graphene planes is a properties of the local geometry which is identical in bulk LiC 6 and in the C 6 LiC 6 bilayer.

In conclusion, the three dimensional character of the Fermi surface of bulk LiC 6 , neglected in the analysis of [1], invalidates the charge-transfer estimate of [1] and the subsequent claim of a graphene-sheet-driven superconducting state.

  Comment on "Electronic Structure of Superconducting KC8 and Nonsuperconducting LiC6 Graphite Intercalation Compounds: Evidence for a Graphene-Sheet-Driven Superconducting State" Matteo Calandra, Claudio Attaccalite, Giovanni Profeta, Francesco Mauri To cite this version: Matteo Calandra, Claudio Attaccalite, Giovanni Profeta, Francesco Mauri. Comment on "Electronic Structure of Superconducting KC8 and Nonsuperconducting LiC6 Graphite Intercalation Compounds: Evidence for a Graphene-Sheet-Driven Superconducting State". Physical Review Letters, 2012, 108, pp.149701. 10.1103/PhysRevLett.108.149701. hal-00692115

FIG. 1 .

 1 FIG. 1. (color online):(a) DFT Fermi surface of Bulk LiC6(Li stacking is AA)[2] and cuts over the kz = 0 (b) and kz = π/c (c) planes within the local density approximation (LDA). In (b) empty blue dots are ARPES data from[1]. (d) DFT Fermi surface (red) of a Li-intercalated graphene bilayer (C6LiC6) compared with ARPES data (blue dots) from ref.[3]. (e) DFT and GW bulk LiC6 electronic-structure. The Fermi level is at zero. IL labels the metal-interlayer state. In all plots we use the LiC6 Brillouin zone where the graphene Dirac point is folded at Γ.


Ran Ni

Frank Smallenburg

Laura Filion

Marjolein Dijkstra

Crystal Nucleation in binary hard-sphere mixtures: The effect of order parameter on the cluster composition

Keywords: nucleation, computer simulation, binary mixtures, colloids, hard spheres PACS numbers:

We study crystal nucleation in a binary mixture of hard spheres and investigate the composition and size of the (non)critical clusters using Monte Carlo simulations. In order to study nucleation of a crystal phase in computer simulations, a one-dimensional order parameter is usually defined to identify the solid phase from the supersaturated fluid phase. We show that the choice of order parameter can strongly influence the composition of noncritical clusters due to the projection of the Gibbs free-energy landscape in the two-dimensional composition plane onto a one-dimensional order parameter. On the other hand, the critical cluster is independent of the choice of the order parameter, due to the geometrical properties of the saddle point in the free-energy landscape, which is invariant under coordinate transformation. We investigate the effect of the order parameter on the cluster composition for nucleation of a substitutional solid solution in a simple toy model of identical hard spheres but tagged with different colors and for nucleation of an interstitial solid solution in a binary hard-sphere mixture with a diameter ratio q = 0.3. In both cases, we find that the composition of noncritical clusters depends on the order parameter choice, but are well explained by the predictions from classical nucleation theory. More importantly, we find that the barrier height and properties of the critical cluster do not depend on the order parameter choice.

The process of nucleation in colloidal systems has attracted significant attention in recent years, both in experimental and simulation studies. The framework with which phenomena like these have been described traditionally is classical nucleation theory (CNT), which is based on the notion that a thermal fluctuation spontaneously generates a small droplet of the thermodynamically stable phase into the bulk of the metastable phase. In CNT as developed by Volmer [1], Becker [2], and Zeldovich [3], the free energy of formation of small nuclei of the new phase in the parent phase is described by using the "capillary approximation", i.e., the free energy to form a cluster of the new phase relative to the homogeneous metastable phase is described by their difference in bulk free energy and a surface free-energy term that is given by that of a planar interface between the two coexisting phases at the same temperature. Thus the droplet is assumed to be separated from the metastable bulk by a sharp step-like interface in CNT. The bulk freeenergy term is proportional to the volume of the droplet and represents the driving force to form the new phase, while the surface free-energy cost to create an interface is proportional to the surface area of the cluster. Hence, small droplets with a large surface-to-volume ratio have a large probability to dissolve, while droplets that exceed a critical size and cross the free-energy barrier, can grow further to form the new stable bulk phase.

CNT has successfully explained simulation results for the nucleation of spherical particles, such as the fluidsolid and gas-liquid nucleation in Lennard-Jones systems * These authors contributed equally to this work. [4][5][6] and crystal nucleation of hard spheres [7,8]. A modified CNT has been used to explain the nucleation of anisotropic clusters of the nematic or solid phase (also called tactoids) from a supersaturated isotropic phase of colloidal hard rods [9][10][11] and the nucleation of 2D assemblies of attractive rods [12,13]. This state of affairs should be contrasted with the case of binary nucleation for which various nucleation theories have been developed that differ substantially in the way they describe the composition of the cluster [14][15][16]. For instance, Reiss assumed the surface tension to be independent of composition [14], while Doyle extended CNT by taking into account a surface tension that depends on the cluster composition [17]. However, more than 20 years later, it was shown by Renninger [18], Wilemski [19,20], and Reiss [21] that Doyle's derivation leads to thermodynamic inconsistencies. A revised thermodynamically consistent classical binary nucleation theory was developed by Wilemski in which the composition of the surface layer and the interior of the cluster could vary independently [19,20]. However, in the case of strong surface enrichment effects, this approach can lead to unphysical negative particle numbers in the critical clusters [22,23]. In addition, it was shown in Ref. [24] that the derivation by Wilemski starts off with the wrong equations, but the resulting equations are correct. Moreover, binary nucleation can be accompanied with huge fractionation effects, i.e., the compositions of the metastable phase and of the phase to be nucleated can differ enormously from the compositions of the two coexisting bulk phases. It is therefore unclear i) how to determine the surface freeenergy term for a cluster, which is in quasi-equilibrium with a metastable parent phase with a composition that is very different from those of the two coexisting bulk 2 phases, ii) whether the interfacial tension depends on composition, curvature, and surface enrichment effects, and finally iii) whether or not one can use the capillary approximation in the first place to describe binary nucleation in systems where fractionation and surface activity of the species are important. To summarize, there is no straightforward generalization to multicomponent systems of classical nucleation theory that is thermodynamically consistent, does not lead to unphysical effects, and can be applied to small nuclei [16,25].

Numerical studies may shed light on this issue, as the nucleation barrier can be determined directly in computer simulations using the umbrella sampling technique [26,27]. In this method, an order parameter is chosen and configuration averages for sequential values of the order parameter are taken. While this makes it possible to measure properties of clusters with specific values for the order parameter, it should be noted that the results can depend on the choice of order parameter. In the present paper, we investigate whether the size and composition of (non)critical clusters can be affected by the order parameter choice employed in simulation studies of multicomponent nucleation. For simplicity, we focus here on crystal nucleation in binary hard-sphere mixtures, where surface activity of the species can be neglected, and we assume the surface tension to be composition independent. The paper is organized as follows. In Sec. I, we describe the general nucleation theorem as derived by Oxtoby and Kashchiev [16], which does not rely on the "capillary approximation" and can even be employed to describe small clusters. Starting from the multicomponent nucleation theorem, it is straightforward to reproduce the usual CNT for binary nucleation, which is the focus in the remainder of the paper. In Sec. II and III, we define the (Landau) free energy as a function of an order parameter, and we describe the order parameter that is employed to study crystal nucleation. Additionally, we discuss the effect of order parameter choice on the nucleation barrier in more detail. We present results for binary nucleation for a simple toy model of hard spheres in Sec. IV, and subsequently, we study the nucleation of an interstitial solid solution in an asymmetric binary hard-sphere mixture in Sec. V.

I. CLASSICAL NUCLEATION THEORY FOR MULTI-COMPONENT SYSTEMS

We study the formation of a multicomponent spherical cluster of the new phase in a supersaturated homogeneous bulk phase α consisting of species i = 1, 2, . . . . We note that the thermodynamic variables corresponding to the metastable phase α are denoted by the subscript α, whereas those corresponding to the new phase do not carry an extra subscript to lighten the notation. We first consider a homogeneous bulk phase α characterized by an entropy S o α , volume V o α , and particle numbers N o i,α Note that the superscripts denote the original bulk phase. The internal energy U o α of the original bulk phase reads

U o α = T o S o α -P o α V o α + µ o i,α N o i,α (1)
with T o the temperature, P o α the bulk pressure, µ o i,α the bulk chemical potential of species i, and the summation runs over all species.

Following the derivation in Refs. [16,25], we now consider a spherical cluster of the new phase with a volume V separated from the original phase by an arbitrarily chosen Gibbs dividing surface. The volume of the interface is set to zero, and the particle number of species i in the cluster is given by N i + N i,s , where N i is the number of particles of species i in a volume V which is homogeneous in the new bulk phase, and N i,s is the surface excess number of particles of species i that corrects for the difference between a step-like interfacial density profile and the actual one. The surface excess number N i,s depends on the choice of dividing surface. The internal energy U of the resulting system is then given by

U = T S α + T S -P α V α -P V + Ψ + µ i,α N i,α + µ i N i + µ i,s N i,s , (2)
where P and S denote the bulk pressure and entropy of the nucleated phase, and µ i and µ i,s are the chemical potentials of species i in the new phase and the surface phase, T is the temperature of the system with the cluster, and Ψ = Ψ({N i }, {N i,s }, V ) is the total surface energy of the spherical cluster. As the volume of the surface layer is zero, the corresponding pressure is not defined. The difference in the appropriate thermodynamic potential as a function of cluster size depends on the quantities that are kept fixed during the nucleation process. If the nucleus is formed at constant temperature and constant total number of particles of each species i, and if we keep the pressure of the original phase fixed, then

T = T o , N i,α + N i + N i,s = N o i,α
, and P o α = P α . The corresponding Gibbs free energy of the initial system G o α and that of the final system G are then given by the Legendre transformation

G o α = U o α -T o S o α + P o α V o α = µ o i,α N o i,α G = U -T S + P o α (V α + V ) = (P o α -P)V + Ψ + µ i,α N i,α + µ i N i + µ i,s N i,s . (3)
If we now assume that the composition of the metastable phase α remains unchanged and we consider the Maxwell relation

∂V α ∂N i,α T,Pα,{N j =i,α } = v i,α = ∂µ i,α ∂P α T,{Ni,α} (4)
with v i,α the partial particle volumes of species i in phase α, we find that at constant pressure, the chemical potential for each species i remains constant µ o i,α = µ i,α . Subsequently, we obtain for the change in Gibbs free energy

-P)V + Ψ + (µ i (P) -µ o i,α (P o α ))N i + (µ i,s -µ o i,α (P o α ))N i,s . (5)
Consequently, the Gibbs free energy ∆G of a growing cluster depends on the number of particles N i and N i,s in the cluster and the surface energy of the cluster. Hence, one can define a free-energy surface in the multidimensional composition plane with a saddle point that corresponds to the critical nucleus [14]. The conditions for the critical cluster read

∂∆G ∂N i V,{N j =i },{N i,s } = 0, ∂∆G ∂N i,s V,{Ni},{N j =i,s } = 0, (6)
∂∆G ∂V {Ni},{Ni,s} = 0.

To recover the chemical and mechanical equilibrium conditions, we use the above conditions as well as the Gibbs-Duhem equation and the Gibbs adsorption equation. The Gibbs-Duhem equation at constant temperature for the nucleated bulk phase is

-V dP + N i dµ i = 0, (7)
and the Gibbs adsorption equation for the surface at constant temperature is

Adγ + N i,s dµ i,s = 0, (8)
where we have employed Ψ = γA. Note that γ denotes the surface free energy per unit area and A is the surface area of the cluster. The resulting equilibrium conditions for all particle species i in the critical cluster, the surface, and the metastable parent phase are then given by

µ * i (P * ) = µ * i,s = µ o i,α (P o α ), (9)
and for the pressure difference inside and outside the droplet we find

P * -P o α = ∂γ * A * ∂V * , ( 10
)
where * denotes quantities associated with a system where a critical cluster is present. Hence, the composition of the critical cluster can be determined from these saddle point conditions. In order to obtain the usual classical nucleation theory for multicomponent systems, we assume a spherical droplet with radius R. Note that the surface area is then A = 4πR 2 . In addition, we use the fact that the volume of a spherical droplet can be expressed in terms of the partial particle volumes v i of species i:

V = 4 3 πR 3 = N i v i . ( 11
)
Combining this with Eq. (10), we arrive at the generalised Laplace equation:

P * -P o α = 2γ * R * + ∂γ * ∂R * , ( 12
)
where the square brackets denote a derivative associated with the displacement of the dividing surface. One can now choose the dividing surface so that

∂γ * ∂R * = 0, (13)
and hence one recovers the usual Laplace equation. This choice for the dividing surface, corresponding to a specific value for R * and γ * , is called the surface of tension. In addition, if we use the Gibbs adsorption isotherm (8) and the Maxwell relation (4) for the bulk phase of the nucleated cluster, we find for the critical cluster

dµ * i,s = dµ * i = v i dP ( 14
)
and

A ∂γ * ∂R * = - N i,s v i ∂P * ∂R * = 0, (15)
which is the condition for a curvature independent surface tension. Since ∂P * /∂R * = 0, Eq. ( 15) implies that the dividing surface has to be chosen such that

N i,s v i = 0, (16)
which is called the equimolar surface, as for onecomponent systems N i,s = 0, i.e. the number of particles in the cluster equals the number of particles in a uniform bulk phase with the same volume. It is generally not possible in a multicomponent system to choose the dividing surface such that N i,s = 0 for all species. Thus, as v i is usually positive, N i,s < 0 for at least one of the species. This may lead to (unphysical) negative particle numbers when N i + N i,s < 0 as noted in Refs. [22,23]. However, as will be discussed in sections IV and V, there are cases in which the assumption N i,s = 0 for all i is valid.

If the nucleated phase is assumed to be incompressible, one can integrate the Gibbs-Duhem equation (7) at constant temperature to arrive at

V(P o α -P) = (µ i (P o α ) -µ i (P))N i , (17)
and using Eq. ( 5), we find

∆G = γA + (µ i (P o α ) -µ o i,α (P o α ))N i + (µ i,s -µ o i,α (P o α ))N i,s . ( 18
* i = - 2γ * v i R * , ( 19
)
where ∆µ

* i = µ * i (P o α ) -µ o i,α (P o α ).
The radius of the critical cluster R * and the barrier height ∆G * read

R * = 2γ * v i |∆µ * i | (20) ∆G * = 4πR * 2 γ * 3 = 16πγ * 3 3(∆µ * i /v i ) 2 . ( 21
)
Using Eq. ( 20) or the Maxwell relation (4), one can show:

v i ∆µ i = v j ∆µ j , ( 22
)
and the radius of the critical cluster R * can be expressed in terms of the bulk composition x i = N i / N i of the critical cluster and v = V / N i :

R * = 2γ * v x i |∆µ * i | . ( 23
)
In order to study multi-component nucleation, MC simulations are often performed in the isobaricisothermal ensemble, in which the number of particles N o 1,α and N o 2,α , the pressure of the original bulk phase P o α , and the temperature T are kept fixed. One of the assumptions of classical nucleation theory is that the composition of the metastable bulk phase remains constant, while nucleating the new phase, see Eq. ( 4). In simulations this can only be achieved if the system is sufficiently large, i.e., the volume of the metastable bulk phase is much larger than that of the nucleating cluster. Especially, for binary (multicomponent) nucleation, where the composition of the stable phase is very different from that of the metastable phase, this can lead to a huge depletion of one of the components in the metastable fluid phase, and therefore a change in composition. In order to circumvent this problem, simulation studies on binary nucleation are often carried out in the semi-grand canonical ensemble [28,29], i.e. the total number of particles

N o α = N o i,α , the chemical potential difference ∆µ o 12,α = µ o 2,α -µ o 1,
α between the two species, the pressure P o α , and the temperature T are kept fixed of the original bulk phase. The corresponding thermodynamic potential is obtained by a Legendre transformation Y (N, ∆µ 12 , P, T ) = G(N, N 2 , P, T ) -N 2 ∆µ 12 (24) Combining Eq. (3) with the conditions that the total number of particles are fixed

N o 1,α + N o 2,α = N 1 + N 2 + N 1,α + N 2,
Y o α = G o α -N o 2,α ∆µ o 12,α = µ o 1,α (N o 1,α + N o 2,α ) Y = G -(N 2,α + N 2 )∆µ 12,α = (P o α -P)V + Φ + µ 1,α (N 1,α + N 2,α ) + µ 1 (N 1 + N 2 ) -∆µ 12 N 2 -∆µ 12,α N 2 , ( 25
)
where we have set the surface excess numbers N i,s to zero. Using the Maxwell equation

∂µ 1 ∂P N,∆µ12,T = ∂V ∂N ∆µ12,P,T = v, ( 26
)
we find that due to constant pressure, the chemical potential of species 1 remains unchanged µ o 1,α = µ 1,α . Hence, we find that the change in free energy due to the formation of a nucleus ∆Y = Y -Y o α equals ∆G as given in Eq. ( 5) and the nucleation barrier can be calculated in the semi-grand canonical ensemble. Similarly, one can show that in any statistical ensemble (grand canonical, canonical, etc. ), the change in the corresponding thermodynamic potential as a function of cluster size is always the same, provided that the metastable parent phase is sufficiently large. A similar result was also obtained by Oxtoby and Bob Evans, who showed that the nucleation free-energy barriers in the isobaricisothermal and grand canonical ensemble are identical, i.e., ∆G = ∆Ω for a one-component system [30].

II. FREE-ENERGY BARRIER

While nucleation is an inherently non-equilibrium process, the assumption of local equilibrium is often made to describe the behavior of the system during the nucleation process. In essence, this assumption states that the nucleus is in quasi-equilibrium with the parent phase for every cluster size. This is approximately true if the time required to reach an equilibrium distribution of clusters is short compared to the time needed to nucleate. After the system crosses the free-energy barrier, the cluster of the new phase grows too rapidly for this assumption to be accurate, but during the nucleation process itself, local equilibrium has proven to be a useful assumption.

In order to compute the free-energy barrier that separates the metastable phase from the stable phase, an order parameter Φ (or reaction coordinate) should be defined that quantifies how much the system has transformed to the new phase. A common order parameter that is employed in nucleation studies is the size of the largest cluster in the system as defined by a certain cluster criterion. In the present paper, we restrict ourselves to binary nucleation. From Eq. ( 18), we find that the Gibbs free energy ∆G of a growing binary cluster depends on the number of particles of species 1 and 2 in the cluster, and hence, one can define a free-energy surface in the (N 1 , N 2 )-plane with a saddle point that corresponds to the critical nucleus [14]. By projecting the phase space of the system onto the (usually) one-dimensional order parameter, one can define the (Landau) Gibbs free energy ∆G(Φ) as a function of this order parameter Φ

β∆G(Φ) = G c -ln P (Φ), (27)
where β = 1/k B T , k B Boltzmann's constant, T the temperature, G c is a normalization constant generally taken to correspond to the free energy of the homogeneous metastable phase, and P (Φ) is the probability of observing an order parameter of value Φ. In a system of N particles, at fixed pressure P , and constant temperature T , the probability P (Φ) is given by:

P (Φ) = dV dr N exp[-β(U (r N ) + P V )]δ(Φ -Φ(r N )) dV dr N exp[-β(U (r N ) + P V )] (28)
with V the volume of the system, U the potential energy, and δ the Kronecker delta function. The order parameter function Φ(r N ) is a function that assigns to each configuration r N of the system a value for the order parameter.

The probability distribution P (Φ) can be obtained from Monte Carlo (MC) simulations via the umbrella sampling technique [26,27]. In this method, an additional external potential U bias is added to the system to bias the sampling towards configurations corresponding to a certain window of order parameter values centered around Φ o . By increasing Φ o sequentially, the entire free-energy barrier as a function of Φ can be sampled. The typical biasing potential used in umbrella sampling simulations is given by:

βU bias (r N ) = k(Φ(r N ) -Φ o ) 2 , ( 29
)
where the constants k and Φ o determine the width and location of the window, and r N are the positions of all N particles in the simulation.

III. ORDER PARAMETER

In order to follow a phase transformation, a cluster criterion is required that is able to identify the new phase from the supersaturated phase. In this paper, we focus on the formation of a solid cluster in a supersaturated fluid phase. In order to study crystal nucleation, the local bond-order parameter is used to differentiate between liquid-like and solid-like particles and a cluster algorithm is employed to identify the solid clusters [5]. In the calculation of the local bond order parameter a list of "neighbours" is determined for each particle. The neighbours of particle i include all particles within a radial distance r c of particle i, and the total number of neighbours is denoted N b (i). A bond orientational order parameter q l,m (i) for each particle is then defined as

q l,m (i) = 1 N b (i) N b (i) j=1 Y l,m (θ i,j , φ i,j ), ( 30
)
where Y l,m (θ, φ) are the spherical harmonics, m ∈ [-l, l] and θ i,j and φ i,j are the polar and azimuthal angles of the center-of-mass distance vector r ij = r j -r i with r i the position vector of particle i. Solid-like particles are identified as particles for which the number of connections per particle ξ(i) is at least ξ c and where

ξ(i) = N b (i) j=1 H(d l (i, j) -d c ), (31)
H is the Heaviside step function, d c is the dot-product cutoff, and

d l (i, j) = l m=-l q l,m (i)q * l,m (j) l m=-l |q l,m (i)| 2 1/2 l m=-l |q l,m (j)| 2 1/2 .
(32) A cluster contains all solid-like particles which have a solid-like neighbour in the same cluster. Thus each particle can be a member of only one cluster.

The parameters contained in this algorithm include the neighbour cutoff r c , the dot-product cutoff d c , the critical value for the number of solid-like neighbours ξ c , and the symmetry index for the bond orientational order parameter l. The hard-sphere crystals considered in this paper are expected to have random hexagonal order, thus the symmetry index is chosen to be 6 in the present study.

This choice of order parameter Φ, defined as the number of solid-like particles in the largest crystalline cluster, has been used to study crystal nucleation in various one-component systems, e.g., Lennard-Jones systems [5], hard-sphere systems [7], and Yukawa systems [31].

On the other hand, for binary systems, a variety of crystal structures can appear in the bulk phase diagram, e.g., substitutionally ordered (superlattice) crystal structures with varying stoichiometries, substitutionally disordered solid solutions, interstitial solid solutions, crystalline phases of species 1 with a dispersed fluid of species 2, etc. Nucleation of a substitutionally disordered solid solution and a crystal with the CsCl structure has been studied in a binary mixture of hard spheres using the total number of particles in the largest crystalline cluster as an order parameter, i.e. Φ = N 1 + N 2 [28]. This order parameter has also been employed in a crystal nucleation study of a substitutionally disordered face-centered cubic crystal and a crystal with the CsCl structure of oppositely charged colloids [32], and nucleation of the NaCl salt crystal from its melt using the symmetry index l = 4 instead of l = 6 for the bond orientational order parameter [33]. However, one can also define other linear combinations of N 1 and N 2 as an order parameter. When the partial particle volumes of the two species are very different, one can employ the volume of the largest crystalline cluster Φ = V = N 1 v 1 + N 2 v 2 as an order parameter. While, if the crystal structure consists of only one species, say species 1, with the other species randomly dispersed, the number of particles of species 1 in the largest crystalline cluster would be more appropriate to use as an order parameter Φ = N 1 . On the other hand, one can also use the stoichiometry n of the AB n superlattice structure to define the order parameter Φ = N 1 + N 2 /n in order to prevent a strong bias towards one of the species. More generally, if the cluster size is measured by the order parameter Φ = N 1 + λN 2 , the sensitivity of the order parameter to particles of species 2 can be tuned via the parameter λ. For λ = 1, this corresponds to the total number of particles in the cluster, while for λ = 0, this corresponds to the number of particles of type 1.

As already mentioned above, the umbrella sampling technique is often employed to determine the probability distribution P (Φ) and the Gibbs free energy ∆G(Φ). To this end, a biasing potential is introduced to sample configurations with certain values for this order parameter Φ. In this paper, we investigate the effect of the choice of order parameter on the properties of the clusters during nucleation in a binary mixture of hard spheres, where we assume that the surface excess numbers of species i are negligible. Using Eq. ( 18), we now write down explicitly the change in Gibbs free energy for binary nucleation

∆G = γA + ∆µ 1 N 1 + ∆µ 2 N 2 , ( 33
)
where

∆µ i = µ i (P o α ) -µ o i,α (P o α )
. The Gibbs free energy ∆G depends on the particle numbers N 1 and N 2 and the composition of the critical cluster can be determined from the saddle point conditions for ∆G. The free-energy surface in the two-dimensional composition plane (N 1 , N 2 ) is projected in umbrella sampling MC simulations onto a one-dimensional order parameter, e.g. Φ = N 1 + λN 2 . Hence, the projected ∆G(Φ) and the averaged (or projected) cluster composition of noncritical clusters both depend on the order parameter. We note that this is not an artifact of the umbrella sampling MC simulations, but merely the projection of a correctly measured equilibrium distribution. To determine the averaged composition of noncritical spherical clusters with radius R as a function of Φ, we can minimize ∆G with respect to N 2 while keeping the order parameter Φ fixed:

∂∆G ∂N 2 Φ = ∆µ 2 -λ∆µ 1 + 2γv 1 R (ω -λ) = 0, ( 34
)
where ω = v 2 /v 1 . If we use the umbrella sampling technique in MC simulations to determine the Gibbs free energy ∆G(Φ) as a function of Φ, one can easily determine the slope of the barrier from the simulations, which is equal to

d∆G dΦ = (∆µ 2 + 2γv 1 R ω) ∂N 2 ∂Φ + (∆µ 1 + 2γv 1 R ) ∂N 1 ∂Φ ( 35
)
with

∂N 1 ∂Φ = 1 -x -λN ∂x ∂Φ 1 -x + λx (36) ∂N 2 ∂Φ = x + N ∂x ∂Φ 1 -x + λx , ( 37
)
where we define the composition x = N 2 /N and N = N 1 + N 2 . Combining Eqs. ( 34) and ( 35) yields

ω∆µ 1 -∆µ 2 = (ω -λ) d∆G dΦ . ( 38
)
We wish to make a few remarks here. First, we recover the Gibbs-Thomson equations for the critical cluster (19) when we set d∆G/dΦ in Eq. ( 35) to zero, and we recover Eq. ( 22) from Eq. ( 38) for critical clusters. Consequently, the barrier height, size and composition of the critical cluster are independent of the choice of λ. This can also be understood from the fact that the saddlepoint in the free-energy landscape is invariant under coordinate transformations. As long as the top of the nucleation barrier corresponds to this saddle point, the average properties of the cluster will be dominated by the configurations around this saddlepoint, regardless of the chosen order parameter. While most reasonable choices of order parameter fulfill this requirement, it is possible to design order parameters that shift the top of the barrier away from the saddle point. In this case, the clusters at the top of the barrier are non-critical clusters, and rates calculated from the resulting free energy barrier are unreliable. It is important to note that a different choice of order parameter can change the height of the nucleation barrier, since the barrier height is determined by the fraction of phase space mapped to the same order parameter value at the top of the barrier. However, this effect should be small, as the probability of finding a cluster at the top of the nucleation barrier is dominated by the probability of being in the saddle point of the free-energy landscape. For noncritical clusters, we clearly find that the slope of the barrier, and hence the composition of the cluster, depends on the choice of order parameter via λ. Below, we study the effect of the choice of order parameter for a simple toy model of hard spheres and for the nucleation of an interstitial solid solution in an asymmetric binary hard-sphere mixture. It is interesting to compare this to past studies investigating one-component systems with higher-dimensional order parameters [34,35]. For the Lennard-Jones system, Moroni et al., have shown that the number of particles in the cluster alone is insufficient to provide a good prediction for the probability a cluster will grow out to a large crystal [34]. Using a two-dimensional order parameter, they observed a strong correlation between the crystallinity and the size of clusters with a 50% probability of growing out. Specifically, clusters with a large amount of face-centered-cubic (fcc) ordering require much smaller sizes to grow out than those with more body-centeredcubic (bcc) ordering. They found that this correlation was not visible in the two-dimensional free-energy landscape, and argued that the shape and structure of a nucleus could determine whether it will grow out. However, we note that the two-dimensional order parameter is still a projection from a higher-dimensional phase space. Thus, the properties of non-critical clusters likely depend on the choice of order parameter as well.

IV. A SUBSTITUTIONAL SOLID SOLUTION

In order to obtain more insight in the effect of order parameter choice on the cluster composition of noncritical clusters, we first investigate binary crystal nucleation in a toy model of hard spheres. Here, we consider a system consisting of two species of hard spheres with identical sizes, but tagged with different colors, say species 1 is red and species 2 is blue. Obviously, the stable solid phase to be nucleated is a substitutional disordered face-centeredcubic (fcc) crystal phase with the red and blue particles randomly distributed on an fcc lattice. Refs. [7,8] showed that the nucleation barriers for pure hard spheres are well-described by the predictions from classical nucleation theory, where because of the condition of the equimolar surface, the surface excess number N s = 0. It is therefore safe to neglect the surface excess numbers for the present model as well. In addition, it is clear that the partial particle volumes v i and volume per particle v are identical, and ω = v 2 /v 1 = 1. Using the Gibbs-Thomson equations for a binary critical cluster (19), we find that the supersaturation ∆µ * 1 = ∆µ * 2 = -2γ * v/R * , and hence the composition of the critical cluster follows straightforwardly from the bulk chemical potentials µ * 1 (P o α ) and µ * 2 (P o α ), which depends on the bulk chemical potentials of the original bulk phase and the supersaturation. As already mentioned above, the composition of noncritical clusters depends on the choice of order parameter, i.e., the projection of the two-dimensional composition plane onto a one-dimensional order parameter Φ. Using Eq. ( 38), we find that for λ = 1, the composition of noncritical cluster is determined by the supersaturation ∆µ 1 = ∆µ 2 and the bulk chemical potentials of the original bulk phase. For λ = 0, we only measure the number of particles of one color, say red, in the cluster. However, a thermodynamic average of all clusters with N 1 red particles also includes all post-critical clusters with many blue particles, and as a result, the order parameter fails to work for λ = 0. For non-zero values of λ, the ensemble of clusters of each size is well-defined, and we can perform umbrella sampling MC simulations to measure the average cluster composition.

In order to keep the composition of the metastable fluid fixed, we perform Monte Carlo simulations on a binary mixture with N = 1000 hard spheres in the semi-grand canonical ensemble. Both species of hard spheres are identical in size with diameter σ, and are either tagged red (species 1) or blue (species 2). The simulations were carried out in a cubic box with periodic boundary conditions and the Metropolis sampling consists of particle displacements and volume changes, and attempts to switch the identity (color) of the particles. The acceptance rule for the identity swap moves is determined by the chemical potential difference ∆µ o 12,α [28,29]. We use the umbrella sampling technique to determine the nucleation barrier ∆Y = ∆G as a function of an order parameter Φ = N 1 + λN 2 , where N 1 (N 2 ) denotes the number of red (blue) solid-like particles in the largest crystalline cluster in the system as determined by the local bond-order parameter and cluster criterion described in Sec. III with cutoff radius r c = 1.3σ, dot-product cutoff d c = 0.7, and number of solid bonds ξ c ≥ 6. We first calculate the nucleation barrier for λ = 1, for which the order parameter Φ is simply the total number of solidlike particles in the largest cluster. We set the reduced pressure P * = P o α σ 3 /k B T = 17, and ∆µ o 12,α = 0, which corresponds on average to an equimolar mixture of red and blue hard spheres for the metastable fluid phase. We plot the resulting nucleation barriers ∆G as a function of Φ in Fig. 1. We note that the nucleation barrier for λ = 1 is equivalent to the nucleation barrier for a pure system of hard spheres [7,8]. In addition, we show the composition of the largest cluster as a function of Φ in Fig. 2. We find that the averaged composition x = N 2 /N = 0.5 as it should be since ∆µ 1 = ∆µ 2 and the bulk chemical potentials of the metastable fluid are equal µ o 1,α = µ o 2,α . Using the binomial coefficients and the measured one-dimensional free-energy barrier, we determine the two-dimensional free-energy landscape ∆G(N 1 , N 2 )/k B T = -ln P (N 1 , N 2 ) from the probability distribution function

P (N 1 , N 2 ) = exp[-∆G(N 1 + N 2 )/k B T ] 2 N N N 1 .
(39) Fig. 3 presents a contour plot of the two-dimensional free-energy landscape β∆G(N 1 , N 2 ) as a function of N 1 and N 2 . Exemplarily, we also plot isolines for the order parameter Φ = N 1 + λN 2 for λ = 1 and 0.5 to show the projection of the two-dimensional composition plane onto a one-dimensional order parameter.

In order to check the effect of order parameter choice in the biasing potential (29) on the nucleation barrier and the composition of the clusters, we also calculate the nucleation barrier for λ = 0.5 at the same reduced pressure. We plot the nucleation barrier in Fig. 1 and the averaged composition of the cluster as a function of Φ in Fig. 2. We clearly find that the barrier height is not affected by the choice of order parameter in the biasing potential in agreement with our predictions in Sec. III, but the critical cluster "size" as measured by Φ, i.e. 79 and 96 for λ = 0.5 and 1, respectively, depends on the order parameter choice as expected. While the barrier height is not significantly affected by the choice of order parameter in the biasing potential, in agreement with our predictions in Sec. III, the critical cluster "size" as measured by Φ, i.e.

79 and 96 for λ = 0.5 and 1, respectively, depends on the order parameter choice as expected. In addition, we determine the theoretical prediction for the cluster composition using Eq. (38). Using the measured slope of the nucleation barrier from Fig. 1, we obtain the chemical potential difference ∆µ 12 (Φ) of species 1 and 2 in the cluster from Eq. (38). Using Eq. ( 39), we find

P (N 1 , N 2 ) ∝ 2 N N ! N 2 !(N -N 2 )! exp[-βN 2 ∆µ 12 (Φ)]
(40) from which we determine the most probable (or averaged) composition x = 1 -exp[-β∆µ 12 (Φ)]. The theoretical prediction for the composition is plotted in Fig. 2. We find good agreement with the measured composition, except for very small cluster sizes, where we do not expect CNT to match our nucleation barriers. For comparison, we also plot the same predictions for the nucleation paths in Fig. 3. We clearly observe that the two nucleation paths cross at the saddle point yielding the same size and composition of the critical cluster for both order parameters, as expected.

Finally, we also determine the composition of the clusters from the steady-state distribution. In systems where the nucleation of the new phase is measured directly, either in experiments or simulations, the measured cluster size distribution corresponds to a steady-state distribution rather than an equilibrium distribution. The steadystate distribution observed during the nucleation process is different from the equilibrium distribution, as clusters that exceed the critical cluster size during the steadystate process will continue to grow further. The steadystate distribution depends both on the free-energy landscape and the dynamics of the system, and includes a flux across the free-energy barrier, whereas the equilibrium distribution can only be determined by preventing the system from nucleating, i.e, constraining the maximum cluster size by e.g. umbrella sampling MC simulations. While the equilibrium and steady-state distributions are in good agreement for small cluster sizes, they disagree strongly for postcritical cluster sizes, i.e., when the system crosses the free-energy barrier. In particular, the equilibrium cluster size distribution shows a minimum corresponding with the maximum in the free-energy barrier, and the steady-state distribution generally decreases continuously (even) beyond the critical cluster size.

We calculate the cluster composition from the steadystate distribution for our binary mixture of hard spheres. To this end, we determine the free energy as a function of cluster size N 1 and N 2 from Eq. (39) using a fit to the free-energy barrier obtained from umbrella sampling MC simulations with λ = 1. The dynamics of the cluster are described by the following rates:

k +,1 N1,N2 = 1 k +,2 N 1 ,N 2 = 1 k -,1 N 1 ,N 2 = exp[-β(G(N 1 -1, N 2 ) -G(N 1 , N 2 ))] k -,2 N1,N2 = exp[-β(G(N 1 , N 2 -1) -G(N 1 , N 2 ))].
Here, k

+(-),i N 1 ,N 2
is the rate associated with adding (removing) a particle of species i to (from) the nucleus consisting of N 1 and N 2 particles. Hence, clusters can only grow or shrink by one particle at a time with a rate determined by the corresponding free-energy difference. In order to determine the steady-state cluster size distribution, we set a limit to the steady-state distribution by defining a maximum cluster size, which exceeds the critical cluster size. As a barrier crossing can be considered as a one-way event, subsequent nucleation events should start again from the metastable fluid phase. To this end, we impose that the addition of an extra particle to a nucleus with this maximum cluster size falls back to size zero. We note that this step is not reversible, and results in slightly modified rates for nuclei with the maximum cluster size and for clusters of zero size. With the exception of these steps, the dynamics obey detailed balance.

In order to determine the steady-state distribution, we set the rate at which clusters of size (N 1 , N 2 ) are created to zero. Hence, the flux with which clusters of size (N 1 , N 2 ) are created should balance the flux with which clusters of this size disappear:

P ss (N 1 , N 2 ) i (k +,i N1,N2 + k -,i N1,N2 ) = P ss (N 1 + 1, N 2 )k -,1 N1+1,N2 + P ss (N 1 -1, N 2 )k +,1 N1-1,N2 + P ss (N 1 , N 2 + 1)k -,2 N 1 ,N 2 +1 + P ss (N 1 , N 2 -1)k +,2 N 1 ,N 2 -1 .
Here, P ss (N 1 , N 2 ) denotes the steady-state cluster size distribution. The equations for cluster size zero and the maximum cluster size are slightly different due to a flux of clusters from maximum to zero cluster size. By solving this set of linear equations numerically, we obtain the steady-state distribution. Subsequently, the average cluster composition can be obtained from the steadystate distribution by averaging over clusters with equal Φ = N 1 + λN 2 . The resulting cluster composition is shown in Fig. 2 for λ = 0.5. Since the two-dimensional steady-state cluster size distribution, which is symmetric in N 1 and N 2 decreases monotonically with cluster size, the resulting projected composition is always lower than 0.5 and matches well with the cluster compositions obtained from umbrella sampling MC simulations and the theoretical prediction, except at small cluster sizes as expected. Moreover, in the limit of large (postcritical) clusters, the cluster growth rate approaches a constant for the current choice of dynamics, resulting in a nearly flat steady-state cluster size distribution and a cluster composition of 0.5.

In conclusion, we have shown using a simple model for a binary mixture of hard spheres that the composition of the critical cluster does not depend on the choice of order parameter, while the composition of noncritical clusters is affected by the order parameter. This is a direct consequence of the projection of the two-dimensional free-energy landscape onto a one-dimensional order parameter, say Φ = N 1 + λN 2 , which influence directly the projected (Landau) ∆G(Φ) and the averaged (or projected) cluster composition. Moreover, as the umbrella We also plot a few isolines for the order parameter Φ = N 1 + λN 2 for λ = 0.5 and 1 (dashed lines), and we plot the nucleation path (solid lines labeled with λ = 1 and λ = 0.5) for the two order parameters that we considered as predicted by (38). The two nucleation paths cross at the saddle point corresponding to the critical cluster size.

sampling method allows us to equilibrate the system for various values of the order parameter, the system can be regarded to be in local equilibrium for each value of the order parameter. The nucleation paths that the system then follows remain close to the minimum free-energy path (see Fig. 3), and thus the height of the nucleation barrier is largely unaffected by the choice of order parameter.

V. AN INTERSTITIAL SOLID SOLUTION

We consider crystal nucleation of an interstitial solid solution in a highly asymmetric binary mixture of large and small hard spheres with size ratio q = σ 2 /σ 1 = 0.3, where σ 1(2) denotes the diameter of species 1 (large spheres) and 2 (small spheres). The interstitial solid solution consists of a face-centered-cubic crystal phase of large spheres with a random occupancy of the octahedral holes by small spheres, and hence the composition of the interstitial solid solution can vary from x = N 2 /N ∈ [0, 1] [36]. As the volume of this solid phase is not largely affected by the density of small spheres, we set the partial particle volume v 2 and ω = v 2 /v 1 to zero. Using Eq. ( 38), we find the following relation if the system is in local equilibrium at fixed order parameter Φ

= N 1 + λN 2 ∆µ 2 = λ d∆G dΦ . ( 41
)
For λ = 0, the order parameter Φ = N 1 measures only the large spheres in the cluster, and the cluster composition of both critical and noncritical clusters is deter- mined by the chemical equilibrium condition for the small spheres in the cluster and the metastable fluid phase, i.e., ∆µ 2 = 0. For λ = 1, when all particles in the clusters are counted by the order parameter Φ = N 1 + N 2 , the composition of precritical clusters will have a higher density of small particles compared to the chemical equilibrium condition for the small particles in the cluster and the metastable fluid phase, as the slope of the nucleation barrier is positive, and similarly postcritical clusters will have a lower density of small particles. For both order parameters, we find that the critical cluster satisfies the Gibbs-Thomson equation (19), and thus for a partial particle volume v 2 = 0, we obtain chemical equilibrium for the small particles in the critical cluster and the fluid phase independent of the order parameter choice.

As the composition and size of the critical cluster and the corresponding barrier height are not affected by the choice of order parameter, we set λ = 0 in order to investigate whether or not we observe diffusive equilibrium for species 2 for all noncritical clusters. To keep the composition of the fluid fixed, it would be convenient to use again Monte Carlo simulations in the semi-grand canonical (N P T -∆µ 12,α ) ensemble. However, the acceptance probability of changing small spheres into large spheres is extremely small, which makes the equilibration time of the simulation prohibitively long, even when we use the augmented semigrand ensemble presented in Ref. [28], where the diameter of the particles is changed gradually in different stages. In order to solve this problem, we determine the free-energy barrier using the umbrella sampling technique in isothermal-isobaric MC simulations, in which the pressure P o α , the temperature T , and the particle numbers N o 1,α and N o 2,α are kept fixed of the original metastable bulk phase. We perform successive simulations for each window, but in such a way that the composition

x o α = N 2,α /(N o 1,α + N o 2,α
) of the metastable fluid phase is on average kept fixed during the growth of the nucleus. To this end, we first measure the instantaneous composition x α of the fluid phase in the initial configuration for the successive umbrella sampling windows centered around a new order parameter value Φ. If the composition of the fluid has changed more than 0.1%, we resize random particles in the fluid phase during an equilibration run until the fluid phase reaches its original composition x o α . We then start the production run to measure the probability distribution function P (Φ) and the corresponding part of the free-energy barrier in a normal isobaric-isothermal MC simulation. We assume that the composition of the fluid phase during MC simulations of a single umbrella sampling window does not change significantly, since the cluster size is approximately constant. In order to determine the composition of the fluid phase, we first determine the largest crystalline cluster in the system by using the local bond-order parameter and cluster criterion as described in Sec. III with cutoff radius r c = 1.1σ 1 , dot-product cutoff d c = 0.7, and number of solid bonds ξ c ≥ 6. The composition of the fluid is defined as

x α = (N o 2,α -N 2 )/(N o 2,α +N o 1,α -N 2 -N 1 )
where N 1 is the number of large spheres in the cluster and N 2 is the number of small spheres which have at least 6 neighbors of large spheres in the cluster within cut-off distance r c = 1.1σ 1 . N o 1,α and N o 2,α denote the total number of large and small spheres in the MC simulation.

In addition, we determine the composition of the solid nucleus x = N 2 /N . In order to avoid surface effects and defects in the crystal structure of the solid nucleus, we determine the fraction of octahedral holes that is occupied by a small sphere in the fcc lattice of the large spheres in the solid cluster. An octahedral hole is defined as a set of 6 large particles, where each particle is a neighbour of 4 other particles in the same set, and the octahedral hole is occupied by a small particle if all 6 large particles are within a cutoff radius of 0.22σ 1 of the center-of-mass of this small sphere.

We first determine the nucleation barrier in a normal N o 1,α N o 2,α P o α T MC simulation using the umbrella sampling technique for system sizes N o α = N o 1,α + N o 2,α = 3000, 6000, and 9000 particles. The initial fluid composition is set to x o α = 0.5 and reduced pressure P * = βP o α σ 3 1 = 25. We plot the Gibbs free energy ∆G/k B T as a function of the number of large spheres N 1 in the largest crystalline cluster in Fig. 4. We observe that the nucleation barrier height and critical cluster size decreases upon increasing system size. This can be explained by a change in the composition of the metastable fluid phase during the growth of a crystalline cluster. In Fig. 5, we plot the composition of the metastable fluid phase as a function of the cluster size N 1 for the various system sizes. We clearly find that the fluid composition changes significantly during the growth of a solid nucleus for smaller system sizes. In order to corroborate this result, we perform umbrella sampling MC simulations in which the composition of the metastable fluid phase is kept fixed in each successive umbrella sampling window using the method as described above. The composition of the fluid phase is indeed kept fixed by this method as shown in Fig. 5. The nucleation barrier as obtained by fixing the composition of the metastable fluid phase is presented in Fig. 4. As the nucleation barrier calculated at fixed fluid composition should correspond to an infinitely large system size, we plot the barrier heights

∆G * /k B T as a function of 1/N o α with N o α = N o 1,α + N o 2,α
. We find that the barrier height depends linearly on 1/N o α within errorbars. Moreover, extrapolating the barrier heights obtained from N o 1,α N o 2,α P o α T MC simulations to the thermodynamic limit, we find that the finite-size corrected barrier height agrees well within errorbars with the barrier height determined from umbrella sampling MC simulations with fixed fluid composition corresponding to an infinitely large system size. In addition, we plot the composition of the solid cluster as a function of cluster size N 1 in Fig. 5, and we find no strong dependence of the cluster composition on system size.

Finally, we determine the composition of (non)critical clusters for the nucleation of the interstitial solid solution for four different fluid compositions x o α = 0.2, 0.5, 0.7 and 0.8 at statepoints well-inside the fluid-solid coexistence region using umbrella sampling MC simulations with fixed fluid composition and system size N o α = 3000. Following Ref. [28], the "supercooling" was kept fixed, i.e., P o α /P * coex = 1.2, where P * coex is the pressure at the bulk fluid-solid coexistence at the corresponding fluid composition. We note however that these statepoints correspond to different values for the supersaturation, and can therefore lead to significantly different barrier heights. We determine the Gibbs free-energy barrier and the cluster composition as a function of cluster size N 1 using umbrella sampling MC simulations, and plot the results in Fig. 6 and7 for the four different fluid compositions. In Fig. 7, the dashed lines indicate the compositions predicted by Eq. ( 41) with λ = 0, i.e., chemical equilibrium for species 2 in the clusters and the metastable fluid phase. For comparison, we also plot the composition of the coexisting solid phase at P o α . We clearly observe that the measured cluster compositions obtained from umbrella sampling MC simulations are in good agreement with the predictions from CNT for cluster sizes larger than 30, which predicts chemical equilibrium for the small spheres in the cluster and the metastable parent phase. If we now take a closer look at the statepoint defined by x o α = 0.2 and P o α /P * coex = 1.2 for the metastable fluid phase, we find from Ref. [36] that the composition of the coexisting fluid and solid phase after full phase separation should be x 0.47 and 0.15, respectively. Interestingly, we find that the composition of the nucleating clusters is much lower (x 0.07) than that of the coexisting bulk crystal phase. Hence, the phase separation is mainly driven by the nucleation of large spheres while maintaining chemical equilibrium for the smaller species throughout the whole system. Only when the chemical potential of the large spheres in the metastable fluid is sufficiently low due to a depletion of large spheres as a result of crystal nucleation and crystal growth, small spheres will diffuse into the crystal phase in order to increase the composition of the solid phase. However, we note that the chemical equilibrium condition for the smaller species only holds for the present order parameter choice Φ = N 1 , whereas any other choice of order parameter would certainly yield different results for the cluster composition.

For highly asymmetric binary hard-sphere mixtures, where the stable solid phase corresponds to a fcc of large spheres with a dispersed fluid of small particles, one would naively expect that the small particles are always in chemical equilibrium during the nucleation process. Hence, in order to study crystal nucleation in highly asymmetric mixtures, one can employ an effective pairwise depletion potential description as described by Bob Evans and coworkers in Ref. [37][38][39] provided that threeand higher-body interactions are negligible and the depletion potentials are determined at fixed chemical potential of the small spheres. Such an effective pair potential approach was employed in a nucleation study in the vicinity of a critical isostructural solid-solid transition in a binary mixture of hard spheres with size ratio q = σ 2 /σ 1 = 0.1, but this study showed according to the authors a breakdown of classical nucleation theory [40]. It would be interesting to investigate whether or not the breakdown is caused by the (false) assumption of chemical equilibrium of small spheres during the nucleation process. For less asymmetric binary hard-sphere mixtures, where the small spheres cannot diffuse freely in the solid cluster, chemical equilibrium of the smaller species is harder to maintain, especially when the nucleated crystal phase has long-range crystalline order for both species as in the case of a superlattice structure where the chemical potentials of the two species are not independent as it is determined x o α = 0.8, P * = 47

x o α = 0.7, P * = 37.5

x o α = 0.5, P * = 25

x o α = 0. by the stoichiometry of the crystal structure. It would be interesting to investigate at which size ratio and pressures this crossover occurs.

VI. CONCLUSIONS

In this paper, we have studied crystal nucleation in a binary mixture of hard spheres and we have investigated what the effect is of the choice of order parameter on the composition and size of both critical and noncritical clusters. We have studied nucleation of a substitutional solid solution in a simple toy model of identical hard spheres but tagged with different colors and we investigate the nucleation of an interstitial solid solution in a binary hard-sphere mixture with a diameter ratio q = 0.3. In order to study nucleation of a crystal phase in computer simulations, a one-dimensional order parameter is usually defined to identify the solid phase from the supersaturated fluid phase. We have shown that the choice of order parameter can strongly influence the composition of noncritical clusters, as the free-energy landscape in the twodimensional composition plane (N 1 , N 2 ) is projected onto a one-dimensional order parameter, say Φ = N 1 +λN 2 , in umbrella sampling MC simulations. This is supported by the good agreement that we found between our results on the composition of noncritical clusters obtained from umbrella sampling MC simulations and the predictions from CNT for the nucleation of a substitutional solid solution in a toy model. While the effect is clearly visible in the case of a binary system, it should occur more generally whenever a higher-dimensional free-energy landscape is projected onto a single order parameter. For the nucleation of an interstitial solid solution in a highly asymmetric hard-sphere system, we found that the composition of noncritical clusters is determined by the chemical equilibrium condition of the small spheres in the cluster and the fluid phase, as the partial particle volume of the small spheres in the solid phase can be neglected. We compared the composition of the noncritical clusters obtained from umbrella sampling MC simulations and the theoretical prediction from CNT, and found again good agreement. More importantly, we find that the barrier height and the composition of the critical cluster are not significantly affected by the choice of order parameter. As a result, critical clusters and the barrier height should be comparable even with different order parameters. 13 phase diagram of pure hard spheres is well-studied and understood by now. Amazingly, arriving at a consensus on this very simple system took more than 40 years, starting in 1957 with the computer simulations of Alder and Wainwright [41] and Wood and Jacobson [42], which were heavily disputed by many hard-nosed theorists at that time, and continued with the free-energy calculations by Hoover and Ree [43] in 1968, and the simulation work by Woodcock [44], Bolhuis, Frenkel, Mau, and Huse [45], and Bruce, Wilding, and Ackland [46] in 1997, and the more mathematical approach by Koch, Radin, and Sadun [47,48] in 2005 on the free-energy difference of the hexagonal-close-packed and face-centered-cubic crystal phase. Yet, a formal mathematical proof for the bestpacked structure of hard spheres is still missing, although the work by Hales in 1998Hales in (published in 2006 [49] [49]) appears to be almost a complete proof. Apparently, computer simulations of hard spheres are extremely tricky, which may have contributed to Bob's antipathy to computer simulations. The situation for binary mixtures of hard spheres is not much better. A major flaw in the historical overview of the "big hard-sphere paper" is the sentence "Thus the phase behavior of binary hard-sphere mixtures is well understood for q > 0.4." where q is the diameter ratio of the small and large spheres. Re-cent computer simulation studies have demonstrated the thermodynamic stability of the so-called Laves phases in the size ratio range of 0.76-0.84 [50,51], and the binary superlattice structure with the NaCl structure in the size ratio range of 0.25-0.42 has often been misidentified by a so-called interstitial solid solution [36,52] of which we studied the nucleation in the present paper. On the other hand, the glass transition and random close packing (see the extensive review by Parisi and Zamponi [53] and the many references therein), nucleation of hard spheres [7,8], and the (best-packed) crystal structures of binary hard-sphere mixtures [54][55][56][57] are still far from conclusive and will certainly keep us off the street for some time. To the best of MD's knowledge, the phase diagrams for q = 0.033, 0.05, 0.1 and 0.2 presented in the "big hard-sphere paper" have hitherto been able to stand the test of time, but only time will provide the final verdict. MD thanks Bob Evans for the warm hospitality, support, constructive criticism, inspiring discussions, and the cosmetic changes to her manuscripts during her stay in Bristol, and wishes Bob and Margaret all the best in the coming years.

We acknowledge financial support of a NWO-VICI grant and the high potential programme of Utrecht University. We would like to thank the referee for his/her careful reading of the manuscript and for the suggestions. Below we address in detail the comments of the referee. All changes to the manuscript have been marked in red.

We hope that the new version of the manuscript can be accepted for publication in Molecular Physics. We would like to thank the referee for his/her careful reading of the manuscript and the valuable comments, which were addressed in the revised manuscript.

Response to the Referee

• The referee inquired if the invariance of the barrier height and critical cluster to the choice of the one-dimensional order parameter is a general result. For the properties of the critical cluster, this is indeed a result of the geometrical properties of the saddle point in the free energy landscape, which is invariant under coordinate transformations. However, it should be noted that the height of the nucleation barrier is not exactly independent of the choice of order parameter. We have removed all reference to a constant barrier height from the manuscript, and have added the following to address this question: "This can also be understood from the fact that the saddlepoint in the free-energy landscape is invariant under coordinate transformations. As long as the top of the nucleation barrier corresponds to this saddle point, the average properties of the cluster will be dominated by the configurations around this saddlepoint, regardless of the chosen order parameter. While most reasonable choices of order parameter fulfill this requirement, it is possible to design order parameters that shift the top of the barrier away from the saddle point. In this case, the clusters at the top of the barrier are non-critical clusters, and rates calculated from the resulting free energy barrier are unreliable. It is important to note that a different choice of order parameter can change the height of the nucleation barrier, since the barrier height is determined by the fraction of phase space mapped to the same order parameter value at the top of the barrier. However, this effect should be small, as the probability of finding a cluster at the top of the nucleation barrier is dominated by the probability of being in the saddle point of the free-energy landscape."

• We thank the referee for pointing out the work of Moroni et al. We have now included a reference to their work in the manuscript, and added the following discussion: "It is interesting to compare this to past studies investigating one-component systems with higher-dimensional order parameters [Moroni et have shown that the number of particles in the cluster alone is insufficient to provide a good prediction for the probability a cluster will grow out to a large crystal. Using a two-dimensional order parameter, they observed a strong correlation between the crystallinity and the size of clusters with a 50% probability of growing out. Specifically, clusters with a large amount of face-centered-cubic (fcc) ordering require much smaller sizes to grow out than those with more body-centered-cubic (bcc) ordering. They found that this correlation was not visible in the two-dimensional freeenergy landscape, and argued that the shape and structure of a nucleus could determine whether it will grow out. However, we note that the twodimensional order parameter is still a projection from a higher-dimensional phase space. Thus, the properties of non-critical clusters likely depend on the choice of order parameter as well."

• The typo on page 5 has been corrected.

  -G o α when a nucleus is formed in the bulk of the original phase: ∆G = (P o α

1 FIG. 1 :

 11 FIG.1: Gibbs free energy ∆G(Φ)/k B T as a function of order parameter Φ = N 1 + λN 2 for a binary mixture of red (species 1) and blue (species 2) hard spheres with equal diameter σ as obtained from umbrella sampling MC simulations at a reduced pressure of P * = P o α σ 3 /kBT = 17 with λ = 1 and λ = 0.5.

1 FIG. 2 :

 12 FIG.2: Composition x = N2/N of the largest crystalline cluster as a function of order parameter Φ = N1+λN2 for a binary mixture of red (species 1) and blue (species 2) hard spheres with equal diameter σ as obtained from umbrella sampling simulations at pressure P * = P o α σ 3 /k B T = 17 with λ = 1 (red circles) and λ = 0.5 (green squares). For comparison, we plot the theoretical prediction (38) using the measured nucleation barrier of Fig.1(black solid line) and the composition determined from a steady-state cluster size distribution for λ = 0.5 (blue dashed line). The critical cluster size is Φ 79 and 96 for λ = 0.5 and 1, respectively.

FIG. 3 :

 3 FIG.3: Contour plot of the two-dimensional free-energy landscape ∆G(N 1 , N 2 )/k B T as a function of N 1 and N 2 . We also plot a few isolines for the order parameter Φ = N 1 + λN 2 for λ = 0.5 and 1 (dashed lines), and we plot the nucleation path (solid lines labeled with λ = 1 and λ = 0.5) for the two order parameters that we considered as predicted by(38). The two nucleation paths cross at the saddle point corresponding to the critical cluster size.

FIG. 4 :FIG. 5 :

 45 FIG. 4: a) Gibbs free energy ∆G/k B T as a function of number of large particles N 1 in the largest crystalline cluster using normal N o 1 N o 2 P T MC simulations and isobaric-isothermal MC simulations with fixed fluid composition xα = 0.5, and pressure P * = βP o α σ 3 1 = 25. b) Free-energy barrier height ∆G * /k B T as a function of 1/N o α .

2 FIG. 7 :

 27 FIG.7: Cluster compositions x = N2/N as a function of cluster size N 1 for four different fluid compositions x o α = 0.2, 0.5, 0.7, and 0.8 for a binary mixture of hard spheres with size ratio 0.3 at supercooling P o α /P * coex = 1.2. The long dashed lines denote the composition predicted by CNT, which corresponds to chemical equilbrium of species 2 in the solid clusters and the metastable fluid phase, while the dotted lines denote the composition of the coexisting bulk crystal phase.

  TMPH-2010-0457 Title: Crystal Nucleation in binary hard-sphere mixtures: The effect of order parameter on the cluster composition Authors: Ran Ni, Frank Smallenburg, Laura Filion, and Marjolein Dijkstra Dear Prof. Hansen, Thank you for your letter December 3rd, 2010 and the forwarded referee report.

3 Page 17 of 17 URL

 317 : http://mc.manuscriptcentral.com/tandf/tmph Molecular Physics

  al, Phys. Rev. Lett. 94 (2005), P. R. ten Wolde et al, Science 277 (1997)]. For the Lennard-Jones system, Moroni et al.,

Acknowledgments

It is a great pleasure to contribute to this Special Issue of Molecular Physics dedicated to Bob Evans in recognition of his outstanding contributions to the field. MD had the privilege to be a post-doctoral fellow in the group of Bob Evans in 1997-1999, where she and her husband René van Roij started to work on the phase behaviour of highly asymmetric binary hard-sphere mixtures under the highly inspiring and most enthusiastic supervision of Bob Evans. This work was carried out in the lively atmosphere of the H.H. Wills Physics laboratory at Bristol University and was published in a series of papers [37][38][39]. We, often jokingly, referred to this work as the "truth", especially after long hours of heated discussions with a couple of beers in the pub. But has this work stood the test of time so far? Considering the dazzling number of papers that are still being published on hard spheres, one might get the impression that the situation is still far from clear. A historical overview on the phase behavior of pure and binary mixtures of hard spheres was already presented in Ref. [39], also referred to as the "big hard-sphere paper" by us. To date the