On the jamming phase diagram for frictionless hard-sphere packings

Abstract

We computer-generated monodisperse and polydisperse frictionless hard-sphere packings of 10⁴ particles with log-normal particle diameter distributions in a wide range of packing densities φ (for monodisperse packings φ = 0.46–0.72). We equilibrated these packings and searched for their inherent structures, which for hard spheres we refer to as closest jammed configurations. We found that the closest jamming densities φ_J for equilibrated packings with initial densities φ ≤ 0.52 are located near the random close packing limit φ_RCP; the available phase space is dominated by basins of attraction that we associate with liquid. φ_RCP depends on the polydispersity and is ∼0.64 for monodisperse packings. For φ > 0.52, φ_J increases with φ; the available phase space is dominated by basins of attraction that we associate with glass. When φ reaches the ideal glass transition density φ_g, φ_J reaches the ideal glass density (the glass close packing limit) φ_GCP, so that the available phase space is dominated at φ_g by the basin of attraction of the ideal glass. For packings with sphere diameter standard deviation σ = 0.1, φ_GCP ≈ 0.655 and φ_g ≈ 0.59. For monodisperse and slightly polydisperse packings, crystallization is superimposed on these processes: it starts at the melting transition density φ_m and ends at the crystallization offset density φ_off. For monodisperse packings, φ_m ≈ 0.54 and φ_off ≈ 0.61. We verified that the results for polydisperse packings are independent of the generation protocol for φ ≤ φ_g.

---

I. Introduction

Frictionless hard-sphere packings represent a useful model for atomic systems, liquids, glasses, and crystals,¹ aside from being a system directly utilized in materials science and chemical engineering.^2,3 This simple yet powerful model exhibits a range of diverse phenomena, including melting and freezing transitions,^1,4–9 the ideal glass transition,^1,7,10–13 the ideal glass or the glass close packing (GCP) limit,^1,14 as well as the random close packing (RCP) limit.^1,14–16

There are several attempts to merge the multitude of these effects into a single picture.^1,10 It is a difficult task, as significant debate on some of the concepts above is underway. The first big challenge is the definition and determination of the RCP limit. For monodisperse particles, there exist at least three estimates for the RCP limit, with distinct densities φ: (i) φ = 0.634–0.636;^2,16–18 (ii) φ ≈ 0.64;^15,19–21 and (iii) φ ≈ 0.65.^22–29 In our previous studies,^14,30 we suggested that φ ≈ 0.64 and φ ≈ 0.65 refer to different phenomena and represent the RCP limit φ_RCP (in the sense of the J-point¹⁵) and a lower bound of the GCP limit φ_GCP,¹ respectively. It implies that random jammed packings can systematically be produced at any density in the range [φ_RCP, φ_GCP].^1,10,31 The definition and determination of the GCP limit and the corresponding ideal glass transition represent the second actively discussed topic,^1,10,32 stemming from research on glasses and colloids. It is debated whether the ideal glass exists and if so, what is the density of the ideal glass transition. It is unclear if there are multiple glassy states and what is the lowest glass transition density.

We believe that the RCP and GCP limits shall be studied together, as a part of the systematic investigation of the phase space structure for hard spheres. At each packing density φ, the phase space has areas corresponding to valid packing configurations. These areas comprise basins of attractions of jammed configurations.^33–35 Some of these basins dominate the available phase space. Thus, one of the characteristics of the phase space at a given φ is the jamming density of these dominant basins of attraction, φ_J. The main objective of this paper is to build a map from φ to φ_J for a wide range of initial densities φ.

With this intention to study the structure of the phase space for hard spheres from first principles, we computer-generated monodisperse and polydisperse frictionless hard-sphere packings of 10⁴ particles (cf.Fig. 1) over a wide range of densities φ (for monodisperse packings φ = 0.46–0.72). Polydisperse particles have log-normal diameter distribution with diameter standard deviation σ from 0.05 to 0.3 in steps of 0.05. Then, we equilibrated these packings to let packing configurations arrive at the basins of attraction of inherent structures^33–35 that dominate the phase space at given densities. Finally, we searched for the inherent structures of these equilibrated configurations. In this paper, we use for hard spheres the term “closest jammed configurations” instead of “inherent structures”. An inherent structure for an arbitrary configuration of hard spheres is a jammed configuration that is the closest one to the initial configuration, hence the term closest jammed configuration.¹⁴

Fig. 1 Closest jammed configuration at a density φ = 0.662 for a random packing of 10⁴ polydisperse spheres. The sphere radii distribution is log-normal and has a standard deviation σ = 0.3. The initial unjammed packing was generated with the force-biased algorithm at a density φ = 0.613.

The paper is structured as follows. Before we present any experimental results, we use Section II to start with definitions that are relevant for the subsequent discussion. These include closest jammed configuration, basin of attraction of a closest jammed configuration, bounding region, bounding surface, and others. We describe the methods that we use to generate packings, to conduct equilibration, and to search for the closest jammed configurations in Section III. Section IV contains the results of packing generation, subsequent equilibration, and searching for the closest jammed configurations. We discuss these results in the same section. Section V presents a summary and conclusions.

II. Definitions

In this section, we briefly provide definitions needed for the rest of the paper. A more elaborate discussion and precise mathematical definitions of most of them can be found in our previous paper.¹⁴

We rely on the phase space packing description introduced by Salsburg and Wood³⁶ and therefore use the terms “limiting polytope”, “hypersurface”, and “hypercylinder” from their paper. Particle velocities are not included in the phase space. Under jamming, we understand collective jamming in packings of frictionless particles.^37–40 A packing is called jammed if at least a subset of its particles is jammed; other particles are referred to as rattlers. In this paper, we do not exclude rattler particles from packings, when calculating packing densities.

We utilize the concept of inherent structures, initially introduced by Stillinger for packings of particles with soft potential.^33–35 Inherent structures for such systems are local potential energy minima in the phase space. The potential energy in hard-sphere packings is replaced by the maximum density that a packing can have at a given phase space point, taken with the minus sign. The maximum density is calculated by fixing particle coordinates and inflating particle radii until the first contact between particles occurs. Inherent structures for hard-sphere packings correspond to jammed configurations.¹⁴ To emphasize that we are investigating hard particles, not particles with soft potential, we use below in this paper the term “closest jammed configuration” instead of “inherent structure”.

