Random Close Packing and the Hard Sphere Percus-Yevick Theory

TL;DR

This paper investigates the Percus-Yevick theory for hard spheres across different dimensions, revealing a critical density where the pair distribution function becomes negative, aligning with the onset of random close packing.

Contribution

It identifies a critical volume fraction where the pair distribution function turns negative, linking it to the random close packing transition across dimensions.

Findings

Pair distribution function becomes negative at a critical density.

Critical densities align with random close packing thresholds.

Supports implications for systems modeled by Percus-Yevick theory.

Abstract

The Percus-Yevick theory for monodisperse hard spheres gives very good results for the pressure and structure factor of the system in a whole range of densities that lie within the liquid phase. However, the equation seems to lead to a very unacceptable result beyond that region. Namely, the Percus-Yevick theory predicts a smooth behavior of the pressure that diverges only when the volume fraction $\eta$ approaches unity. Thus, within the theory there seems to be no indication for the termination of the liquid phase and the transition to a solid or to a glass. In the present article we study the Percus-Yevick hard sphere pair distribution function, $g_2(r)$, for various spatial dimensions. We find that beyond a certain critical volume fraction $\eta_c$, the pair distribution function, $g_2(r)$, which should be positive definite, becomes negative at some distances. We also present an intriguing observation that the critical $\eta_c$ values we find are consistent with volume fractions where onsets of random close packing (or maximally random jammed states) are reported in the literature for various dimensions. That observation is supported by an intuitive argument. This work may have important implications for other systems for which a Percus-Yevick theory exists.

Tables

Table 1. Table 1: A summary of the values η c ​ ( d ) subscript 𝜂 𝑐 𝑑 \eta_{c}(d) for dimensions in the range 3 ≤ d ≤ 9 3 𝑑 9 3\leq d\leq 9 .

d	3	4	5	6	7	8	9
${\eta }_c\left(d\right)$	$0.613$	$0.467$	$0.367$	$0.230$	$0.207$	$0.087$	$0.112$

Equations

$$\rho (\mathbf{r})=\sum _{i=1}^N\delta (\mathbf{r}-{\mathbf{r}}_i),$$

$$g_2(\mathbf{r})=\frac{1}{{\bar{\rho }}^2}⟨\rho (0)\rho (\mathbf{r})-\bar{\rho }\delta (\mathbf{r})⟩$$

$${\rho }_{\mathbf{q}}=\frac{1}{\sqrt{N}}\int d\mathbf{r}\rho (\mathbf{r})e^{-i\mathbf{q}\cdot \mathbf{r}},$$

$$\frac{\partial P^I}{\partial t}=\frac{k_{B}T}{\sqrt{N}}\left[\sum _{\mathbf{k},\mathbf{\ell }}\mathbf{k}\cdot \mathbf{\ell }\frac{\partial }{\partial {\rho }_{\mathbf{k}}}{\rho }_{\mathbf{k}+\mathbf{\ell }}\frac{\partial }{\partial {\rho }_{\mathbf{\ell }}}-\sqrt{N}\sum _{\mathbf{k}}{\mathbf{k}}^2\frac{\partial }{\partial {\rho }_{\mathbf{k}}}{\rho }_{\mathbf{k}}\right]P^I\{\rho \}=0,$$

$$\left\{\sum _{\mathbf{k}}{\mathbf{k}}^2\frac{\partial }{\partial {\rho }_{\mathbf{k}}}\left[\frac{\partial }{\partial {\rho }_{-\mathbf{k}}}+{\rho }_{\mathbf{k}}+{\lambda }_{\mathbf{k}}{\rho }_{\mathbf{k}}\right]-\frac{1}{\sqrt{N}}\sum _{\mathbf{k},\mathbf{\ell }}\mathbf{k}\cdot \mathbf{\ell }\frac{\partial }{\partial {\rho }_{\mathbf{k}}}{\rho }_{\mathbf{k}+\mathbf{\ell }}\frac{\partial }{\partial {\rho }_{\mathbf{\ell }}}\right\}{P}_S^I\{\rho \}=0.$$

