Glassy dynamics of a binary Voronoi fluid: A mode-coupling analysis
Céline Ruscher, Simone Ciarella, Chengjie Luo, Liesbeth M. C. Janssen, Jean Farago, Jörg Baschnagel
aa r X i v : . [ c ond - m a t . d i s - nn ] A ug Glassy dynamics of a binary Voronoi fluid: Amode-coupling analysis
C Ruscher , , S Ciarella , C Luo ‡ , L M C Janssen , JFarago , J Baschnagel Universit´e de Strasbourg, Institut Charles Sadron, CNRS–UPR22, 23 rue duLoess, BP 84047, 67034 Strasbourg Cedex 2, France Department of Physics and Astronomy and Quantum Matter Institute,University of British Columbia, Vancouver BC V6T 1Z1, Canada Theory of Polymers and Soft Matter, Department of Applied Physics,Eindhoven University of Technology, P.O. Box 513, 5600MB Eindhoven, TheNetherlandsE-mail: [email protected]
Abstract.
The binary Voronoi mixture is a fluid model whose interactions arederived from the Voronoi-Laguerre tessellation of the configurations of the system.The resulting interactions are local and many-body. Here we perform molecular-dynamics (MD) simulations of an equimolar mixture that is weakly polydisperseand additive. For the first time we study the structural relaxation of this mixturein the supercooled-liquid regime. From the simulations we determine the time-and temperature-dependent coherent and incoherent scattering functions for alarge range of wave vectors, as well as the mean-square displacements of bothparticle species. We perform a detailed analysis of the dynamics by comparing theMD results with the first-principles-based idealized mode-coupling theory (MCT).To this end, we employ two approaches: fits to the asymptotic predictions of thetheory, and fit-parameter-free binary MCT calculations based on static-structure-factor input from the simulations. We find that many-body interactions of theVoronoi mixture do not lead to strong qualitative differences relative to similaranalyses carried out for simple liquids with pair-wise interactions. For instance,the fits give an exponent parameter λ ≈ .
746 comparable to typical valuesfound for simple liquids, the wavevector dependence of the Kohlrausch relaxationtime is in good qualitative agreement with literature results for polydisperse hardspheres, and the MCT calculations based on static input overestimate the criticaltemperature, albeit only by a factor of about 1.2. This overestimation appears tobe weak relative to other well-studied supercooled-liquid models such as the binaryKob–Andersen Lennard-Jones mixture. Overall, the agreement between MCTand simulation suggests that it is possible to predict several microscopic dynamicproperties with qualitative, and in some cases near-quantitative, accuracy basedsolely on static two-point structural correlations, even though the system itself isinherently governed by many-body interactions.
Keywords : Voronoi liquid, binary mixture, glass transition, molecular-dynamicssimulations, mode-coupling theory
Submitted to:
J. Phys.: Condens. Matter ‡ S Ciarella and C Luo have contributed equally to this work. lassy dynamics of a binary Voronoi fluid
1. Introduction
Disordered materials, such as polymers, metallic alloys,and polydisperse colloidal suspensions, are of hugepractical interest as they can be designed with specificmechanical, optical, or thermal properties. At lowdensity or high temperature these materials are in theliquid state. Provided the liquid can be supercooledwithout undergoing structural ordering, the dynamicsstrongly slows down upon cooling or increasing thedensity, shifting the time scale for viscous flow toever longer times. Ultimately, the glass transition isreached, beyond which structural relaxation can nolonger occur on experimental time scales. Such systemsare then in a nonequilibrium solid-like state where theyexhibit mechanical rigidity but, contrary to crystallinematerials, they lack any long-range order. Developinga microscopic understanding of the nature of the glassystate and the glass transition is one of the challengingproblems in condensed matter physics [1, 2].In the dense liquid phase, relaxation takes placethrough cooperative rearrangements of (groups of)neighboring particles. That is, in order for a particleto move, its neighbors must also move, and hencethe local particle environment plays an importantrole in the dynamics. A way to access informationon the neighborhood of a particle is to apply aVoronoi tessellation [3]. Voronoi tessellation is ageometric partitioning of space into contiguous cellswhose volume can be thought of as the zone of influenceof a given particle. Such tessellations have beenextensively used for glass-forming [4, 5] or granularsystems [6, 7, 8], mostly as a tool to define free volumeor to obtain local geometric information. For instance,Morse and Corwin [6, 7] emphasize the geometricnature of the jamming transition in granular systemsby showing that a large set of geometrical observables(surface area, aspect ratio, standard deviation of thevolume, etc.) extracted from Voronoi tessellationshows a marked signature at the jamming point.A similar observation was made by Rieser et al [8] who found a strong signature of jamming in aquantity related to the relative free volume of theparticles. These observations highlight the importanceof Voronoi tessellation to get a deeper level of structuralinformation which is either not contained or toostrongly averaged in the usual two-point correlationfunctions, like the radial distribution function or staticstructure factor that are known to vary only weakly onapproach to the glass transition.Voronoi tessellation offers more than only adiagnostic tool, however; it also provides the basis fora new class of complex liquid models. During thepast five years, two new models have emerged whoseinteractions are intrinsically many-body and derivedfrom the inherent geometrical properties of the Voronoi tessellation: the “Voronoi liquid” introduced by someof us [9] and the “self-propelled Voronoi (SPV) model”proposed by Bi et al [10]. The SPV model aimsat describing cell motility and cell-cell interactions inconfluent tissues. One major achievement of the SPVmodel has been the identification of a structural orderparameter, the shape index, which depends only theperimeter and area of the Voronoi cell, and whichidentifies, for given single-cell motility and persistencetime, a liquid-to-solid transition reminiscent of theglass transition. The SPV model has also found usein understanding collective cell phenomena such asthe epithelial-to-mesenchymal transition—a key stepin the propagation of cancer cells [11]—and in thedesign of a new generation of bioinspired materials,such as tunable photonic fluids [12]. Moreover, themodel can also shed light on the glass transitionin active matter [13]. Recently, Sussman et al [14] studied a passive version of the SPV model.Their findings differ from the usual phenomenologyof glass formers as they observed a sub-Arrheniusbehavior of the relaxation time—i.e. an anomalousfragility which thus far has been found in only a verylimited number of systems [15]—and a high density ofcollective low-frequency vibrational modes associatedwith low-temperature energy minima. The specificmany-body nature of the interactions is at the core ofthis anomalous dynamics, meaning that going beyondthe usual pairwise potentials could broaden the viewon the “stylist facts” [16] commonly associated withthe glass transition phenomenology. This has alsobeen a main motivation for the introduction of theVoronoi liquid. From its conception the Voronoi liquidis a passive fluid model. In many ways, perhapssurprisingly, the model behaves like an ordinary simpleliquid regarding traditional structural and dynamicalcorrelation functions [9]. However, the Voronoi fluidalso has some unique specificities. The most strikingone is arguably that the potential energy of an N + 1 particle system, E p ( r , · · · , r N , r N +1 ), becomesequal to that of N particle system E p ( r , · · · , r N ),if r N +1 → r N [9]. In this sense, the potential is“hypersoft”. This property does not compromise thestability of the liquid because the interactions arelocally repulsive and the superposition of two particleshas a finite energy cost. Hypersoftness has, however,an impact in situations where the dynamics is slow,e.g. for the crystalline phase. At low temperature themonodisperse Voronoi liquid forms a bcc crystal thatis “plastic” in that the particles can diffuse freely inthe solid without destroying the crystalline structure[17]. A further striking property of the Voronoi liquidis an anomalous scaling of the sound attenuation rate( ∝ q instead of ∝ q with q being the modulus ofthe wave vector) at mesoscopic scales and a shift of lassy dynamics of a binary Voronoi fluid q -values withrespect to a standard Lennard-Jones (LJ) system [18].This specific behavior can be attributed to a very weakresistance to shear deformations at high frequency.For the Voronoi liquid the product of the infinite-frequency shear modulus ( G ∞ ) and the isothermalcompressibility ( χ T ) is much smaller than 1, whereas G ∞ χ T ∼
2. Model and simulation method
Consider a system of N point particles at positions r j ( j = 1 , . . . , N ) in a three-dimensional volume V .To each configuration { r j } i = j,...,N one can associate aVoronoi tessellation, a space-filling partitioning of V into N cells assigning one cell to each particle. Thecell of particle j is defined as the region of space beingcloser to j than to any other particle in the system.The cell has a volume v j and a centroid at position g j . Since g j does in general not coincide with r j , wecan introduce the “geometric polarization” τ j of a cellby τ j = v j ( g j − r j ). Analysis of a supercooled liquidof short polymer chains revealed that the geometricpolarization τ j ( t ) at time t is correlated to the totalinteraction force F j ( t ) on particle j and τ j obeys the conservation law P Nj =1 τ j = , analogous to F j [5]. The equilibrium properties of the geometricpolarization are thus reminiscent of those of a force.This observation motivated us to introduce a newmodel for a liquid—the Voronoi liquid—where the forceon particle j is taken proportional to τ j [9]: F j = γ τ j = γ Z v j d r r . (1)Here the constant γ is a parameter of the model and r denotes the vector from the particle position r j tothe boundary of its Voronoi cell. The force F j can bewritten as F j = − ∇ j E p with [9] E p ( r , · · · , r N ) = N X j =1 " γ Z v j d r r . (2)This defines the potential energy of the monodisperseVoronoi liquid as the sum of the interaction energiesof all particles. The interaction energy of a particle ispositive and determined by its nearest-neighbor shell:it is local, many-body, and soft in the sense that theenergy cost for particle overlap is finite.Thermodynamic, structural and dynamic proper-ties [9, 18] have been studied for the monodispersesystem. Upon cooling the liquid becomes metastableand eventually crystallizes in a bcc structure [17]. Forthe study of glasses the tendency of structural order-ing has to be suppressed. This can be achieved byusing systems with multiple components of differentsizes (and interaction energies) [20, 21]. Therefore, weintroduce size dispersity into our model by choosing theVoronoi-Laguerre generalization of the Voronoi tessala-tion [19]. The Voronoi-Laguerre tessellation assigns a“natural radius” R j ( >
0) to every particle j , whichenters the construction of its cell, and has the advan-tage of preserving the defining features of the Voronoiliquid. For the Voronoi-Laguerre tessellation we stillhave P j τ j = and F j = γ τ j = − ∇ j E p with thefollowing generalization of the potential energy [19]: E p ( r , · · · , r N ) = N X j =1 " γ Z v j d r (cid:0) r − R j + R (cid:1) . (3)Here R = P Nj =1 R j /N is mean-square natural radiusaveraged over the polydispersity. Equation (3) reducesto (2) in the monodisperse case.Since E p is defined in terms of v j , the relevantlength scale of the Voronoi liquid is given by v / ,where v = P Nj =1 v j /N = V /N is the average volumeper particle, and the temperature scale by γv / /k B with k B being the Boltzmann constant. As in previouswork [9, 18, 19] we choose the density v − = 1 and take γ = 1000 so that the temperatures of the liquid phaseare in the range T ∼ k B = 1).Here we examine the simplest representative ofa polydisperse system, a binary mixture of N A large lassy dynamics of a binary Voronoi fluid R A and N B = N − N A small particlesof radius R B < R A . The mixture is characterizedby the number concentration of small particles x B = N B /N and the size ratio R B /R A . The choice ofthese parameters is motivated by theoretical work onbinary hard-sphere mixtures [22], suggesting that thepropensity to form a glass is enhanced, relative tothe monodisperse system, for size ratios ∼ . x B ∼ . R B /R A = 0 .
83 and x B = x A = 0 . ξ = R − R (recall that the length scale v / = 1).As pointed out in [19], ξ needs to be smaller than 1 toavoid unphysical situations where small particles aresituated outside their Voronoi cell. On the other hand, ξ has to be large enough to suppress crystallization.From continuous cooling runs at finite rates, it wasfound that the binary mixture forms glasses for 0 . . ξ . .
36 [17]. Here, to probe the glassy regime,we choose ξ = (0 . ≃ . R A = 0 . R B = 0 . ξ is the relevant parameter. The physicalproperties of the mixture are not changed when varying R A and R B but keeping ξ the same. Our choices forthe natural radii, however, turn out to be physicallymeaningful. The partial pair-distribution functions ofthe A particles, g AA ( r ), and the B particles, g BB ( r ),show a first maximum at r AA = 1 . ≈ R A and r BB = 0 . ≈ R B [19]. Therefore, R A and R B canbe thought of as the radii of soft repulsive particles.Moreover, the cross pair-distribution function of A andB, g AB ( r ), peaks at r AB = 1 . ≈ ( r AA + r BB ) / ≈ ( R A + R B ), suggesting that the studied binary Voronoimixture is intrinsically additive [19]. We performed molecular-dynamics (MD) simulationswith a modified version of the LAMMPS code [23],enabling the computation of the geometric polarization τ j , and therefore of F j , by using the Voro++ library[24]. The system contained N = 1000 particles ofmass m = 1 in a cubic box of linear dimension L withperiodic boundary conditions. Since v = V /N = 1,this implies L = 10. Thus, the smallest accessiblewave vector ( q ) has the modulus q min = 2 π/L = 0 . γv / , the mass m and the lengthscale v / , the characteristic time scale of the Voronoiliquid is τ voro = p m/γv = p / ≈ .
03 with γ = 1000, m = 1 and v = 1. The time step ( δt ) of theMD simulation has to be smaller than τ voro . We used δt = 0 . ≈ . τ voro when integrating the classicalequations of motion by the velocity-Verlet algorithm.In the following, all times are measured in units of τ voro . The simulations were carried out in the canonicalensemble with the Nos´e-Hoover thermostat (using adamping parameter of T damp = 0 . . ≤ T ≤
2. This intervalranges from the regime of the normal to the moderatelysupercooled liquid (above the critical temperature ofMCT T c = 0 . T , the temperaturewas instantaneously decreased by a small step to T 84, corresponding to about 10 times the α relaxation time at that temperature. 3. Static structure factors Previous work studied the thermodynamics, the stresstensor, and structural properties of the binary Voronoiliquid [19]. Due to their importance for mode-couplingtheory we here revisit the discussion of the staticstructure factors. The collective static structure factor S ( q ) = 1 N h ρ ( q ) ρ ( − q ) i (4)is defined in terms of the coherent density fluctuationsfor wave vector q , ρ ( q ) = N X j =1 exp (i q · r j ) (for q = 0) , (5)where h . . . i denotes the canonical average and r j is theposition of particle j . For a spatially homogeneous andisotropic system, the structure factor depends only onthe modulus of the wave vector, q = | q | . Figure 1(a)presents S ( q ) for four temperatures in the investigatedinterval 0 . ≤ T ≤ 2. We see that the collectivestructure of the Voronoi mixture is typical of a densedisordered system. In the limit q → S ( q ) is smallbecause the fluctuations of the particle number relativeto the average h N i (= N ) are weak in a dense system, S ( q → 0) = h N i − h N i h N i ≪ . (6)With increasing q , S ( q ) increases toward a maximumthat occurs around q ∗ = 6 . 85. The correspondinglength scale 2 π/q ∗ ∼ S ( q ∗ )comes from the amorphous packing in neighbor shellsaround a particle. Upon cooling the packing becomes lassy dynamics of a binary Voronoi fluid q S ( q ) T=2.00T=1.00T=0.85T=0.83 q* (a) q S t r u c t u r e f a c t o r s q S(q) S AA (q)+S BB (q) 2 S AB (q) S mono (q) S(q) (b) k B T ρχ T q S t r u c t u r e f a c t o r s S nn (q)S cc (q)S nc (q) A x B S cc (0) = 0.1895 < x A x B S nc (0) = - δ S cc (0)S nn (0) = k B T ρχ T + δ S cc (0) (c) Figure 1. (a) Collective static structure factor S ( q ) as afunction of the modulus of the wave vector q for T = 0 . · · · · · · ), 0.85 (——), 1 (– – –) and 2 (— · —). S ( q ) has amaximum near q ∗ = 6 . 85 which is indicated by an arrow. (b) S ( q ) versus q at T = 0 . 85 and its decomposition into partialstructure factors according to (9): S AA ( q ) + S BB ( q ) (– – –)and 2 S AB ( q ) (— · —). The horizontal dotted line indicatesthe value of k B T ρχ T (= 0 . 006 [19]) with ρ = 1 /v the particledensity and χ T the isothermal compressibility. Inset: Zoom forsmall q focusing on S ( q ) (——). The dashed line presents thecollective structure factor S mono ( q ) of the monodisperse Voronoiliquid at T = 1 . 05 which has a comparable compressibility asthe binary mixture. k B T ρχ T of the mixture from the mainfigure is shown as a horizontal dotted line. (c) Bhatia–Thorntonstructure factors versus q at T = 0 . S nn ( q ) = S ( q ) ( × ), S cc ( q )( ◦ ) and S nc ( q ) ( ♦ ). The horizontal dashed lines represent thelarge- q limits: S nn ( q → ∞ ) = 1, S cc ( q → ∞ ) = x A x B and S nc ( q → ∞ ) = 0. The horizontal full lines indicate the limits for q → tighter, which is reflected by the increase of the heightand the decrease of the width of S ( q ) near q ∗ .Further insight can be obtained from the partialstatic structure factors S αβ ( q ) = 1 N h ρ α ( q ) ρ β ( − q ) i ( α, β = A , B) (7)defined by the partial density fluctuations ρ α ( q ) = N α X j α =1 exp (i q · r j α ) , (8)where r j α is the position of particle j α of species α . As ρ A ( q ) + ρ B ( q ) = ρ ( q ) ≡ ρ n ( q ), the collective structurefactor can be expressed as S ( q ) = S nn ( q ) = S AA ( q ) + S BB ( q ) + 2 S AB ( q ) . (9)While S αβ ( q ) characterize spatial correlations betweenlike and unlike particles, S ( q ) describes number-number (nn) correlations (whence the notation S = S nn ). Equation (9) is not the only physically significantlinear combination of the partial structure factors.Since composition (or concentration) fluctuations ρ c ( q )are defined by ρ c ( q ) = ρ A ( q ) − x A ρ ( q ) = x B ρ A ( q ) − x A ρ B ( q ), the structure factor S cc ( q ) = 1 N h ρ c ( q ) ρ c ( − q ) i = x S AA ( q ) + x S BB ( q ) − x A x B S AB ( q ) (10)represents composition-composition (cc) correlations,and the structure factor between ρ n and ρ c , S nc ( q ) = 1 N h ρ n ( q ) ρ c ( − q ) i = x B S AA ( q ) − x A S BB ( q ) + ( x B − x A ) S AB ( q ) , (11)describes number-composition (nc) correlations. Thestructure factors S nn ( q ), S cc ( q ) and S nc ( q ) are oftenreferred to as Bhatia–Thornton structure factors [25].They have been studied extensively for metallic alloys[26, 27] or colloidal suspensions [28].Figure 1(b) compares S ( q ) with S AA ( q ) + S BB ( q )and 2 S AB ( q ) at T = 0 . 85. In the limit q → ∞ , thesystem behaves like an ideal mixture with vanishingcorrelations. This implies S AB ( q → ∞ ) = 0 as wellas S AA ( q → ∞ ) = x A and S BB ( q → ∞ ) = x B . Forlarge q , say q & 20, the behavior of S ( q ) is thereforedominated by correlations between like particles. Thesum S AA ( q ) + S BB ( q ) is positive for all q , whereas S AB ( q ) oscillates around 0 and remains negative for q < q ∗ , displaying a minimum at q ≈ 4. These negativevalues indicate that long-range AB correlations aresuppressed in the Voronoi mixture, a feature also foundin other binary systems [27, 29, 30]. The minimum of S AB ( q ) at q ≈ S AA ( q ) + S BB ( q ), leading to a dip in S ( q ) at q ≈ S ( q ) increases again toward a plateau as q → 0. Such a dip is not observed for the monodisperseVoronoi liquid. Here S mono ( q ) continuously decreases lassy dynamics of a binary Voronoi fluid S mono ( q → 0) = k B T ρχ T with χ T being the isothermal compressibility(cf inset in figure 1(b)). For the mixture, however, wesee that S ( q → 0) adopts a value larger than k B T ρχ T .For binary mixtures this deviation between S ( q → 0) and k B T ρχ T is expected from the work of Bhatiaand Thornton [25] and also from the Kirkwood–Bufftheory for multicomponent solutions [31]. For q → S ( q → 0) = k B T ρχ T + δ S cc ( q → , (12) S cc ( q → 0) = N k B T ( ∂ G/∂x ) p,T,N , (13) S nc ( q → 0) = − δ S cc ( q → , (14)where G is the Gibbs free energy, p the pressure and δ = ρ ( v A − v B ) (15)is a dilatation factor given by the partial molar volumes v A = ( ∂V /∂N A ) p,T,N B and v B = ( ∂V /∂N B ) p,T,N A .Equation (12) shows that in a mixture fluctuations ofthe total particle number, S ( q → ∂ G/∂x ) p,T,N > 0, we have in general S ( q → > k B T ρχ T , as seen in the inset of figure 1(b).This implies that the second term in the right-hand-side of (12) does not vanish, in particular δ = 0 or v A = v B . The molar volumes of the two species canbe calculated via the Kirkwood–Buff theory from thepartial structure factors in the limit q → v A = v x A S BB (0) − x B S AB (0) x S BB (0) + x S AA (0) − x A x B S AB (0) , (16) v B = v x B S AA (0) − x A S AB (0) x S BB (0) + x S AA (0) − x A x B S AB (0) , (17)where S αβ (0) is an abbreviation for S αβ ( q → S cc ( q ) and S nc ( q ), as obtained from (10) and (11),together with S ( q ) at T = 0 . 85. We find that S cc ( q )is positive for all q . For large q , S cc ( q ) oscillatesaround x A x B (= 0 . S cc ( q → 0) = 0 . q limit. The ratio Φ( x A , T ) = x A (1 − x A ) /S cc ( q → 0) enters the definition of theinterdiffusion coefficient of the mixture (cf (61)) [26].In systems that favor mixing, as the binary Voronoimixture [19], one has Φ > x A = 1 / , T = 0 . ≃ . v A = 1 . 266 and v B = 0 . δ = 0 . 532 at T = 0 . 85. If we also take k B T ρχ T = 0 . S cc ( q → 0) from figure 1(c), the values of S ( q → 0) and S nc ( q → 0) can be computed.These results are shown as horizontal full lines infigure 1(c). As can be seen, we find excellent agreementbetween simulation and theoretical expectation. 4. Mode-coupling theory The idealized mode-coupling theory (MCT) and itsapplication to simple and molecular glass formers aredescribed in detail in a monograph [33], as well asin several review papers [34, 35]. Specifically forbinary mixtures, MCT is also discussed in severalpublications, see e.g. [22, 27, 29, 32, 36, 37, 38]. Belowwe first recapitulate the main equations and generalpredictions of binary MCT, and subsequently apply thetheory to our Voronoi mixture. MCT assumes that the slow dynamics of glass-formingliquids results from the relaxation of collective densityfluctuations. For binary mixtures the central dynamiccorrelation functions are therefore the partial dynamicstructure factors S αβ ( q, t ) = 1 N h ρ α ( q , t ) ρ β ( − q ) i ( α, β = A , B) (18)where ρ α ( q , t ) = N α X j α =1 exp [i q · r j α ( t )] (19)and r j α ( t ) is the position of particle j α of species α attime t . Let us combine these functions in a 2 × S ( q, t ) with ( S ( q, t )) αβ = S αβ ( q, t ). By means of theZwanzig–Mori projection operator formalism, an exactequation of motion for S ( q, t ) is derived [36, 38] ∂ t S ( q, t ) + J ( q ) S − ( q ) S ( q, t )+ J ( q ) Z t d t ′ M ( q, t − t ′ ) ∂ t ′ S ( q, t ′ ) = . (20)The matrix J αβ ( q ) = q ( k B T /m α ) δ αβ is given in termsof the square of the thermal velocities, k B T /m α , where m α is the mass of a particle of species α ; this matrixdescribes inertial effects. Outside the initial inertialregime the dynamics is determined by the memorykernels M αβ ( q, t ). These kernels are fluctuating force-correlation functions, reflecting many-body interactioneffects. MCT writes M ( q, t ) as a sum of two terms: M ( q, t ) = M reg ( q, t ) + M MCT ( q, t ). The “regular”term M reg ( q, t ) is supposed to describe the normalliquid-state dynamics; it decays on short time scalesand is not responsible for slow glassy dynamics. Wemodel the regular term as a Markovian process withfriction constant ν : M reg ( q, t ) = J − ( q ) νδ ( t ) [39]. Theslow dynamics is encapsulated in the second term.For M MCT ( q, t ) the theory assumes that the dominant lassy dynamics of a binary Voronoi fluid M MCT ( q, t ) = F [ S ( t ) , S ( t )]( q ) , (21)where the components of the mode-coupling functional F are given by F αβ [ S ( t ) , S ( t )]( q )= 12 q ρx A x B Z d k (2 π ) X α ′ β ′ α ′′ β ′′ V αα ′ α ′′ ( q , k ) × V ββ ′ β ′′ ( q , k ) S α ′ β ′ ( k, t ) S α ′′ β ′′ ( | q − k | , t ) . (22)Here V αα ′ α ′′ are the coupling vertices V αα ′ α ′′ ( q , k ) = q · k q c αα ′ ( k ) δ αα ′′ + q · ( q − k ) q c αα ′′ ( | q − k | ) δ αα ′ , (23)which only depend on the equilibrium structure of thesystem via the matrix of direct correlation functions, c αβ ( q ), defined in terms of S ( q ) by the Ornstein–Zernike equation ρc αβ ( q ) = δ αβ x A − (cid:0) S − ( q ) (cid:1) αβ . (24)In writing (23) we assume that static triple correlationscan be treated by the convolution approximation. Thisapproximation has been justified for the Kob–Andersenmixture [40] and we suppose that it also holds for theVoronoi mixture.Equation (18) to (24) establish a link between theequilibrium structure and dynamics of a glass former.This opens the possibility to predict the temperaturedependence of the dynamics based on static inputfrom simulations and to compare these predictionsagainst the simulated relaxation behavior. Here wewill carry out such a comparison for the Voronoimixture. Similar comparisons have been performed fora variety of different models, including binary [36, 38]and polydisperse hard-sphere systems [39], the Kob–Andersen Lennard-Jones mixture [37, 41, 42], metallicglasses [27, 32], strong liquids [40, 43], orthoterphenyl[44, 45], and polymer melts [15, 46, 47, 48]. To describe the single-particle dynamics, MCTconsiders the correlation function of the tagged-particledensity, i.e. the incoherent intermediate scatteringfunction φ s ,α ( q, t ) = 1 N α N α X j =1 h exp { i q · [ r j ( t ) − r j (0)] }i , (25) of species α (= A , B). By the Zwanzig–Mori projectionoperator formalism the following equation of motion isobtained m α q k B T ∂ t φ s ,α ( q, t ) + φ s ,α ( q, t )+ Z t d t ′ M s ,α ( q, t − t ′ ) ∂ t ′ φ s ,α ( q, t ′ ) = 0 . (26)As for the coherent density fluctuations, the memorykernel is approximated by a sum of a regular part,modeled as a damped Markovian process M s ,α reg ( q, t ) = m α ν/ ( q k B T ) δ ( t ) [39], and an MCT contribution. Theexpression for the latter reads [32, 39] M s ,α MCT ( q, t ) = F s ,α [ S ( t ) , φ s ,α ( t )]( q )= ρq Z d k (2 π ) X α ′ β ′ (cid:18) q · k q (cid:19) c αα ′ ( k ) c αβ ′ ( k ) × S α ′ β ′ ( k, t ) φ s ,α ( | q − k | , t ) . (27)The solution of (26) requires not only static input, butalso the collective S ( q, t ) which needs to be determinedfrom (22). Using bipolar coordinates and the rotational symmetryof the system, the three-dimensional integral over k in (22) and (27) is written as a double integral over k = | k | and p = | q − k | . Then, q is discretized byintroducing a finite, equally spaced grid of M points q = q + ˆ q ∆ q with ˆ q = 0 , , . . . , M − 1. This allows usto replace the double integral by Riemann sums Z ∞ d k Z q + k | q − k | d p → (∆ q ) M − X ˆ k =0 min[ M − , ˆ q +ˆ k ] X ˆ p = | ˆ q − ˆ k | . (28)Following commonly made choices [39] we took M =300, ∆ q = 0 . . . . ) and q = 0 . . . . )so that 0 . . q . . t = 10 − , and ∆ t was subsequently doubledfor every 32 new points. The friction constant of theregular kernel was set to ν = 1. MCT makes a number of “universal” predictions. Theyare universal in the sense that they do not depend onthe details of the static input, but are mathematicalconsequences of the form of the MCT equations[33]. Here we summarize those results which will beimportant for the analysis of the MD simulations. lassy dynamics of a binary Voronoi fluid F ( q ) = lim t →∞ S ( q, t ) ,f s ,α ( q ) = lim t →∞ φ s ,α ( q, t ) . (29)By means of the Laplace transform and the final valuetheorem, one can show that F ( q ) and f s ,α ( q ) obey theequations F ( q ) = S ( q ) − (cid:0) S − ( q ) + F [ F , F ]( q ) (cid:1) − , (30) f s ,α ( q )1 − f s ,α ( q ) = F s ,α [ F , f s ,α ]( q ) . (31)These equations are defined by the static structure;neither the inertia matrix J nor the regular kernel M reg enter. Therefore, the solutions are independent of themicroscopic dynamics. Equations (30) and (31) can besolved by an iteration procedure [50].The solutions of (30) display bifurcations. Forstructural glasses usually the A bifurcation is relevant[33]. If T is the control variable, the bifurcationoccurs at a critical temperature T c (depending oncomposition and particle size ratio [22]). For T > T c ,one has F ( q ) = . This behavior corresponds to anergodic liquid where density correlations decay to 0for t → ∞ . For T ≤ T c , the long-time limits aregiven by a (nondegenerate symmetric) positive-definitematrix F ( q ). Since density correlations no longer decayto zero, MCT describes an amorphous solid, i.e. anonergodic ideal glass. Accordingly, the corresponding F ( q ) are called “nonergodicity parameters”. The glasstransition point T c can be identified with the highesttemperature at which the system is a glass, i.e. at which F jumps from to some finite F c . In the generic case,the tagged-particle dynamics is strongly coupled to thecollective dynamics and undergoes a glass transition.This implies that the solution of (31) also jumps fromzero to a finite value f sc ,α ( q ) at T c . Along with thefinite value F c of the nonergodicity parameter, thecorresponding stability matrix C c of (30), defined by C c [ H ( q )] = (cid:2) S c ( q ) − F c ( q ) (cid:3)(cid:2) F c [ F c , H ]( q )+ F c [ H , F c ]( q ) (cid:3)(cid:2) S c ( q ) − F c ( q ) (cid:3) , (32)has a unique right eigenvector H ( q ) with eigenvalue E = 1. Here the superscript “c” means that all staticinput is evaluated at T c . The normalization factors of H ( q ) are determined by the convention b H ( q ) : H ( q ) = 1 , b H ( q ) : (cid:8) H ( q )[ S c ( q ) − F c ( q )] − H ( q ) (cid:9) = 1 , (33)where b H ( q ) is the left eigenvector of C c with E = 1and the double-dot operator includes integration over q . Close to T c the solutions of (20) and (26) showthat S ( q, t ) and φ s ,α ( q, t ) stay close to a plateau given by F c ( q ) and f sc ,α ( q ) for an intermediate time interval.This time interval is called “ β relaxation regime” inMCT. Whereas the β relaxation exists both above andbelow T c (below T c , F c and f sc ,α are replaced by the T dependent long-time limits (29)), a decay of S ( q, t )from the plateau to zero—i.e. the α relaxation—canonly occur in the liquid phase for T > T c . The β andthe α process are characterized by two time scales: the β relaxation time t σ , t σ = t | σ | / a (for T → T ± c ) , (34)and the α relaxation time t ′ σ , t ′ σ = t ( − σ ) γ , γ = 12 a + 12 b (for T → T +c ) . (35)Here t represents a system-specific microscopic timescale and σ is the “separation parameter” quantifyingthe distance to the critical point where the bifurcationoccurs. Close to T c the separation parameter can beexpressed as σ = Cε, ε = T c − TT c (36)with C being a constant. MCT refers to a as “criticalexponent” and to b as “von Schweidler exponent”.They are connected to one another by the “exponentparameter” λ , λ = Γ (1 − a ) Γ (1 − a ) = Γ (1 + b ) Γ (1 + 2 b ) (1 / ≤ λ < , (37)where Γ is the Gamma function. The parameter λ is a static quantity that can be calculated from theequilibrium structure of the glass former at T c by λ = b H ( q ) : (cid:8)(cid:2) S c ( q ) − F c ( q ) (cid:3) F c [ H , H ]( q ) (cid:2) S c ( q ) − F c ( q ) (cid:3)(cid:9) . (38)Since 1 / ≤ λ < A bifurcation [33], (37)gives 0 < a < . < b ≤ 1, and so γ > . T c , the ratio t ′ σ /t σ increases. The smaller T − T c , the more separatedare the β and α relaxation regimes. MCT thereforepredicts a two-step relaxation. The intermediate timeinterval of the β regime is defined by t ≪ t ≤ t ′ σ .This interval comprises t ∼ t σ where S ( q, t ) ∼ F c or φ s ,α ( q, t ) ∼ f sc ,α ( q ). The α regime begins for t > t σ and leads to S ( q, t ) → or φ s ,α ( q, t ) → t ≫ t ′ σ .Both regimes overlap for t σ ≤ t ≤ t ′ σ . The latter timeinterval is called late β or early α process in MCT.For both the α and β process, MCT makesdetailed predictions [33, 51], many of which have beentested in fits to experimental and simulation data (forreviews of these tests see e.g. [33, 35, 48, 52, 53, 54,55, 56]). In the following, we will also perform suchfits for the binary Voronoi mixture. This analysis will lassy dynamics of a binary Voronoi fluid φ ( q, t ) = S AA ( q, t ) + S BB ( q, t ) + 2 S AB ( q, t ) S ( q ) , (39)the incoherent scattering functions φ s ,α ( q, t ) and themean-square displacements (MSDs), g ,α ( t ) = 1 N α N α X j =1 D [ r j ( t ) − r j (0)] E , (40)of species α = A , B. The MSD is related to φ s ,α ( q, t )by g ,α ( t ) = lim q → − φ s ,α ( q, t )] /q . Therefore, wesummarize below the MCT predictions pertinent forthis analysis. Predictions for the β regime In the β regime MCTpredicts a “factorization theorem” according to whichall density correlators (and all quantities coupling tothem) can be expressed as a sum of the nonergodicityparameter and a correction term that exhibits afactorization into a wavevector-dependent and a time-dependent part [33, 51]: φ ( q, t ) = f c ( q ) + h ( q ) G ( t ) , (41) φ s ,α ( q, t ) = f sc ,α ( q ) + h s ,α ( q ) G ( t ) . (42)The nonergodicity parameters, f c ( q ) = lim t →∞ φ ( q, t )and f sc ,α ( q ) = lim t →∞ φ s ,α ( q, t ), and the “criticalamplitudes”, h ( q ) and h s ,α ( q ), are evaluated at T c and are thus independent of T . The temperaturedependence resides in the “ β correlator” G ( t ) which,for T → T +c , is given by G ( t ) = p | σ | g ( t/t σ ) t ≫ t σ −→ − B ( λ ) (cid:18) tt ′ σ (cid:19) b . (43)Here B ( λ ) is a T -independent constant and G ( t ) ∼− ( t/t ′ σ ) b is the so-called von Schweidler law whichholds for t σ ≪ t ≪ t ′ σ . Both the factorization theoremand the von Schweidler law are MCT results in leadingorder of | σ | . Second order corrections to (41) and (42)are also known [57, 58]: φ ( q, t ) = f c ( q ) − h ( q ) B ( λ ) (cid:18) tt ′ σ (cid:19) b + h ( q ) B ( λ ) B ( q ) (cid:18) tt ′ σ (cid:19) b , (44) φ s ,α ( q, t ) = f sc ,α ( q ) − h s ,α ( q ) B ( λ ) (cid:18) tt ′ σ (cid:19) b + h s ,α ( q ) B ( λ ) B s ,α ( q ) (cid:18) tt ′ σ (cid:19) b . (45)The q dependence of the correction amplitudes B ( q )and B s ,α ( q ) implies a violation of the factorizationtheorem. Both amplitudes are again evaluated at T c ; inthe β regime the T dependence therefore solely resultsfrom the time scale t ′ σ . Predictions for the α regime In the α regime MCTpredicts that the density correlators are described by T -independent master curves (for T → T +c ): φ ( q, t ) = e φ ( q, t/t ′ σ ) , φ s ,α ( q, t ) = e φ s ,α ( q, t/t ′ σ ) , (46)which have the following limits for t → φ ( q, t → 0) = f c ( q ) , φ s ,α ( q, t → 0) = f sc ,α ( q ) . (47)Equation (46) implies a time-temperature superposi-tion principle (TTSP): For fixed q , φ ( q, t ) and φ s ,α ( q, t )collapse for different T onto master curves when rescal-ing t by some relaxation time that is proportional to t ′ σ . For instance, we can choose φ ( q, t ) at the peak po-sition q ∗ of S ( q ) to define the relaxation time τ q ∗ bythe criterion φ ( q ∗ , τ q ∗ ) = const. Then, we have τ q ∗ = C q ∗ t ′ σ , (48)where the T -independent prefactor C q ∗ is determinedby the constant used in the definition φ ( q ∗ , τ q ∗ ) =const.For t ≪ t ′ σ , (46) recovers the von Schweidler law.This justifies the statement made above that the late β and early α process overlap for t σ ≪ t ≪ t ′ σ .Moreover, model calculations within MCT reveal thatthe α master curves are stretched. As for experimentalor simulation data, this stretched relaxation can befitted well by a Kohlrausch–Williams–Watts (KWW)function. For φ ( q, t ) the KWW function reads φ ( q, t ) ≃ A ( q ) exp " − (cid:18) tτ K ( q ) (cid:19) β K ( q ) ( t ≥ t σ ) , (49)where A ( q ) is an amplitude, τ K ( q ) the relaxation timeand β K ( q ) ≤ α process, except inthe special limit of large q . In this limit, it was proved[59] that there is a time interval t/t ′ σ ≪ t K q /t ′ σ ≤ α process obeyslim q →∞ φ ( q, t ) = f c ( q ) exp " − Γ( q ) (cid:18) tt ′ σ (cid:19) b , (50)with Γ( q ) ∝ q . This implieslim q →∞ A ( q ) = f c ( q ) , lim q →∞ β K ( q ) = b, lim q →∞ τ K ( q ) ∝ t ′ σ q /b . (51)Equations analogous to (50) and (51) also hold for theincoherent scattering functions φ s ,α ( q, t ). 5. Results Our MCT analysis of the Voronoi mixture starts witha test of the TTSP. To rescale the time axis we follow lassy dynamics of a binary Voronoi fluid -5 -4 -3 -2 -1 t/ τ q* φ ( q * ,t ) T=0.88T=0.86T=0.85T=0.84T=0.83 τ q* Arrhenius τ q* Figure 2. Test of the TTSP for 0 . ≤ T ≤ . φ ( q ∗ , t )as a function of t/τ q ∗ where q ∗ = 6 . 85 is the peak positionof the first maximum of S ( q ). τ q ∗ is defined by the condition φ ( q ∗ , t = τ q ∗ ) = 0 . T = 0 . 88, the temperature abovewhich the TTSP is violated. Inset: Arrhenius plot of τ q ∗ versus1 /T . The circles represent the simulation data. The solid lineshows a fit to τ q ∗ ( T ) = τ ∞ q ∗ exp( E A /k B T ) with τ ∞ q ∗ = 0 . E A /k B = 5 . (48) and define the α relaxation time as the time when φ ( q ∗ , t ) has decayed to 10% of its initial value, i.e. φ ( q ∗ , t = τ q ∗ ) = 0 . 1. The threshold of 0.1 is arbitrary,but convenient: The choice ensures that the densitycorrelator is small enough to be well in the α regime,but still sufficiently above the noise level so that thestatistical accuracy of the data remains satisfactory.Figure 2 shows φ ( q ∗ , t ) as a function of t/τ q ∗ for 0 . ≤ T ≤ . 88. This interval corresponds to the regimeof the supercooled liquid where a super-Arrheniusincrease of τ q ∗ with decreasing T is observed (cf inset offigure 2). For these temperatures we find that φ ( q ∗ , t )decays in two steps, developing an intermediate timeinterval where φ ( q ∗ , t ) plateaus. This time intervalextends upon cooling, and the second relaxation stepaway from the plateau toward zero obeys the TTSPfor T . . 88. These observations are in qualitativeagreement with MCT, suggesting to focus on T . . t and t ( t > t )in the β regime and calculate the ratio R ( q, t ) = φ ( q, t ) − φ ( q, t ) φ ( q, t ) − φ ( q, t ) = G ( t ) − G ( t ) G ( t ) − G ( t )= φ s ,α ( q, t ) − φ s ,α ( q, t ) φ s ,α ( q, t ) − φ s ,α ( q, t ) = R s ,α ( q, t ) , (52)where α = A , B. This equation shows that R ( q, t ) and t -3-2-10123 R q ( t ) q=3.00q=6.85q=9.55q=11.75q=15.35 (a) t t t -3-2-10123 R q R s , A ( t ) R q R s , B ( t ) R(q*,t)A: q=5.5A: q=12.5B: q=5.5B: q=12.5 t -101 R(q*,t) b=0.5652b=0.6172 (b) t t Figure 3. Test of the factorization theorem at T = 0 . t = 6 and t = 16. Panel(a): R ( q, t ) for 3 ≤ q ≤ . 35. Panel (b): R s ,α ( q, t ) for q = 5 . . R ( q, t = 6) = R s ,α ( q, t = 6) = 1 and R ( q, t = 16) = R s ,α ( t = 16) = 0. Thecoordinate positions of t and t are indicated by a plus sign.The red filled circles in panel (b) reproduce R ( q ∗ , t ) from panel(a) to illustrate that R ( q, t ) and R s ,α ( q, t ) collapse onto the samemaster curve. In both panels, the dashed black line presents (53)with b = 0 . − . ≤ R ( q ∗ , t ) ≤ . b = 0 . b = 0 . R s ,α ( q, t ) are independent of q and superimpose on thesame curve in the time window where the factorizationtheorem holds. Using furthermore (43), R ( q, t ) and R s ,α ( q, t ) are given by R ( q, t ) = R s ,α ( q, t ) = t b − t b t b − t b . (53)Equations (52) and (53) are predicted to hold closeto T c . In the following, we therefore focus on a lowtemperature, T = 0 . T = 0 . 84 with the choice t = 6 and t = 16. Wesee that there is a time interval comprising t and t where R ( q, t ) and R s ,α ( q, t ) are indeed independent of q (cf top and bottom panel of figure 3) and collapse ontothe same master curve (bottom panel). The mastercurve tends to persist for t < t in the case of R s ,α ( q, t )and also for R ( q, t ) if q ≥ . 85, while for q ≤ lassy dynamics of a binary Voronoi fluid q -dependent way. This finding is expected fromMCT which predicts an ordering rule [57, 58]: Since thesecond-order corrections to the factorization theoremhave the same q dependent amplitudes both for theearly-time and long-time corrections, correlators thatlie, for example, above the factorization theorem forshort times must also lie above it for long times.Therefore, if we number the correlators in the order inwhich they enter the collapse regime, this numberingis preserved when the correlators leave the regime [57].This prediction has been observed in many simulations[38, 39, 61, 65, 66, 67, 46]. Figure 3 suggests that italso holds for our Voronoi mixture.Finally, the black dashed lines in figure 3(a) andfigure 3(b) indicate that the master curve is welldescribed by (53) with b = 0 . b = 0 . b values in section 5.3).Therefore, in the present case, we see that (53) doesnot allow to determine b precisely: Equation (53) isan asymptotic result for T close to T c . Apparently, T = 0 . 84 is still too far above T c so that the timeinterval over which the factorization holds—about adecade in figure 3—is too narrow. From this analysiswe conclude that, albeit a precise determination of b via (53) may be difficult in practice, it certainly allowsto obtain bounds for b that can serve as valuable inputto guide the fits via (44). We turn to such fits in thenext section. β regime We examine the dynamics in the supercooled regimeby fitting the asymptotic MCT predictions, (44) and(45), to our simulation data. To this end, we write (44)in the following form: φ ( q, t ) = f c ( q ) − h fit ( q ) (cid:18) tt ′ σ (cid:19) b + h fit ( q ) B fit ( q ) (cid:18) tt ′ σ (cid:19) b . (54)The fit constants h fit ( q ) and B fit ( q ) are related to theamplitudes h ( q ) and B ( q ) of (44) by h fit ( q ) = B ( λ ) h ( q ) , B fit ( q ) = B ( λ ) B ( q ) . (55)The same equations are also valid for φ s ,α ( q, t ) aftersubstituting f c ( q ) → f sc ( q ), h ( q ) → h s ,α ( q ) and B ( q ) → B s ,α ( q ). Five fit parameters are involved in (54). Fourof them are independent of T , namely f c ( q ), h fit ( q ), B fit ( q ) and b . One parameter, the α time t ′ σ , dependson T . To carry out the fits it is judicious to work atlow temperature. Guided by the tests of the TTSPand of the factorization theorem, we begin the analysiswith the coherent scattering function at q = q ∗ and T = 0 . 84 because inspection of figure 2 suggests thatthe plateau, i.e. f c ( q ∗ ), is large and so the late β process is pronounced. This allows us to determine b . Fixing b and performing the fits for different T gives t ′ σ ( T ). Keeping then b and t ′ σ ( T ) constant, thewavevector dependence of f c ( q ), h fit ( q ) and B fit ( q ) canfinally be determined. In practice, we utilize again T = 0 . 84 for the latter fits.It is known that information from the α relaxationis crucial to guide the fit in the β regime [69]. The five-parameter fit is thus subjected to two constraints:(i) The nonergodicity parameter f c ( q ) is the initialvalue of the α master curve [cf (47)], implying that e φ ( q, t/t ′ σ ) < f c ( q ) for t/t ′ σ > 0. This imposes a lowerbound on f c ( q ). The fit result for f c ( q ) cannot besmaller than the value of φ ( q, t ) at the shortest timewhere the TTSP still holds. To verify this constraintfigure 2 serves as a guideline.(ii) Equation (44) is invariant under the rescaling h ( q ) → ℓ b h ( q ), B ( q ) → ℓ b B ( q ) and t ′ σ → ℓt ′ σ where ℓ isa constant scale factor [36]. Thus, the same fit resultcan be obtained for a small (small ℓ ) or large (large ℓ ) α time t ′ σ , provided the amplitudes are rescaledaccordingly. To guide the fit, we make use of the fullymicroscopic MCT calculations based on static input.Early work on binary soft-sphere mixtures showed that0 . < h ( q ) < . q ∗ / . q . q ∗ [70]. When fittingthe data one can constrain h ( q ) to lie within thesebounds. In the present case, we take advantage of theMCT calculations using the static structure factors ofour simulations. These calculations provide h ( q ) andwe adjust the constant ℓ such that the fit result matchesthe theoretical h ( q ).A final technical aspect is to choose the timeinterval [ t min , t max ] where the fit is carried out becausethe latter can have a significant influence on the fit[71, 72, 73, 69]. Certainly, t min should be largerthan the time associated with the initial relaxation( t min & 1, cf figure 4), whereas t max ( > t min ) maynot be taken too large to ensure that the second ordercorrection in (54) remains small in comparison to thevon Schweidler law. Preliminary tests at T = 0 . t min = 10 , t max = 500] satisfiesthese requirements (e.g. we find that | B fit ( q ∗ ) | ( t/t ′ σ ) b ∼ . β analysis in the following. lassy dynamics of a binary Voronoi fluid T c λ a b γ B Table 1. MCT parameters obtained from fits of the MD datato the asymptotic MCT predictions. T c is the average valueof the critical temperature determined in figure 5. The fit of φ ( q = q ∗ , t ) to (54) gives the von Schweidler exponent b fromwhich λ , a and γ are calculated by (35) and (37). The constant B = B ( λ ) appearing in (55) is obtained by interpolation of thedata in Table 3 of [74]. -2 -1 t φ ( q ,t ) T=0.88T=0.86T=0.85T=0.84T=0.83 t’ σ (0.88) t’ σ (0.83)f c (q) Figure 4. Time evolution of the coherent scattering function φ ( q, t ) for q = q ∗ = 6 . 85 and 0 . ≤ T ≤ . 88 (full lines). Theblack dashed lines correspond to the result of the fit to (54)carried out for the interval 10 ≤ t ≤ f c ( q ) = 0 . α relaxation time( t ′ σ ) at T = 0 . 83 and T = 0 . 88, respectively. T τ - / γ τ q * ABcoh T ( t ’ σ ) - / γ Figure 5. Rectification plot of the α relaxation times for theA particles 1 / ( τ s , A q ∗ ) /γ (crosses), the B particles 1 / ( τ s , B q ∗ ) /γ (triangles) and the coherent dynamics 1 / ( τ q ∗ ) /γ (circles). Therelaxation times are defined as the time when the incoherent orcoherent intermediate scattering functions take a value of 0.1. γ is given in table 1. The black full lines are linear extrapolationsto zero giving T c = 0 . 799 (A particles), T c = 0 . 799 (B particles)and T c = 0 . 801 (coherent dynamcis). Inset: Rectification plot ofthe MCT α relaxation time t ′ σ (squares). The black full line isan extrapolation to zero giving T c = 0 . Figure 4 depicts the simulation results for φ ( q ∗ , t ) inthe temperature interval 0 . ≤ T ≤ . 88 (full lines).The dashed lines present the fits to (54). The fits yielda good description of the MD data, over about twodecades in time at T = 0 . 88 and extending to aboutthree decades at T = 0 . 83. The fits extend to fairlyshort times; they begin to describe the MD data afterthe first relaxation step for t & 4. The shape of thisfirst step depends on the microscopic dynamics of thesimulation method—e.g. Newtonian or Langevin-based[61, 75, 76]. For Newtonian MD simulations, as in ourcase, the first step masks the early β relaxation towardthe nonergodicity parameter [39, 55]. Due to thisreason, we based our analysis on the MCT predictionsfor the late β process, as many other works have doneas well [36, 39, 45, 64, 66, 67, 77].From the fits to (54) we find b = 0 . λ = 0 . λ are found forhard spheres [57, 58, 39] and various binary mixtures[27, 29, 36, 37] (see also [66] for an overview of λ and further MCT parameters for simple and polymericliquids). In this respect, our binary Voronoi mixture iscomparable to other glass-forming systems.The fit also provides t ′ σ ( T ). Following (35) a plotof 1 / ( t ′ σ ) /γ against T , with γ calculated from a and b via (35), should give a straight line that extrapolatesto 0 at T c . By virtue of (48), the same behavioris expected for the α relaxation times defined by φ ( q ∗ , τ q ∗ ) = 0 . φ s ,α ( q ∗ , τ s ,αq ∗ ) = 0 . 1. Figure 5 teststhese expectations. For all relaxation times we findstraight lines extrapolating to almost the same valueof T c . From these results we calculate the average T c =0 . 798 given in table 1. This T c is in excellent agreementwith the independent estimate T c = 0 . 804 determinedfrom the vanishing of the negative directions associatedwith saddle points of the potential energy landscape[17]. As described in section 4.4, the critical tempera-ture and the MCT parameters can also be predictedby MCT calculations in a fit-parameter-free mannerbased on the static input of the system. More specifi-cally, once the critical temperature is determined, thelong time limit of density correlation functions F c ( q )and the matrix form of the critical amplitude H ( q ) canbe obtained by solving (30), (32), and (33). The non-ergodicity parameters and the critical amplitude whichwill be used to compare to the simulation results arerelated to the components of F c ( q ) and H ( q ) via f c ( q ) = F cAA ( q ) + F cBB ( q ) + 2 F cAB ( q ) S ( q ) ,h ( q ) = H AA ( q ) + H BB ( q ) + 2 H AB ( q ) S ( q ) , (56) lassy dynamics of a binary Voronoi fluid T c E λ a b γ Table 2. Impact of the precision of the MCT critical point location ( T c ) on the value of λ and the exponents a , b and γ , going fromthe most precise (first row, with an eigenvalue E of the stability matrix C close to 1) to the least precise (last row). The parametershave been obtained from binary MCT using the (interpolated) static structure factors as input. with S ( q ) being given by (9). The exponent parameter λ and the exponents a , b can be obtained via (38) and(37). Since all these quantities are determined at T c , itis vital to accurately predict the critical temperature.To determine T c , we used linear interpolations for thepartial static structure factors between T = 0 . 97 and0 . 98. The results are summarized in table 2. Fromthe bottom to the top the precision of T c increases.The best estimate for T c is T c = 0 . 979 245, leading to λ = 0 . ϕ c = 0 . 52 and λ = 0 . 758 [78]. Later, the estimate of the criticalpoint was refined to ϕ c = 0 . 515 912 13(1), leading to λ = 0 . 723 [57]. Reference [57] points out that thishigh accuracy of ϕ c is necessary to reproduce the slowdynamics over many orders of magnitude within MCT.A similar sensitivity of λ on T c is also reported forthe Kob–Andersen binary mixture. Using static inputfrom simulations the first predictions were T c = 0 . λ = 0 . 708 [29], whereas later work suggested amore precise estimate of T c = 0 . 951 5 and along withthat a different value of ( γ = 2 . 46 corresponding to) λ = 0 . 735 [37].When comparing the results of table 1 and table 2two differences can be noted. First, T MDc = 0 . Figure 6. Time dependence of φ ( q, t ) for various q . The fulllines depict the MD data for q = 6 . , . , . 35 at T = 0 . ε = ( T c − T ) /T c = − . × − . For q = 6 . 85 the dottedline reproduces the fit result to (54) from figure 4. The dashedlines present the results of the MCT calculations based on thestatic input for q = 6 . , . , . 29 at T = 0 . ε = − . × − . The MCT results are shifted along the t axisso as to optimize the overlap with the MD data for 5 . t . . × for q = 6 . 89, 1 . × for q = 11 . × for q = 15 . fits. A smaller value of b implies more stretchingof the α relaxation. This is illustrated in figure 6which compares the MD results for φ ( q, t ) at T = 0 . q (full lines) with the MCT calculations(dashed lines). The MCT calculations correspond toa temperature very close to T MCTc and are thereforegood proxies for the α master curve at the wave vectorsconsidered. The MCT curves are shifted along thetime axis so as to optimize the overlap with the MDdata for t ∼ 10, that is in the time window shownin the inset of figure 3 where a distinction between b MCT = 0 . b MD = 0 . q = 6 . 85 (dottedline). Figure 6 also shows that the MD data at longtimes lie above the MCT calculations and are thusmore stretched (for q = 15 . 35 this is not visible onthe scale of the figure). The fit based on the MD datamodels this enhanced stretching by a smaller value ofthe von Schweidler exponent. lassy dynamics of a binary Voronoi fluid b MCT decreases with increasing preci-sion of T MCTc , table 2 indicates that the decrease isfairly weak. It is thus unlikely that further improve-ment of T c might make b MCT converge to b MD . Doesthis mean that MCT cannot account for the enhancedstretching of the α relaxation? Not necessarily, al-beit (presumably) not within the idealized MCT inthe present case. Extensions of the theory need toround off the ideal glass transition and account foractivated processes. Recent efforts in this directioninvolve the inclusion of activated events at the sin-gle particle level [79, 80, 81, 82], the implementationof spatially heterogeneous relaxation by consideringthe distance to T c as a spatially fluctuating variable[83], or a hierarchical framework systematically shiftingthe factorization approximation to high-order dynamicmulti-point correlations [84, 85, 86, 87]. The latter ap-proach, referred to as generalized mode-coupling the-ory (GMCT), has recently been examined numericallyfor Percus-Yevick (PY) hard spheres by performing ex-plicitly wavenumber- and time-dependent calculationsup to sixth order [87]. The results indicate that theinclusion of more levels in the GMCT hierarchy leadsto a systematic increase of γ (and of the predicted crit-ical packing fraction ϕ c ). Due to (35) an increase of γ implies a decrease of b and so more stretching (cftable 1 in [87]). Qualitatively, we can thus expect thatthe inclusion of more levels in the GMCT hierarchywould probably lead to a better approximation of theactivated dynamics. From the fits to (54) we obtain the q dependence of f c ( q ), h ( q ), B ( q ) and of their incoherent counterparts.The nonergodicity parameters, f c ( q ) and f sc ,α ( q ), andthe critical amplitude, h ( q ), were also calculated bybinary MCT based on the simulated static input.Figure 7 to figure 10 show the results.As seen in figure 7(a), the nonergodicity parame-ters from the fits (symbols) and the MCT calculations(lines) are in semiquantitative agreement. For q & q ∗ the MCT calculations tend to lie below the fit results.This trend is evident for the incoherent scattering andalso visible for q & 10 in the coherent scattering. Asimilar underestimation was observed for polydispersehard spheres and rationalized as follows [39]: MCT pre-dicts structural arrest at T MCTc > T MDc . As the glassstiffens with decreasing T , one can expect f c ( q ) fromthe fits to be larger than from the MCT calculations.This argument is corroborated by the GMCT analy-sis of the PY hard sphere system, which finds f c ( q ) toincrease with increasing ϕ c [87].For q < q ∗ the agreement between the fit results q f c ( q ) f sc ( q ) MD: f sc,A MD: f sc,B MD: f c MCT: f sc,A MCT: f sc,B MCT: f c S(q)/10 (a) q f c ( q ) q f sc (q) AB S(q)/10 (b) Figure 7. Panel (a): q dependence of f sc , A ( q ), f sc , B ( q ) and f c ( q ). The symbols (labeled “MD”) depict the results of thefit to (54). The full lines (labeled “MCT”) correspond to theresults of the MCT calculations based on static input. Panel(b):The symbols in the main figure and in the inset reproduce thefit results from panel (a). The full line in the main figure shows S ( q ) / [ S cc ( q ) S ( q )] related to composition fluctuations [cf (57)]. S nc ( q ), S cc ( q ) and S ( q ) are obtained from the partial structurefactors at T = 0 . 85, cf section 3. The dashed lines in the insetpresent the Gaussian approximation (58) with r sc , A = 0 . r sc , B = 0 . f sc , A ( q )and f sc , B ( q ) to q → 0. In both panels the dotted line shows S ( q )divided by 10 for comparison. and MCT calculations improves with decreasing wavevector. For small q , f c ( q ) strongly increases and tendsto a value of about 0.9 in the q → q dependence for small q and f c ( q → ≈ . q → f c ( q → 0) = lim q → (cid:20) S ( q ) S cc ( q ) S ( q ) (cid:21) . (57)The ratio S ( q → /S cc ( q → 0) corresponds to the lassy dynamics of a binary Voronoi fluid δ S cc ( q → T = 0 . 85 (cf figure 1(c)) we canestimate the term in the square brackets of (57). Thefull line in figure 7(b) presents the result. We find goodagreement with f c ( q ) from the fits for q < 5, therebyconfirming (57). For q & q ∗ , on the other hand, f c ( q ) isin phase with S ( q ) (cf dotted line in figure 7(b)) and thecontribution due to composition fluctuations decreasesin amplitude with increasing q . This suggests thatcomposition fluctuations do not play a prominent rolefor q & q ∗ , a conclusion that resonates with the findingsof [39] and the MCT predictions in [70]. Figure 7(b)thus indicates that a crossover between a composition-fluctuation dominated small- q regime and a packingdominated large- q regime occurs at q ≈ 5, leading to aminimum in f c ( q ) at q ≈ q ≈ F c αβ ( q ) determining f c ( q ) via (56). From figure 8(a)we see that the components of like particles, F cAA ( q )and F cBB ( q ), are always positive (cf inset), and so istheir sum (crosses). By contrast, the component ofunlike particles, F cAB ( q ), becomes negative for q < q ∗ ,leading to a shallow minimum at q ≈ F c ( q ) = F cAA ( q ) + F cBB ( q ) + 2 F cAB ( q ) of (56)is calculated (full line). The depth of the minimumis amplified after division by S ( q ) (dotted line). Asseen in figure 8(a), S ( q ) is similar in magnitude to F c ( q ) for q → q ≈ q ∗ , while S ( q ) > F c ( q )for q ≈ 5. This gives rise to values near 1 for thenormalized nonergodicity parameter f c ( q ) for q → q ≈ q ∗ , and explains the pronounced minimum at q ≈ f sc ,α ( q ) = 1 − ( qr sc ,α ) for q → r sc ,α is the “Lindemann localization length” of species α .Fitting this relation for q . r sc , A = 0 . ≈ . × (2 R A ) and r sc , B =0 . ≈ . × (2 R B ) where R A and R B are thenatural radii of the Voronoi mixture (cf section 2). Ifwe take R A and R B as approximations for the particleradii, we see that the localization lengths are on theorder of 10% of the particle diameters, as suggested byMCT [22, 58, 70]. Moreover, MCT predicts that theGaussian approximation, f sc ,α ( q ) = exp( − q r ,α ) ( α, β = A , B) , (58)gives a reasonable description of the q dependence ofthe nonergodicity parameter. The inset in figure 7(b)confirms this expectation.Figure 9(a) displays the critical amplitude h ( q )for the coherent scattering. The circles correspond tothe fit results, the full line to the MCT calculations.Recall from section 5.2 that the fits involve a constant, q -0.500.511.522.53 F c F α β ( q ) F c F AA +F c F BB F c F AB F c q -0.4-0.20.00.20.40.60.81.0 F c F αβ (q) AABBAB S(q) (a) q H α β ( q ) H AA +H BB H AB H q -0.20.00.20.40.6 H αβ (q) AABBAB S(q)/5 (b) Figure 8. Results from the fully microscopic MCT calculationsfor the nonergodicity parameters F c αβ ( q ) (panel (a)) and thecritical amplitudes H αβ ( q ) (panel (b)) at T c = 0 . 979 245(cf table 1). Both the nonergodicity parameters and criticalamplitudes are not normalized by S ( q ). The insets in panel(a) and panel (b) show the partial components of like (AA,BB) and unlike (AB) particles. The main figure in panel(a) presents F cAA ( q ) + F cBB ( q ) (crosses), F cAB ( q ) (circles), and F c ( q ) = F cAA ( q ) + F cBB ( q ) + 2 F cAB ( q ) (full line). The dottedline indicates S ( q ) at T c . The main figure in panel (b) depicts H AA ( q ) + H BB ( q ) (crosses), H AB ( q ) (circles), and H ( q ) = H AA ( q )+ H BB ( q )+2 H AB ( q ) (full line). The dotted line indicates S ( q ) / T c for comparison. but arbitrary, scale factor ℓ : To fix this factor weadjust ℓ so that the fitted h ( q ) closely matches the h ( q ) from the MCT calculations (here we took ℓ =0 . q dependence can be bettercompared. For q > q ∗ figure 9(a) shows that fitsand MCT agree well with each other, albeit theagreement is a bit worse than for f c ( q ) (cf squaresand dashed line). The MCT calculations indicatethat h ( q ) oscillates in phase with f c ( q ) for q > q ∗ ,whereas it is in antiphase with f c ( q ) for q ≤ q ∗ . For q > q ∗ the fitted h ( q ) has the same q dependenceas the MCT calculations. On the other hand, for2 . q . q ∗ —that is, in the regime where composition lassy dynamics of a binary Voronoi fluid q f c ( q ) h ( q ) MD: f c MCT: f c MD: hMCT: h S(q)/10 (a) q h ( q ) h s ( q ) h s,A h s,B h S(q)/10 (b) Figure 9. Panel (a): Critical amplitude h ( q ) versus q . Thecircles (labeled “MD”) depict the results of the fit to (54). Thefull line (labeled “MCT”) corresponds to the result for h ( q ) fromthe MCT calculations based on static input. For comparison f c ( q ) from figure 7 is reproduced (squares: fit results to (54),dashed line: MCT calculations). The dotted line shows S ( q )divided by 10. Panel (b): q dependence of h s , A ( q ), h s , B ( q )and h ( q ). The symbols depict the results of the fit to (54).The full lines present the Gaussian approximation (59) with h Amsd = 0 . h Bmsd = 0 . h s ,α ( q ) to q . r ,α fixed atthe values given in figure 7. The dotted line shows S ( q ) dividedby 10 for comparison. fluctuations become important—qualitative differencesoccur. The fit results exhibit an oscillation, while thecalculations rather predict a weak shoulder at q ≈ q , in qualitativeagreement with other MCT studies [70]. The presenceof the shoulder can be understood from the interplayof the partial components H αβ ( q ) determining h ( q )via (56). Figure 8(b) shows that, similar to thenonergodicity parameters, the components of likeparticles, H AA ( q ) and H BB ( q ), are always positive,whereas the component of unlike particles, H AB ( q ),becomes negative for q < q ∗ (cf inset). However,contrary to the nonergodicity parameters, the sum H AA ( q ) + H BB ( q ) (crosses) increases steeply in therange q ∼ 4. This increase cannot be outweighedby H AB ( q ) so that the numerator H ( q ) = H AA ( q ) + H BB ( q )+ 2 H AB ( q ) of (56) plateaus for q ≈ S ( q ). To verify our fit results for 2 . q . q ∗ weattempted to impose in (54) the value of h ( q ) fromthe MCT calculations, whereas f c ( q ) and b were fixedat the values found before from the fits. A fit ofcomparable quality is only obtained, if we allow t ′ σ to depend on q , which is not acceptable within MCT.Therefore, it seems we cannot achieve better agreementbetween fits and MCT calculations for 2 . q . q ∗ . Forsmaller q , however, the fit and MCT results appear toconverge again toward one another. Both approachessuggest that h ( q ) . . q < 1. This value is muchsmaller than in the one-component PY hard-spheresystem [57] and may be attributed to a composition-fluctuation effect [39, 70].For the critical amplitudes h s ,α ( q ) of the species α = A , B no MCT calculations are currently availablefor the Voronoi mixture. The corresponding fit resultsare shown in figure 9(b). For q & 10 we find that h s , A ( q ) and h s , B ( q ) bracket h ( q ), as it is also observedfor the nonergodicity parameters in figure 7(a). MCTcalculations for the one-component PY hard-spheresystem suggest that h s ( q ) vanishes in the limits q → q → ∞ and has a maximum near the second peakof S ( q ) [58]. Similar behavior is found here for theA and B particles. In particular, figure 9(b) showsthat h s ,α ( q ) → q → 0. In this limit, the criticalamplitude is supposed to be well described by theGaussian approximation [58], h s ,α ( q ) = h msd ,α q exp( − q r ,α ) , (59)where h msd ,α are constants. We fix r ,α to the valuesfrom figure 7(b) and fit (59) for q . h msd ,α . This gives h msd , A =0 . h msd , B = 0 . q < B ( q ) and B s ,α ( q ), change sign. This isexpected from the literature on MCT [57, 58, 70].However, comparison of these literature results andthe data in figure 10 also reveals some differences, for q . q ∗ . From MCT calculations for binary mixtures[70] one expects B ( q ) to be in phase with h ( q ) for q . q ∗ , to be negative at q ∗ and to tend to a smallpositive value for q → 0. Figure 10(a) does not supportthis expectation. Moreover, for the tagged-particledynamics the Gaussian approximation should becomevalid in the limit q → 0, predicting that B s ,α ( q ) islarger than the constant B s ,α = lim q → B s ,αq [58]. Theconstant B s ,α enters the long-time correction to thevon Schweidler law for the mean-square displacement(MSD) of species α , cf (60). We determined B s ,α from the MSD and the results are shown as horizontaldashed lines in figure 10(b). While B s , B q > B s , B0 , this lassy dynamics of a binary Voronoi fluid q -0.50.00.51.0 h ( q ) B ( q ) hB S(q)/3 - 0.6 (a) q -3-2-101 B ( q ) B s ( q ) B s,A B s,B B S(q)/2 - 3.25 (b) B B s,A B B s,B Figure 10. Panel (a): Long-time correction amplitude B ( q )versus q . The circles depict the results of the fit to (54). Thesquares reproduce the fit results for the critical amplitude h ( q )from figure 9. The dotted line shows S ( q ) / − . q dependence of B s , A ( q ), B s , B ( q ) and B ( q ). Thehorizontal dotted lines indicate the limit B s ,α = lim q → B s ,α ( q )obtained from fits of (60) to the MSD of species α , B s , A0 = − . B s , B0 = − . S ( q ) / − . 25 for comparison. is not the case for the A particles. As pointed out in[36], the determination of the correction amplitudes isimpeded for data which cannot be chosen close enoughto T c [36]. In part, the here described differences maybe attributed to such uncertainties. Figure 11(a) shows the MSD of the A particles, g , A ( t ),and of the B particles, g , B ( t ), at T = 0 . 84. Forboth species the MSD starts from the ballistic regime(3 T t ). Outside this regime, the small (B) particlesalways move much farther than the large (A) particlesin a given time. For t > . -1 t -3 -2 -1 g , A ( t ) g , B ( t ) AB A t6D B tt’ σ sc,A r sc,B r (a) T=0.84 -4 -3 -2 -1 t/t’ σ -3 -2 -1 g , A ( t ) T=0.88T=0.86T=0.85T=0.84T=0.83 -2 -1 | ε | -1 D α t ’ σ sc,A r (b) BA~| ε | -1 Figure 11. Panel (a): Time dependence of the mean-squaredisplacements (MSDs) at T = 0 . 84 for the A particles g , A ( t )(blue full line) and the B particles g , B ( t ) (orange full line). Theballistic motion (3 T t ) at short times and the diffusive motion(6 D α t ) at long times are indicated by dotted lines. The diffusioncoefficients are D A = 1 . × − and D B = 5 . × − . Thedashed lines present the fit results to (60) where only the long-time correction coefficients B s ,α were adjusted, yielding B s , A0 = − . B s , B0 = − . r sc , A = 0 . r sc , B = 0 . h Amsd = 0 . h Bmsd = 0 . r ,α ,are indicated by horizontal dotted lines. The vertical dottedline shows the value of the MCT α time scale t ′ σ (= 359) at T = 0 . 84. Panel(b): Test of the TTSP for the A particles byplotting g , A ( t ) versus t/t ′ σ . Deviations are visible for T = 0 . T = 0 . 88 (dashed line). The inset shows D α t ′ σ as a function of | ε | = | ( T c − T ) /T c | with T c = 0 . / | ε | . the plateau MCT predicts the following relation [58] g ,α ( t ) = 6 r ,α + 6 h α msd B (cid:18) tt ′ σ (cid:19) b − h α msd B B s ,α (cid:18) tt ′ σ (cid:19) b , (60)with the localization lengths r sc ,α [(58)], the criticalamplitudes h α msd [(59)], and the long-time correctioncoefficients B s ,α . Equation (60) is a consequence of(45), since g ,α ( t ) = lim q → − φ s ,α ( q, t )] /q . Whencomparing (60) only the long-time corrections needto be fitted; all other parameters are taken from lassy dynamics of a binary Voronoi fluid -2 -1 t φ ( q ,t ) q=6.85q=11.75q=15.35 KWW Figure 12. Plot of φ ( q, t ) for q = 6 . , . , . 35 at T = 0 . the previous analysis. Figure 11(a) shows that (60)describes the MSD over approximately four decades intime for both species before the crossover to diffusionoccurs at late times. In this long-time regime, g ,α ( t ) =6 D α t with D α being the self-diffusion coefficient ofspecies α .From (46) it follows that the MSD should obey theTTSP when plotting g ,α ( t ) against t/t ′ σ . Figure 11(b)tests this prediction for the A particles in the T intervalwhere φ ( q ∗ , t ) obeys the TTSP (cf figure 2). We seethat the TTSP holds for the MSD only in a narrowertemperature interval (for T = 0 . 84, 0.85 and 0.86),whereas deviations occur for higher and lower T . Thisis highlighted in the inset which plots D α t ′ σ against | ε | = ( T − T c ) /T c . The product D α t ′ σ is not constantover the whole interval 0 . ≤ T ≤ . 88, but appears toincrease as 1 / | ε | . With (35) and (36) this would implya fractional Stokes-Einstein relation [88] D ∼ / ( t ′ α ) ξ with exponent ξ = ( γ − /γ ≈ . α relaxation The KWW function (49) is often used as a convenientparameterization of the α process in experiments andsimulations [53, 54, 89]. When fitting the α relaxationwith (49) similar caveats as discussed for the late β analysis (cf section 5.2) apply: The parameters A ( q ), τ K ( q ) and β K ( q ) are sensitive to the choice of the timeinterval employed for the fit [39, 90, 91], in particularthe stretching exponent appears to be plagued by thiseffect [52, 92]. To guide the KWW fits we here drawupon the asymptotic MCT results from section 4.4 andsubject the fits to two constraints. First, since (49) isa model for the α process, we require A ( q ) ≤ f c ( q ).Second, the early α process should be excluded fromthe fit because β K ( q ) = b for finite q and so theshort-time expansion of (49) cannot agree with the von Schweidler law (43) [62]. Different strategies to copewith this problem have been proposed (see [39, 52, 92]and references therein). One possibility is to focus onthe late α process only [93, 94] by restricting the fit totimes for which φ ( q, t ) is smaller than f c ( q ) by somefactor x cut < 1. We varied x cut in the interval [0 . , . x cut = 0 . φ ( q, t ) at T = 0 . 84 and three wave vectors. Asdesired, the KWW function (dotted lines) providesa good description of the final relaxation and barelyoverlaps with the early β process (von Schweidler law,dotted lines) for q = 6 . 85 and 11.75. For q = 15 . q & 15 corresponds to the asymptotic large- q regimewhere we may expect (50) and (51) to hold. Analysisof the q dependence of the stretching exponents andrelaxation times can test this expectation.Figure 13(a) shows the results for the stretchingexponents and figure 13(b) for the relaxation times.For q & q ∗ the stretching exponent β K ( q ), obtainedfrom φ ( q, t ), is roughly in phase with S ( q ) and tendsto the von Schweidler exponent b for large q . Thesame large- q asymptote is also found for β K ,α , thestretching exponents of φ s ,α ( t ). Along with that, therelaxation times τ K ( q ) and τ K ,α ( q ) for coherent andincoherent scattering also converge to the same large- q asymptote which is proportional to 1 /q /b . Thesefindings agree with the MCT predictions (50) and (51).However, a reservation has to be mentioned: Fromfigure 13(a) it seems as if the limit lim q →∞ β K ( q ) = b isapproached from below. However, according to theory[90, 87], the limit should be approached from above.Such an approach has been seen in several simulations[36, 45, 48, 62, 77, 95]. Certainly, data with highaccuracy at long times are needed to verify (51), sincethe amplitude of the α process becomes small at large q (cf figure 7). This may be a prime source of uncertaintyin the present analysis.In the hydrodynamic limit we expect all scatteringfunctions to decay as single exponentials: φ s ,α ( q, t ) ∝ exp( − q D α t ) due to self diffusion and φ ( q, t ) ∝ exp( − q D int t ) due to interdiffusion, with D int beingthe interdiffusion coefficient [70]. Therefore, β K ( q → 0) = β K ,α ( q → 0) = 1 and τ K ( q ) ∼ τ K ,α ( q ) ∼ /q for q → 0. For q < q ∗ we see from figure 13(a) that thestretching exponents increase toward 1 with decreasing q , but clearly the linear dimension of the simulationbox is still too small so that the hydrodynamic limitis not reached for the smallest accessible q values. Bythe same token, we cannot expect τ K ( q ) or τ K ,α ( q )to attain the hydrodynamic limit. Still, figure 13(b)shows that τ K ,α ( q ) tend to the expected behavior, lassy dynamics of a binary Voronoi fluid q β K ( q ) β K , α ( q ) β K,A β K,B β K S(q)/10 (a) b q τ K ( q ) τ K , α ( q ) τ K,A τ K,B τ K B q (b) A q ~1/q Figure 13. Panel (a): q dependence of the KWW stretchingexponents β K , A ( q ) (crosses), β K , B ( q ) (triangles) and β K ( q )(circles). The horizontal dashed line indicates the value of vonSchweidler exponent b = 0 . β regime.The dotted line shows S ( q ) divided by 10 for comparison. Panel(b): Log-log plot of the KWW relaxation times at T = 0 . q : τ K , A ( q ) (crosses), τ K , B ( q ) (triangles) and τ K ( q )(circles). For the tagged-particle dynamics the dashed linespresent the behavior 1 / ( D α q ) expected for q → D α taken from figure 11. The fullline indicates the MCT prediction ∼ /q /b for large q with b = 0 . τ K ,α ( q ) = 1 /D α q , for q → D int (this would bepossible via an Einstein relation similar to the one forthe self-diffusion coefficients [26]). However, [26] sug-gests that the following linear combination of the self-diffusion coefficients, known as the “Darken equation”, D int = x A x B S cc ( q → (cid:0) x A D B + x B D A (cid:1) , (61)represents a good approximation even in the super-cooled regime. We estimate D int from the data shownin figure 1 and figure 11. The result (1 /D int q ) is in-cluded as a dashed line in figure 14. This figure com-pares the Voronoi mixture to the polydisperse hard-sphere-like model studied in [39] in order to assess towhat extent the q dependence of τ K ( q ) is model spe-cific. For a better comparison we superimpose the dataat one point, q max and τ K ( q max ), where q max ≈ q ∗ for -1 q/q max -1 τ K ( q ) / τ K ( q m a x ) Voronoipolydisperse HS int q ~1/q Figure 14. Comparison of the binary Voronoi mixture (circles)and the polydisperse hard-sphere-like model (crosses) of [39]:The figure shows KWW relaxation times τ K ( q ). The data of thehard-sphere model were digitized from the upper panel of figure 7in [39]. The axes are scaled by q max ≈ q ∗ (Voronoi: q max = 7,Hard spheres: q max = 7 . 1) and τ K ( q max ) (Voronoi: τ K ( q max ) =677 . 74, Hard spheres: τ K ( q max ) = 0 . ∼ /q /b with b = 0 . b = 0 . 53 from [39]. The dashed line shows the hydrodynamicbehavior 1 /D int q where the interdiffusion coefficient D int wasestimated for the Voronoi mixture from the Darken equation(61). both models. We see that the relaxation times for bothmodels are in good qualitative agreement. For large q they are compatible with the scaling ∼ /q /b with amodel-specific von Schweidler exponent and for small q they tend to the hydrodynamic behavior. For the q regime near q max ≈ q ∗ the agreement is even semiquan-titative. In particular, the drop of τ K ( q ≈ . q max ) byan order of magnitude relative to τ K ( q max ) is presentfor both models. This drop is accompanied by a lowamplitude of the α process (cf figure 7 and figure 5 in[39]) and a pronounced stretching of the KWW func-tion (cf figure 13 and figure 8 in [39]). These fea-tures therefore appear to be independent of the modeland rather characteristic of the collective dynamics inmulticomponent systems on length scales where thecrossover between large-scale composition fluctuationsand local-scale liquid-like packing constraints occurs.Figure 14 also suggests that the hard-sphere-likemodel reaches the hydrodynamic limit ( ∼ /q ) earlierthan the Voronoi mixture. A slow convergence to thehydrodynamic limit was also observed for the soundattenuation in the monodisperse Voronoi liquid andcould be traced back to the fact that the product ofthe infinite frequency shear modulus ( G ∞ ) and theisothermal compressibility ( χ T ) is exceptionally small(compared Lennard-Jones systems) [18]. It would beworthwhile to explore whether a similar mechanismalso protracts the crossover to the hydrodynamic limitfor the interdiffusion process in the binary Voronoimixture. lassy dynamics of a binary Voronoi fluid 6. Summary and discussion The Voronoi liquid is a fluid model whose interactionsare local, many-body and soft [9, 18]. Here we study ageneralization of the Voronoi liquid to binary mixtures.Our mixture is equimolar, weakly polydisperse andadditive. This binary Voronoi mixture is a relativelynew model. Up to now, only its thermodynamic andstructural properties, from the normal liquid to thesupercooled state, have been investigated [19]. Withthe present work we extend the analysis to dynamicproperties. The focus of our analysis is a comparison ofMD results for the incoherent and coherent scatteringfunctions with the idealized MCT. Overall, we findthat the glassy dynamics of the binary Voronoi fluidconforms to the same qualitative phenomenology asthat of simple liquids, albeit with a few subtleties.As in every multicomponent system, the binaryVoronoi mixture exhibits transport processes related tocomposition fluctuations. In the hydrodynamic limit,these processes are described by the interdiffusion ofthe two particle species. The idealized MCT obeys thishydrodynamic limit and makes a number of predictions[70]. For q → φ ( q, t )is determined by the ratio S ( q → /S cc ( q → 0) ofthe Bhatia–Thornton structure factors, φ ( q, t ) decaysexponentially and the corresponding relaxation time isgiven by 1 /D int q . Although the systems simulated arestill too small to fully realize the hydrodynamic limit,figures 7, 13 and 14 reveal that our simulation resultsapproach the predicted behavior with decreasing q . Inthis small- q regime the α process of φ ( q, t ) is dominatedby transport processes due to composition fluctuations.A hallmark of glassy slowing down is the super-Arrhenius increase of the local relaxation times withdecreasing T . Figure 2 provides an example for τ q ∗ .MCT attributes this slowing down to the nonlinearcoupling between dynamic density fluctuations, whichamplifies weak structural changes of the dense packingin the neighbor shells of the liquid (“cage effect”).As a consequence, the α process of φ ( q, t ) exhibitsthe fingerprint of S ( q ) for q & q ∗ . We find evidencefor this in-phase modulation with S ( q ) for f c ( q )(figure 7), β K ( q ) (figure 13) and τ K ( q ) (figure 14).Therefore, at intermediate q a crossover exists betweenthe composition-fluctuation dominated small- q regimeand the cage-effect dominated large- q regime. Thiscrossover occurs in the range q ≈ . q ∗ , not only for theVoronoi mixture but also for polydisperse hard spheres(figure 14). Here the amplitude of the α process is weakand τ K ( q ) is about an order of magnitude smaller than τ K ( q ∗ ), while the decay of φ ( q, t ) is strongly stretched.We compare our MD simulations with two MCTapproaches, with fits to the asymptotic predictionsvalid for T & T c and with MCT calculations usingthe partial static structure factors from the simulations as input to compute the dynamics. Fits to theasymptotic predictions have been carried out formany experimental and simulated systems in the past[33, 53], including binary Lennard-Jones [60, 61, 94]and hard-sphere mixtures [36] or metallic alloys [27].Compared to these studies, we get similar results forthe Voronoi mixture, despite its more complicatedmany-body potential. The MCT α time ( t ′ σ ) is stronglycoupled to the α relaxation times of the coherentand incoherent scattering functions at q ∗ (cf figure 5),allowing for a consistent extrapolation from all of theserelaxation times to estimate T c (= 0 . T & T c we find evidence for the space-time factorizationin the β regime (figure 3) and the TTSP in the α regime (figure 2) from the scattering functions at finitewave vectors. On the other hand, time-temperaturesuperposition by scaling time with t ′ σ appears tobecome violated for q → 0, as shown for the MSDin figure 11, implying a decoupling of the α relaxationtime and self-diffusion. It could be that single-particlehopping processes are responsible for this decoupling[37, 42, 79, 80, 96]. Investigations in this directionfor the Voronoi mixture, following e.g. the lines of[97, 98, 99], would be interesting.The binary MCT calculations based on staticinput give very good agreement for f c ( q ) (figure 7),whereas the agreement is worse for h ( q ), in particularin the regime of the crossover between compositionfluctuations and cage effect (figure 9). We note thatour MCT calculations have used only the partial staticstructure factors, i.e. two-point correlation functions,as input, even though the fluid itself contains many-body interactions by construction. In this regard,it may be considered striking that some of theMCT predictions are in such good agreement withsimulation. Indeed, our work suggests that even fora complex fluid such as the Voronoi mixture, oneof the simplest measures of structure (i.e. S αβ ( q ))already constitutes a major portion of the relevantstructural information needed to predict the dynamics.Nonetheless, discrepancies in e.g. the prediction for h ( q ) highlight the need for more refined theory.Currently, the origin of these discrepancies is unclear.To resolve this issue, it would be worthwhile to carryout the comparison between MCT and simulation forthe partial dynamic structure factors S αβ ( q, t ) becausethey are the primary correlators calculated by thetheory (cf section 4). Such a comparison would allowone to identify whether the observed differences in h ( q )stem from one particle species (A or B), or from theinterplay between them. Unfortunately, S αβ ( q, t ) wasnot determined in the present simulations, but work inthis direction is planned for the future.The MCT calculations also illustrate the very highprecision required of T c to get convergent results for λ lassy dynamics of a binary Voronoi fluid λ only settles if T c is accurate to the fifth orsixth digit after the decimal point. Still, the final valueis not so satisfying when compared to the results fromthe asymptotic analysis (cf table 1). The α process ismore stretched than predicted by MCT (figure 6). Thisdifference could be related to the overestimation of T c ( ≃ . α process increases. Itmight therefore be worthwhile to extend the GMCT tobinary mixtures, as studied here. Acknowledgments Financial support by the ANR LatexDry projectgrant ANR-18-CE06-0001 of the French AgenceNationale de la Recherche, the Canada First ResearchExcellence Fund, Quantum Materials and FutureTechnologies Program and by the Dutch ResearchCouncil (NWO) through a START-UP grant isgratefully acknowledged. We are indebted to MFuchs (Konstanz), H Meyer, A N Semenov (bothStrasbourg) for very helpful discussions and to OBenzerara (Strasbourg) for valuable technical supportwith the MD simulations. The simulations were madepossible by a generous grant of computer time on theHPC cluster of the University of Strasbourg. [1] Cavagna A 2009 Phys. Rep. Rev. Mod. Phys. Spatial Tessellations: Concepts and Applications ofVoronoi Diagrams (Wiley)[4] Starr F W, Sastry S, Douglas J F and Glotzer S C 2002 Phys. Rev. Lett. Eur.Phys. E Phys. Rev. Lett. Soft Matter Phys. Rev. Lett. EPL Phys. Rev. X Proceedings of the NationalAcademy of Sciences Proceedings of the NationalAcademy of Sciences J. Phys. Condens. Matter EPL Proc.Natl. Acad. Sci. USA J. Chem. Phys. The Voronoi liquid : a new modelto probe the glass transition J. Chem. Phys. Phys. Rev. E Phys. Rev. X Phys. Rev.X Phys. Rev. E Comput. Phys. Chaos [25] Bhatia A B and Thornton D E 1970 Phys. Rev. B Phys. Rev. B Phys. Rev. B Mol. Phys. Phys. Rev. E Phys. Rev. E Molecular Theory of Solutions (Oxford:Oxford University Press)[32] Kuhn P, Horbach J, Kargl F, Meyer A and Voigtmann T2014 Phys. Rev. B Complex Dynamics of Glass-FormingLiquids: A Mode-Coupling Theory (Oxford: OxfordUniversity Press)[34] Reichman D and Charbonneau P 2005 J. Stat. Mech.Theor. Exp. P05013[35] Janssen L M C 2018 Front. Phys. Phys. Rev. E Phys. Rev. E Phys. Rev. E (4) 041503[39] Weysser F, Puertas A M, Fuchs M and Voigtmann T 2010 Phys. Rev. E Phys. Rev. Lett. J. Non-Cryst.Solids Phys. Rev. E Europhys. Lett. Phys. Rev. E Phys. Rev. E Eur. Phys. E Phys. Rev. E J. Phys.: Condens. Matter J. Phys.:Condens. Matter J. Stat. Phys. Proceedings of the Les Houches Summer School ofTheoretical Physics, Les Houches 1989, Session LI edHansen J P, Levesque D and Zinn-Justin J (Amsterdam:North-Holland) pp 287–503[52] Baschnagel J and Varnik F 2005 J. Phys.: Condens.Matter R851[53] G¨otze W 1999 J. Phys.: Condens. Matter A1[54] G¨otze W and Sj¨ogren L 1992 Rep. Prog. Phys. lassy dynamics of a binary Voronoi fluid [55] Kob W 2003 Supercooled liquids, the glass transition, andcomputer simulations Slow relaxations and nonequi-librium dynamics in condensed matter ed BarratJ L, Feigelmann M, Kurchan J and Dalibard J (LesUlis/Berlin: EDP Sciences/Springer) pp 201–269[56] Kob W 1999 J. Phys.: Condens. Matter R85–R115[57] Franosch T, Fuchs M, G¨otze W, Mayr M R and Singh A P1997 Phys. Rev. E Phys. Rev. E J. Non-Cryst. Solids Phys. Rev. E Eur. Phys. J. B Phys. Rev.E J. Phys.: Condens. Matter Phys. Rev. E J. Chem.Phys. J. Phys.: Condens. Matter Phys.Rev. E Eur. Phys. J. E. J. Phys.: Condens.Matter A261[70] Fuchs M and Latz A 1993 Physica A Phys. Rev. E Phys. Rev. E Phys. Rev. E J. Phys.: Condens. Matter Phys. Rev. Lett. J. Phys.: Condens. Matter Phys. Rev. E Phys. Rev. A Phys. Rev. E J. Phys.:Condens. Matter J. Chem. Phys. J. Chem. Phys. Europhys. Lett. Phys. Rev. Lett. Phys. Rev. Lett. J. Stat.Mech. arXiv:1909.0042 [88] Parmar A D S, Sengupta S and Sastry S 2017 Phys. Rev.Lett. J. Appl. Phys. Phys. Rev. A Phys. Rev. E , 5625(1999)[92] Aichele M and Baschnagel J 2001 Eur. Phys. J. E Phys. Rev. E Phys. Rev. E Phys. Rev.E PNAS Soft Matter J. Stat. Mech. Theory Exp. Phys. Rev. X Phys.Rev. E90