Each potential energy minimum can be associated with a corresponding basin of attraction, i.e., the energy minimum of a given basin of attraction is the termination point of energy minimization—the steepest descent procedure—for any configuration in this basin of attraction. The bounding region of a given jammed configuration at a given density is by definition the intersection of this configuration's basin of attraction with the available phase space (when contact hypercylinders for the given density are excluded from the phase space). We define bounding surfaces as surfaces of bounding regions. A bounding region is closed if the bounding surface is fully formed by hypercylinder surfaces.

We define the GCP limit (φ_GCP) for sufficiently polydisperse packings as the highest possible jamming density of these packings.¹⁴ For monodisperse and slightly polydisperse packings, it is the highest jamming density that can be achieved if crystallization is artificially suppressed.¹⁴ If crystallization is allowed, the GCP limit is revealed by an entropy minimum.^22,30 We define the RCP limit φ_RCP as the upper bound of the J-segment.^15,41,42 Similarly, the lowest typical (LT) jamming density φ_LT is the lower bound of the J-segment.

We distinguish between typical basins of attraction and untypical ones. Basins of attraction with jamming densities in the range [φ_LT, φ_RCP] are typical by definition; the others are untypical. With increasing number of particles in the packings, φ_RCP is almost unchanged and φ_LT increases.¹⁵ Thus, we may estimate φ_RCP in the thermodynamic limit by the upper boundary of the J-segment for sufficiently large finite packings. In the present paper, we assume that the J-segment converges in the thermodynamic limit to a single value φ_RCP,¹⁵ though this question is still discussed.⁴¹ Let φ_max be the highest possible packing density for a given sphere diameter distribution: it is the crystalline density (∼0.74) for monodisperse packings and φ_GCP for sufficiently polydisperse packings. Then, untypical basins of attraction have jamming densities in the range [φ_L, φ_LT) ∪ (φ_RCP, φ_max], where φ_L is the lowest jamming density, which for monodisperse packings equals at least (the density of tunneled crystals).^34,43–46

Typical and untypical basins of attraction have just been defined for Poisson packings, i.e., when the entire phase space is available or, in other words, when initial packing densities are zero. Under typical closest jamming densities for non-zero initial packing densities we understand the closest jamming densities that will be almost always found for packings created at a given density using a given algorithm, if it starts generation at a Poisson configuration. If the jamming density of a bounding region lies in the range [φ_LT, φ_RCP], we refer to this region and all the configurations in this region as liquid. Earlier, we talked about the part of the phase space that is available at a given density. Parts of the available phase space may be completely separated with contact hypercylinders and may have different properties. Thus, it is important to talk about the part of the phase space that is achievable from a given type of configuration. If a bounding region covers almost the entire part of the phase space that is achievable from a given type of configuration (at a given density), we call this region dominant. Liquid, typical, and dominant basins of attraction coincide for zero initial packing density.

III. Methods

A. Packing generation

Particles in our polydisperse packings have log-normal radii distributions with standard deviations σ from 0.05 to 0.3 in steps of 0.05 (the particle mean diameter is normalized to unity). All packings are generated in a fully periodic cubic box and contain 10⁴ particles (cf.Fig. 1). Packings are created in a wide range of compression rates using the force-biased (FB) protocol.^47,48 This protocol is a modification of the Jodrey–Tory algorithm.^49,50 The FB algorithm starts from a random distribution of particle centers in a simulation box. Each particle is supplied with an inner diameter chosen to be proportional to the desired particle diameter and to make particles in the closest pair touch each other with their inner diameter shells. Alternatively, a single inner diameter ratio can be specified for the entire packing as the ratio of inner diameters to the desired particle diameters. Similarly, a packing is supplied with an outer diameter ratio, initially larger than unity. The initial outer diameter ratio is chosen to ensure that the total volume of the particles equals the box volume. Particles are also supplied with elastic potential of the third order by overlap distance,⁴⁸ which is cut-off at the outer particle shell. It is now possible to compute the forces between each pair of particles, as well as the net forces for each particle. The algorithm is iterative and each iteration proceeds as follows: (i) determine a net force for each particle; (ii) displace all the particles by distances proportional to the particles' net forces and in the direction of the net forces; (iii) decrease the outer diameter ratio according to a specified contraction rate; and (iv) update the inner diameter ratio so that the inner diameter shells for the pair of closest particles touch each other. Though the inner diameter ratio may decrease through the iterations, its value has an increasing trend. The algorithm terminates when the outer diameter ratio is equal to the inner diameter ratio. The lower the outer diameter ratio contraction rate, the denser is the final configuration. The source code used in this paper is available under the MIT free software license.⁵¹

B. Packing equilibration

Salsburg and Wood³⁶ derived an equation of state for hard spheres, p = 1 + 1/[(φ_CJ/φ)^1/d − 1], where p is the estimated reduced pressure, φ is the current packing density, φ_CJ is the (closest) jamming density for the polytope where the given packing configuration resides, and d is the dimensionality of the system.⁴ The derivation of Salsburg and Wood³⁶ assumes that the pressure is stationary and packings are in equilibrium. But all the packings produced in computer simulations or experiments are intrinsically out of equilibrium,^30,52 because the generation process is non-stationary by definition, especially for fast compressions. The pressure measured in the course of packing generation should therefore not be used for the estimation of φ_CJ or for the tracking of jamming; instead, packings should be preliminarily equilibrated, i.e., exposed to molecular dynamics simulation with zero compression rate until the pressure is stationary. Equilibration moves the packing to bounding regions that dominate the part of the phase space achievable from an initial configuration.

We equilibrate the packings by performing sets of 2 × 10⁷ collisions with zero compression rate in a loop until the relative difference of reduced pressures in the last two sets is less than 10⁻⁴, so the pressure can be regarded as stationary. More precisely, to measure the pressure during 2 × 10⁷ collisions, we average pressures for 100 sub-sets of 2 × 10⁵ collisions, which amounts to 20 collisions in a sub-set per particle. We use our own implementation⁵¹ of the Lubachevsky–Stillinger (LS) packing generation algorithm^53,54 to carry out the equilibration.

C. Searching for the closest jammed configurations

