{'fraction_non_alphanumeric': 0.030776699029126213, 'fraction_numerical': 0.007766990291262136, 'mean_word_length': 4.549078763064325, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 8, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Polydispersity induced solid-solid transitions in model colloids

Specialized Monte Carlo simulation techniques and moment free energy method calculations, capable of treating fractionation exactly, are deployed to study the crystalline phase behaviour of an assembly of spherical particles described by a top-hat"parent"distribution of particle sizes. An increase in either the overall density or the degree of polydispersity is shown to generate a succession of phase transitions in which the system demixes into an ever greater number of face-centred cubic"daughter"phases. Each of these phases is strongly fractionated: it contains a much narrower distribution of particle sizes than is present in the system overall. Certain of the demixing transitions are found to be nearly continuous, accompanied by fluctuations in local particle size correlated over many lattice spacings. We explore possible factors controlling the stability of the phases and the character of the demixing transitions.


I. INTRODUCTION
Hard spherical particles can be packed to fill maximally just over 74% of space, in the face centred cubic (fcc) structure. 1 For systems in thermal equilibrium such as a colloidal suspension, this structure remains preferred 2,3 for packing fractions down to about 55% 4 where melting occurs. But what is the thermodynamically optimal structure for spherical colloids that are "polydisperse", i.e. have a spread of diameters? Polydispersity should act to destabilize a colloidal crystal because of the difficulty of accommodating a range of particle sizes within a single lattice structure; but despite sustained attention spanning over three decades (see e.g. [5][6][7][8][9][10][11][12][13] ), there is a lack of consensus as to what stable structures arise instead.
Attempts to address this matter have focused on the use of analytical theory and simulation to predict the fate of a single crystal in the dense regime (above typical fluid densities) when the degree of polydispersity becomes large. Broadly speaking, two incompatible proposals have emerged: either the system demixes into multiple coexisting crystalline phases 7,8 or, alternatively, crystalline phases disappear altogether, 10 the crystal being replaced by an equilibrium glassy phase. 11 Ideally, of course, one should like to settle the matter as to which (if either) of these scenarios is correct by simply performing an experiment with a suitable suspension of colloids. But the inhibition of diffusion in crystalline phases is expected to render solid-solid demixing transitions largely unobservable on experimental timescales, even if they are thermodynamically favoured 14 . This would -on the face of it -appear to render our central question moot. However, one should recognize that even when equilibrium is itself unattainable in practical situations, independent knowledge of the stable state represents an important baseline for interpreting dynamical properties of colloidal systems (such as crystallization kinetics 15 ) which can be understood in terms of the topology of the free energy surface. 16 There are also suggestions 17 that the equilibrium phase diagram sheds light e.g. on the ability of a glassy phase to crystallize. The question as to the nature of the true stable state is thus of more than merely academic interest.
In our view, the disparity in the predictions of previous theoretical and computational work is traceable to the fact that when considering phase separation, little or no account was taken of "fractionation", i.e. the phenomenon whereby the distribution of the particle diameters, σ, can vary from one coexisting phase to another. 18,19 Indeed it is now well established that fractionation can radically alter the qualitative features of phase behaviour in polydisperse systems compared to their monodisperse counterparts (see 20 for a review). Accordingly, it is essential to fully incorporate its effects if one hopes to describe the equilibrium phase behaviour of polydisperse systems correctly. 21 To quantify fractionation 22 one simply counts, for a certain phase (labeled α), the number density of particles having diameters in the range σ . . . σ + dσ. This serves to define a density distribution ρ (α) (σ). However, in real colloidal suspensions, one has the constraint that the overall distribution of sizes (across all phases) has a form fixed by the synthesis of the suspension. This gives rise to a generalized lever rule: with ρ (0) (σ) the "parent" density distribution, ρ (α) (σ) the "daughter" distributions, and λ (α) the fractional volume occupied by phase α (so that α λ (α) = 1). Since the form of the parent is fixed, only its scale is free to vary, e.g. by dilution with solvent, and one writes ρ (0) (σ) = n (0) f (σ), where n (0) is the total number density and f (σ) is a prescribed normalized shape function. The degree of polydispersity, δ, is then defined as the standard deviation of the parent distribution f (σ), expressed in units of its mean. Fractionation greatly complicates the task of determining the phase behaviour of polydisperse systems compared to their monodisperse counterparts. To illustrate this, consider (for a given colloidal system) increasing n (0) from an initially low value, i.e. following a "dilution line" through the phase diagram of the system. For sufficiently large n (0) the system typically encounters a coexistence region of the phase diagram, which is entered at a "cloud" point 20 value of n (0) . At and beyond this density the system separates into differently fractionated daughter phases. However, as a consequence of fractionation both the daughter distributions ρ (α) (σ) themselves and their associated fractional volumes λ (α) depend nonlinearly on n (0) . Thus in order to quantify the phase behaviour one is faced with the challenge of determining the daughter phase properties for all values of n (0) within the coexistence region -a situation which contrasts with the monodisperse case where the coexistence densities are independent of the total density, whilst the fractional volumes depend linearly on it.
One theoretical technique that does take fractionation into account exactly (within the context of a mean field framework) is the moment free energy (MFE) method. Previous work using this method by one of us 23 predicts that, for polydisperse spheres, increasing δ or n (0) within the solid region leads to a succession of phase transitions in which the system demixes into an ever greater number of differently-fractionated daughter phases. Each daughter phase contains a narrower distribution of particle diameters than the parent. This MFE calculation thus provides clear evidence for the scenario of multiple coexisting solids. But it uses approximate free energy expressions, which for solids are derived from those of binary mixtures and implicitly already assume that all solids are fcc. Independent confirmation of its predictions is then highly desirable, but has hitherto been lacking. The purpose of the present report is therefore to provide a definite answer to the question of the nature of the equilibrium phase behaviour via state-of-the-art Monte Carlo (MC) simulations, and to compare with MFE calculations; both fully provide for fractionation and employ a fixed parent size distribution.
Our paper is organized as follows. In Sec. II we introduce our model systems: size disperse hard spheres (which we have studied by the MFE method), and soft spheres (which we have studied by MC simulation). Sec. III provides a brief description of both the MFE method and the bespoke MC techniques required for dealing with fixed polydispersity and fractionation. Thereafter, in Sec. IV, we report our observations concerning the phase behaviour of these models, the central finding being that the original MFE calculations are indeed correct: as δ and/or n (0) are increased, a succession of transitions occurs in which the system demixes into first two, then three, then four fractionated coexisting fcc crystalline phases. We analyse the observations to arrive at a qualitative picture of when a crystalline phase will become unstable to phase separation. In Sec. V, we investigate in detail the character of these phase transitions, finding that some are strongly first order, whilst others are quasi-continuous. To quantify the differences, we introduce and measure a susceptibility that probes particle size fluctuations. The associated correlation length at the near continuous transitions is found to extend over several crystal lattice spacings. Prompted by these observations, we use MFE calculations to study in detail how the shape of the size distribution affects the tendency of a given solid phase to exhibit a near continuous demixing transition. This leads us to a simple approximate criterion that quantifies this tendency. Finally, in Sec. VI we summarize and discuss the significance of our results, and indicate some issues for future work.

II. MODELS
The systems that we shall consider in this work are assemblies of spheres interacting either by a repulsive soft sphere potential (as considered by simulation) or a hard sphere potential (as studied in our MFE calculations). The soft sphere interaction potential between two particles i and j with position vectors r i and r j and diameters σ i and σ j is given by with particle separation r ij = |r i −r j | and interaction radius σ ij = (σ i +σ j )/2. The choice of this potential rather than hard spheres is made on pragmatic grounds; in our isobaric SGCE simulations (to be reported below), any MC contraction of the simulation box that leads to an infinitesimal overlap of two hard spheres would always be rejected, so (particularly at high densities) we can expect higher MC acceptance rates using a "softer" potential. In common with hard spheres, the monodisperse version of our model freezes into an fcc crystalline structure, [24][25][26] and temperature only plays the role of a scale: the thermodynamic state depends not on n (0) and T separately but only on the combination n (0) ( /k B T ) 1/4 . Phase diagrams for different T then scale exactly onto one another, and we can fix /k B T = 1. In all cases we consider parent size distributions of the top-hat form: Here the width parameter c controls the degree of polydispersity δ = c/ √ 3. In the following we use the mean particle diameterσ as our unit of length.

III. METHODOLOGIES
A. Analytical calculations: the moment free energy method Calculating analytically the phase behaviour of polydisperse systems is a challenging problem. 27 This is because for each of the infinitely many different particle sizes σ one has a separate conserved density ρ(σ). Effectively one thus has to study the thermodynamics of an infinite mixture, where e.g. from the Gibbs rule there is no upper limit on the number of phases that can occur.
The moment free energy (MFE) method 27-30 is designed to get around this issue by effectively projecting the infinite mixture problem down to that for a finite mixture of "quasi-species". This is possible when the free energy density has a so-called truncatable form, where the excess part f ex depends on a finite number of moments of the density distribution, This truncatable structure obtains for a large number of models of mean field type. Importantly for our purposes, it is also found in accurate free energy expressions for polydisperse hard spheres, with the simple weight functions w i (σ) = σ i (i = 0, 1, 2, 3). Specifically, we use the free energy developed by Bartlett 31 on the basis of the simulation data of Kranendonk et al. 32 for binary mixtures. As mentioned above, this effectively presupposes fcc structures for all solids, so that validation e.g. by simulations, as provided in this paper, is important.
The MFE method provides a way of expressing the ideal contribution to the free energy from Eq. (4), which depends on the complete shape of the density distribution, in terms of the moment densities ρ i . The result is the moment free energy. The key feature of the method is that if one then treats the quasi-species densities ρ i as if they were densities of ordinary particle species, and calculates phase equilibria accordingly, the results for cloud points are fully exact. Within coexistence regions, the method can be extended by including additional moments. 33 Their weight functions can be chosen adaptively, and using the resulting approximation as an initialization, 34 the exact phase equilibrium conditions can then be solved numerically, even if e.g. three or four daughter phases are present. Overall, the MFE approach is therefore the method of choice for our current investigation. We do not give further details of the numerical implementation here as these are set out in full in Ref. The appropriate ensemble for determining phase behaviour in dense assemblies of polydisperse particles is the isobaric variant of the semi-grand canonical ensemble (SGCE). 21 Within this ensemble, the particle number N , pressure p, temperature T , and a distribution of chemical potential differencesμ(σ) are all prescribed, while the system volume V , the energy, and the form of the instantaneous density distribution ρ(σ) all fluctuate. 35 The fluctuations in ρ(σ) are linked to the volume fluctuations by the relation V ρ(σ)dσ = N . Importantly, they permit the sampling of many realizations of the polydisperse disorder, thus ameliorating finite-size effects. Moreover, in conjunction with volume fluctuations, they facilitate separation into differently fractionated phases. Coexistence of two or more phases is signalled by a multimodal form for the distribution of some order parameter such as the overall number density n = N/V or the volume fraction η, which for a phase with density distribution ρ(σ) can be written as η = dσρ(σ)(π/6)σ 3 .
Operationally, the sole difference between the isobaric semi-grand canonical ensemble and the constant-NpT ensemble 36 is that one implements MC updates that select a particle at random and attempt to change its diameter σ to σ by a random amount σ −σ drawn from a zeromean uniform distribution. This proposal is accepted or rejected with a Metropolis probability controlled by the change in the internal energy and chemical potential: 35 where ∆Φ is the internal energy change associated with the resizing operation and β = 1/(k B T ).
For SGCE simulations of a polydisperse system at some given N and T , it is necessary to first determine the pressure p and distribution of chemical potential differencesμ(σ) such that a suitably defined ensembleaveraged density distribution matches the prescribed parent ρ (0) (σ) = n (0) f (σ). Unfortunately, this task is complicated by the fact p andμ(σ) are unknown functionals of the parent. 37 To solve this problem -and hence determine correct coexistence properties -we shall employ a version of a scheme originally proposed in the context of grand canonical ensemble studies of polydisperse phase coexistence 38 and later extended to the SGCE, 21,39 the latter implementation of which we now summarize.
The strategy is as follows. For a given choice of n (0) and temperature T , one tunes p,μ(σ) and the λ (α) iteratively within a histogram reweighting (HR) framework, 40 such as to simultaneously satisfy both a generalized lever rule and equality of the probabilities of occurrence of the phases, i.e.
with E as defined in Eq. (8) below. In the first of these constraints, Eq. (6a), the ensemble averaged daughter density distributions ρ (α) (σ) are assigned by averaging only over configurations belonging to the respective phase, distinguishable via the multimodal character of the order parameter distribution p(n). The deviation of the weighted sum of the daughter distributions ρ(σ) ≡ α λ (α) ρ (α) (σ) from the target n (0) f (σ) is conveniently quantified by a "cost" value: In the second constraint, Eq. (6b), provides a measure of the extent to which the probability of each phase occuring, p (α) , is equal for each of the m coexisting phases. Imposing this equality ensures that finite-size errors in coexistence parameters are exponentially small in the system volume. 38,41 The iterative determination of p,μ(σ) and λ (α) such as to satisfy Eqs. (6a) and (6b) proceeds thus: 1. Guess initial values of the fractional volumes λ (α) corresponding to the chosen value of n (0) . Usually if one starts near a cloud point, the fractional volume of the incipient phase will be close to zero.
2. Tune the pressure p (within the HR scheme) such as to minimize ∆.

4.
Measure the corresponding value of E.
In step 3 the minimization of ∆ with respect to variations inμ(σ) is most readily achieved 42 using the following simple iterative scheme forμ(σ): for iteration k → k + 1. This update is applied simultaneously to all entries in the histogram ofμ(σ), and thereafter the distribution is shifted so thatμ(σ 0 ) = 0, where σ 0 is the chosen reference size. The quantity 0 < a < 1 appearing in Eq. (9) is a damping factor, the value of which may be tuned to optimize the rate of convergence. Note that (as described in 38 ) it is important that one minimizes ∆ and E to a very high precision in order to ensure that the finite-size effects are exponentially small in the system size. Typically we iterated until both were less than 10 −12 .
The values of λ (α) and p resulting from the application of the above procedure are the desired fractional volumes and pressure corresponding to the nominated value of n (0) . As mentioned above, daughter phase densities and volume fractions are obtainable by monitoring the multimodal nature of the order parameter distribution p(n), which allows configurational properties to be assigned to a given phase. 39

IV. PHASE DIAGRAM AND SOLID STABILITY
We consider first the overall phase diagram of the soft sphere system as studied in our simulations for a system of N = 256 particles. Fig. 1a shows (empty symbols) the boundaries of the fluid-solid (FS) coexistence region at low densities. These boundaries are the cloud curves coming from the low and high density regimes, respectively, and were previously determined by us using MC phase-switch techniques. 21 Our focus in this paper is on the solid region at higher densities. Here a comprehensive exploration of the (n (0) -δ) plane is impractical because of the relatively high computational cost of our specialized simulation technique. But we can understand important qualitative features by following the dashed trajectory included in Fig. 1a 43 . Along this path, we monitored the state of the system via the probability distribution of the fluctuating total number density p(n), which serves as an order parameter for phase changes. Starting from the fcc solid cloud point at δ = 6.3%, we initially increased n (0) in a stepwise fashion (filled circles) to n (0) = 1.45, and then switched to increasing δ at constant n (0) as a potentially faster route to demixing. Indeed, at δ ≈ 8% there was a smooth change in p(n) from single to double peaked; an example of the double peaked form is shown in Fig. 2a. The two associated phases were identified as being fcc solids. As is physically reasonable, the higher density solid (HDS) daughter phase contains a surplus of the smaller particles while the lower density solid (LDS) phase has more of the larger particles; see Fig. 3 below.
Continuing to higher δ eventually led to spontaneous melting of the system at δ = 13.7%, implying that the limit of metastability with respect to a fluid-solid-solid (FSS) coexistence had been overstepped, as is indeed predicted by our MFE calculations (see Fig. 1b). We therefore backtracked slightly into the solid-solid (SS) region, embarking on a new trajectory with increasing n (0) at constant δ = 13.5%. This produced a third peak in p(n) at n (0) ≈ 1.475. The corresponding intermediate density solid (IDS) was again found to be isostructural with the other two, with dominant particle sizes between those in the HDS and LDS. Finally, increasing the overall density to n (0) ≈ 1.68 we observed that the central IDS peak in p(n) split rather smoothly into two peaks, yielding a four peaked structure (Fig. 2b). All four solids were again identified as having an fcc structure.
We next compare to our theoretical MFE calculations. These used the same parent size distribution (3) but, as explained above, the analysis was performed for hard spheres. The reason is that no suitable polydisperse model free energies are available for the soft repulsive potential (2). Nevertheless, the qualitative physics should be the same. Indeed, taking a comparable path (see Appendix A) through the calculated phase diagram, we find the same features as in the simulations, as shown in Fig. 1b. Quantitatively, the fluid-solid coexistence region is narrower, and transitions to multiple solids occur at lower n (0) and δ, presumably because with a hard repulsion, a crystal can accommodate above average-sized particles less easily.
A key feature of the phase diagram is the absence of glassy phases. The frustration that could otherwise en-  gender such phases is avoided precisely by fractionation.
To illustrate this, we show in Fig. 3 the density distributions for two coexisting solids, at the state points marked by the circles in Fig. 1. The figure also shows the parent density distribution. It is likely that if a single solid were forced to have this size distribution at the density considered, it would indeed assume a disordered, glassy structure. Our results show that at equilibrium, this is avoided by effectively splitting the range of particle sizes among two phases, allowing each phase to remain crystalline on account of its now narrower range of particle size variation. This scenario is then broadly in line with that proposed by Bartlett,7 but the split in sizes is not "sharp" in the sense that particles of a given size would be found exclusively in one phase or the other. Such a sharp split would require infinite differences between phases of the relevant size-dependent chemical potentials. Apart from the general phenomenon of fractionation, Fig. 3 also demonstrates good agreement between the simulation results and the MFE predictions, with e.g. the cross-  ing point between the three density distributions in both cases located somewhat to the right of the parental mean.
As the density or parent polydispersity of a system in the SS regime are increased further, the polydispersity in the two daughter phases becomes unfavorably large. At this point a third solid appears that takes up the middle of the size distribution, producing three daughters whose size distributions are again sufficiently narrow. This is illustrated in Fig. 4. At the transition from this SSS regime to four solids (SSSS), we then see a process that is qualitatively similar to the S-SS transition: the middle (IDS) phase splits into two phases, each again with a narrower size distribution (Fig. 5). It is worth emphasizing that also in these more complicated fractionation scenarios, the agreement between simulations for soft spheres and theory for hard spheres remains good.
A natural question to ask about the results so far is: what determines the stability of solid phases, i.e. when do new solids appear? Intuitively one would expect that there should be a certain threshold in polydispersity beyond which a given single solid phase would become thermodynamically unfavourable. This threshold should then depend on how dense the phase is: a denser solid can accommodate less variation in particle sizes. To test this idea quantitatively, we plot along the path through our phase diagram the polydispersity δ versus the volume fraction η of all coexisting phases. The results are shown in Fig. 6. One sees that the coexisting phases do indeed cluster around a line in the (η, δ) plane, although the clustering is clearly tighter for the MFE (hard sphere) theory. In the plot for the latter case we also show the S-SS phase boundary from Fig. 1b as a dashed line. Recall that this is the boundary as it applies to solids with a top hat size distribution. Most of the coexisting phases that we find lie inside this phase boundary, implying that with their smoother size distributions they can tolerate a somewhat larger amount of polydispersity. In summary, while the general picture of a line in the volume fractionpolydispersity plane where solids become unstable holds true, this line is broadened into a transition region by its dependence on the shape of the size distribution. Motivated by this finding we also experimented with other measures of polydispersity to see whether they would reduce this dependence on distribution shape. In particular, we considered δ 2n = [ 1 2 (σ − σ ) 2n ] 1/(2n) / σ where the averages are over particle sizes σ and σ randomly drawn from the relevant size distribution. For n = 1 this gives the conventional δ; for n → ∞ it becomes the difference between the largest and the smallest particle size present, normalized by the mean size. While one may imagine the latter quantity to be the most relevant one for determining crystal stability, we found in practice that the clustering in the (η, δ 2n ) plane becomes worse for larger n, with the most easily interpretable results being the ones shown above for n = 1. As a final comment on Fig. 6 it is worth highlighting that the phase with the highest density (HDS, shown by squares) in fact always has the smallest volume fraction among the daughter phases for a given parent. This again reflects the strong fractionation effects: as illustrated in Figs. 3-5, the HDS phase contains the smallest particles, and this reduces its volume fraction η to the point where it is smaller than for all other phases. This trend is true throughout, i.e. the ordering of the daughter phases by density n is always the reverse of the ordering by volume fraction η.

A. Order parameter distributions and fractional volumes
In this section we discuss the nature of the transitions as our system of polydisperse spheres fractionates into an increasing number of solids. Our focus will be on the rather surprising finding that these transitions can be nearly continuous in character.
Initial evidence for this claim is provided by Fig. 2 above. This shows the distributions p(n) of the fluctuating number density in the MC simulations, with each peak corresponding to one of the solid phases. One sees in Fig. 2a, for the S-SS transition, that the initial sin- gle peak splits smoothly into two nearby peaks which then rapidly move outwards towards more clearly separated densities. This contrasts with what one would have expected for a first order transition, where a new peak appears at some finite distance from the initial peak and gradually acquires more and more weight. Such a scenario is found, along our particular path through the phase diagram, for the SS-SSS transition (data not shown). The SSS-SSSS transition, on the other hand, is again nearly continuous, like the S-SS transition. This can be seen in Fig. 2b, where the middle peak splits smoothly into two new peaks which move apart and form the IDS1 and IDS2 phases. Further evidence for nearly continuous transitions to multiple solids is provided by the variation of the fractional phase volumes λ (α) , shown in Fig. 7. One observes that at the S-SS transition, the fractional volume occupied by the new phase has a strongly nonlinear variation with the parent polydispersity. In fact, looking at the simulation results (Fig. 7a), where we cannot get reliable data close to the transition, one would guess that the fractional volume of the new phase has a discontinuous onset, as is typical of phase transitions which are continuous in the thermodynamic sense. The difficulty in obtaining data close to the transition in simulations stems from the fact that in a finite-sized system the crit- ical density distribution p(n) has two peaks, so that one has to proceed some way into the two phase region before one can be sure that peaks observed in p(n) indicate genuine phase coexistence. Looking at the right half of Fig. 7, the behaviour at the SS-SSS transition is rather different, with the fractional volume taken up by the new phase increasing smoothly from zero in an almost linear fashion. This is in line with expectations for a first order transition. The SSS-SSSS transition, on the other hand, again shows nearly continuous behaviour. As for the S-SS transition, phase coexistence cannot be determined unambiguously from the simulation data for our finite systems, and the data outside of the resulting gap are again suggestive of a jump in the new fractional volume at the transition. The MFE calculations show that there is no real jump, rather a strongly nonlinear increase from zero, so that the transition is close to but not fully critical. Taken together, the above observations of the behaviour of the density distribution p(n) in the simulations, and of the variation of the fractional phase volumes, provide strong evidence that demixing transitions to multiple solids can be near critical. Along our specific path through the phase diagram, it is the S-SS and SSS-SSSS transitions that are of this type. To investigate this issue in more detail, we now turn to characterizing the near critical properties at the level of single phases, via appropriate correlation functions of particle size fluctuations.

B. Correlations in size fluctuations
To define a measure of how strongly spatially correlated size fluctuations are in our solids, we consider first a grand canonical setting for a single phase in a fixed volume V , and with imposed chemical potentials µ(σ).
The fluctuating density distribution ρ(σ) has ensemble average ρ(σ) . If we define moment densities ρ n = dσρ(σ)σ n , then the normalized ensemble average size distribution is ρ(σ) / ρ 0 . Its variance Σ = ρ 2 / ρ 0 − ( ρ 1 / ρ 0 ) 2 sets the scale for any particle size fluctuations. To define our correlation measure χ, we measure the mean particle size in any configuration, which is ρ 1 /ρ 0 , and construct its variance across the ensemble. This is then normalized by Σ and multiplied by system volume V to get a quantity with the dimension of a volume: In the thermodynamic limit of large V , ρ 0 and ρ 1 have small fluctuations so one can expand ∆( Abbreviating the ensemble-averaged mean size asσ = ρ 1 / ρ 0 , this gives 11) or in terms of the fluctuating density distribution The denominator here could also be written as ρ 2 ρ 0 − ρ 1 2 . To motivate further the above definition of our measure of correlations χ, one can express it via correlation functions of the full spatially-resolved density ρ(r, σ). The fluctuations of the latter can be expressed in terms of the pair correlation function g σσ (r) between particles of sizes σ and σ as 45 So the numerator of (12) is, using ∆ρ(σ) = V −1 dr∆ρ(r, σ), This shows that our definition of χ is physically reasonable: it is the volume integral of a correlation function that measures the spatial correlations of fluctuations in particle size away from the ensemble mean. We will therefore also refer to χ as the size fluctuation susceptibility. Note that the trivial first term above makes a contribution of 1/ ρ 0 to χ. This is the unit volume per particle and of order unity in the density range we are considering. We will see below that it is negligible compared to the main contribution from the correlation function integral. Some care is needed when relating the susceptibility χ as defined above to a length scale ξ for the spatial correlations of size fluctuations. Away from criticality, and in d spatial dimensions, then since the correlation function being integrated decays on a spatial scale of ξ, one estimates χ ∼ ξ d . This is the identification we made previously. 44 At criticality, on the other hand, the correlation function appearing above will have a spatial power law decay with |r| −d+2−η up to the cutoff, and hence the susceptibility scales as χ ∼ ξ 2−η where η is the standard critical exponent (and not, as elsewhere in the paper, the volume fraction).
One can show that, for large systems, the size fluctuations we are considering are the same in all reasonable ensembles, for example a semi-grand canonical ensemble where particle number N is fixed and the volume V can fluctuate. In this case the factor V in (10) is replaced by N/ ρ 0 to give 13) and this is the method we use to extract χ from simulation data. For the theoretical calculations, we employ (11) and extract the (co-)variances of the fluctuations of the moment densities ρ 0 , ρ 1 from the appropriate curvature matrix of the moment free energy. 20 Note finally that in the context of experiments on colloids a canonical ensemble, with fixed particle number N , volume V and parent size distribution, would be the most natural description. For a single phase, the mean size is then fixed and no size fluctuations occur. But χ can still be defined in terms of the pair correlation function g σσ (r) as described above, provided the spatial integration over r is cut off at some distance much larger than the correlation length but much smaller than the system size. This eliminates the contribution from the nonzero values g σσ (r) − 1 = O(1/N ) that remain at larger r when the total particle number N is fixed. 45 Once several phases appear, each phase has fluctuating particle numbers and volume, but one can check that the size fluctuations in each phase, and hence the size fluctuation susceptibility, are as would be calculated for single phases in the grand canonical ensemble.
Having defined how we will quantify the strength of correlations in spatial particle size fluctuations, we show in Fig. 8 results for χ along the vertical and final horizontal paths through the phase diagrams of Fig. 1. One observes that χ grows large near the transitions to two and four solids, confirming their near continuous character. In the latter case, the splitting of the middle peak seen earlier in p(n) suggests that the new solids arise out of the IDS phase, and this is consistent with large fluctuations occurring (see Fig. 8) only in this phase and not the HDS or LDS. The MFE predictions are, again, in good qualitative accord with the simulation data.
To summarize our observations in this section, the behaviour of fractional phase volumes and of the order parameter distributions p(n) suggested that phase transitions to multiple solid phases can be near critical in nature; for our path through the phase diagram this applies to the S-SS and SSS-SSSS transitions. We proposed the size fluctuation susceptibility χ as a quantitative measure of the range of correlations in the spatial fluctuations of particle sizes. Results for this from both simulations and MFE calculations then demonstrated that these transitions are indeed close to critical, being characterized by values of χ far above the unit volume per particle. That such critical or near critical transitions from one to several solids might occur is plausible given that S-SS critical points are observed also in simulations of binary hard sphere mixtures. 32 The free energy expression for polydisperse hard spheres that we use in our MFE calculations was devised by Bartlett 31 on the basis of free energies fitted to these binary mixture simulations. The polydisperse system must then "inherit" the existence of critical points, though not in any trivial way. For example, the polydisperse system has many more degrees of freedom for fluctuations in its size distribution, and one can show from this that spinodal densities are always lower in the polydisperse than in the corresponding binary case.
It is worth stressing that even in transitions involving multiple solid phases (SSS-SSSS), criticality is essentially a single-phase property. Indeed, we have found above quite distinct values of χ in the three coexisting solid phases before the transition to four solids. In the simulations, we have only a single phase in the simulation box for most of the time, which emphasizes further that χ is determined from the properties of this single phase. To be more precise, there is an effect of the presence of other phases: the lever rule forces the density distributions of all phases to add up to the parent, and this provides constraints on the chemical potentials µ(σ). Once we know these chemical potentials, however, we can determine χ individually for every phase, independently of the others.
We next ask what features of a given size distribution make it undergo a critical or near critical transition to multiple solids. Having recognized that criticality is a single-phase property, we focus in this enterprise on the cloud point for the transition S-SS from a single solid to two fractionated solids.