$$S^{HS}(\mathbf{k})=\frac{1}{1+{\lambda }_{\mathbf{k}}},$$

$$\lambda (r)=\int d\mathbf{k}{\lambda }_{\mathbf{k}}{e}^{i\mathbf{k}\cdot \mathbf{r}}=0\text{for }r>R,$$

$$g_2(r)=1+\frac{1}{N}\sum _{\mathbf{k}\ne 0}\left(S^{HS}(\mathbf{k})-1\right)e^{i\mathbf{k}\cdot \mathbf{r}}=0\text{for }r<R.$$

$$P_{\mathrm{CS}}=\bar{\rho }{k}_{B}T\frac{1+\eta +{\eta }^2-{\eta }^3}{(1-\eta {)}^3},$$

$$\eta =\frac{\pi R^3}{6}\bar{\rho },$$

$$g_2(r)=1+\sum _{n\ge 1}{\eta }^n{g}_2^{(n)}(r).$$

$$g_2(r)\simeq \frac{{\sum }_{n=0}^N{\eta }^n{u}_n(r)}{{\sum }_{n=0}^M{\eta }^n{d}_n(r)}.$$

$${\eta }_{\mathrm{MRJ}}=\frac{c_1+c_{2}d}{2^d},$$

$${\eta }_{\mathrm{GCP}}\propto \frac{d\cdot \ln d}{2^d}.$$

$$h(\mathbf{r})=c(\mathbf{r})+\bar{\rho }\underset{0}{\overset{\infty }{\int }}h({\mathbf{r}}^{\prime })c(\left|\mathbf{r}-{\mathbf{r}}^{\prime }\right|)d{\mathbf{r}}^{\prime },$$

$$\eta =\bar{\rho }{V}_d\left(\frac{R}{2}\right)={\left(\frac{\pi }{4}\right)}^{d/2}\frac{\bar{\rho }{R}^d}{\Gamma \left(\frac{d+2}{2}\right)},$$

$$\tilde{c}(k)=(2\pi {)}^{d/2}{k}^{-\frac{d-2}{2}}\sum _{i=0}^d{c}_i(\eta )\underset{0}{\overset{1}{\int }}{u}^{i+d/2}{J}_{\frac{d-2}{2}}(ku)du,$$

Full text

I Introduction

The hard sphere model provides a canonical minimalistic model that captures the main ingredient in the description of simple liquids, namely the strong short-range repulsion between atoms in the liquid. As in other systems in equilibrium statistical physics, the model is used to obtain macroscopic observables from the microscopic description of the system. In the case of the hard sphere model the goals are the equation of state, the liquid structure factor and a description of the solidification of the liquid in terms of the average particle density, $\bar{\rho}$ and $R$, the range of the hard sphere interaction (namely the diameter of the hard spheres). A wide arsenal of methods has been applied over the years to the hard sphere problem with considerable success. Monte-Carlo and Molecular Dynamics simulations have been applied to that model as early as the 1950s Wood1957 ; Alder1957 ; Alder1959 and extended much later. Experimental studies, first with ball-bearings Bernal1960 ; Bernal1964 and later with colloids Vladeeren1995 , helped elucidate properties of very dense packings which where not attainable in simulations. This gave rise to the concept of the random close packing (RCP) density, defined as the maximal density among amorphous packings. Improved computation technologies went hand in hand with hard-sphere simulations. For example, the hard sphere system is one of the first systems to be simulated on the Small Web Computing (SWC) platform in recent years Bishop2005 . These important numerical efforts resulted in obtaining the phase diagram of the system, including crystallization and a super dense rotation invariant phase in three dimensions Hoover1968 ; Frenkel1984 ; Bolhuis1996 ; Gruhn2001 . The most trusted analytic tool applied successfully to the hard sphere problem is the virial expansion, which is based in turn on the cluster expansion Mayer1941 ; Morita1961 ; Dominicis1962 ; Dominicis1963 ; Stell1964 ; Hansen2006 ; Clisby2005 ; Clisby2006 ; Zhang2014 .