To search for the closest jammed configurations, we do not follow the definition of these configurations through the steepest descent energy minimization, but modify the LS algorithm instead. We run the LS algorithm with a high compression rate of 10, until the non-equilibrium reduced pressure reaches a conventional high value of 10¹², then decrease the compression rate by a factor of two and run the LS algorithm again, until the pressure is high enough again (10¹²). We repeat this procedure until the compression rate is ≤10⁻⁴. The Boltzmann constant and masses of all particles are set to unity; the temperature is set to 0.2. Fast compressions at the beginning of the search make the initial bounding region collapse as much as possible and at the same time retain the configuration point in this bounding region. Slow compressions at the end of the search allow arriving at a truly jammed configuration. The details and validation of this algorithm are provided in our previous paper.¹⁴ There we ensured that the distribution of closest jamming densities for packings before equilibration is independent of a particular set of algorithm parameters in a wide range of the latter. We did not repeat this validation for equilibrated packings, as packings before equilibration in the previous paper covered the entire range of polydispersities as well as the entire range of initial and final packing densities considered in the present paper. Searching for the closest jammed configurations after equilibration produces dominant jamming densities.

IV. Results and discussion

In this section, we present the results of the packing equilibration and of searching for the closest jammed configurations of the equilibrated packings. We give an overview of the data in Subsection A (Data overview), analyze the data in Subsection B (Data analysis), test our conclusions for independence from the packing generation protocol in Subsection C (Protocol independence), and finally provide a schematic diagram with the phase space structure in Subsection D (Schematic phase space structure). At the end of the section, in Subsection E (Applicability of liquid equations of state), we check if it is possible to recover the properties of the phase space through comparison of the reduced pressure to predictions from liquid equations of state. To ease the reading of this section, we provide with Table 1 an overview of the symbols used below. Some of them have already been introduced, others will be defined later.

Table 1 Important symbols used in the text

Symbol	Brief description	Key figures and tables	Values for σ = 0
σ	Standard deviation of the log-normal particle radii distributions
γ	Compression rate for initial packing generation	X-axis in Fig. 2a and 7
φ	Initial packing density after force-biased generation	Y-axis in Fig. 2a, X-axis in Fig. 2b–d
φ CJ	Closest jamming density of a packing	Y-axis in Fig. 2b
φ J	Dominant jamming density of a packing	Y-axis in Fig. 2d
φ LT	Lowest typical jamming density
φ RCP	Random-close packing limit (J-point^15)	Left sides of Fig. 2b–d	∼0.64
φ GCP	Glass close packing limit	Right sides of Fig. 2b–d; Table 2	∼0.65
φ HCP	Crystalline packing density for monodisperse packings (FCC or HCP crystals)		∼0.74
φ max	Highest packing density: φHCP for monodisperse packings, φGCP for sufficiently polydisperse packings		∼0.74
φ L	Lowest possible jamming density, at least 2/3×φHCP for monodisperse packings (density of tunnelled crystals^43)		∼0.49
φ m	Melting transition density (onset of crystallization)	Fig. 2c and d	∼0.545
φ off	Offset of crystallization	Fig. 2c and d	∼0.61
φ f	Freezing transition density	Fig. 4	∼0.494
φ g	Ideal glass transition density. φJ(φg) = φGCP	Fig. 2c and d; Table 2	∼0.585
φ MCT	Density at which available phase space becomes relatively disjoint. φJ(φ < φMCT) = φRCP	Fig. 2c and d	∼0.52

---

A. Data overview

The dependence of initial packing densities φ on the inverse compression rate γ⁻¹ for packings produced with the FB algorithm is shown in Fig. 2a. Even before equilibration, we searched for the closest jamming densities φ_CJ for these packings, as described in Section III C; Fig. 2b depicts for these packings the closest jamming densities φ_CJvs. initial packing densities φ. Then, we equilibrated the initial packings in Fig. 2a and calculated closest jamming density estimates φ_JE from the equation of state by Salsburg and Wood, as described in Section III B. Closest jamming density estimates after equilibration φ_JEvs. initial packing densities φ are shown in Fig. 2c. Finally, we searched for the actual closest jammed configurations of the equilibrated packings (dominant jammed configurations). Dominant jamming densities φ_Jvs. initial packing densities φ are shown in Fig. 2d. The same four plots, built vs. γ⁻¹, can be found in the Appendix.

Fig. 2 (a) Initial packing density φ vs. inverse compression rate γ⁻¹. (b) Closest jamming density before equilibration φ_CJvs. initial packing density φ. (c) Closest jamming density estimate after equilibration φ_JEvs. initial packing density φ. (d) Closest jamming density after equilibration φ_Jvs. initial packing density φ. All the packings were generated with the force-biased (FB) algorithm. Colours for the different standard deviations σ of the log-normal particle radii distributions are depicted in the legends.

We did not average the data in Fig. 2; each point in these figures corresponds to a single packing. To guide the eye, points have been connected by straight lines. Averaging assumes that fluctuations in the data will disappear in the thermodynamic limit. This question is still unresolved.^15,41 If fluctuations disappear for infinite packings, averaging also assumes that the noise stemming from finite-size effects is symmetrical around the true value. This is not the case for the closest jamming densities of Poisson packings: the J-segment decreases with the increase of the number of particles, but its upper boundary φ_RCP is almost unchanged, only the lower boundary φ_LT is moving upwards.¹⁵ Thus, averaging would not produce meaningful results; taking the upper boundary of the J-segment instead will give a better estimate of φ_RCP in the thermodynamic limit. Additionally, averaging would remove the information about the exact boundaries of jamming intervals for finite packings.

The RCP densities φ_RCP are included as horizontal lines on the left side of the plots in Fig. 2. Similarly, the GCP densities φ_GCP are shown as horizontal lines on the right side of these plots. We determined the RCP limits as the upper boundaries of the horizontal parts of the plots in Fig. 2b, because the horizontal parts of the plots correspond to the closest jamming densities of Poisson packings.¹⁴ We determined the GCP limits for σ ≥ 0.05 by asymptotically expanding the plots φ(γ⁻¹) (Fig. 2a) and φ_CJ(γ⁻¹) (Fig. 7b in the Appendix) into infinite generation time or zero compression rate. φ_GCP for monodisperse packings was taken at the value of a structural transition φ ≈ 0.65.^14,22–30 More details on the determination of φ_RCP and φ_GCP can be found in our previous paper.¹⁴

