{'fraction_non_alphanumeric': 0.03593083314619358, 'fraction_numerical': 0.012575791601167752, 'mean_word_length': 4.474083179676295, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
A theory of amorphous packings of binary mixtures of hard spheres

We extend our theory of amorphous packings of hard spheres to binary mixtures and more generally to multicomponent systems. The theory is based on the assumption that amorphous packings produced by typical experimental or numerical protocols can be identified with the infinite pressure limit of long lived metastable glassy states. We test this assumption against numerical and experimental data and show that the theory correctly reproduces the variation with mixture composition of structural observables, such as the total packing fraction and the partial coordination numbers.

We extend our theory of amorphous packings of hard spheres to binary mixtures and more generally to multicomponent systems. The theory is based on the assumption that amorphous packings produced by typical experimental or numerical protocols can be identified with the infinite pressure limit of long lived metastable glassy states. We test this assumption against numerical and experimental data and show that the theory correctly reproduces the variation with mixture composition of structural observables, such as the total packing fraction and the partial coordination numbers.
Amorphous packings of hard spheres are ubiquitous in physics: they have been used as models for liquids, glasses, colloidal systems, granular systems, and powders. They are also related to important problems in mathematics and information theory, such as digitalization of signals, error correcting codes, and optimization problems. Moreover, the structure and density (or porosity) of amorphous multicomponent packings is important in many branches of science and technology, ranging from oil extraction to storage of grains in silos.
Despite being empirically studied since at least sixty years, amorphous packings still lack a precise mathematical definition, due to the intrinsic difficulty of quantifying "randomness" [1]. Indeed, even if a sphere packing is a purely geometrical object, in practice dense amorphous packings always result from rather complicated dynamical protocols: for instance, spheres can be thrown at random in a box that is subsequently shaken to achieve compactification [2], or they can be deposited onto a random seed cluster [3]. In numerical simulations, one starts from a random distribution of small spheres and inflates them until a jammed state is reached [4,5]; alternatively, one starts from large overlapping spheres and reduces the diameter in order to eliminate the overlaps [6][7][8][9]. In principle, each of these dynamical prescriptions produces an ensemble of final packings that depends on the details of the procedure used. Still, very remarkably, if the presence of crystalline regions is avoided, the structural properties of amorphous packings turn out to be very similar. This observation led to the proposal that "typical" amorphous packings should have common structure and density; the latter has been denoted Random Close Packing (RCP) density. The definition of RCP has been intensively debated in the last few years, in connection with the progresses of numerical simulations [1,10].
Nevertheless, the empirical evidence, that amorphous packings produced according to very different protocols have common structural properties, is striking and call for an explanation. This is all the more true for binary or multicomponent mixtures, where in addition to the usual structural observables, such as the structure factor, one can investigate other quantities such as the coordination between spheres of different type, and study their varia-tion with the composition of the mixture.
In earlier attempts to build statistical models of packings, only the main geometrical factors, such as the relative size and abundance of the different components, were taken into account [11][12][13]. More precisely, these models focus on a random sphere in the packing and its first neighbors, completely neglecting spatial correlations beside the first shell and all the global geometric constraints. This already accounts for the main qualitative structural properties of random packings. However, in order to obtain a quantitative description, some free parameters have to be introduced and adjusted to match with experimental data.
To go beyond these simple models, many authors proposed that random packings of hard spheres can be thought as the infinite pressure limit of hard sphere glasses [14][15][16][17][18][19][20]. This is very intuitive since a glass is a solid state in which particles vibrate around amorphous reference positions, and vibrations are reduced on increasing pressure. A typical algorithm attempting to create a random packing starts at low density and compresses the system at a given rate. During this evolution, when the density is high enough, relaxation becomes more and more difficult until at some point the system is stuck into a glass state [19,20]; at this point further compression will only reduce the amplitude of the vibrations. In a nutshell, this is why amorphous packings can be identified with glasses at infinite pressure.
The main advantage of this identification, if it holds, is that a glass is a metastable state that has a very long lifetime; therefore, its properties can be studied using concepts of equilibrium statistical mechanics. In this way a complicated dynamical problem (solving the equations of motion for a given protocol) is reduced to a much more simple equilibrium problem. In [18,20] it was shown, in the case of monodisperse packings, that this strategy is very effective since it allows to compute structural properties of random packings directly from the Hamiltonian of the system, without free parameters and in a controlled statistical mechanics framework. Note that the existence of an equilibrium glass transition in hard sphere systems (or in other words the existence of glasses with infinite life time) has been questioned [5,21]. Although very inter-esting, this problem is not relevant for the present discussion since we are only interested in long-lived metastable glasses that trap dynamical algorithms. At present it is very well established by numerical simulations [22,23] that for system sizes of N < ∼ 10 4 particles and on the time scales of typical algorithms, metastable glassy states exist, at least in d ≥ 3. This is enough to compare with most of the currently available numerical and experimental data. Finally, the relation of this approach to special packings such as the MRJ state [1] and the J-point [8] has been discussed in detail in [20].
The aim of this paper is to extend the theory of [18,20] to binary mixtures. This allows to compare quantitatively the predictions of the theory and the results of numerical simulations. We will focus in particular on the variation of density and local connectivity as a function of mixture composition. These results constitute, in our opinion, a stringent test of the assumption that random packings reached by standard algorithms can be identified with infinite pressure metastable glasses.
Theory -The equilibrium statistical mechanics computation of the properties of the glass is based on standard liquid theory [24] and on the replica method [17,25] that has been developed in the context of spin glass theory [26,27]. For monodisperse hard spheres, it has been described in great detail in [20]. The extension to multicomponent systems is straightforward following [28]; details are given in Appendix.
Here we just recall some features of this approach, based on simple physical considerations. The basic assumption of the method are that i) crystallization and phase separation are strongly suppressed by kinetic effects so that the liquid can be safely followed at high density, and ii) that at sufficiently high density, the liquid is a superposition of a collection of amorphous metastable states. Namely, in the liquid, the system spends some time inside one of these states, and sometimes undergoes a rearrangement that leads to a different state [29]. Each state is characterized by its vibrational entropy per particle, denoted by s, and the number of such states is assumed to be exponential in N , so defining a configurational entropy Σ(ϕ, s) = N −1 log N (ϕ, s), being N (ϕ, s) the number of states having entropy s at density ϕ. On increasing the density, the liquid is trapped for longer and longer times into a metastable state, until at some point the transition time becomes so long that for all practical purposes the system is stuck into one state: it then becomes a glass. To compute the properties of the glassy states, the central problem is to compute the function Σ(ϕ, s). This can be done by means of a simple replica method introduced by Monasson [27]. One introduces m copies of each particle, constrained to be close enough, in such a way that they must be in the same metastable state. Then, the total entropy of the system of m copies is given by S(m, ϕ|s) = ms + Σ(ϕ, s); the first term gives the entropy of the m copies in a state of entropy s, while the second term is due to the multiplicity of possible states. The total entropy of the system at fixed density Ref. [30]. Lines are predictions from theory, obtained fixing Σj = 1.7. Note that the large r-small η region could not be explored, since for such very asymmetric mixtures the large spheres form a rigid structure while small spheres are able to move through the pores and are not jammed [12,13].
ϕ is obtained by maximizing over s, i.e.
The knowledge of S(m, ϕ) allows to reconstruct the curve Σ(ϕ, s) for a given density by a parametric plot of Eqs. (2) by varying m. The function Σ(ϕ, s) gives access to the internal entropy and the number of metastable glassy states; from this one can compute their equation of state, i.e. the pressure as a function of the density; in particular, for each set of glassy states of given configurational entropy Σ j , one can compute the density ϕ j (jamming density) at which their pressure diverges. Since ϕ j turns out to depend (slightly) on Σ j , a prediction of the theory is that different glasses will jam at different density: amorphous packings can be found in a finite (but small) interval of density [19,20].
Results for binary mixtures -The details of the computation of the function S(m, ϕ) for a general multicomponent mixture, based on [20], can be found in Appendix.
Here we consider a binary mixture of two types of spheres µ = A, B in a volume V , with different diameter D µ and density ρ µ = N µ /V . We define r = D A /D B > 1 the diameter ratio and x = N A /N B the concentration ratio;  Open symbols are experimental data from Ref. [31]. Note that in the lower right panel a different scale is used for z ls .
Once an equation of state for the liquid has been chosen, the jamming packing fraction ϕ j is given in term of Σ j by the solution of Σ j (ϕ) = Σ j . The average coordination numbers at ϕ j are denoted z µν (ϕ j ), but we checked that their variations with ϕ j are negligible. We used in d = 3 the equation of state proposed in [32], using the Carnahan-Starling equation for the monodisperse system [24]. The latter, as well as Σ j (ϕ) and z µν (ϕ j ), are given respectively in Eqs. (20), (16), (19), in Appendix.
Numerical simulations -We produced jammed packings of binary mixtures of N = 1000 hard spheres using the code developed by Donev et al. [33,34]. In this algorithm spheres are compressed uniformly by increasing their diameter at a rate dD/dt = 2γ while event-driven molecular dynamics is performed at the same time. In order to obtain a perfectly jammed final packing, the later stages of compression must be performed very slowly. On the other hand, at low density slow compression is a waste of time, since the dynamics of the system is very fast. Following [35], we find a good compromise by performing a four stages compression: starting from random configurations at ϕ = 0.1, i) the first stage is a relatively fast compression (γ = 10 −2 ) up to a reduced pressure p = βP/ρ = 10 2 ; then we compress at ii) γ = 10 −3 up to p = 10 3 ; iii) γ = 10 −4 up to p = 10 9 ; iv) γ = 10 −5 up to p = 10 12 . The first stage terminates at a density ϕ ∼ 0.6, and is fast enough to avoid crystallization and phase separation. During the following stages the system is already dense enough to stay close to the amorphous structure reached during the first stage. Still, little rear-rangements (involving many particles) are possible and allow to reach a collectively jammed final state [35]. In the final configurations, we observe a huge gap between contacting (typical gap ∼ 10 −11 D) and non-contacting (typical gap > ∼ 10 −6 D) particles. We then say that two particles are in contact whenever the gap is smaller than 10 −8 D. Typically, a small fraction ( < ∼ 5% of the total) of rattlers, i.e. particles having less then 4 contacts, is present. Once these are removed, the configuration is isostatic (the total number of contacts is 6N ) within 1% accuracy.
Comparison of theory and numerical/experimental data -During the four stages of compression, the pressure initially follows the liquid equation of state up to some density close to the glass transition density [20,36,37]. Above this density, pressure increases faster and diverges on approaching jamming at ϕ j . The exact point where this happens depends on compression rate. This is a nice confirmation of a prediction of the theory, that different glassy states jam at different density; it was already observed in [5,35] and recently discussed in great detail in [23,38]. Within the theory ϕ j is related to Σ j , the value of configurational entropy at which the system falls out of equilibrium; hence there is one free parameter, Σ j , that depends on the compression protocol. The equation of state of the glass obtained numerically with our protocol corresponds within our theory to Σ j ∼ 1.5. We decided to use the value Σ j = 1.7 that gives the best fit to the numerical data. This is consistent with previous observations, that the configurational entropy is close to 1 when the system falls out of equilibrium [39]. A detailed discussion of the behavior of pressure can be found in Appendix.
In figure 1, we report the jamming density for different mixtures, putting together our numerical results and experimental data from Ref. [30], and the theoretical results. Note that a single "fitting" parameter Σ j , that is strongly constrained, allows to describe different sets of independent numerical and experimental data. The prediction of our theory are qualitatively similar to previous ones [12,13], but the quantitative agreement is much better. Interestingly, a similar qualitative behavior for the glass transition density has been predicted in [36,37]; although there is no a priori reason why the jamming and glass transition density should be related [23], it is reasonable to expect that they show similar trends [37].
Finally, in figure 2 we report the average partial contact numbers for different mixtures. These values have been obtained by removing the rattlers from the packing. As discussed above, the total coordination is close to the isostatic value z = 6. This is also a non trivial prediction of the theory, see Appendix. As it can be seen from figure 2, the computed values agree very well with the outcome of the numerical simulation, at least for r not too large. Some discrepancies are observed in the contacts of the large particles for large r. We produced packings of N = 10 4 and checked that these are not finite size effects. Also, inspection of the configurations seem to exclude the presence of phase separation. However, for these values of r and x, a large fraction of rattlers (∼ 10%) is present within the small particles. This might affect the determination of the partial contacts. It would be interesting to check if better results are obtained us-ing different algorithms. Experimental data from [31] are also reported in the right panel of figure 2.
Conclusions -In this paper we have extended our theory of amorphous packings to binary mixtures and have tested it against numerical and experimental data. In particular we have shown that the theory correctly predicts the variation of total density (or porosity) and local coordination with mixture composition. We have also shown that the behavior of pressure during compression follows the predictions of the theory. A striking prediction of the theory is that different compression procedures lead to different final densities, which seems to be confirmed by numerical data, see also [23,38]. Note that the only free parameter is the value of Σ j , that is still strongly constrained (it must be close to 1). It affects slightly the values of density (by varying Σ j in the reasonable range one can change ϕ j of ∼ 10%, see figure 3 in Appendix) and does not affect at all the curves in figure 2 for the local coordinations. As stated in the introduction, we believe that these results constitute a stringent test of the idea that amorphous packings can be considered as the infinite pressure limit of metastable glassy states. Note that our results have no implications on the existence of an ideal glass transition [5,21]. Indeed, we only used numerical data obtained using fast compressions. It is possible that using much slower compressions physics changes dramatically and the glass transition is avoided. Still, these time scales are out of reach of current algorithms and experimental protocols.
Acknowledgments -We wish to thank A. Donev for his invaluable help in using his code, and L. Berthier for many important suggestions and for providing equilibrated configurations at high density. Here we discuss in detail the behavior of pressure during compression in our numerical simulations. In figure 3 we report the evolution of the inverse reduced pressure during the first and second stages of compression (Numerical 1) for a mixture with x = 1 and r = 1.4. We observe, as in similar studies [5,23,35,38], that the pressure follows the liquid equation of state up to some density, that for this system is around ϕ ∼ 0.56. At this value of density, the relaxation time of the liquid becomes long enough that liquid relaxation is effectively frozen on the compression time scale used in this work and the system falls out of equilibrium. Although we did not measure relaxation time directly, this has been done in related works that confirmed the correctness of this statement [23,37,38]. Above this density, pressure increases faster and diverges on approaching jamming around ϕ j ∼ 0.65. The numerical equation of state is compared with that of a glass state corresponding to Σ j = 1.5. This is consistent with previous observations, that the configurational entropy is close to 1 when the system falls out of equilibrium [39]. In the same plot, we report the curve (Numerical 2) obtained starting the first stage from a carefully equilibrated liquid configuration of the same mixture at ϕ = 0.58, kindly provided by L. Berthier (see [23] for details on how this configuration was produced and equilibration was checked). In this case, since the relaxation time of the liquid at that density is already very long compared to our compression rate, the system falls immediately out of equilibrium and the pressure increases fast until jamming occurs at a higher density compared to the previous case. In our interpretation, this corresponds to a glassy state with lower Σ j (compare with the theoretical curve for Σ j ∼ 1.2). This is a nice confirmation of a prediction of the theory, that different glassy states jam at different density. Finally, we report, for the same system, the numerically extrapolated value of the ideal glass transition pressure, p K = 34.4, see Ref. [23] for details of the careful extrapolation procedure. Again, this corresponds well (within 10%) to the computed value p K = 31.8 from the theory, see figure 3. Note that this coincidence does not prove the existence of the Kauzmann transition, since on much larger time and length scales than the ones explored in these simulations, a crossover to a non mean-field behavior might happen. Still, the coincidence shows that, on the length and time scales explored by current numerical simulations, the full mean field phenomenology is observed, including the apparent extrapolation of Kauzmann transition. Whether this transition really exists remains a major open point of the field, that however does not affect the results presented here.