The other two interesting analytic approaches are the Hyper-Netted-Chain (HNC) approximation Rowlinson1965 and the Percus-Yevick (PY) equation Percus1958 for the structure factor of the hard sphere system. The most appealing to us is the PY equation and that is for several reasons. First the equation has been given exact analytic solutions in odd dimensions $d\leq 7$ Thiele1963 ; Wertheim1963 ; Wertheim1964 ; Baxter1968a ; Leutheusser1984 ; Rohrmann2007 ; Robles2007 ; Robles2004 (where $d$ is the dimension of the system). In fact, an exact analytic solution can be obtained in principle for any odd dimension but it involves solving a polynomial equation of degree $2^{(d-3)/2}$. Thus, the highest dimension for which a strict analytic solution in closed form exists is $7$, due to the Abel-Ruffini theorem. However, thanks to the existence of this developed analytic structure it is possible to obtain semi-analytic results for higher odd dimensions with a simple numerical computation Robles2004 ; Rohrmann2008 . More recently, a systematic analytic method of solution, based on the virial expansion for the PY equation, has been obtained for general dimensions including the even ones AddaBedia2008a ; AddaBedia2008b .

The paper is organized as follows: In section II we present an alternative derivation of the Percus-Yevick equation. In section III we show that unlike the common lore the PY theory breaks down at a critical density, which is intrinsic to the theory. In section IV we discuss the physical meaning of the reported observations and in section V we summarize and provide future perspective.

II An alternative derivation of the Percus-Yevick equation

The PY equation is usually seen as a certain diagrammatic approximation or closure scheme of the full problem Hansen2006 . It was shown, however, that the PY equation for the hard sphere system can be given a very simple and intuitive meaning Schwartz1999 ; Edwards2003 , which elucidates the assumptions underlying it, as well as possible extensions and refinements. For the benefit of the readers we reproduce here a short derivation of the PY equation. Consider the particle number density,

[TABLE]

where ${{\bf{r}}_{i}}$ is the location of particle $i$ which is one of $N$ identical particles enclosed in a cubic container of linear size $L$ and periodic boundary conditions. The pair distribution function,

[TABLE]

yields the $d$ dimensional distribution to find a particle at ${\bf{r}}$ given the existence of another particle at ${\bf{0}}$. The hard sphere system is then viewed as an ideal gas with a pair distribution function which is constrained to vanish for $\left|{\bf{r}}\right|<R$ (where $R$ is the diameter of the hard spheres). To see how it works we have to transform from particle coordinates to collective coordinates Schwartz2002 ; Edwards2003 as described shortly in the following. The natural collective coordinates are the Fourier components of the density,

[TABLE]

for ${\bf{q}}\neq 0$ and with discrete Fourier modes ${\bf{q}}_{\ell}=\frac{2\pi n_{\ell}}{L}$, where $n_{\ell}$ is an integer. The ideal gas Fokker-Planck equation for the distribution of the $N$ free particles is translated into a functional Fokker-Planck equation for the probability to obtain a given configuration of the density, $P^{I}\{\rho\}$ Edwards2003 , which reads at equilibrium,

[TABLE]

where $\rho_{0}=\sqrt{N}$ is not a dynamical variable and $T$ is the absolute temperature.

To make the pair distribution function vanish within the hard sphere range we introduce into Eq. (4) a Lagrange multiplier function $\lambda_{\bf{k}}$, which is a Fourier transform of a yet unknown function $\lambda(\bf{r})$, that vanishes outside the hard sphere interaction range. The last requirement reflects the fact that the pair distribution function is constrained only within that range. The equation now reads,

[TABLE]

Note that by making $\lambda$ depend only on the absolute value of $\bf{k}$ ,we assume implicitly a rotation invariant phase. Also, the subscript $S$ on the density distribution function denotes steady state, which is not necessarily the equilibrium state. The reason is that at densities corresponding to the equilibrium solid phase we constrain our system to have a spherically symmetric structure factor. Thus the state described by the solution of Eq. (5) and the constraint on the pair distribution function is a steady state metastable state rather than an equilibrium state. Multiplication of Eq. (5) by $\frac{1}{2}\rho_{k}\rho_{-k}$ and functional integration by parts yields the structure factor in terms of the unknown Lagrange multiplier $\lambda_{\bf{k}}$, namely