B. Data analysis

The closest jamming densities of initial packings (without preliminary equilibration), Fig. 2b, belong to the interval [φ_LT, φ_RCP] for φ ≤ 0.61 in the case of monodisperse packings and for φ ≤ 0.63 for packings with σ = 0.3. We may say that for these initial densities liquid bounding regions are typical for the FB algorithm. When we talk about available phase space below in this section, we will always mean available and achievable from the configurations typical for the FB algorithm (which are liquid for φ < [0.61, 0.63], the exact value depending on σ). In the remainder of the paper, we focus on Fig. 2d.

Spontaneous crystallization, estimation of φ_m and φ_off. For monodisperse packings, bounding regions of configurations with crystalline inclusions dominate the available phase space in the density range φ ≈ 0.54–0.61 (σ = 0, red line in Fig. 2d). It is manifested by a sudden departure of the φ_J(φ) plot up to almost crystalline densities. This interval conforms to other studies of spontaneous crystallization and crystal nucleation.^5–8 These studies show that the crystal nucleation rate in monodisperse hard sphere packings is negligible for φ < 0.54, then grows rapidly at φ ≈ 0.54 and reaches a plateau, and then rapidly decreases at φ ≈ 0.61.⁵ The density ∼0.54 is interpreted as the melting transition density φ_m (dedicated studies produce the value φ_m ≈ 0.545 (ref. 1, 4, and 9) for monodisperse packings). The density φ ≈ 0.61 was earlier interpreted as the (ideal) glass transition density φ_g for monodisperse hard spheres,⁵⁵ which is usually detected by a rapid decrease of compressibility or self-diffusivity. Recent accurate studies show that the ideal glass transition occurs at φ_g ≈ 0.585, while at the same time crystallization is prevented after φ ≈ 0.61.^6–8 Thus, we refer to φ ≈ 0.61 as the density of the offset of crystallization φ_off. The interval of densities where bounding regions of configurations with crystalline inclusions dominate the phase space of packings with σ = 0.05 is φ ≈ 0.56–0.58, which can be seen in a sudden departure of the dominant jamming densities in the plot φ_J(φ) (cyan line in Fig. 2d). This range is also consistent with other results (φ = 0.56–0.59).^7,56 Our data report that crystallization becomes impossible for a certain σ in the range (0.05, 0.1). Nucleation studies determine that crystallization in packings with Gaussian diameter distribution becomes impossible for σ ≈ 0.07,^7,56 which conforms to our results (though we investigate log-normal sphere diameter distributions, for such a small polydispersity the two distributions almost coincide).

Change in the phase space structure at φ_MCT ≈ 0.52. For φ ≤ 0.52, the available phase space is strongly connected and dominated by liquid bounding regions for all the packing types. Thus, the closest jamming densities after equilibration are always obtained in the range of liquid jamming densities, φ_J ∈ [φ_LT, φ_RCP]. Starting from a certain density φ ≈ 0.52, the interval of dominant jamming densities moves upward, and starting from a slightly higher characteristic density the dominant jamming density φ_J is always higher than φ_RCP. It means that none of the liquid bounding regions participate in dominating the available phase space any longer. Packings with slight polydispersity σ = 0.05 clearly demonstrate that the onset of crystallization does not coincide in general with this characteristic density associated with the changes in the structure of the phase space. Changes in the phase space at φ ≈ 0.52, independent of the particle size distribution, are predicted under certain assumptions by the mode-coupling theory;^7,13 thus, we denote this transition density with φ_MCT. We assume that bounding regions become relatively disjoint at this density (i.e., the fraction of the “wormholes” in bounding surfaces becomes low).¹ We show later, in Subsection D, that the available phase space truly splits into disconnected portions and becomes non-ergodic only at the ideal glass transition density φ_g, which for monodisperse packings equals ∼0.585.

Glass close packing limit, estimation of φ_GCP. For packings with σ ≥ 0.1 and φ ≥ φ_MCT, the jamming density of the dominant bounding regions increases and reaches maxima at certain densities depending on σ. We associate these maxima in Fig. 2d with the GCP limits φ_GCP (and corresponding packing configurations with the ideal glass), as the GCP limits for packings where crystallization is impossible are by definition the highest jamming densities of these packings. For packings with σ = 0 and σ = 0.05 crystallization effects are superimposed on the φ_J(φ) plots in Fig. 2d. For monodisperse packings crystallization is no longer possible at φ_off ≈ 0.61. To the right of this density the dominant jamming density φ_J is ∼0.65, which also conforms to the density of the onset of crystalline inclusions in jammed packings and the density where the entropy is minimal φ_GCP ≈ 0.65.^14,22–30 For polydisperse packings with σ = 0.05, crystallization is superimposed on the φ_J(φ) plot in the range φ ≈ 0.56–0.58, but it does not cover the local maximum φ_J ≈ 0.652 at φ ≈ 0.586; so we attribute φ_J ≈ 0.652 to the GCP limit of particles with σ = 0.05. The estimates of the GCP limits compare very well with the GCP limit estimates from our previous paper,¹⁴ where we extrapolated the closest jamming densities of computer-generated packings with σ ≥ 0.05 to infinite generation time. The values from our previous paper¹⁴ are displayed as the lines of corresponding color to the right of Fig. 2d. The differences between the two estimates are <10⁻³ for all the packing types.