II. REPLICA THEORY FOR MULTICOMPONENT MIXTURES
Here we show in detail how to generalize the computation of [20] to multicomponent mixtures. We will not repeat the discussion of [20] but we only explain how to modify it for mixtures. Reading Appendix B and C and section VII of [20] is necessary to follow the discussion.
We consider a multicomponent system in a (large) volume V with partial densities ρ µ = N µ /N . The total density is ρ = µ ρ µ and we define x µ = ρ µ /ρ. Spheres of type µ have diameter D µ and we define D µν = (D µ + D ν )/2. We denote by Ω d = 2π d/2 /Γ(d/2) the d-dimensional solid angle, and by V d (D) = π d/2 D d /Γ(1 + d/2) the volume of a d-dimensional sphere of diameter D. The hard sphere potential φ µν (r) is infinite for r < D µν and zero otherwise: we denote χ µν (r) = exp[−φ µν (r)]. g µν (r), as usual, is the µ-ν pair correlation function [24]. We will use the shorthand notations g µν ≡ g µν (D µν ) for the contact values of g µν (r), and V µν d = V d (D µν ) for the volume of a sphere of diameter D µν .
As discussed in [28], we assume that in the replicated liquid molecules are built of particles of the same type. This amounts to assume that in a glassy state particles of different type cannot easily exchange. However, as the diffusion constant is always finite in a glass, this assumption is wrong, since particles can exchange also within a state. This is all the more true in the liquid phase before the glass transition. Therefore, in order to correct this error, we will subtract the mixing entropy from the configurational entropy we will compute in the following, in order to obtain the correct physical result.
According to [28], to each molecule of the replicated liquid we can attach a label µ = 1, · · · , n according to the type of particles that build the molecule. We denote coordinates in a molecule byx = (x 1 , · · · , x m ). We have then to compute the free energy functional of a liquid made of n species of molecules, with a single molecule density ρ µ (x) and an interactionχ µν (x,ȳ) = m a=1 χ µν (x a − y a ).