[TABLE]

where $\lambda_{\bf{k}}$ has to obey the two following conditions,

[TABLE]

and

[TABLE]

It turns out that these two conditions are in fact the hard sphere PY equation. In classical liquid theory a quantity termed direct correlation function is used extensively and is traditionally denoted by $c(r)$. In our language, the Lagrange multiplier function $\lambda(r)$ is simply $-\bar{\rho}c(r)$ (see Appendix A for more details on the notation used here).

Although Eq. (5) includes terms tri-linear in the operators $\rho_{k}$ and $\frac{\partial}{\partial\rho_{k}}$ , the solution Eq. (6) is exact , because of the specific form of those terms. The next step beyond PY can be affected by constraining the triplet distribution function to vanish whenever any pair of the triplet is closer than the hard sphere interaction range. The equation for the distribution in that case will include tri-linear terms of a different nature, which will prevent an exact solution for the structure factor in terms of the two point Lagrange multiplier function. The equation in that case belongs, however, to a wide family of stochastic nonlinear systems, described by a functional Fokker-Planck equation that has been treated successfully by the Self-Consistent Expansion (SCE) Schwartz2008 ; Schwartz1992 ; Schwartz1998 ; Katzav1999 ; Katzav2003 ; Katzav2007 ; Edwards2002 . Thus, the PY equation is not the last word, as it can be systematically improved to include higher order correlation functions. The interesting thing is that in spite of its simplicity, the equation of state the PY approximation produces is in very good agreement with simulations for liquids Hansen2006 . In fact, since the PY is only an approximation it produces two very good but different equations of state, depending on the route of derivation. When a proper weighted average of the two is constructed the really excellent Carnahan–Starling (CS) equation of state Carnahan1969 is obtained

[TABLE]

where $P_{\rm{CS}}$ is the pressure and the volume fraction $\eta$ is given (in three dimensions) by

[TABLE]

with similar expressions in other dimensions (see Appendix A). Recall that $R$ is the range of the hard sphere interaction, namely, the diameter, and not the radius of a single sphere.

Since the CS equation of state holds for volume fractions below crystallization, the fact that it holds also above crystallization seems to be irrelevant. The reason is that the PY approximation assumes invariance under rotation and the emergence of a crystalline structure is just due to the fact that the free energy associated with the solid is lower than the one associated with the rotation invariant phase. It is interesting to note, however, that the hard sphere system possesses a metastable super dense rotation invariant phase which is disordered. Actually, the pressure in that phase is well described by the CS equation of state up to $\eta=0.57$ for monodisperse hard spheres, while for polydisperse hard spheres the related BMCSL approach Boublik1970 ; Mansoori1971 extends way beyond that Berthier2009 ; Berthier2016 . This super dense branch should, however, have terminated at random close packing, where the pressure is expected to diverge Parisi2010 . Furthermore, that branch as predicted by PY continues into nonphysical volume fractions, even above the crystalline close packing. The main trouble with PY is therefore that there seems to be no intrinsic indication within the PY theory that something goes wrong at higher volume fractions. The first message of the present article is that, contrary to the above statements, an intrinsic indication for the failure of the theory at a certain density does exist in PY.

III The breakdown of PY at $\eta_{c}$

Consider the pair distribution function $g_{2}({\bf{r}})$ defined above in Eq. (2). $g_{2}({\bf{r}})$ is obtained, within the PY approximation, in the following way. The exact solution in odd $d$ dimensions provides the so-called direct correlation function $c(r)\equiv-\lambda(r)/\bar{\rho}$, for $r<R$. As it happens, those are polynomials of degree $d$ in $r$ with coefficients that are functions of the volume fraction $\eta$. Since for $r>R$, the direct correlation vanishes, obtaining the corresponding Fourier transform $\lambda_{k}$’s is a straightforward analytic calculation. The last step to obtain $g_{2}(r)$ is to use Eq. (6) for the structure factor and finally use Eq. (8) to obtain $g_{2}(r)$ in the limit of infinite volume by numerical integration.