Ideal glass transition, estimation of φ_g. The initial packing density φ for which the dominant jamming density reaches its maximum is called the density of the ideal glass transition φ_g. It is usually measured through the jump in compressibility and divergent alpha-relaxation time.^1,4,7,10 We estimated φ_g from Fig. 2d in the following way: (i) first, we took the data points φ_J(φ) from Fig. 2d in the vicinity of the expected ideal glass transitions (i.e., in the interval ±0.02 around the global maxima of the φ_J(φ) plots); (ii) then, we selected local maxima from these data points (because local maxima represent upper boundaries of jamming intervals at each φ); and (iii) fitted these local maxima with third-order polynomials and found the positions of maxima for these polynomials. We consider these positions of maxima as the estimates for φ_g, which are depicted in Fig. 3. We fitted φ_g for σ ≥ 0.05 with the third-order polynomial and display it in Fig. 3 as well. The extrapolation of this polynomial to monodisperse packings gives a value φ_g = 0.585. We assume that if crystallization is artificially suppressed in monodisperse packings (e.g., by pinning a certain fraction of particles⁵⁷), the φ_J(φ) plot will look similar to those for polydisperse packings with σ ≥ 0.1, reaching its maximum value at φ = φ_g ≈ 0.585 with φ_J(φ_g) = φ_GCP ≈ 0.65. It would explain why crystalline inclusions appear in generated packings (prior to equilibration) only at the initial density φ = φ_GCP ≈ 0.65:^14,22–30 if φ_GCP ≈ 0.65 is the highest (jamming) density for monodisperse packings with suppressed crystallization, the only way to produce denser packings—for protocols that try to avoid crystallization as long as possible—is to introduce crystalline inclusions in the packing structure. We provide the values for φ_g and φ_GCP obtained so far in Table 2. Our data comply well with predictions from the mode-coupling theory,¹³ simulations and experimental observations of divergent relaxation times: φ_g = 0.582 is reported for packings of 10³ particles with uniform distribution having σ = 0.082;¹²φ_g = 0.59 is reported for a binary 50 : 50 mixture of 10³ particles with diameter ratio 1.4;¹¹φ_g = 0.585 and 0.586 are reported for packings of 2 × 10³ particles with Gaussian distributions of σ = 0.07 and 0.085, respectively.⁷ We display these values in Fig. 3 as well. For binary mixtures, σ is computed as the standard deviation of corresponding discrete probabilities of encountering a specific sphere radius.

Fig. 3 Ideal glass transition density φ_gvs. sphere radii standard deviation σ of the log-normal radii distributions (○). Values of φ_g are summarized in Table 2. Estimates for φ_g from other works are also presented: 10³ particles with uniform distribution (+);¹² binary 50 : 50 mixture of 10³ particles with diameter ratio 1.4 (×);¹¹ and 2 × 10³ particles with Gaussian distribution (□).⁷

Table 2 Ideal glass transition densities φ_g and corresponding ideal glass densities φ_GCP for packings of hard frictionless spheres with log-normal sphere diameter distributions having different diameter standard deviations σ (as indicated). φ_gvs. σ is plotted in Fig. 3

	σ
	0	0.05	0.1	0.15	0.2	0.25	0.3
φ g	0.585	0.586	0.59	0.595	0.602	0.61	0.622
φ GCP	0.65	0.6518	0.6549	0.6596	0.6645	0.6711	0.6779

---

We can now call the basins of attraction, the bounding regions, and all the configurations with φ_J ∈ (φ_RCP, φ_GCP] as glassy, for both monodisperse and polydisperse packings. This definition conforms to the fact that for σ ≥ 0.1 and φ > φ_MCT glassy bounding regions dominate the available phase space instead of liquid bounding regions. For monodisperse packings it also conforms to the fact that crystalline inclusions appear in jammed packings only for φ > φ_GCP.^{22,24–26,29,30} In our previous paper,¹⁴ we referred to the states in closed bounding regions as glassy (following Parisi and Zamponi¹ and Brambilla et al.¹¹) because for such states there is only one jammed configuration in the achievable phase space and the configurational entropy is zero. This definition is unsuitable in the light of the current results, as the states in closed liquid bounding regions (that dominate the phase space for φ ≤ φ_MCT) would be called glassy as well.

C. Protocol independence

To check if our results are protocol-dependent, we did the following: we (i) took the densest jammed packings obtained in Fig. 2 (the FCC crystal for monodisperse packings, the densest partially crystallized packing for σ = 0.05, and the ideal glass packings for σ ≥ 0.1); (ii) proportionally reduced the radii of the particles in these packings to produce unjammed packings in the entire range of densities starting from φ = 0.4 (we call these packings “diluted densest packings”, they represent a completely different packing generation protocol); and (iii) repeated the procedure for Fig. 2d, i.e., we equilibrated these packings and searched for the closest jammed configurations φ_J for the equilibrated packings. Fig. 4 depicts the plot φ_J(φ) for these diluted densest packings. The horizontal lines to the left and to the right of the figure represent the RCP and GCP limits from Fig. 2d.

Fig. 4 Closest jamming density after equilibration φ_Jvs. initial packing density φ. Colors for the different standard deviations σ of the log-normal particle radii distributions are depicted in the legend. All the initial packings were created by proportional decrease of the radii of the densest jammed packings of a corresponding σ (represented as the global maxima of the plots in Fig. 2d).

Monodisperse packings exhibit a well-known freezing transition at φ_f ≈ 0.5 (the value obtained in dedicated studies is ∼0.494).^1,4,9 Polydisperse packings with σ = 0.05 exhibit a freezing transition at φ_f ≈ 0.52. For σ ≤ 0.05 and φ > φ_f, the achievable phase space is dominated by the bounding regions of the densest configurations. For σ ≤ 0.05 and φ ≤ φ_f, the achievable phase space is dominated by liquid basins of attraction. The plots in Fig. 2d for σ ≤ 0.05 differ from Fig. 4 for φ ∈ [φ_f, φ_m] ∪ [φ_off, φ_max]. For these φ, the bounding regions for the densest configurations dominate the achievable phase space in Fig. 4 but are not achievable from the initial configurations in Fig. 2d.

For polydisperse packings with σ ≥ 0.1 and φ ≤ φ_g, the plots in Fig. 4 are qualitatively and quantitatively similar to the ones in Fig. 2d (see the RCP and GCP lines from Fig. 2d in Fig. 4). It means that the conclusions reached in Subsection B (Data analysis) are protocol-independent for these packings; the dominant jamming density does not depend on the initial packing configuration, which determines the achievable portion of the phase space. Thus, we may draw our conclusions in Subsection C for σ ≥ 0.1 and φ ≤ φ_g for the entire available phase space. For φ > φ_g, the plots differ: the dominant jamming density depends on the initial packing configuration. For diluted densest packings with σ ≥ 0.1 and φ > φ_g, the achievable phase space is dominated by the bounding regions of the ideal glass. It is sometimes argued that the available phase space is fully connected for φ ≤ φ_g and splits into disconnected portions at φ_g.¹

