Local order and crystallization of dense polydisperse hard spheres
LLocal order and crystallization of dense polydisperse hard spheres
Daniele Coslovich, Misaki Ozawa, and Ludovic Berthier
Laboratoire Charles Coulomb, Universit´e de Montpellier, CNRS, Montpellier 34095, France
Computer simulations give precious insight into the microscopic behavior of supercooled liquidsand glasses, but their typical time scales are orders of magnitude shorter than the experimentallyrelevant ones. We recently closed this gap for a class of models of size polydisperse fluids, which wesuccessfully equilibrate beyond laboratory time scales by means of the swap Monte Carlo algorithm.In this contribution, we study the interplay between compositional and geometric local orders ina model of polydisperse hard spheres equilibrated with this algorithm. Local compositional orderhas a weak state dependence, while local geometric order associated to icosahedral arrangementsgrows more markedly but only at very high density. We quantify the correlation lengths and thedegree of sphericity associated to icosahedral structures and compare these results to those for theWahnstr¨om Lennard-Jones mixture. Finally, we analyze the structure of very dense samples thatpartially crystallized following a pattern incompatible with conventional fractionation scenarios.The crystal structure has the symmetry of aluminum diboride and involves a subset of small andlarge particles with size ratio approximately equal to 0.5.
I. INTRODUCTION
Computer simulations play an important role inthe study of amorphous materials, since they provideparticle-scale resolution of their structural and dynam-ical properties. There is however a huge gap between thetimescales accessible in conventional simulations and inexperiments on molecular and polymeric liquids. Despitethe continuous increase of computing power, simulationstimescales are still about 8 orders of magnitudes shorterthan experimental ones. Numerical analysis is thereforelimited to studies of moderately supercooled liquids orpoorly annealed glasses [1, 2].Recently, we have developed a very efficient simulationsetup by applying the swap Monte Carlo algorithm [3–5] to realistic models of polydisperse particles [6–8].Some of these models can be equilibrated even beyondthe timescales accessible in the laboratory. Thanks tothis simulation approach, it becomes possible to scru-tinize under experimentally relevant conditions severaloutstanding issues concerning glass formation, such asthe entropy crisis [8], the kinetic stability of ultrastableglasses [9, 10], jamming [11], and the Gardner transi-tion [12–14]. These aspects are central in the currentdebate on the thermodynamic and dynamical propertiesof amorphous materials.Local structure is another feature that may provideimportant insight into the thermodynamic and dynamicbehavior of glass-formers [15, 16]. Multi-component mix-tures are characterized by local “compositional” order,which emerges due to preferential interactions betweendifferent chemical species [17]. Systems with continu-ous polydispersity might have even more complex formsof compositional ordering [18–20]. A large body of ex-perimental and simulation studies further demonstratedthat simple glass-formers, such as colloids [21], metallicglasses [22] as well as simple simulation models [23–25],display a tendency to form locally favored structures astemperature decreases or density increases. The symme-try of these local structures is often incompatible with the one of the underlying crystalline ground state [26],either because of compositional [27] or geometric frus-tration [28]. In other cases, however, the preferred localorder is the crystalline one [15], but the latter competeswith an alternate local structure. The influence of this“geometric” local order on the dynamics of supercooledliquids [23, 29] and on their rheological properties [30–32]has been the focus of several numerical studies. How-ever, these studies were limited to the moderately super-cooled regime and the ultimate role of local structure inthe overall picture of glass formation is still under de-bate [16, 33, 34].One outstanding issue of the local structure descriptionis that the spatial correlations associated to the geometricorder are fairly small. The correlation lengths associatedto locally favored structures remain small in the range oftemperature and density accessible to conventional simu-lations [35, 36]. This behavior contrasts with the appar-ent increase of dynamic correlations, as measured fromtime-dependent multi-point functions [1]. These discrep-ancies might be attributed to model dependence [29] orto the existence of different dynamic regimes not coveredby standard simulations [37], but also raise some doubtsabout the physical relevance of local geometric order inthe process of glass formation.In this paper, we address these issues by analyzing apolydisperse fluid equilibrated very deeply with the swapMonte Carlo algorithm. We carefully analyze the role ofcompositional fluctuations and identify the preferred ge-ometric motif of the system. We find that local composi-tional order increases smoothly with increasing density.On the other hand, the geometric order associated to thepreferred icosahedral order starts growing markedly onlyat sufficiently large volume fractions. We extract thecorrelation lengths associated to icosahedral structuresand compare these results with a representative binaryLennard-Jones mixture. We further elucidate the inter-play between compositional and geometric order in thepolydisperse system. Finally, we analyze a partially crys-tallized sample that we obtained during long swap Monte a r X i v : . [ c ond - m a t . s t a t - m ec h ] M a r Carlo simulation and rule out a conventional fractiona-tion scenario for the model at hand. Overall, our resultsshow that size polydisperse systems represent good glass-formers that are difficult to crystallize over a broad dy-namical range, and are characterized by only weak staticcompositional fluctuations. By comparison with earliermodels of glass-formers, they appear to contain much lesslocal order at equivalent degree of supercooling.This paper is organized as follows. In Section II wepresent the numerical methods we use. In Section III wepresent the results, which we organise into compositionalorder (Section III A), geometric order (Section III B), fol-lowed by an analysis of the crystal structure occasionallyfound in long simulations (Section III C). We concludethe paper in Section IV.
II. METHODS
We study systems composed of N polydisperse addi-tive hard spheres of diameter σ in three dimensions. Thediameter distribution is the same as in Ref [7], P ( σ ) = Aσ − , σ min ≤ σ ≤ σ max with σ min /σ max = 0 . A is a normalization constant. We use the average di-ameter σ = N (cid:80) Ni =1 σ i as the unit of length. In thefollowing we mostly focus on samples of N = 8000 parti-cles, but we also carried out simulations for N = 64000to check for finite size effects. The simulations were per-formed using the swap Monte Carlo algorithm [3, 5, 38]using the same setup as in [7, 8]. This simulation ap-proach is extremely efficient and allows one to equili-brate the fluid at least as deeply as conventional lab-oratory experiments on molecular liquids [6]. We notethat this is enabled by the combined optimization ofboth the Monte Carlo algorithm and the model param-eters, which must be chosen such that the system isrobust enough against crystallization or phase separa-tion [6]. Typical reference volume fractions of the sys-tem are onset of two-step relaxation ( φ onset (cid:39) .
56) andmode coupling crossover ( φ mct (cid:39) . φ mct [39, 40]. The initial configurations were prepared byfast Monte Carlo compressions of a low density fluid [41],which was subsequently equilibrated at the target pack-ing fraction. We have checked that the configurationsanalyzed in the following correspond to an equilibrium,disordered fluid, by carefully monitoring possible signs ofcrystallization or phase separation [18, 20] using the samestructural tools described in [7, 8]. One smaller sampleof N = 1000 particles showed clear signs of partial crys-tallizationous during long simulations at a high volumefraction ( φ = 0 . S w ( k ) = 1 N (cid:104) δρ w ( k ) δρ w ( − k ) (cid:105) , (1) δρ w ( k ) = ρ w ( k ) − (cid:104) ρ w ( k ) (cid:105) , (2) where (cid:104)· · · (cid:105) is the statistical average and ρ w ( k ) is theFourier transform of a weighted microscopic density ρ w ( k ) = (cid:88) j w j exp ( − i k · r j ) . (3)Here, the field w j is a generic particle property and entersas a weight in the calculation of the structure factor. Inthe following, we will consider various fields w j .The structure of simple mixtures and polydisperse par-ticle systems is often characterized by some preferred lo-cal arrangements, also known as locally favored struc-tures [16]. To identify this kind of geometric local or-der we perform a radical Voronoi tessellation using thevoro++ package [42]. In this construction, the totalvolume is partitioned into cells surrounding each par-ticle in the system. Cells are then classified accordingto their signature ( n , n , n , . . . ), where n q is the num-ber of faces of the cell with q vertices. Icosahedral localstructures correspond to cells with the (0 , ,
12) signa-ture. In the following, we will further distinguish be-tween particles that are at the center of an icosahedralstructure and icosahedral structures as a whole [23]. Acluster of neighboring icosahedral centers will be called“backbone”, while a cluster of neighboring icosahedralstructures will be called “domain”, see Section III B.
III. RESULTSA. Compositional order
In simple mixtures of particles, compositional order re-flects the tendency of particles to coordinate according totheir chemical species. A typical example is the presenceof chemical short range ordering in metallic alloys, whichreflects the tendency of particles to be surrounded byneighbors of a different chemical species [43]. This effectcan be quantified by computing partial structure factorsand their linear combinations [17]. This amounts to as-signing a weight w i equal to 1 or 0 depending on whetherparticle i belongs to a given species or not.In a polydisperse system, the relevant microscopicweight associated to compositional order varies contin-uously. Since our system is size-dispersed, we are actu-ally interested in the spatial structure associated to thediameter field. The simplest two-point correlation func-tion that captures the local fluctuations of σ i is the di-ameter structure factor S σ ( k ), defined by setting w i = σ i in Eq. (3). In Figure 1, we show S σ ( k ) for several vol-ume fractions ranging from the moderately dilute regime( φ = 0 .
5) to highly packed configurations ( φ = 0 . S σ ( k )strongly resembles the one of the total structure factor S ( k ) [11], shown in Figure 1(c) at the largest volumefraction.Since the diameter field by itself is weakly coupled tothe local structure, we consider instead the structure fac- (a) S σ ( k ) k φ =0.500 φ =0.550 φ =0.587 φ =0.597 φ =0.608 φ =0.618 φ =0.629 φ =0.635 φ =0.6400.000.020.040.060.08 0 5 10 (b) S δ σ ( k ) k0.00.51.01.5 0 5 10 (c) S ( k ) k 0.000.020.040.060.08 0.4 0.45 0.5 0.55 0.6 0.65 (d) MinMax S δ σ m i n ( k ) , S δ σ m a x ( k ) φ FIG. 1. (a) Structure factor of the diameter field S σ ( k ). (b) Structure factor of the diameter fluctuations field S δσ ( k ). (c)Total structure factor of the density field S ( k ) at φ = 0 .
64. (d) Absolute minimum and absolute maximum values (specified bythe arrows in (b)) of S δσ ( k ) as a function of φ . tor S δσ ( k ) associated to the fluctuating part of the diam-eter field. To this end we use w i = δσ i = σ i − σ , where σ is the average particle diameter. We note that S δσ ( k )is related to S σ ( k ) in a non-trivial way because of thepresence of cross terms S δσ ( k ) = S σ ( k ) + σ S ( k ) − σN Re[ (cid:104) δρ σ ( k ) δρ ( − k ) (cid:105) ] , (4)where Re[( · · · )] is the real part of ( · · · ) and ρ ( k ) is theFourier transform of the microscopic density with w j = 1.We expect S δσ ( k ) to capture local composition fluctua-tions better than S σ ( k ).The structure factor S δσ ( k ) is shown in Figure 1(b).Like the full diameter structure factor, S δσ ( k ) shows onlymild changes as a function φ . In contrast to S σ ( k ),however, S δσ ( k ) presents a more complex pattern anda marked suppression around wave numbers k ∗ ≈
7, cor-responding to typical length scales of particles of interme-diate sizes. Superficially, this dip might indicate an anti-correlation between diameter fluctuations over lengths oforder 2 π/k ∗ . This, in turn, suggests the presence of localcompositional order involving particles of different sizes,similar to the chemical ordering known in simple binarymixtures [17]. This dip gets more pronounced as φ in-creases, but its depth varies smoothly, see Figure 1(d).In addition, the smooth evolution of these structure fac-tors does not seem to correlate with the evolution of theglassy behavior of the system. Note finally that the fluc-tuations observed at the smallest wave number compat-ible with the simulation cell are not systematic and aredue to statistical noise.In polydisperse hard spheres, fluctuations of local vol-ume fraction do not occur only via variation of the localnumber density, but can also be mediated by size dis-persity. An appropriate correlation function to capture these fluctuations is the local volume structure factor S v ( k ), obtained by setting w i = v i = 4 πR i /
3, where R i = σ i /
2. We computed S v ( k ) and found that its over-all shape is similar to that of S σ ( k ), except at small k (not shown). In this regime, S v ( k ) behaves asymp-totically as the spectral density χ ( k ) [44], for which w i = πk (sin( kR i ) − ( kR i ) cos( kR i )) in three dimensions.This latter quantity, which provides direct insight intohyperuniform behavior in jammed packings [44, 45], alsoevolves smoothly by increasing density in equilibriumpolydisperse hard spheres [11]. Thus, we conclude thatfluctuations of both the local diameter and of the localvolume fraction evolve gradually with volume fractionand they reveal only weak compositional order. This con-clusion is consistent with the very smooth evolution thatwe have observed of partial structure factors and paircorrelation functions obtained by discretizing the parti-cle size distribution into discrete families (not shown).The overall conclusion of this analysis of compositionalorder is that the present system is a good glass-formerthat presents very weak fluctuations of the composition,and remains amorphous and well-mixed even in deep su-percooled states.We now look for an order parameter that detects moreprecisely fluctuations of local compositional order in theneighborhood of a given particle. Specifically, we con-sider fluctuations of the diameter within the first shellof neighbors. The neighbors of a particle are obtainedfrom the radical Voronoi tessellation, see Section II. Weintroduce the average neighbor diameter ˜ σ ( i ) of the i -thparticle ˜ σ ( i ) = 1 n b ( i ) n b ( i ) (cid:88) j =1 σ j , (5) . . . . . . . . . . Diameter0.80.91.01.11.21.3 A v e r age ne i ghbo r s d i a m e t e r φ =0.525 . . . . . . . . . . Diameter0.80.91.01.11.21.3 A v e r age ne i ghbo r s d i a m e t e r φ =0.597 . . . . . . . . . . Diameter0.80.91.01.11.21.3 A v e r age ne i ghbo r s d i a m e t e r φ =0.618 . . . . . . . . . . Diameter0.80.91.01.11.21.3 A v e r age ne i ghbo r s d i a m e t e r φ =0.629 . . . . . . . . . . Diameter0.80.91.01.11.21.3 A v e r age ne i ghbo r s d i a m e t e r φ =0.635 . . . . . . . . . . Diameter0.80.91.01.11.21.3 A v e r age ne i ghbo r s d i a m e t e r φ =0.640 FIG. 2. Conditional probability distribution P (˜ σ | σ ) = P (˜ σ, σ ) /P ( σ ) for several volume fractions, as indicated by the labels. where the sum runs over the n b ( i ) neighbors of the i -thparticle. We then compute the joint probability distribu-tion for ˜ σ and σ , P (˜ σ, σ ). To account for the polydisper-sity of the system, we actually focus on the conditionalprobability distribution P (˜ σ | σ ) = P (˜ σ, σ ) /P ( σ ). Thisdistribution is shown in Figure 2 for various volume frac-tions. At small and intermediate volume fractions, theaverage neighbor diameter is essentially independent of σ . This confirms that, at least for not too dense condi-tions, smaller particles tend to be surrounded on averageby larger ones and vice versa. This chemical local or-dering is however only apparent, as it simply means thateach particle (small and big) feels the same mean-fieldenvironment.For volume fraction above φ ≈ .
6, however, the dis-tribution P (˜ σ | σ ) presents an additional feature at inter-mediate values of σ . In the range 0 . < σ < .
1, weclearly see a spot displaying an excess of positive correla-tion between σ and ˜ σ , which becomes more marked withincreasing φ . We found that a similar excess correlationis also visible when the central particle is included in thedefinition of ˜ σ ( i ) in Eq. (5) (not shown). A possible in-terpretation of this excess correlation is that the systempresents local arrangements that involve particles of simi-lar sizes, forming more regular and symmetric structures. We will show in the next section that this feature is due tothe appearance of icosahedral structures, which providethe most regular local arrangements at high density. B. Geometric order
Recent numerical and experimental studies provideevidence of preferred geometric motifs in simple glass-formers [21, 23, 25]. These motifs include local icosahe-dral structures, which are the structural building blockof some metallic glass-formers [22], but also polytetra-hedral structures [46] or compositionally frustrated lo-cal crystalline structures [27, 47]. One important anddelicate question is to what extent these locally favoredmotifs correlate over larger length scales. Malins etal. [35, 36] have analyzed the structure factors associ-ated to locally favored structures in two Lennard-Jonesmixtures. The results of [35] indicate that icosahedraldomains are weakly correlated, even though the concen-tration of icosahedral structures is high enough that thedomains percolate. Recent works [34, 48] have furtheremphasized that static correlations are essentially localin the temperature regime accessible to conventional sim-ulations. In this section, we analyze the preferred local ( , , ) ( , , , ) ( , , , ) ( , , ) ( , , , ) ( , , , ) ( , , ) ( , , , ) (a) P e r c en t age φ =0.59 φ =0.64 φ =0.648 (crystal) 0 10 20 30 40 50 60 0.4 0.5 0.6 (b) P e r c en t age o f ( , , ) c e ll s φ BackboneDomains
FIG. 3. (a) Percentage of most frequent Voronoi cells for twoselected volume fractions, and for the sample which partiallycrystallized. (b) Fraction of particles forming the backbone oficosahedra, i.e. at the center of a (0,0,12) cell, (empty points)and belonging to icosahedra domains, i.e. at the center or onthe vertices of a (0,0,12) cell, (filled points) as a function ofvolume fraction. order of the polydisperse model and quantify its spatialextent over a very broad range of supercooling.In Figure 3 we show the fraction of most frequentVoronoi cells found around the mode-coupling crossover( φ mct ≈ .
6) and at the largest volume fraction φ = 0 . N = 1000)that partially crystallized during our swap MC simula-tions at a slightly larger volume fraction ( φ = 0 . , , φ mct and amount to about 10%of the total number of cells at large volume fractions. Wefind that at the largest φ about 60% of particles are in-volved in icosahedral domains, i.e. either being at thecenter or at the vertices of a (0 , ,
12) cell. The fractionof particles at the center of a (0,0,12) cell and that ofthe particles involved in icosahedral domains is shown inFigure 3(b) as a function of φ . Both quantities increasesteadily with increasing φ .To quantify the spatial correlations associated to icosa-hedral order, we introduce a microscopic field w i whichequals 1 if the i -th particle belongs to an icosahedralstructure and 0 otherwise. We further distinguish be-tween icosahedral backbones and icosahedral domains.For the former, w i = 1 only if the i -th particle is at thecenter of a (0 , ,
12) cell. For the latter, w i = 1 if the i -th particle is the center of a (0 , ,
12) cell or at the ver-tices of a (0 , ,
12) cell. In contrast to [35], we normalizethe corresponding structure factors S b ( k ) and S d ( k ) bythe average number of particles forming icosahedral back-bones and domains, respectively, and not by N . This isdone to remove the trivial part of the state dependenceof the correlation functions.The resulting structure factors are shown in Figure 4. (a) D o m a i n s t r u c t u r e f a c t o r , S d ( k ) k φ =0.640 φ =0.635 φ =0.629 φ =0.618 φ =0.608 φ =0.597 φ =0.587 φ =0.550 φ =0.500 φ =0.450 0.5 1 1.5 2 2.5 3 3.5 0 2 4 6 8 10 (b) B a ck bone s t r u c t u r e f a c t o r , S b ( k ) k FIG. 4. Structure factors of (a) the icosahedral domains S d ( k )and (b) the icosahedral backbone S b ( k ) for several volumefractions. (a) (b) FIG. 5. (a) Particles in icosahedral domains for a configura-tion at φ = 0 .
64. Red particles are at the center of a (0,0,12)cell, orange ones are at the vertices. (b) Icosahedral back-bones of the configuration in (a). The white bonds connectneighbors centers of (0,0,12) cells.
The domain structure factor S d ( k ) shows a peak at k = 0at any density (including the unstructured, non-glassyfluid at moderate density), whose height decreases withincreasing φ . This indicates that the icosahedral do-mains are only weakly correlated. The presence of apeak around k = 0 thus merely reflects the icosahedralform factor but not a nontrivial large-scale correlation.By contrast, correlations in the icosahedral backboneare nearly absent at low density and start to increasesmarkedly beyond φ ≈ .
60. The icosahedral backbonethus reveals subtle but nontrivial changes in the structureof the fluid. Typical snapshots of icosahedral domainsand backbones at high density are shown in Figure 5(a)and (b), respectively.To determine the correlation length associated to icosa-hedral domains and backbones, we fitted to the low k portion of the structure factors S d ( k ) and of S b ( k ), re-spectively, using the Ornstein-Zernicke function, S α ( k ) = S α (0)1 + ( ξ α k ) , (6)where α = d or b. We restricted our fits to k < . (a) Polydisperse HSMCT ξ φ Gyration radiusBackboneDomainPTS 0.00.20.40.60.81.01.2 0.8 1.2 1.6 (b)
Wahnström LJMCT ξ FIG. 6. Length scales of icosahedral order. Domain correla-tion length ξ d , backbone correlation length ξ b and backbonegyration radius R g as a function of (a) φ in the polydispersesystem and (b) 1 /T in the Wahnstr¨om mixture. Open circlesin panel (a) represent the scaled PTS length ξ PTS / are consistent with those obtained by manually rescal-ing the structure factors so as to optimize data collapseat small k . In Figure 6 we show the variation of thecorrelation lengths ξ α as a function of φ for both icosa-hedral domains and backbones. The domain correlationlength remains approximately constant around 1 inter-particle distance throughout the studied range of volumefraction. This confirms that the peak observed around k = 0 in S d ( k ) has a trivial origin. By construction,the backbone correlation length ξ b is smaller than ξ d ,but it increases markedly (by about a factor 3) uponincreasing φ beyond the mode-coupling crossover den-sity. The maximal value reached, ξ d ∼ .
5, remains how-ever quite modest and this small growth of the corre-lation length should be contrasted with the striking vi-sual impression provided by the snapshot in Figure 5(a)where the configuration appears full of icosahedral struc-tures. For comparison, we also include the static point-to-set [49, 50] correlation length obtained in Ref. [8],scaled to roughly match the backbone length around theMCT crossover. We see that the relative increase of ξ P T S over the studied range of densities is qualitatively simi-lar to the one of the icosahedral backbone. It would beinteresting to further investigate the connection betweenorder agnostic correlations, such as point-to-set correla-tions, and locally favored structures, as already suggestedin Refs. [29, 51, 52].A different approach to characterize the extent of icosa-hedral order is to compute the gyration radius R g = M M (cid:88) j =1 ( r j − r b ) / , (7)where r b = M (cid:80) Mj =1 r j . M is the number of connectedicosahedral particles, see [35]. The results for the back- . . . All(0,0,12) A s phe r i c i t y φ (a) All(0,0,12) (b) . . . . . . Diameter0.80.91.01.11.21.31.4 A v e r age ne i ghbo r s d i a m e t e r (a) FIG. 7. (a) Conditional probability distribution P ico (˜ σ | σ ) foricosahedra centers at φ = 0 .
64. The straight line is a guidefor the eyes. (b) Asphericity of all Voronoi cells (circles) and(0 , ,
12) cells (squares) as a function of φ . bone gyration radius are included in Figure 6. R g in-creases by increasing φ following the trend of the corre-lation function. We also found that the domain gyrationradius tends to substantially overestimate the correla-tion length, consistent with the results of [35]. We notethat R g is not well defined as soon as the cluster perco-lates through the system and therefore we do not showthese results here. Moreover, percolation of icosahedraldomains has no obvious connection with the (swap) dy-namics of the system, which evolves smoohtly throughoutthe studied temperature regime [8].It is interesting to compare the behavior of the modelat hand with one of the Wahnstr¨om Lennard-Jones mix-ture [53], which is a binary glass-former displaying afairly strong icosahedral ordering [23]. We computed thebackbone correlation length and gyration radius for thismodel, using the same parameters, density and units asin [23, 35, 53]. Both quantities increase markedly by de-creasing temperature already above the mode-couplingcrossover, see Fig 6(b). The domain correlation lengthalso increases slightly at sufficiently low temperature, butthe absolute values of all these lengths remain small, be-cause geometric frustration is strong in this system [54].The behavior of this Lennard-Jones mixture is thus qual-itatively similar to the one of the polydisperse system,but the latter differs for two main reasons. First, thestructure of the polydisperse system remains highly dis-ordered in the moderately supercooled regime and onlystarts to develop some geometric order well beyond thecrossover volume fraction φ mct . Thus, the increase of thelocal order is shifted to a considerably deeper degree ofsupercooling compared to the Wahnstr¨om Lennard-Jonesmixture. Second, we recently showed [11] that icosahe-dral structures are actually more distorted than in simplemixtures due to local compositional fluctuations.We now shed some light on the interplay between geo-metric and compositional order. In Figure 2 we noticedthe emergence of an excess correlation in the P (˜ σ | σ ) dis-tribution, which becomes increasingly visible at largerdensities. We argue that this correlation is due to thegrowing icosahedral order. In Figure 7 we show the dis-tribution P ico (˜ σ | σ ) restricted to particles at the center ofan icosahedron, evaluated at the largest volume fraction.The clear correlation between σ and ˜ σ corresponds nicelyto the feature observed in the full distribution P (˜ σ | σ ).A plausible interpretation of the excess correlation in P (˜ σ | σ ) is that icosahedra tend to be more regular andspherical than other structures. At sufficiently high den-sity, a high degree of sphericity also likely implies thatparticles involved in the local structure have similar sizes,thus a correlation between σ and ˜ σ . Note that this trendby itself does not imply fractionation, which should leadto a much sharper structural change.To confirm that icosahedra are indeed the most regularstructures in the model, we measure the asphericity ofthe Voronoi cell by computing the normalized standarddeviation of the distances from the i -th particle s i = 1˜ r i (cid:118)(cid:117)(cid:117)(cid:116) n b ( i ) n b ( i ) (cid:88) j =1 ( r ij − ˜ r i ) , (8)where ˜ r i = ( (cid:80) n b ( i ) j =1 r ij ) /n b ( i ) and r ij are respectively theaverage nearest neighbors distance from the i -th particleand the distance between particles i and j . We expectthis measure to be closely related to other measures ofregularity (i.e., tetrahedricity) previously introduced inthe study of simple particulate systems [55, 56]. First,we observe from Figure 7(b) that the average asphericityof the Voronoi cells decreases as the system gets denser,as expected. We find that icosahedra are appreciablymore regular than the other structures. By restrictingourselves to the most spherical structures, we find thatthe proportion of icosahedra is significantly higher thanin the bulk. Specifically, we computed the Voronoi cellstatistics for the 2% most spherical particles at φ = 0 . C. Partial crystallization
The structural analysis carried out so far concernsequilibrium, disordered fluid states. However, like any su-percooled fluid, the model at hand is thermodynamicallymetastable and sufficiently long simulations with theswap Monte Carlo algorithm may trigger a fluctuationtowards the ground state. In systems with sufficientlyhigh polydispersity, reaching the crystalline ground state P ( σ ) σ CrystalDisorderedOverall (a) (b)
FIG. 8. (a) Snapshot of a partially crystallized sample at φ = 0 . P ( σ ) inthe crystalline region (red curve with shaded area), in the dis-ordered region (black curve), and in the bulk (dashed curve). may involve fractionation into families of particles char-acterized by similar diameters [19]. However, this processmay take extremely long times, and alternate crystalliza-tion scenarios, not involving fractionation, have also beenobserved [57, 58].In this section, we focus on a sample of N = 1000 parti-cles at a volume fraction φ = 0 . φ > . φ = 0 . − .
645 [8]),and after long simulations times, see below. These sam-ples were excluded of all published analysis of dense su-percooled liquid states [7, 8]. In the following, we pro-vide some quantitative insight into the structure of thepartially crystallized sample and highlight the differencescompared to the normal fluid states. Note that from theviewpoint of glass transition studies, these crystallizationevents are only problematic when they occur on timescales comparable to the structural relaxation time τ α ,which controls the equilibration of density fluctuations.In a related study we have shown that it is possible to al-leviate this problem by introducing non-additive pair in-teractions, which help to suppress crystallization eventseven further [6].Visual inspection of the particles configurations in Fig-ure 8(a) immediately shows that the system phase sep-arates into a disordered and a crystalline region. Wefound that particles in the crystalline region are clearlyassociated to (0 , ,
6) and (0 , , ,
8) Voronoi cells. Theproportions of these cells can thus be used as a marker of (a) (b)
FIG. 9. Structure of a representative (0,0,12,8) cell foundin the crystalline region. In panel (a) spheres are drawn toscale. In panel (b) spheres are scaled to half of their size andbonds are added between neighboring particles to highlightthe hexagonal symmetry of the crystal. In (b) small and largeparticles are shown as yellow and red spheres, respectively. the system instability, since they increase markedly uponpartial crystallization, see Figure 3. Particles located atthe center of (0 , ,
6) and (0 , , ,
8) cells are the small-est and the largest particles, respectively, and are high-lighted accordingly in Figure 8(a). Thus, while the disor-dered region has a local polydispersity similar to the oneof the homogeneous fluid, the crystalline one comprisesonly a subset of the particles, and is completely devoidof particles of intermediate sizes. This is demonstratedin Figure 8(b), where we compare the overall diameterdistribution P ( σ ) to the one measured in the crystallineregion, i.e. for particles at the center of either (0 , ,
6) or(0 , , , , with small and large particles playing therole of B and Al atoms, respectively. The crystal struc-ture has an hexagonal symmetry and is formed by in-terleaved layers of small and large particles. The typicalshape of the first coordination shell around a large parti-cle of the crystal is illustrated in Figure 9, where we showthe structure of a (0,0,12,8) Voronoi cell. The typicalsize ratio γ ≈ . γ = 0 .
58) of AlB -forming binary hard colloids [59] and lies in the stabilityrange (0 . ≤ γ ≤ .
59) expected from theoretical stud-ies of binary hard spheres [60, 61]. Finally, we note thatthe Voronoi tessellation of the AlB lattice comprises in-deed only (0 , , ,
8) and (0 , ,
6) cells, centered aroundAl and B atoms respectively, see e.g. [62].To investigate the structure of the partially crystal-lized sample more quantitatively, we compute the distri-bution of the compositional order parameter (˜ σ, σ ). Incontrast to the equilibrium fluid states studied in the pre-vious sections, the partially crystallized sample displaysa complex distribution P (˜ σ, σ ) characterized by multiplespots. Two of these spots, marked by arrows in Figure 10,are clearly associated with the crystalline region of thesample and involve the smallest and the biggest particlesin the sample. Note that, to enhance visualization, weshow here the full distribution P (˜ σ, σ ) instead of P (˜ σ | σ ). . . . . . . . . . . Diameter0.80.91.01.11.21.3 A v e r age ne i ghbo r s d i a m e t e r (0,3,6) (0,0,12,8) FIG. 10. Joint probability density P (˜ σ, σ ) for the partiallycrystallized sample at φ = 0 . By computing the Voronoi statistics restricted to specificranges of σ and ˜ σ , we confirmed that the crystalline spotsindicated by arrows in Figure 10 correspond to (0 , , , , ,
8) Voronoi cells. Note, however, that someof the smallest particles have values of ˜ σ comparable tothe ones found in fluid states. Visual inspection of theparticle configurations indicates that these small parti-cles populate the disordered portion of the sample.The distribution P (˜ σ, σ ) displays two additional spotslocated at intermediate values of σ and characterized by adistinct positive correlation. One of these spots, between0 . < σ < .
1, is associated to icosahedral structuresand is similar to the one found in the fluid samples. Bycontrast, no clear structural signature stands out in thesecond spot between 1 . < σ < .
35 [63]. Surprisingly,we found that particles that contribute to the two cen-tral spots of P (˜ σ, σ ) are not spatially segregated fromone another and are characterized by very similar localpolydispersities [64]. Thus, the existence of a positivecorrelation between σ and ˜ σ does not imply per se frac-tionation and may be attributed instead to a subtle localgeometric ordering. Overall, our results confirm that, inpractice, crystallization and phase separation in polydis-perse hard spheres follow a more complex pattern thanpredicted by existing theoretical models [19].Finally, we estimated the crystallization time by per-forming additional simulations of independent samplesof 1000 particles at several volume fractions. Crystal-lization was detected by inspecting the evolution of thepercentage of (0,0,12,8) cells, which typically fluctuatesbetween 0.5% and 2% for fluid states and exceeds 5% inpartly crystallized samples. For φ = 0 .
635 and φ = 0 . τ α and 200 τ α , respectively, where τ α is the structural relaxation time measured from the selfintermediate scattering function [8]. Thus, in this den-sity regime and for this system size, our swap MonteCarlo simulations can safely probe the structure of themetastable equilibrium fluid. For φ = 0 . τ x = 6 × Monte Carlo steps, correspond-ing to about 30 τ α . A conservative upper bound to thecrystallization rate 1 / ( τ x V ) is thus 3 × − in reducedunits [65]. We emphasize that this value may still bestrongly affected by finite size effects induced by the pe-riodic boundary conditions. IV. CONCLUSIONS
We characterized the local structure of a fluid of poly-disperse hard spheres over a wide range of volume frac-tions, where conventional computer simulations fail toequilibrate. We showed that local compositional orderincreases smoothly with increasing the volume fraction,with little correlation with the glassy evolution of the sys-tem. Concomitantly, local geometric order associated toicosahedral particle arrangements grows steadily at verylarge volume fractions. We extracted correlation lengthsassociated to icosahedral structures using weighted struc-ture factors and the gyration radius. These correlationlengths increase appreciably only at large packing frac-tions and their absolute values remain small over the en-tire glassy regime we could probe. It is interesting tocompare this behavior with the results for the Wahn-str¨om Lennard-Jones mixture, which is a glass-formingmodel displaying a large amount of icosahedral struc-tures. The measured growth of icosahedral length scalesis qualitatively similar in the two models, but the glassyregime explored in the Wahnstr¨om mixture is much nar-rower. Thus, in the regime probed by conventional sim-ulations, the local structure of the polydisperse modelappears highly disordered, as is the case in simple bi-nary mixtures of hard [34] or quasi-hard spheres [29].This indicates that the role of local structure is not onlysystem-dependent [29], but also highly state-dependent.Finally, we characterized the structure of a partiallycrystalline sample of N = 1000 particles obtained dur-ing long simulations at a packing fraction φ = 0 . crystal comprising only the smallest and largestparticles of the sample. We emphasize that crystalliza- tion only occurs at a packing fraction that lies beyondthe estimated laboratory glass transition, in a regimewhere crystal growth, driven by physically realistic dy-namics, would be extremely slow on observational timescales [66]. The fact that swap Monte Carlo simulationsinvolve non-physical moves, which accelerate sampling ofconfiguration space, makes it difficult to infer the abso-lute glass-forming ability of the models using ordinarydynamics, but we expect that the relative trends acrosssystems, see Ref. [6], will be preserved. Investigations totackle these issues are currently under way. Comparingthe glass-forming ability of these models with experimen-tal systems represents a major challenge that is left forfuture investigations.Continuously polydisperse systems can be regardedas an extreme case of multicomponent mixtures. Atfirst glance, these systems may appear peculiar [67],for instance because of their formally infinite mixingentropy [68, 69]. However, their glassy phenomenol-ogy strongly resembles the one of conventional glass-formers [8]. Moreover, our detailed structural character-ization demonstrates that the local structure of polydis-perse system at hand shows qualitatively similar featuresas other representative glass-formers such as colloidal [21]and metallic glasses [22], which also display growingicosahedral order. Thus, overall, our results confirm thatcontinuous polydisperse systems can be regarded as goodmodels to study the glass transition. Whether more gen-eral classes of glass-formers, such as molecular or poly-meric liquids, display or not a pronounced local geometricorder remains an open question to be addressed in futurenumerical studies. ACKNOWLEDGMENTS
We thank G. Tarjus for enduring discussions on therole of icosahedral order and its length scales. Theresearch leading to these results has received fundingfrom the European Research Council under the EuropeanUnions Seventh Framework Programme (No. FP7/2007-2013)/ERC Grant Agreement No.306845. This workwas supported by a grant from the Simons Founda-tion (No. 454933, Ludovic Berthier). Data relevant tothis work have been archived and can be accessed athttps://doi.org/10.5281/zenodo.1183325. [1] L. Berthier and G. Biroli, Rev. Mod. Phys. , 587(2011).[2] K. Binder and W. Kob, Glassy materials and disor-dered solids: An introduction to their statistical mechan-ics (World Scientific, 2011).[3] D. Gazzillo and G. Pastore, Chem. Phys. Lett. , 388(1989).[4] D. Frenkel and B. Smit,
Understanding molecular simu-lation: from algorithms to applications , Vol. 1 (Academic press, 2001).[5] T. S. Grigera and G. Parisi, Phys. Rev. E , 045102(2001).[6] A. Ninarello, L. Berthier, and D. Coslovich, Phys. Rev.X , 021039 (2017).[7] L. Berthier, D. Coslovich, A. Ninarello, and M. Ozawa,Phys. Rev. Lett. , 238002 (2016).[8] L. Berthier, P. Charbonneau, D. Coslovich, A. Ninarello,M. Ozawa, and S. Yaida, Proc. Natl. Acad. Sci. , , 36003 (2017).[10] L. Berthier, P. Charbonneau, E. Flenner, and F. Zam-poni, Phys. Rev. Lett. , 188002 (2017).[11] M. Ozawa, L. Berthier, and D. Coslovich, SciPostPhysics , 027 (2017).[12] L. Berthier, P. Charbonneau, Y. Jin, G. Parisi,B. Seoane, and F. Zamponi, Proc. Natl. Acad. Sci USA , 8397 (2016).[13] Y. Jin and H. Yoshino, Nat. Comm. , 14935 (2016).[14] C. Scalliet, L. Berthier, and F. Zamponi, Phys. Rev.Lett. , 205501 (2017).[15] H. Tanaka, Eur. Phys. J. E , 113 (2012).[16] C. P. Royall and S. R. Williams, Phys. Rep. , 1(2015).[17] A. Bhatia and D. Thornton, Phys. Rev. B , 3004 (1970).[18] N. B. Wilding and P. Sollich, EPL (Europhysics Letters) , 219 (2004).[19] P. Sollich and N. B. Wilding, Phys. Rev. Lett. ,118302 (2010).[20] N. B. Wilding and P. Sollich, J. Chem. Phys. , 224102(2010).[21] M. Leocmach and H. Tanaka, Nat. Comm. , 974 (2012).[22] A. Hirata, L. Kang, T. Fujita, B. Klumov, K. Matsue,M. Kotani, A. Yavari, and M. Chen, Science , 376(2013).[23] D. Coslovich and G. Pastore, J. Chem. Phys. , 124504(2007).[24] R. Soklaski, V. Tran, Z. Nussinov, K. Kelton, andL. Yang, Philosophical Magazine , 1212 (2016).[25] C. P. Royall and W. Kob, J. Stat. Mech. , 024001(2017).[26] W. Lechner and C. Dellago, J. Chem. Phys. , 114707(2008).[27] P. Crowther, F. Turci, and C. P. Royall, J. Chem. Phys. , 044503 (2015).[28] G. Tarjus, S. A. Kivelson, Z. Nussinov, and P. Viot, J.Phys.: Condens. Matt. , R1143 (2005).[29] G. M. Hocky, D. Coslovich, A. Ikeda, and D. R. Reich-man, Phys. Rev. Lett. , 157801 (2014).[30] R. Pinney, T. B. Liverpool, and C. P. Royall, J. Chem.Phys. , 234501 (2016).[31] S. Feng, L. Qi, L. Wang, S. Pan, M. Ma, X. Zhang, G. Li,and R. Liu, Acta Mater. , 236 (2015).[32] J. Ding, S. Patinet, M. L. Falk, Y. Cheng, and E. Ma,Proc. Natl. Acad. Sci. , 14052 (2014).[33] G. Tarjus, Dynamical Heterogeneities in Glasses, Col-loids, and Granular Media , 39 (2011).[34] P. Charbonneau and G. Tarjus, Phys. Rev. E , 042305(2013).[35] A. Malins, J. Eggers, C. P. Royall, S. R. Williams, andH. Tanaka, J. Chem. Phys. , 12A535 (2013).[36] A. Malins, J. Eggers, H. Tanaka, and C. P. Royall, Fara-day Discussions , 405 (2014).[37] C. P. Royall, A. Malins, A. J. Dunleavy, and R. Pinney,J. Non-Cryst. Solids , 34 (2015).[38] P. Sindzingre, C. Massobrio, G. Ciccotti, and D. Frenkel,Chem. Phys. , 213 (1989).[39] G. Brambilla, D. El Masri, M. Pierno, L. Berthier,L. Cipelletti, G. Petekidis, and A. B. Schofield, Phys.Rev. Lett. , 085703 (2009). [40] E. Zaccarelli, S. M. Liddle, and W. C. Poon, Soft Matter , 324 (2015).[41] L. Berthier and T. A. Witten, Phys. Rev. E , 021502(2009).[42] C. H. Rycroft, Chaos , 041111 (2009).[43] J. H. Na, M. D. Demetriou, M. Floyd, A. Hoff, G. R.Garrett, and W. L. Johnson, Proc. Natl. Acad. Sci. ,9031 (2014).[44] Y. Wu, P. Olsson, and S. Teitel, Phys. Rev. E , 052206(2015).[45] C. E. Zachary, Y. Jiao, and S. Torquato, Phys. Rev. E , 051308 (2011).[46] A. V. Anikeenko and N. N. Medvedev, Phys. Rev. Lett. , 235504 (2007).[47] R. G. Della Valle, D. Gazzillo, R. Frattini, and G. Pas-tore, Phys. Rev. B , 12625 (1994).[48] M. Wyart and M. E. Cates, Phys. Rev. Lett. , 195501(2017).[49] J.-P. Bouchaud and G. Biroli, J. Chem. Phys. , 7347(2004).[50] G. Biroli, J.-P. Bouchaud, A. Cavagna, T. S. Grigera,and P. Verrocchio, Nature Physics , 771775 (2008).[51] J. Russo and H. Tanaka, Proc. Nat. Acad. Sci. , 6920(2015).[52] S. Yaida, L. Berthier, P. Charbonneau, and G. Tarjus,Phys. Rev. E , 032605 (2016).[53] G. Wahnstrom, Phys. Rev. A , 3752 (1991).[54] F. Turci, G. Tarjus, and C. P. Royall, Physical ReviewLetters , 215501 (2017).[55] W. P. Krekelberg, V. Ganesan, and T. M. Truskett, J.Chem. Phys. , 214502 (2006).[56] A. Anikeenko, N. Medvedev, and T. Aste, Phys. Rev. E , 031101 (2008).[57] L. Fern´andez, V. Martin-Mayor, and P. Verrocchio,Phys. Rev. Lett. , 085702 (2007).[58] L. Fern´andez, V. Mart´ın-Mayor, B. Seoane, and P. Ver-rocchio, Phys. Rev. E , 021501 (2010).[59] P. Bartlett, R. H. Ottewill, and P. N. Pusey, Phys. Rev.Lett. , 38013804 (1992).[60] X. Cottin and P. A. Monson, J. Chem. Phys. ,33543360 (1995).[61] L. Filion and M. Dijkstra, Phys. Rev. E , 046714(2009).[62] A. Travesset, Phys. Rev. Lett. , 115701 (2017).[63] The most frequent signatures in this spot are (0 , , , , , , σ and of their reespective neighbors.[65] The time unit is given by one Monte Carlo step, compris-ing N attempts to either displace a particle or swap theidentities of a pair ( i, j ) of particles such that | σ i − σ j | < . , 51124 (2009).[67] V. Lubchenko and P. G. Wolynes, arXiv preprintarXiv:1710.02857 (2017).[68] D. Frenkel, Mol. Phys. , 2325 (2014).[69] M. Ozawa and L. Berthier, J. Chem. Phys.146