We begin with the one-dimensional case. In Fig. 1 we present the pair distribution function for three different volume fractions in one dimension. Of particular interest is the high volume fraction graph. The apparent peaks are related to the short range order in the system but it is clear enough that nothing spectacular happens as the peaks are broadened and reduced in height as a function of the distance.

We continue with Fig. 2, where we present the corresponding pair distribution function in three dimensions.

The low volume fraction graph shows no interesting features, but the intermediate and high density graphs show clearly short range order, which is manifested by the oscillations of the pair distribution function. The alert reader may have already detected a serious problem in the $\eta=0.65$ case. The pair distribution function, $g_{2}(r)$ as defined by Eq. (2), is by definition nonnegative, while in Fig. 2 the pair distribution function is negative in a certain region of $r$. Since we also compare our numerical integration with an exact representation of $g_{2}(r)$ Smith1970 ; Kelly2016 which is available in three dimensions up to $r\leq 5R$ this negativity cannot be attributed to an artifact of the approximate numerical integration needed to obtain the pair distribution function. More quantitatively, we can look for the lowest volume fraction for which a negative part appears. This is actually the point $r_{c}$ where both the function and its derivative become zero simultaneously, i.e. $g_{2}(r_{c})=g^{\prime}_{2}(r_{c})=0$. Using the analytical result of Refs. Smith1970 ; Kelly2016 we find that the lowest volume fraction for which a negative part appears is $\eta_{c}(d=3)\simeq 0.612574..$. Similar regions of negative values of the colloid-colloid pair distribution function were obtained for colloids immersed in a fluid Henderson1986 .

The inevitable conclusion thus is that the PY approximation breaks down intrinsically at high enough volume fraction. Namely, in contrast to the reasoning based on the continuity of the equation of state across physically impossible densities that was discussed above, the PY calculation itself indicates that something must go wrong, by giving negative values to a function that is nonnegative definite. This is obviously good news since it sets an internal limit, $\eta_{c}$, on the applicability of the fluid equation of state (or rotationally invariant phase) at high volume fractions.

At this point it is natural to ask whether this $\eta_{c}$ has any physical meaning beyond being an intrinsic upper bound on the theory? The first thing to check is whether it happens in higher dimensions as well. We have obtained the pair distribution function for dimensions $3<d\leq 9$. For the odd dimensions we used the exact solution, with a numerical solution of the appropriate polynomial equation when required, as well as a numerical Fourier transform to obtain the pair distribution function - see Refs. Robles2004 ; Rohrmann2008 and the appendices for more details (in these dimensions there is not any direct analytical representation of the pair distribution function like the one available in $d=3$ in Refs. Smith1970 ; Kelly2016 ). It turns out that similar to three dimensions, for all odd dimensions in the range $5\leq d\leq 9$, the pair distribution function becomes negative at some range of $r/R$. and above some critical value of the volume fraction, $\eta_{c}(d)$. The results are summarized in Table 1 below, as well as graphically in Fig. 4 (red circles). The results of our numerical integration are supported by the exact 3D result Smith1970 ; Kelly2016 as well as intrinsically by comparison with improved approximate integration.

For even dimensions we use the method and results reported in previous work AddaBedia2008a ; AddaBedia2008b that provides the pair distribution function as a power series in the volume fraction $\eta$ as

[TABLE]