These results comply with the general agreement that the available phase space is ergodic at φ ≤ φ_g.⁷ For packings that allow crystallization and for which φ_off > φ_g, one may assume that ergodicity breaks at φ_off, when the φ_J(φ) plot becomes protocol-dependent. It was a general assumption in the colloidal literature as well, but recent careful studies of diffusion dynamics show that ergodicity for such packings breaks already at φ_g.⁷

D. Schematic phase space structure

Finally, we present a schematic image with the phase space structure, where we incorporate all the results obtained so far (Fig. 5). We assume that packings that allow crystallization (σ ≤ 0.07) are initially generated in liquid bounding regions for φ ∈ [φ_f, φ_off] (as it occurs in Fig. 2d); thus, melting transition and crystallization offset are exhibited for such packings, but not a freezing transition (as it occurs in Fig. 4). Symbols for different characteristic densities are displayed to the left and below the image. We provide the corresponding values for monodisperse packings to the right and above the image. All the characteristic densities have already been introduced; we only mention that under φ_g for σ < 0.05 we understand the ideal glass transition densities from extrapolating the φ_g(σ) plot to σ < 0.05, as done in Fig. 3 (which corresponds to divergent alpha-relaxation time). If crystallization is impossible, the protocol-dependent (non-ergodic) region for the φ_J(φ) plot is φ > φ_g (as depicted at the top of Fig. 5). For packings that allow crystallization and for which φ_off > φ_g, we follow Zaccarelli et al.⁷ and assume that ergodicity and protocol independence break at φ_g as well.

Fig. 5 Schematic jamming phase diagram for frictionless hard-sphere packings, depicting dominant jamming density φ_Jvs. initial packing density φ. Opaque areas represent jamming densities that do not dominate the phase space. Red dashed line and red opaque area refer to packings which allow crystallization. φ_L is the lowest jamming density, φ_RCP is the random close packing limit (the J-point), φ_GCP is the glass close packing limit, φ_max is the highest jamming density, φ_MCT is the density of transition from liquid to glass, φ_m is the melting transition density (onset of crystallization), φ_g is the ideal glass transition density for packings where crystallization is impossible, and φ_off is the density of the offset of crystallization.

We assume in Fig. 5 that in the thermodynamic limit the J-segment [φ_LT, φ_RCP] converges to a J-point (φ_RCP);¹⁵ we also assume that the dominant jamming density for any given initial density converges to a point in the thermodynamic limit as well. Opaque areas represent the jamming densities that do not dominate the phase space, but whose basins of attraction are available. Opaque areas are protocol-independent. The red dashed line and the red opaque area refer to packings that allow crystallization.

The procedure that we used to produce Fig. 2d can be extended to particles with soft potential. In that case, it is necessary to use the particle number density ρ instead of volume density φ. For each number density ρ and temperature T, one shall sample the phase space according to the usual equilibrium distribution of states in the canonical ensemble. For each sampled configuration, one shall perform the steepest descent in the potential energy landscape (infinitely fast quenching) and track the value of the obtained potential energy minimum U. Presumably, these values will be distributed in one or several narrow intervals around some dominant minima Ũ(ρ, T). There may be more than one of them at a given ρ and T if the phase space is not ergodic. Particle interactions are usually pair-wise; typical model potentials for interactions between soft particles are Gaussian core potential^58,59 and inverse-power potential.^60,61 An example for the rapid quenching of several equilibrated monodisperse two-dimensional packings with Gaussian core potential can be found in Stillinger and Weber.⁵⁹ The potential energy of interaction between two particles in inverse-power law systems is ε(D/r)ⁿ, where r is the distance between the two particles and ε, D, and n are parameters of the potential. Inverse-power law systems have a convenient scaling property. Specifically, it can be shown that all dimensionless excess thermodynamic properties of these systems depend only on a dimensionless reduced number density ρ_* ≡ ρD^d(ε/kT)^d/n, where d is the packing dimensionality and k is the Boltzmann constant.^61,62 For polydisperse packings, D depends on a given particle pair, but it can be split into a dimensionless pair-dependent part and a dimensional pair-agnostic part, of which the latter shall be used for calculating ρ_*. Thus, the plot Ũ(ρ, T) will be transformed into Ũ_*(ρ_*), where Ũ_* is a dominant dimensionless potential energy minimum. We assume that this plot will resemble Fig. 5, only with inverted Y-axis, as far as for hard spheres we defined Ũ ≡ −φ_J. Previous papers on particles with soft potentials typically focused on the distinction between liquid and crystalline phases.^58,60,61 We believe that finding dominant potential energy minima and building Ũ(ρ, T) plots may shed light on other important properties of such systems, including the ideal glass transition and the distinction between liquid and glass phases.

E. Applicability of liquid equations of state

At the end of this paper, we investigate for which density ranges liquid equations of state are applicable to the hard-sphere packings under study. We compared reduced pressures p(φ) in equilibrated packings (i.e., before densification) from Fig. 4 (diluted densest packings) with values analytically predicted by equations of state for hard spheres. There are many liquid equations of state for polydisperse packings.⁶⁴ We confirm that eqn (9)–(13) from Ogarko and Luding⁶⁴ produce very similar results, as well as eqn (6) from Mansoori et al.⁶⁵ We use the simplest of these equations of state, the one of Boublik–Carnahan–Starling–Mansoori (eqn (4) in Boublik⁶³ or eqn (9) in Ogarko and Luding⁶⁴): , where and ; 〈rⁱ〉 is the ith raw moment of the distribution of particle radii r. We calculated relative differences δ between experimental and theoretically predicted reduced pressures, δ = |p(φ) − p_cs(φ)|/p_cs(φ), and present in Fig. 6 the δ(φ) plots for different σ.

Fig. 6 Relative difference δ between experimental and theoretically predicted reduced pressures vs. initial packing density φ. δ = |p(φ) − p_cs(φ)|/p_cs(φ), p_cs(φ) is the reduced pressure from the Carnahan–Starling equation of state for polydisperse hard spheres.^63,64 Experimental pressures p(φ) were measured for equilibrated packings from Subsection IV C, cf.Fig. 4. Colors for the different standard deviations σ of the log-normal particle radii distributions are depicted in the legend.