C. Predicting criticality
The question of determining whether the S-SS transition from a parent with a given size distribution is close to critical can be cast in quantitative terms as follows: how large is the size fluctuation susceptibility χ at the S-SS cloud point? We investigate this using MFE calculations for the hard sphere case; precise simulation studies would inevitably require finite-size scaling to larger system sizes than we can access using our computational resources.
As explained in Sec. III A, the free energy expression that we use for polydisperse hard spheres have excess contributions that depend only on the moment densities ρ i with i = 0, 1, 2, 3, defined by the weight functions w i (σ) = σ i . One can then show in generality (see e.g. 28 ) that the criterion for a spinodal, where a phase becomes unstable to local density fluctuations, involves these moments as well as those defined by the secondorder weight functions w i (σ)w j (σ) = σ i+j , giving in our case moments up to ρ 6 . For a given particle size distribution, all the ratios ρ 1 /ρ 0 , . . . , ρ 6 /ρ 0 are fixed and the density ρ 0 ≡ n at the spinodal can be found from the spinodal criterion. The additional condition for a spinodal point also to be a critical point involves in addition the third-order weight functions w i (σ)w j (σ)w k (σ) = σ i+j+k , which produce moments up to ρ 9 . Inserting the spinodal density, the exact critical point condition resulting from our model free energies is then some function of ρ 1 /ρ 0 , . . . , ρ 9 /ρ 0 . These are the 1st to 9th moments of the normalized size distribution, and so whether a parent phase with a given size distribution will exhibit a critical S-SS transition or not depends only on these moments.
Unfortunately, because the solid free energies we use are derived from fits to simulation data, 31 the critical point condition that results is far too complicated to allow for any analytical progress. We therefore proceed initially by solving the condition numerically for a range of parent size distributions of interest. The first case to consider is evidently the top hat distribution studied throughout the paper so far. Here there is only a single parameter to vary, namely the polydispersity δ. We solve for each δ the spinodal condition to find the spinodal density, and then evaluate the critical point condition at this density. It turns out that there is indeed a critical point in the phase diagram, at (n (0) = 1.1669, δ = 0.0472). It is marked in Fig. 1b, and lies close to the path through the phase diagram that we have considered above. This rationalizes why the S-SS transition along this path is near critical, with a large value of χ: at the critical point itself, we would have found χ diverging to infinity at the transition.
The situation with regard to the shape of the parent distribution is not trivial, however. For example, in previous work we considered both triangular and Schulz dis-tributions, 46 and found no critical points on the S-SS cloud curve in the physically relevant ranges of density and polydispersity. To get more insight, we consider next families of parent distributions where we can tune both the width, as measured by δ, and the shape. Generalizing from the top hat case studied above, we look first at "slanted top hat" parents where the size distribution is f (σ) = A + Bσ in some interval σ − ≤ σ ≤ σ + , and zero otherwise. We adjust A, B, σ − and σ + so that f (σ) is normalized, has mean 1 as before, and the desired value of δ. This leaves one degree of freedom, which we express via the slant ratio R = f (σ + )/f (σ − ), with R = 1 giving back the simple top hat distribution. Proceeding as for the top hat parent, we can now determine numerically for fixed slant ratio R the critical value (if any) of δ, or vice versa. In the resulting Fig. 9 we observe that whether or not there are critical points for a given parent shape depends on R: for R below around 0.81, no critical points appear; for slightly larger values, two critical points can exist in the phase diagram, and for values of R around unity and above we generically find one critical point.
So far we have asked what marks out parent size distributions that have critical S-SS transitions, which corresponds to χ = ∞ at the cloud point. Here we digress slightly to ask how χ at the cloud point then varies as we move away from the critical parent shape. Data from MFE calculations are shown for this in Fig. 10, where we consider parents with fixed density n (0) = 1.133 as on the vertical path in Fig. 1b. For given slant ratio R we find the polydispersity δ at the cloud point; see the bottom inset of Fig. 10. The main plot displays the resulting cloud point value of the susceptibility χ against R. It is seen to diverge as a critical value R = R c is approached, and indeed by solving the critical point criterion for the given parent density we find a single such critical value, The top inset of Fig. 10 plots the susceptibility χ versus the distance from the critical parent shape. The data are consistent with a divergence as χ ∼ |R − R c | −2 , except for the points nearest R c either side, where our numerics become unreliable. That exponent value may seem surprising at first: for our mean field free energy, the susceptibility for models in the Ising universality diverges as χ ∼ |T − T c | −γ with γ = 1. But R smoothly changes the parent shape, and the latter is analogous to the Ising magnetization m. We should then identify |R − R c | with m and this leads to the scaling χ ∼ m −γ/β ∼ |R − R c | −γ/β . For our mean field free energy this gives an exponent value γ/β = 1/(1/2) = 2, exactly as observed. A more accurate theory which captures the non-mean field Ising singularities would then be expected to give in d = 3 the exponent γ/β = δ −1 ≈ 3.9.
The bottom inset of Fig. 10 shows the value of the polydispersity δ at the cloud point against the slant ratio R, for the same fixed parent density as in the main plot. The variation in δ is very small, between around 0.052 and 0.056, even though the parent shape changes quite dramatically from R = 1/4 to R = 4. This is in line with the expectation that δ is the main aspect of the size distribution that determines solid stability.
Returning now to the question of what determines whether a given particle size distribution will produce a critical S-SS transition, we broaden our investigation to a wider class of distributions, namely the Beta distributions. These are of the form f (σ) ∝ (σ −σ − ) a (σ + −σ) b in some interval σ − ≤ σ ≤ σ + , and zero otherwise. The values of the smallest and largest sizes σ − and σ + and the proportionality coefficient are again adjusted to make f (σ) normalized with unit mean and standard deviation δ. The advantage of Beta distributions is that with their two shape parameters, a and b, they are more flexible than e.g. the slanted top hat parents from above. In particular, they can interpolate from distributions with fairly sharp cutoffs at the extreme sizes -for low a and b, where in particular a = b = 0 gives back a top hat distribution -to ones with almost Gaussian shape (large a and b) where the cutoffs are in the far tails of the distribution. We can now proceed as above and solve the MFE critical point criterion to find out what shape parameters a and b produce critical S-SS transitions. We do this at fixed density, taking again n (0) = 1.133. The spinodal condition then fixes δ for given a and b, and the critical point gives one additional condition, so that we get a line of critical points in the a, b plane as shown in Fig. 11. What is noticeable is that the critical points lie significantly away from the line a = b where the parent density distribution is symmetric. This is also clearly visible in Fig. 12a, with the critical size distributions having distinct peaks to the right of the mean. One is led to ask whether there are other quantities, related but not identical to particle diameter, that would have more symmetric distributions. An obvious choice is the particle volume, which is proportional to v = σ 3 . As Fig. 12b shows, the distributions f v (v) = f (σ)/(3σ 2 ) are indeed much more nearly symmetric at criticality. This suggests that deviations from such symmetry, as measured by the skew indicate deviations from criticality, and conversely s v = 0 might be a reasonable approximate way of identifying critical size distributions. We have included the line in the (a, b) plane that results when we solve this condition (at the same values of δ as previously) in Fig. 11. The agreement with the critical line calculated directly from the MFE criticality condition is qualitatively quite good.
In particular, the criterion s v = 0 captures the fact that the critical size distributions are asymmetric when ex- Corresponding distributions of v = σ 3 , which is proportional to particle volume; these distributions are much more nearly symmetric. pressed in terms of particle size, with a > b throughout.
We have also calculated s v along the path through the phase diagram in Fig. 1 for top hat parents, and show the results in Fig. 13. One observes that the parent phase has relatively low s v at the S-SS transition, in agreement with the large values of the size fluctuation susceptibility χ. Likewise, in the SSS-SSSS transition, the phase that exhibits large χ and splits in a near critical fashion into two solids also has small s v .
Further support for the use of s v = 0 as an approximate criterion for criticality comes from the fact that s v can be written in terms of moment densities of σ as which entails exactly the moment densities ρ 0 , . . . , ρ 9 (though not all of them) that we would expect from the general discussion above. Nevertheless the criterion s v = 0 clearly remains approximate: for the slanted top hat parents, the results in Fig. 9 show that here the agree- ment with the full criticality criterion is less good. In particular, from s v = 0 we would predict that there are no critical points for slant ratio R < 1, whereas in fact critical size distributions exist down to R ≈ 0.81. The question of whether there is a more accurate yet still simple criterion for S-SS criticality remains open.