In practice, the expansion functions, $g_{2}^{(n)}(r)$ up to $n=13$ for $d=4,6$ and $8$ are available numerically from Ref. AddaBedia2008b . These series work very well for small volume fractions. However, in the current work we are interested in fairly high volume fractions and, in particular, in identifying the lowest volume fraction for which $g_{2}(r)$ develops a negative part. Note that generically $g_{2}(r)$ is a decreasing function, exhibiting oscillations that become more and more pronounced as the density rises. Based on this observation (and on the odd dimensional cases discussed above) the first negative part should appear at the first minimum of $g_{2}(r)$, which is obtained in the interval $1<r/R<2$. The technical difficulty we encounter is that the radius of convergence of the series Eq. (11) is not large, and scales as $2^{-d}$ as the dimension, $d$, grows (see Ref. AddaBedia2008b for a more complete discussion). In particular, for the densities that are of interest the series does not converge, and we need to use some method to re-sum it or analytically-continue it. One such popular method is the Padé approximation Bender1999 . We look at various Padé approximants of $g_{2}(r)$, which are composed of a polynomial of order $N$ in $\eta$ divided by a polynomial of order $M$ in $\eta$, of the general form

[TABLE]

The coefficients of the two polynomials are chosen in such a way, that the expansion of the ratio of the two polynomials, recovers the series $g_{2}(r)=\sum\limits_{n\geq 0}{\eta^{n}}g_{2}^{(n)}(r)$ up to order $(N+M)$ in $\eta$. Since there are in principle many ways to choose $N$ and $M$, we mapped all the options up to order $N+M=13$ and looked at the density for which the first zero crossing occurs. We considered only the Padé approximants for which no spurious pole appears inside the interval $1<r/R<2$, i.e., no spontaneous divergence appears where we expect no real physical divergence to occur (this is a well-known artifact of the Padé method). The various results in $d=4$ are presented in Fig. 3, and they lead to the following estimate of the largest volume fraction $\eta_{c}(d=4)\simeq 0.467\pm 0.013$. A similar analysis has been performed for $d=6$ and $8$ and the results are summarized in Table 1, as well as graphically in Fig. 4 (blue squares).

IV Discussion

The next question to consider is whether the appearance of $\eta_{c}(d)$ carries more physical meaning than the obvious one, namely, the inadequacy of the PY equation at high volume fraction. If indeed it carries any physical meaning it must be related to the termination of the super dense, rotation invariant steady-state, metastable phase. What is that termination about? The termination point should be the point where the pressure diverges. A natural candidate is the volume fraction of random close packing (RCP). That concept is a bit vague, however, as evident from the wide range of RCP volume fraction obtained in three dimensions (the values of RCP volume fraction in three dimensions, obtained by numerous authors Bernal1960 ; Robles1998 ; Scott1969 ; Finney1970 ; Berryman1983 ; Visscher1972 ; Tobochnik1988 are spread between $0.6$ Visscher1972 and $0.68$ Tobochnik1988 depending on the method of derivation). Due to this unsatisfactory situation, Torquato, Truskett and Debenedetti criticized the validity of the concept of RCP altogether and introduced instead the concept of MRJ Torquato2000 ; Torquato2010 , which is defined as the maximal random packing among all jammed configurations. The MRJ volume fraction is about $0.64$ well within the range of RCP volume fractions obtained by others. Skoge et al. Skoge2006 give the MRJ volume fraction in four, five and six dimensions and suggest also a fit for the MRJ density as a function of dimension for $3\leq d\leq 6$,

[TABLE]

where $c_{1}=-2.72$ and $c_{2}=2.56$.

More recently, Parisi and Zamponi Parisi2006 , using Replica theory, predicted a rich phase space that contains some transition densities describing the sphere packing in the super dense regime. In particular, they predicted the $J$-point, $\phi_{J}$, the point where the system jams, and the Glass Close Packing (GCP) density, $\phi_{\rm{GCP}}$, which is the point where the pressure diverges, and beyond which there is no disordered packing. Therefore, the termination point we are talking about should be compared to the GCP density, whose large $d$ dependence is given by Parisi2006

[TABLE]

We compare these results to the terminal volume fraction for which the PY pair distribution function becomes first negative $\eta_{c}(d)$ in Fig. 4. Note the GCP density based on Eq. (14) leaves the proportionality coefficient undetermined, and we fitted it to the data in Table 1 for the sake of comparison. We also tried to fit our results using the functional form given by Eq. (13), giving rise to the estimated values $\hat{c}_{1}=-5.397$ and $\hat{c}_{2}=3.385$. As can be seen, the fit based on Eq. (13) (solid line) and the fit based on Eq. (14) (dashed line) are almost indistinguishable within the dimensions under discussion.