Naturally, the reduced pressure in all the equilibrated packings at low densities is very close to the theoretical prediction; pressures differ by no more than 1% (δ ≤ 0.01). Pressure in packings that allow crystallization (σ < 0.1) starts to deviate from theoretical predictions at φ_f. Later on, we will discuss only packings with σ ≥ 0.1, because crystallization effects are not superimposed on their δ(φ) plots. Surprisingly, theoretical predictions for p(φ) for these packings remain valid with the same high accuracy δ ≤ 0.01 even for φ > φ_MCT, until φ reaches certain values depending on σ (Fig. 6). More precisely, δ is 0.01 at φ ≈ 0.581, 0.587, 0.591, 0.596, and 0.602 for σ = 0.1, 0.15, 0.2, 0.25, and 0.3, respectively. Selecting a different threshold than δ = 0.01 to consider deviations in reduced pressures as high leads to slightly different results. For example, δ is 0.06 at φ ≈ 0.588, 0.596, 0.602, 0.612, 0.622, for σ = 0.1, 0.15, 0.2, 0.25, and 0.3, respectively (Fig. 6). These values of φ are very close to φ_g from Table 2 for the corresponding σ. We confirm that the reduced pressure for equilibrated packings from Fig. 2c with σ ≥ 0.1 exposes behavior that is qualitatively and quantitatively similar to the one in Fig. 6 (data not shown).

V. Summary and conclusions

We computer-generated monodisperse and polydisperse frictionless hard-sphere packings of 10⁴ particles with log-normal particle diameter distributions in a wide range of densities φ (for monodisperse packings φ = 0.46–0.72). Then we equilibrated these packings and searched for their closest jammed configurations (inherent structures of hard spheres).

We found that the available phase space is dominated at φ ≤ 0.52 by liquid bounding regions with jamming densities φ_J ∈ [φ_LT, φ_RCP] (φ_LT ≈ 0.635 and φ_RCP ≈ 0.64 for monodisperse packings). At φ ≈ 0.52, independent of the particle radii distribution, the structure of the available phase space changes: bounding regions become relatively disjoint (i.e., the fraction of the “wormholes” in bounding surfaces becomes low). The value for this density, also independent of the particle radii distribution, is predicted under certain assumptions by the mode-coupling theory. Thus, we refer to this transition density as φ_MCT ≈ 0.52.

For φ > φ_MCT, the dominant jamming densities φ_J increase with φ and the available phase space is dominated by basins of attraction that we call glassy. When φ reaches the ideal glass transition density φ_g, φ_J reaches the ideal glass density (the glass close packing limit) φ_GCP, so that the available phase space is dominated at φ_g by the basin of attraction of the ideal glass. φ_g and φ_GCP depend on the particle size distribution. For packings with sphere diameter standard deviation σ = 0.1, φ_GCP ≈ 0.655 and φ_g ≈ 0.59. For monodisperse and slightly polydisperse packings, crystallization is superimposed on these processes: it starts at the melting transition density φ_m and ends at the crystallization offset density φ_off. For monodisperse packings, φ_m ≈ 0.54 and φ_off ≈ 0.61. If we extrapolate the ideal glass transition densities φ_g for polydisperse packings to σ = 0 (monodisperse packings), we obtain φ_g ≈ 0.585, in agreement with experiments and simulations on divergent alpha-relaxation time and the jump in compressibility. We verified that the results for packings with σ ≥ 0.1 and φ ≤ φ_g are independent of a packing generation protocol.

We also discovered that the reduced pressure in equilibrated packings complies with liquid equations of state for hard polydisperse spheres for φ > φ_MCT. Thus, the comparison with liquid equations of state is not sensitive enough to reveal the changes in the structure of the available phase space at φ_MCT, which are detected by dominant jamming densities.

VI. Appendix

A. Densities vs. inverse compression rate

In Fig. 7, we provide the same plots as in Fig. 2, but built vs. the inverse compression rate γ⁻¹, not vs. initial packing densities φ (as done in Fig. 2b–d). Fig. 7b shows the φ_CJ(γ⁻¹) plots that we used to estimate the GCP limits by asymptotic expansion to infinite generation time or zero compression rate. The plots in Fig. 7b have a structure predicted by Parisi and Zamponi¹ (Fig. 2a in that paper).

Fig. 7 (a) Initial packing densities φ vs. the inverse compression rate γ⁻¹. (b) Closest jamming densities before equilibration φ_CJvs. the inverse compression rate γ⁻¹. (c) Closest jamming density estimates after equilibration φ_JEvs. the inverse compression rate γ⁻¹. (d) Closest jamming densities after equilibration φ_Jvs. the inverse compression rate γ⁻¹. All the packings were generated with the force-biased (FB) algorithm. Colours for the different standard deviations σ of the log-normal particle radii distributions are depicted in the legend.

References