Effective liquid
We assume that the vibrations of the m copies are described by a Gaussian distribution, corresponding to harmonic vibrations: The widths A µ are variational parameters and we will maximize the entropy with respect to them at the end. We choose a replica (say replica 1) as a reference and consider the vibrations of the other m − 1 particles around the reference one. We want to construct an expansion assuming that A µ is small. For A µ = 0, the m − 1 copies coincide with the reference one, and S(m, ϕ) is given by the entropy S(ϕ) of the non-replicated liquid plus the entropy of m − 1 harmonic oscillators of spring constant A µ (averaged over the concentration of the different species): The harmonic part of the entropy is computed straightforwardly from the term µ dxρ µ (x)[1 − log ρ µ (x)] in the free energy functional, see section V of [20], where the function S harm (m, A) is also defined. A first order approximation is obtained by considering the effective two-body interaction induced on the particles of replica 1 by the coupling to the m − 1 copies. This can be justified on the basis of a diagrammatic expansion following the derivation in Appendix B of [20], see also [15]. It is possible to show [24] that the diagrammatic expansion of the free energy functional for a multicomponent liquid is the same as the one of a simple liquid, provided a label µ i is attached to each vertex i, and in such a labeled diagram ρ µi (x i ) is placed on each vertex, χ µiµj (x i − x j ) on each link, and a sum over µ i is performed in addition to the integration over x i . The same is true for the molecular liquid. Hence, the treatment of diagrams performed in Appendix B of [20] can be repeated exactly for a mixture, taking into account the presence of the additional labels. This is straightforward and leads to the introduction of effective potentials that now depend on the type of particles involved. In particular, one obtains: A first order approximation to S(m, ϕ) is then obtained by substituting to the entropy S(ϕ) of the simple hard sphere mixture the free energy of a liquid of particles interacting via the potential φ ef f µν (r): It is now evident that we can obtain better approximations of the true function S(m, ϕ) by considering also the three body interactions induced on particles of the replica 1, and so on. Here, for simplicity, we will limit ourselves to consider the two-body interactions. To the same approximation, the pair correlation function g µν (r) of the system is given by the correlation of the effective liquid of particles interacting via φ µν (r). For small A µ , φ ef f µν is close to φ µν and the free energy in Eq. (6) can be computed in perturbation theory around the normal liquid [24], which is taken as an input for the calculation of the properties of the glass.
To compute the previous expression, all the considerations of Appendix C of [20] can be repeated. Instead of a single function Q(r), one has Q µν (r), and in the two Gaussians one has A µ and A ν instead of A. At the very beginning of the discussion, then, one finds that it is enough to replace 2A → A µ + A ν and follow the same steps. We obtain where Q 0 (m) is defined in Appendix C of [20].
Putting all together, we find the final result This expression generalizes the result of section VII of [20] to multicomponent mixtures; it has to be optimized over all the A µ to get the replicated entropy S(m, ϕ).

Replicated free energy for binary mixtures
The optimization with respect to A µ leads to n coupled equations that are not easy to solve in general. Therefore in the following we focus on the case of binary mixtures that is of interest here.
In the limit ϕ → ϕ j , one has A µ ∝ m → 0 and the integral of Q µν (r) can be easily evaluated [20]. The result is z µν = Ω d ρ ν D d−1 µν g µν 2(a µ + a ν ) (18) where a µ = lim m→0 A µ Q0(m) 1−m 2 . In the explicit case of binary mixtures, note that the equation (13) for δ does not depend on m. Using also Eq. (14), we obtain      