VI. DISCUSSION AND FUTURE WORK
In summary we have deployed tailored Monte Carlo simulation methods and moment free energy calculations to provide conclusive evidence that dense polydisperse spheres at equilibrium demix into coexisting fcc phases, with more phases appearing as the spread of diameters and the number density increase. Up to four coexisting phase were tracked, each of which contained a narrower distribution of particle sizes than is present in the system overall. Interestingly it was observed that for our systems the S-SS and the SSS-SSSS transitions are quasi-critical, characterised by a large correlation length for fluctuations in local particle size. By contrast the SS-SSS tran-sition was found to be strongly first order. To rationalize these observations, we investigated the features of the parental size distribution that control the character of solid demixing transitions. It was found that small skew in the parent distribution of particle volumes (s v ≈ 0) correlates well with the existence of a quasi-continuous transition, at least for one class of parental distribution shapes.
Whilst our results settle the matter of the true equilibrium behaviour, they leave open the question as to the extent to which this behaviour will be observable in experimental studies of polydisperse systems. Initial indications from recent experiments on colloid-polymer mixtures are that solid-solid demixing does not occur on the timescale of weeks. 47 Thus the best opportunity to see evidence may be to focus on regions of the phase diagrams where polydisperse solid(s) coexist with a fluid that can transport particles to their preferred solid phase. Additionally it would be interesting to try to manufacture a distribution of particle sizes that has s v ≈ 0 and then look for an increase in particle size fluctuations in the single solid region, even if the full transition itself is not seen.
As regards the questions that our results pose for further simulation and theoretical work, an interesting matter is that of the fate of the regions of multiple solid coexistence at high volume fraction. As the polydispersity δ is reduced, it seems clear that all the transition lines to multiple solids (S-SS, SS-SSS, SSS-SSSS etc) must converge on (but never quite reach) the monodisperse close packed limit at δ = 0, η = π/ √ 18 ≈ 74% since the close packed crystal will be unstable to any finite degree of polydispersity. We can also consider what happens if we fix the polydispersity δ > 0 and increase the parent volume fraction. The number of fractionated solids will increase without bound as the pressure increases, until at some volume fraction the pressure diverges and the system cannot be compressed further. The locus of these points in the phase diagram forms the infinite pressure line. Also this line must, as δ is decreased to zero, approach the monodisperse close packed limit η ≈ 74%.
An intriguing question is whether criticality can play a role in the approach to the close packed limit along the S-SS boundary. As we have seen, it is primarily the parent shape that controls the nature of the S-SS transition. Thus there may exist parent forms for which S-SS demixing is critical at or very near to the close packed limit, and it would be interesting to see whether a simple characterization of such narrow critical parent size distribution forms can be found.
Further solid phases are likely to arise in the phase diagram at values of δ beyond those that we have explored. For example one could imagine that at very large δ (for which the system separates into multiple coexisting phases) the smallest particles, rather than forming their own fcc phase, might instead secrete themselves in the interstitials of the fcc solid formed by the largest particles, thus potentially permitting the volume fraction to exceed 74%. Indeed this could be a mechanism whereby a unimodal parental distribution might produce familiar substitutionally-ordered phases, such as CsCl, or exotic phases such as AB 2 and AB 13 that can appear in binary colloid mixtures. 48 Investigating this question could -in principle -be tackled by simulation, but is probably beyond the present capabilities of the MFE calculations which are based on free energies that are reliable only for small to moderate δ.