1. G. Parisi and F. Zamponi, Rev. Mod. Phys., 2010, 82, 789–845.
2. S. Khirevich, A. Höltzel and U. Tallarek, Comm. Comput. Phys., 2013, 13, 801–822.
3. S. Khirevich, A. Höltzel, A. Daneyko, A. Seidel-Morgenstern and U. Tallarek, J. Chromatogr. A, 2011, 1218, 6489–6497.
4. M. Skoge, A. Donev, F. H. Stillinger and S. Torquato, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2006, 74, 041127.
5. L. Filion, M. Hermes, R. Ni and M. Dijkstra, J. Chem. Phys., 2010, 133, 4115.
6. E. Sanz, C. Valeriani, E. Zaccarelli, W. C. K. Poon, P. N. Pusey and M. E. Cates, Phys. Rev. Lett., 2011, 106, 215701.
7. E. Zaccarelli, C. Valeriani, E. Sanz, W. C. K. Poon, M. E. Cates and P. N. Pusey, Phys. Rev. Lett., 2009, 103, 135704.
8. C. Valeriani, E. Sanz, E. Zaccarelli, W. C. K. Poon, M. E. Cates and P. N. Pusey, J. Phys.: Condens. Matter, 2011, 23, 194117.
9. W. G. Hoover and F. H. Ree, J. Chem. Phys., 1968, 49, 3609–3617.
10. L. Berthier and T. A. Witten, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2009, 80, 021502.
11. G. Brambilla, D. El Masri, M. Pierno, L. Berthier, L. Cipelletti, G. Petekidis and A. B. Schofield, Phys. Rev. Lett., 2009, 102, 085703.
12. G. Pérez-Ángel, L. E. Sánchez-Díaz, P. E. Ramírez-González, R. Juárez-Maldonado, A. Vizcarra-Rendón and M. Medina-Noyola, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2011, 83, 060501.
13. T. Voigtmann, A. M. Puertas and M. Fuchs, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2004, 70, 061506.
14. V. Baranau and U. Tallarek, Soft Matter, 2014, 10, 3826–3841.
15. C. S. O'Hern, L. E. Silbert, A. J. Liu and S. R. Nagel, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2003, 68, 011306.
16. C. Song, P. Wang and H. A. Makse, Nature, 2008, 453, 629–632.
17. A. Zinchenko, J. Comput. Phys., 1994, 114, 298–307.
18. G. D. Scott and D. M. Kilgour, J. Phys. D: Appl. Phys., 1969, 2, 863–866.
19. Y. Jin and H. A. Makse, Phys. A, 2010, 389, 5362–5379.
20. J. D. Bernal and J. Mason, Nature, 1960, 188, 910–911.
21. J. G. Berryman, Phys. Rev. A, 1983, 27, 1053–1061.
22. A. V. Anikeenko, N. N. Medvedev and T. Aste, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2008, 77, 031101.
23. A. V. Anikeenko and N. N. Medvedev, Phys. Rev. Lett., 2007, 98, 235504.
24. S. C. Kapfer, W. Mickel, K. Mecke and G. E. Schröder-Turk, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2012, 85, 030301.
25. B. A. Klumov, S. A. Khrapak and G. E. Morfill, Phys. Rev. B: Condens. Matter Mater. Phys., 2011, 83, 184105.
26. M. Bargieł and E. M. Tory, Adv. Powder Technol., 2001, 12, 533–557.
27. T. Aste and T. Di Matteo, Eur. Phys. J. B, 2008, 64, 511–517.
28. T. Aste and A. Coniglio, Europhys. Lett., 2004, 67, 165–171.
29. K. Lochmann, A. Anikeenko, A. Elsner, N. Medvedev and D. Stoyan, Eur. Phys. J. B, 2006, 53, 67–76.
30. V. Baranau, D. Hlushkou, S. Khirevich and U. Tallarek, Soft Matter, 2013, 9, 3361–3372.
31. P. Chaudhuri, L. Berthier and S. Sastry, Phys. Rev. Lett., 2010, 104, 165701.
32. L. Berthier and G. Biroli, Rev. Mod. Phys., 2011, 83, 587–645.
33. F. H. Stillinger, Science, 1995, 267, 1935–1939.
34. S. Torquato and Y. Jiao, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2010, 82, 061302.
35. F. H. Stillinger, E. A. DiMarzio and R. L. Kornegay, J. Chem. Phys., 1964, 40, 1564–1576.
36. Z. W. Salsburg and W. W. Wood, J. Chem. Phys., 1962, 37, 798–804.
37. S. Torquato and F. H. Stillinger, Rev. Mod. Phys., 2010, 82, 2633–2672.
38. S. Torquato, T. M. Truskett and P. G. Debenedetti, Phys. Rev. Lett., 2000, 84, 2064–2067.
39. A. Donev, J. Appl. Phys., 2004, 95, 989–999.
40. A. Donev, S. Torquato, F. H. Stillinger and R. Connelly, J. Comput. Phys., 2004, 197, 139–166.
41. M. Pica Ciamarra, M. Nicodemi and A. Coniglio, Soft Matter, 2010, 6, 2871–2874.
42. R. Ni, M. A. C. Stuart and M. Dijkstra, Nat. Commun., 2013, 4, 2704.
43. S. Torquato and F. H. Stillinger, J. Appl. Phys., 2007, 102, 093511.
44. Y. Jiao, F. H. Stillinger and S. Torquato, J. Appl. Phys., 2011, 109, 013508.
45. M. Pica Ciamarra and A. Coniglio, Phys. Rev. Lett., 2008, 101, 128001.
46. M. Pica Ciamarra, P. Richard, M. Schröter and B. P. Tighe, Soft Matter, 2012, 8, 9731–9737.
47. J. Mościński, M. Bargieł, Z. A. Rycerz and P. W. M. Jacobs, Mol. Simul., 1989, 3, 201–212.
48. A. Bezrukov, M. Bargieł and D. Stoyan, Part. Part. Syst. Charact., 2002, 19, 111–118.
49. W. S. Jodrey and E. M. Tory, Phys. Rev. A, 1985, 32, 2347.
50. M. Bargieł and J. Mościński, Comput. Phys. Commun., 1991, 64, 183–192.
51. https://code.google.com/p/packing-generation, 2013.
52. M. C. Vargas and G. Pérez-Ángel, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2013, 87, 042313.
53. B. D. Lubachevsky and F. H. Stillinger, J. Stat. Phys., 1990, 60, 561–583.
54. B. D. Lubachevsky, J. Comput. Phys., 1991, 94, 255–283.
55. W. van Megen and S. M. Underwood, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 1994, 49, 4206–4220.
56. P. N. Pusey, E. Zaccarelli, C. Valeriani, E. Sanz, W. C. K. Poon and M. E. Cates, Philos. Trans. R. Soc., A, 2009, 367, 4993–5011.
57. R. L. Jack and L. Berthier, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2012, 85, 021120.
58. F. H. Stillinger, J. Chem. Phys., 1976, 65, 3968–3974.
59. F. H. Stillinger and T. A. Weber, Phys. Rev. A, 1982, 25, 978–989.
60. S. Prestipino, F. Saija and P. V. Giaquinta, J. Chem. Phys., 2005, 123, 144110.
61. C. N. Likos, B. M. Mladek, D. Gottwald and G. Kahl, J. Chem. Phys., 2007, 126, 224502.
62. J. D. Weeks, Phys. Rev. B: Condens. Matter Mater. Phys., 1981, 24, 1530–1535.
63. T. Boublík, J. Chem. Phys., 1970, 53, 471–472.
64. V. Ogarko and S. Luding, J. Chem. Phys., 2012, 136, 124508.
65. G. A. Mansoori, N. F. Carnahan, K. E. Starling and T. W. Leland, J. Chem. Phys., 1971, 54, 1523–1525.

---