The statement that the terminal volume fraction identified by the Percus-Yevick theory has a physical meaning seems to be justified by the observation that $\eta_{c}(d)$ is close in value and behaves as a function of dimension similarly to the theoretical predictions for MRJ/GCP. This observation is, of course, not a proof, but we can still present some supporting arguments. First, some support is gained by the one dimensional result where the PY pair distribution function remains positive for all volume fractions consistent with the fact that the real one-dimensional system of hard spheres never crosses to an RCP or MRJ state as the volume fraction is increased all the way up to $\eta=1$, where it crystallizes. Second, one can gain further insight into why $\eta_{c}(d)$ may be related to RCP and MRJ. Recall that PY provides a rotation invariant steady-state solution, which above the crystallization density describes a meta-stable disordered state. The current work shows that it ceases to exist at $\eta_{c}(d)$, which is consistent with the definition of RCP, namely the largest density for which a disordered phase exists. From another perspective, in a frozen system we expect each particle to sit in a cage formed by other particles. Vanishing $g_{2}(r)$ at a certain distance indicates the existence of such a cage. To understand that, recall that $g_{2}({\bf r})$ is the probability density to find a particle at $\bf r$. Therefore, if $g_{2}(r)$ vanishes at some distance $r_{c}$, no particle can cross that point, which means that the particles at $r<r_{c}$ are blocked within this cage. This argument is supported by simulation results Estrada2011 in which a freezing transition is accompanied by a strong decrease in the first minimum of the pair distribution function. Although these results are associated with freezing into an ordered solid, we expect this to be also the case when the frozen state is metastable and disordered.

The derivation of the Percus-Yevick equation presented above highlights its underlying assumptions. The main ingredients are the approximate treatment of the hard sphere interaction by constraining the pair distribution function and the assumption regarding the existence of a steady-state spherically invariant state. In that sense this approximation does not assume a fluid phase as such, and can, in principle, describe a glassy phase as well. Previous attempts to relate the RCP density and the glass transition include experimental effort using colloids Vladeeren1995 , and theoretical efforts that use the Rational Function Approximation (closely related to PY) Robles1998 ; Yuste1998 . However, the theoretical efforts are not based on an intrinsic criterion of the theory, but rather on an assumption of the existence of a glassy phase and a comparison to an equation of state of that phase developed by Torquato Torquato1995 ; Rintoul1996 . This approach has the advantage of determining the glass transition density $\eta_{g}$, which is not available in our approach. Basically, PY does not take into account the dynamics and hence it is clear that one would need further input if it is to be used to describe the glass transition. We already explained that PY cannot capture any nontrivial higher order correlation or response functions, and definitely not time-dependent quantities, such as $\chi_{4}$, which are extensively used in the glass community Toninelli2005 ; Berthier2011a ; Berthier2011b to characterize the glass transition. This observation motivates future effort to take into account the dynamics (e.g. in the Fokker-Planck Eq. (4) along the lines of Ref. Schwartz1999 ) as well as go beyond the present approximation (via constraining the triplet distribution function and using the self-consistent expansion), to gain more insight into the glassy state. Actually, this may be particularly useful since previtrification arising from ergodicity breaking Brito2009 is already present way below the glass transition in the supercooled phase.

V Summary and future perspective

To summarize, in this paper we show that unlike the common lore, the Percus-Yevick theory for monodisperse hard spheres provides an intrinsic indication for its limitation in the regime of high densities. More specifically, the positivity of the pair-correlation function $g_{2}({\bf r})$ in $d$ dimensions is violated at a certain volume fraction which we denote $\eta_{c}(d)$, and thus beyond it the PY theory is no longer consistent. It turns out that this phenomenon occurs in all dimensions in the range $3\leq d\leq 9$, suggesting that it should hold also beyond. A comparison of $\eta_{c}(d)$ to the various results and predictions for the RCP, GCP or MRJ volume fractions shows that they are close and behave similarly as a function of dimension. This observation suggests that the terminal volume fraction in the PY theory actually indicates the largest density for which a spherically invariant state can exist, even if the solid phase is already preferred at this point, and without being able to distinguish a fluid from a glass. In that sense, the PY theory can provide a simple approximation or indication for the existence of a termination density. Further investment using more elaborate methods such as replica theory could then reveal the full picture, which may contain glassy phases.

We hope this work will motivate other researchers to check this phenomenon in many other systems described by a Percus-Yevick theory. A few examples are - hard spheres in curved space Lishchuk2006 ; Lishchuk2009 ; Sausset2009 ; Tarjus2011 and hard spheres experiencing more complicated interactions such as sticky hard spheres Sengers2000 ; Baxter1968b or square-well fluids Mulero2008 . Other important direction are systems composed of polydisperse or mixture of hard spheres Hansen2006 ; Sengers2000 ; Mulero2008 ; Brouwers2006 ; Lebowitz1964 ; Gray2011 ; Yuste1998 , various charged hard sphere fluids Sengers2000 such as the hard sphere Yukawa fluid Chialvo1996 , ionic liquids Hansen2006 ; Mulero2008 , polarizable fluids Gray2011 , and even fluids of non-spherical shapes such as ellipsoids Letz1999 ; Khordad2009 , spherocylinders Martinez2003 , and chain-like molecules Sengers2000 ; Mulero2008 .

Appendix A The Percus-Yevick approximation in terms of the direct correlation function

The Percus-Yevick approximation derived in section II and summarized by Eqs. (7) and (8) is usually written in terms of the direct correlation function $c(r)$, which is related to the Lagrange multiplier introduced in Eq. (5) by $c({\mathbf{r}})=-\lambda({\mathbf{r}})/\bar{\rho}$. The direct correlation,$c({\mathbf{r}})$ is determined by the so-called Ornstein-Zernike equation,

[TABLE]

where $h({\mathbf{r}})=g_{2}({\mathbf{r}})-1$ is called the total correlation function. Note that this equation is equivalent to Eq. (8) by using Eqs. (6) and (7). Here and in the following, we take $r$ in units of $R$ , the diameter of the hypersphere to be unity, and thus In $d$ dimensions we have for the volume fraction $\eta$

[TABLE]

where $V_{d}(r)$ is the volume of a $d$-dimensional hypersphere of radius $r$, and $\Gamma(x)$ is the Euler gamma function. This equation generalizes Eq. (10) to any dimension.

In odd dimensions, a highly non trivial result Leutheusser1984 ; Rohrmann2007 ; Robles2007 ; Robles2004 is that the direct correlation function $c(r)$ within the PY approximation turns out to be a polynomial of degree $d$, namely

$c(r)=\theta\left(1-r/R\right)\sum\limits_{i=0}^{d}{c_{i}(\eta)(r/R)^{i}}$, where $\theta(x)$ is the Heaviside function. Therefore, obtaining the corresponding Fourier components $\tilde{c}(k)$ is a straightforward analytical calculation

[TABLE]

where $J_{\nu}(x)$ is the Bessel function order $\nu$. From this the structure factor is obtained via $S(k)=1/\left({1-\bar{\rho}\tilde{c}(k)}\right)$, and the radial distribution function $g_{2}(r)$ can be obtained using Eq. (8). Note that in the dimensions discussed in the appendices, namely for $d>3$, there is no direct analytical representation of $g_{2}(r)$ as the one available in three dimensions Smith1970 ; Kelly2016 , and therefore there is no alternative to performing a numerical inverse Fourier transform.

Appendix B The pair distribution functions in five, seven and nine dimensions

In this part we show the function $g_{2}(r)$ in five (Figs. B1 and B2), seven (Figs. B3 and B4) and nine (Figs. B5 and B6) dimensions similarly to what has been presented in the text for one and three dimensions. As can be seen in all cases there is a critical volume fraction $\eta_{c}(d)$ at which $g_{2}(r)$ starts to develop a negative part, which marks the termination density of applicability of the PY theory.