# Exploring the jamming transition over a wide range of critical densities

Misaki Ozawa, Ludovic Berthier, Daniele Coslovich

► To cite this version:

Misaki Ozawa, Ludovic Berthier, Daniele Coslovich. Exploring the jamming transition over a wide range of critical densities. SciPost Physics, 2017, 3 (4), pp.027. 10.21468/SciPostPhys.3.4.027. hal-01685133

HAL Id: hal-01685133

https://hal.science/hal-01685133

Submitted on 31 May 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Exploring the jamming transition over a wide range of critical densities

Misaki Ozawa, Ludovic Berthier and Daniele Coslovich

L2C, Université de Montpellier, CNRS, Montpellier 34095, France

## Abstract

We numerically study the jamming transition of frictionless polydisperse spheres in three dimensions. We use an efficient thermalisation algorithm for the equilibrium hard sphere fluid and generate amorphous jammed packings over a range of critical jamming densities that is about three times broader than in previous studies. This allows us to reexamine a wide range of structural properties characterizing the jamming transition. Both isostaticity and the critical behavior of the pair correlation function hold over the entire range of jamming densities. At intermediate length scales, we find a weak, smooth increase of bond orientational order. By contrast, distorted icosahedral structures grow rapidly with increasing the volume fraction in both fluid and jammed states. Surprisingly, at large scale we observe that denser jammed states show stronger deviations from hyperuniformity, suggesting that the enhanced amorphous ordering inherited from the equilibrium fluid competes with, rather than enhances, hyperuniformity. Finally, finite size fluctuations of the critical jamming density are considerably suppressed in the denser jammed states, indicating an important change in the topography of the potential energy landscape. By considerably stretching the amplitude of the critical “J-line”, our work disentangles physical properties at the contact scale that are associated with jamming criticality, from those occurring at larger length scales, which have a different nature.

<figure><img src="image_2.png" /></figure>

Copyright M. Ozawa *et al.*

This work is licensed under the Creative Commons Attribution 4.0 International License.

Published by the SciPost Foundation.

Received 30-05-2017

Accepted 29-09-2017

Published 12-10-2017

doi:10.21468/SciPostPhys.3.4.027

<figure><img src="image_3.png" /></figure>

Check for
updates

## Contents

1 **Introduction** 2

2 **Numerical methods** 4

2.1 The model 4

2.2 Equilibration of very dense fluid states 4

2.3 Compression towards the jamming transition 5

3 **Extending the range of jamming densities** 5

4 **Structure at the contact scale: Isostaticity and pair correlation** 7

4.1 Isostaticity 7

4.2 Pair correlation function at contact 8
<table><tbody><tr><td><b>5</b></td><td><b>Structure at the local scale</b></td><td><b>9</b></td></tr><tr><td></td><td>5.1 Bond-orientational order</td><td>9</td></tr><tr><td></td><td>5.2 Locally favored structure</td><td>11</td></tr><tr><td></td><td>5.3 Spatial organization of rattlers</td><td>14</td></tr><tr><td><b>6</b></td><td><b>Structure at the large scale</b></td><td><b>15</b></td></tr><tr><td></td><td>6.1 Hyperuniformity</td><td>15</td></tr><tr><td></td><td>6.2 Finite-size fluctuations of the critical density of jamming</td><td>18</td></tr><tr><td><b>7</b></td><td><b>Discussion and conclusions</b></td><td><b>21</b></td></tr><tr><td></td><td><b>References</b></td><td><b>23</b></td></tr></tbody></table>

# 1 Introduction

Granular materials such as grains, pinballs, and large colloids flow when some external forces are applied. However, when the volume fraction of these particles is increased above a certain value, particle motion is no longer allowed and these systems become amorphous solids. This mechanical transition corresponds to a jamming transition because it occurs in the absence of thermal fluctuations. The critical volume fraction of this critical “J-point” is $\phi_J$, and the associated critical properties have been intensively studied in the last few decades [1–4]. In particular, structural properties of jammed states have been examined over a wide range of length scales from microscopic to macroscopic, providing a starting point for understanding the mechanical and rheological properties of jammed systems [3, 5–8].

In practice, the structural properties of jammed systems can be sorted out by the typical length scale or wave number $k$ that is being probed. First, the contact scale, corresponding to $\delta \to 0$, where $\delta$ is the typical gap between particles for $\phi < \phi_J$, or to the typical overlap between particles for $\phi > \phi_J$. The corresponding wave number is $k \sim 1/\delta \to \infty$. A remarkable property at the contact scale is isostaticity, which implies that the average number of contacts per particle, $Z$, becomes twice the number of spatial dimensions exactly at $\phi = \phi_J$. In addition, a critical power law behavior characterizes the pair distribution $g(r)$ near contact at $\phi_J$ [3, 5, 9, 10]. The isostatic nature of the system is at the core of theoretical descriptions of the jamming transition [3, 4, 11]. Second, one can consider distances at the neighbor scale, corresponding to $k \sim 2\pi/\sigma$, where $\sigma$ is the averaged particle diameter. The major goal at this scale is to properly characterize the local geometry of amorphous configurations [12–14], to possibly define appropriate order parameters for the jamming transition [15, 16]. The role of local crystalline and icosahedral order in monodisperse packings has been discussed before [2, 17–21]. In particular, suitable modifications of bond orientational order parameters have revealed local crystalline order in (non-isostatic) jammed packings [18, 19], while icosahedral order is of more limited extent [20, 21]. Finally, we can analyse the large scale corresponding to $k \to 0$, where fluctuations over the entire sample are considered. It has been reported that systems prepared at the jamming transition have unexpected density fluctuations at large scales, possibly corresponding to hyperuniform behavior. Physically, hyperuniformity implies that the amplitude of volume fraction fluctuations at wave number $k$ vanishes as $k \to 0$ [22, 23]. However, several recent works have reported numerical evidence that hyperuniformity is not strictly obeyed at the jamming transition [24–26].

The structural properties mentioned above have been studied at “the” jamming transition.
However, theoretical studies have established that $\phi_J$ is not uniquely determined but is in fact strongly protocol dependent [15,27–41]. This implies that the jamming critical point, J-point, actually corresponds to a line of critical points, thus forming a “J-line” [31]. A simple way of exploring the J-line is to perform compressions of hard sphere configurations using different compression rates. A more controlled method consists in first thermalising the hard sphere fluid at finite temperature at some volume fraction and then performing a rapid compression towards jamming [31, 36, 42]. By varying the volume fraction of the parent fluid, a finite range of jamming volume fractions can be explored, while easily keeping crystallisation under control.

Although the existence of a J-line is well accepted, much less is known about how structural properties of jammed systems may change along the J-line. Some previous studies reported that the structural properties at the contact and neighbor scales hardly change [30, 31, 33, 43] and thus one might conclude that the structural properties are essentially invariant along the J-line. Scaling theories of the jamming transition are based on this assumption [3, 4, 11]. On the contrary, some other studies pointed out that there are tiny but systematic structural changes [34, 36], in particular at the neighbor scale. To our knowledge, the evolution of the large scale structural properties along the J-line has not been studied in detail. Analyzing the evolution of structure along the J-line is a numerical challenge, because changing the value of $\phi_J$ requires changing the typical preparation timescale over orders of magnitude, and even large numerical efforts may yield relatively minute changes to the value of $\phi_J$. Therefore, previous studies have accessed finite, but quantitatively modest, variations of $\phi_J$ along the J-line.

In this paper, we develop a numerical strategy that allows us to stretch the extension of the J-line of frictionless hard spheres considerably. As a result, we can characterize the variation of structural properties with $\phi_J$ over an unprecedentedly broad range of volume fractions, $\phi_J \approx 0.65–0.70$. This range is at least 3 times wider than in any previous study of frictionless jammed packings [31, 40]. The decisive factor allowing the present analysis is the recent development of an efficient thermalisation algorithm for polydisperse hard sphere fluids. We have recently shown that this approach allows the equilibration of very dense fluid states [44], bypassing any alternative method by many orders of magnitude. Here, we use these deeply thermalised fluid configurations as starting point for rapid compressions towards jammed states.

We find that at the contact scale, both isostaticity and the same critical behavior of the radial distribution function hold over the entire J-line, thus corroborating previous results [31, 33]. To characterize the structure at the neighbor scale, we apply analysis tools appropriate to polydisperse packings, based on the detection of locally favored structures [45, 46]. Our key finding is the detection of distorted local icosahedral structures that become increasingly numerous as $\phi_J$ is increased, so that about 80 % of the particles form such structures in our densest packings. This suggests that this local structural motif might be a key geometric factor allowing the production of very dense jammed packings [20]. At the large scale, we find that a nearly hyperuniform behavior [47] for the smallest $\phi_J$, but deviations from hyperuniformity increase rapidly as $\phi_J$ is increased, suggesting that optimized packings display stronger volume fraction fluctuations at large scale. On the other hand, finite size effects of the mean value of $\phi_J$ are considerably suppressed when $\phi_J$ is large. Overall, our work disentangles contact scale properties that are deeply connected to the criticality of the jamming transition, to structural properties at larger scale which evolve significantly along the J-line, and have therefore a distinct physical origin.

This paper is organized as follows. In Sec. 2, we describe the details of the simulation methods. The resulting extended range of jamming densities is discussed in Sec. 3. We discuss the physics at contact scale (isostaticity, pair correlation) in Sec. 4, the physics at the local scale
(bond orientational order, locally favored structures, rattlers) in Sec. 5, and at the large scale (hyperuniformity, global finite-size effects) in Sec. 6. We conclude our paper and discuss our results in Sec. 7.

# 2 Numerical methods

## 2.1 The model

We employ the standard model of additive frictionless hard spheres in three dimensions [48]. The pair interaction is zero for non-overlapping particles, infinite otherwise. We use a continuous size polydispersity, where the particle diameter $\sigma$ is distributed according to $f(\sigma) = A\sigma^{-3}$, $\sigma \in [\sigma_{\min}, \sigma_{\max}]$, where A is a normalization constant. Following previous work [44], we use the size polydispersity $\Delta = \sqrt{\overline{\sigma^2} - \overline{\sigma^2}/\overline{\sigma}^2} = 23\%$, corresponding to $\sigma_{\min}/\sigma_{\max} = 0.4492$, where $\overline{\cdot \cdot \cdot} = \int d\sigma f(\sigma)(\cdots)$. We use $\overline{\sigma}$ as the unit of length. $f(\sigma)$ used in this study is shown in Fig. 6(a). We simulate systems composed of $N$ particles in a cubic cell with periodic boundary conditions. We mainly use the system sizes $N = 1000$ and $8000$, but we also perform selected simulations for $N = 150, 300, 600, 2000, 4000$ to systematically investigate the finite-size effects described in Sec. 6. The state of the hard sphere system is uniquely characterized by the volume fraction $\phi = \pi N \overline{\sigma^3}/(6V)$, where $V$ is the volume of the system. For the fluid state, we measure the reduced pressure $p = P/(\rho k_B T)$, where $\rho = N/V$, $k_B$ and $T$ are the number density, Boltzmann constant and temperature, respectively. We set $k_B$ and $T$ to unity. The pressure $P$ is calculated from the contact value of the pair correlation function properly scaled for a polydisperse system [49]. The fluid state has a finite $p$, whereas the jammed state corresponds to $p \to \infty$.

Note that the functional form of $f(\sigma)$ with $\Delta = 23\%$ is chosen so that the system is fairly robust against crystallization and fractionation at extremely high densities [50]. With smaller $\Delta$, the system easily crystallizes when using the efficient swap Monte Carlo method described below. On the other hand, as demonstrated in previous studies, large values of $\Delta$ or different forms of $f(\sigma)$ lead to fractionation at sufficiently high density [51,52]. By contrast, the model parameters employed in this work are optimized to avoid such instabilities and thus enable us to explore the J-line over an unprecedented range of packing fractions.

## 2.2 Equilibration of very dense fluid states

To obtain equilibrium fluid configurations, we perform Monte Carlo (MC) simulations which combine traditional translational particle displacements and non-local particle swaps [50, 53-59]. Translational displacements are drawn from a cube of linear side 0.115, and a trial displacement is accepted if it does not create an overlap between particles [48]. For a trial swap move, we randomly pick a pair $(i, j)$ of particles with $|\sigma_i - \sigma_j| < a$ (we choose $a = 0.2$) and attempt to exchange their diameters [50]. The swap is accepted if it does not create an overlap. We perform translational moves with probability 0.8, and swap moves with probability 0.2 [50]. We have previously established [44] that this swap Monte Carlo setting is extremely efficient to thermalise hard sphere fluid states up to very large volume fractions, $\phi_{\text{fluid}} \approx 0.655$. We have checked that all our equilibrium configurations are taken in the fluid state, carefully monitoring possible signs of demixing or crystallization using the same structural tools described below to also characterize jammed states. Therefore, the parent fluid configurations used to produce jammed states all belong to the (metastable) equilibrium fluid branch, and we do not artificially vary the jamming density by introducing partially crystallized or demixed states. Instead, we explore the equilibrium fluid branch over a broad range of densities.
## 2.3 Compression towards the jamming transition

Having prepared equilibrium fluid configurations of hard spheres, we use the non-equilibrium compression algorithm introduced in Refs. [60, 61] to reach the jamming transition. Briefly, this algorithm replaces the hard sphere potential with a soft repulsive harmonic potential, given by

$$v_{ij}(r_{ij}) = \frac{\epsilon}{2} [1 - (r_{ij}/\sigma_{ij})]^2 \theta(1 - r_{ij}/\sigma_{ij}), \quad (1)$$

where $\epsilon$ is the energy scale, $r_{ij}$ is the distance between the spheres i and j, $\sigma_{ij} = (\sigma_i + \sigma_j)/2$ is the distance of the two spheres at contact, and $\theta(x)$ is the Heaviside step function. We use $\epsilon$ as the unit of energy. The idea of the algorithm is to alternate instantaneous compression steps and energy minimization to iteratively converge to the jamming density.

During a compression step, we increase the volume fraction of the system $\phi$ by $\delta\phi = 5 \times 10^{-4}$. This compression introduces finite overlaps between particles, such that the potential energy $U = \sum_{i<j} v_{ij}(r_{ij})$ becomes finite. These overlaps are then eliminated by performing an energy minimization. To this end, we apply the conjugate gradient minimization method [62] for $U$.

If the system still has a finite $U$ after minimization, we decrease $\delta\phi$ by factor of 2 and use a series of decompression and compression steps until the overlaps are eliminated. This process is interpreted as pulling the system from trivial local potential energy minima. We stop the algorithm when $\delta\phi < 1.0 \times 10^{-6}$, and the resulting system is essentially a hard sphere jammed packing. However, we find that the algorithm explained so far [61] produces a fraction of slightly hypostatic packings. Thus, we perform an additional process introduced in Ref. [31] to get more accurate isostatic jammed packings. We simulate sequential compression and minimization with $\delta\phi = 1.0 \times 10^{-5}$ until $U/N > 1.0 \times 10^{-6}$. Then, we decompress the system with $\delta\phi = 1.0 \times 10^{-6}$ until $U/N < 10^{-16}$, which determine the jamming transition point [31].

Note that the jamming algorithm employed in this work is based on purely isotropic compressions and decompressions. Thus, stability against shear deformation is not guaranteed [63], but we expect that further optimization against shear would not change our conclusions.

## 3 Extending the range of jamming densities

We present the equation of state of the system in Fig. 1(a). Equilibrium fluid configurations at finite pressure $p$ are compressed using the non-equilibrium algorithm described in Sec. 2.3 towards jammed states at $p \to \infty$. Figure 1(a) demonstrates that the higher the volume fraction of the parent fluid, $\phi_{fluid}$, the higher the jamming transition volume fraction, $\phi_J$. Thus, enhanced thermalization is the key to extend the J-line.

In Fig. 1(b), we show $\phi_J$ as a function of $\phi_{fluid}$, varying $\phi_{fluid}$ over a very broad range. Below $\phi_{fluid} \lesssim 0.53$, the observed $\phi_J$ is almost independent of $\phi_{fluid}$, as suggested by our horizontal dashed line. Also, we show the result obtained for a compression starting from a Poisson distributed system of harmonic soft spheres with $\phi_{fluid} = 0.3$. We find that $\phi_J$ from the Poisson distribution takes the same value as from dilute hard sphere fluid configurations, $\phi_J \approx 0.655$. Thus, we confirm that the protocol dependence of $\phi_J$ is essentially absent at $\phi_{fluid} \lesssim 0.53$. This value for $\phi_J$ is consistent with recent numerical results for a similar continuously polydisperse system [65].

As expected, $\phi_J$ starts to depart from the plateau value $\phi_J \approx 0.655$ when $\phi_{fluid} \gtrsim 0.53$, and it then monotonically increases with increasing $\phi_{fluid}$. The largest value we obtain is $\phi_J \approx 0.7$. A qualitatively similar behavior was observed before in a binary hard sphere mixture [31, 36].
<figure><img src="image_4.png" /></figure>

Figure 1: (a) Equation of state of the system for two different system sizes. Equilibrium fluid states are shown as empty points. A modified version of the Carnahan-Starling equation of state for polydisperse systems is drawn as a solid line [64]. Jammed states (filled points) obtained by the compression algorithm are located at $p \to \infty$. The dashed arrows connect an equilibrium parent fluid state to the corresponding jammed state. (b) Evolution of the averaged jamming transition volume fraction $\phi_J$ obtained by the compression from an equilibrium parent fluid with volume fraction $\phi_{\text{fluid}}$. We also report $\phi_J$ starting from a Poisson distributed configuration with $\phi_{\text{fluid}} = 0.3$. The horizontal dashed line is $\phi_J$ generated from a Poisson distribution for $N = 8000$. The vertical bars correspond to the width of the J-line for this work (thick bar) which is about 3 times larger than previous works using a binary mixture [31] (thin bar).

The variation of the jamming density with the initial fluid density is reminiscent of the variation of the energy of inherent structures with initial temperature for glass-forming materials with a continuous pair interaction [36, 66, 67]. Note that the finite size effect of $\phi_J$ between $N = 1000$ and $8000$ is small at low $\phi_{\text{fluid}}$, and it essentially vanishes in the large $\phi_{\text{fluid}}$ regime, at least in this graphical representation. The finite size effect of $\phi_J$ will be more systematically examined in Sec. 6. We draw as vertical bars the widths of the J-line obtained in the present study (thick bar) and in a previous study of a binary mixture [31] (thin bar). The extension of the J-line achieved in this work is about 3 times larger than previous works [31, 40]. Thus, we have indeed succeeded in stretching the J-line considerably. In the next sections, we analyze the structural properties of the obtained packings along this stretched J-line.
<figure><img src="image_5.png" /><figcaption>Figure 2: (a): The averaged contact number Z as a function of $\phi_J$ remains close to the isostatic value Z = 6 along the entire J-line. (b): The fraction of rattlers, $N_r/N$, increases substantially as a function of $\phi_J$.</figcaption></figure>

# 4 Structure at the contact scale: Isostaticity and pair correlation

## 4.1 Isostaticity

We first examine structural properties at the contact scale. The relevant length scale for $\phi \le \phi_J$ is the gap between particles, $x-1$, where $x = r_{ij}/\sigma_{ij}$. A central quantity describing jammed states at this contact scale is the averaged contact number, Z. In practice, two particles, i and j, are considered in contact if $r_{ij} \le (1+a)\sigma_{ij}$, where $a = 1 \times 10^{-5}$ and $5 \times 10^{-6}$ for N = 1000 and 8000, respectively [33]. With this definition, Z is given by $Z = N_c/(N - N_r)$, where $N_r$ is the number of the rattlers, $N_c$ is the number of contact pairs among $N - N_r$ particles which make a contact network. We define a rattler as a particle whose contact number is smaller than $d+1$, where d is the spatial dimension. Note that a determination of rattlers has to be performed iteratively for a given configuration until all rattlers are removed. The isostatic packing is defined by the condition $Z = 2d$ (thus, $Z = 6$ in three dimensions), reflecting the fact that the system is mechanically marginally stable [68].

In Fig. 2(a), we show Z as a function of $\phi_J$ over the entire range of jamming densities that we were able to explore. We find that the isostatic condition, $Z = 6$, indeed holds over the entire J-line obtained by our protocols. In Fig. 2(b) we show the fraction of rattlers, $N_r/N$. It is around 6% at the lowest jamming density, $\phi_J \approx 0.655$, which is comparable with previous
<figure><img src="image_6.png" /></figure>

Figure 3: Radial distribution function $\rho g(x)$ near contact, where $\rho = N/V$ and $x = r_{ij}/\sigma_{ij}$. The closed and open symbols are obtained with and without rattlers, respectively. The black straight lines characterize a power law behavior, $g(x) \propto (x-1)^{-\gamma}$, with $\gamma = 0.5$ (the upper line) and $\gamma = 0.41269$ (the lower line), respectively. When rattlers are removed, a universal power law with the value predicted by the mean-field theory [70] is observed for all $\phi_J$ values. The gray shaded region marks the regime, $x - 1 \lesssim 10^{-3}$, where critical scaling with $\gamma = 0.41269$ holds.

results using energy minimization in three dimensions for binary [5,31] and polydisperse [69] systems. Interestingly, the fraction of rattlers increases steadily with $\phi_J$. Note that a slight increase of the fraction of rattlers along the J-line was also reported previously in a binary mixture [31,43]. This growth of the number of rattlers might appear counterintuitive at first sight, because rattlers tend to occupy a larger volume [43] and increasing their number should decrease the efficiency of the packing, in contradiction with the results shown in Fig. 2(b). The issue of the increase of the number of rattlers will be further discussed in Sec. 5.3.

## 4.2 Pair correlation function at contact

It is known that the radial distribution function shows a power law critical behavior near contact, $g(x) \propto (x-1)^{-\gamma}$, for isostatic packings. We show $\rho g(x)$, where $\rho = N/V$, as a function of the gap, $x-1$, in Fig. 3 for several $\phi_J$ values along the J-line. Note that we multiply $g(x)$ by $\rho$ to remove the trivial effect of the density change. The data for all $\phi_J$'s follow a power law with an exponent $\gamma = 0.5$, when we compute $g(x)$ from all particles in a given configuration, in particular including rattlers [71]. Instead, when the rattlers are removed and only the particles participating in the contact network are taken into account [72], the $g(x)$ for all $\phi_J$'s now follow a distinct power law which is compatible with the value $\gamma = 0.41269$ predicted by the mean-field theory of the jamming transition [70]. Thus, a careful treatment of the rattlers is essential to assess the critical behavior of the pair correlation function at contact [10,72].

From these observations, we are able to confirm that the critical behavior in $g(x)$ holds over the entire J-line and is therefore universal. This is not a trivial observation. For instance, it may have happened that this property only holds at the lowest end of the J-line, which is the only point where jamming criticality and mean-field predictions had been analyzed before [73]. Notice also that establishing this result was not easy. We found appreciable deviations from a power law behavior for smaller systems at low values of the argument $x-1$ (not shown). On
the other hand, we also find that by increasing $\phi_J$ the power law regime is entered for lower values of $x-1$, see Fig. 3. Therefore, observing a power law for large $\phi_J$ is difficult, as the critical regime becomes less easily observed for a given system size as $\phi_J$ is increased.

We expect that other critical behaviors such as a power law distribution of the contact forces are also observed over the entire J-line since these critical behaviors are all directly connected to one another [68, 70]. However, this direction is beyond the scope of this paper, since it would require a more important computational effort to prepare the packings exactly at the jamming transition [73].

# 5 Structure at the local scale

Our results so far confirm that the combined use of the swap Monte Carlo technique and poly-disperse spheres stretch substantially the J-line, while maintaining isostaticity, $Z = 6$. This already indicates that the system does not have crystalline order, which would produce $Z > 6$ (hyperstatic packings), as explicitly demonstrated in Ref. [33]. For monodisperse spheres, packing fractions comparable and even larger to the ones observed here can be attained. However, these packings are not isostatic and there is evidence that they accumulate local crystalline order beyond a threshold density [18]. Thus, the next question is: How do jammed configurations of polydisperse spheres attain such large densities while remaining amorphous? In this section, we address this issue by investigating the subtle evolution of the geometry of local packing at the neighbor scale, $k \sim 2\pi/\bar{\sigma}$. In particular, we confirm a mild increase of bond orientational order upon increasing density (Sec. 5.1) and reveal the emergence of distorted icosahedral local structures (Sec. 5.2). Above the onset density, these structural features are present in both fluid and jammed states.

## 5.1 Bond-orientational order

Bond-orientational order (BOO) parameters [74] have been used extensively to characterize the local structure of jammed packings [15, 21, 36]. The idea is to compute rotational invariants from a multiple expansion of the nearest neighbors distance distribution, and compare the results to the values observed for known crystal structures [74]. Inspection of the combined distribution of BOO parameters of different orders allows one to disentangle structures characterized by different local symmetries [75]. For a given particle $i$, the local BOO parameter is defined as

$$Q_{l,i} = \sqrt{\frac{4\pi}{2l+1} \sum_{m=-l}^{l} \left| \frac{1}{n_b(i)} \sum_{j=1}^{n_b(i)} Y_{l,m}(\mathbf{r}_{ij}) \right|^2}, \quad (2)$$

where $n_b(i)$ is the number of nearest neighbors of the $i$-th particle, $Y_{l,m}(\mathbf{r}_{ij})$ is the spherical harmonic of degree $l$ and order $m$, and $\mathbf{r}_{ij}$ is the vector distance between particle $i$ and $j$. The average BOO parameters, $Q_l = (1/N) \sum_i Q_{l,i}$, can then be used to characterize the local structure of a packing.

It is well known that the values of local BOO parameters in dense, disordered packings depend sensitively on the definition of the neighbors network surrounding each particle [19]. In disordered packings, neighbors are often defined as particles whose distance contributes to the first peak of the radial distribution function; alternatively, the nearest neighbors network is obtained from a Voronoi tessellation of the particles' coordinates [76]. In order to cure the sensitivity of BOO parameters to the details of neighbors network, Mickel et al. [19] have proposed to weight each bond entering the calculation of $Q_{l,i}$ by the area of the corresponding face of the surrounding Voronoi cell. The resulting weighted BOO parameters are then defined
<figure><img src="image_7.png" /><figcaption>Figure 4: (a, b): Bond orientational order parameter $Q_6^w$ for the fluid and jammed states as a function of their volume fraction (a), or the volume fraction of the parent fluid, $\phi_{fluid}$, (b). (c, d): Similar plots for the fraction of icosahedral structures characterized by the (0, 0, 12) Voronoi signature. The vertical lines are a guide for the eye.</figcaption></figure>

as

$$Q_{l,i}^w = \sqrt{\frac{4\pi}{2l+1} \sum_{m=-l}^{l} \left| \sum_{j=1}^{n_b(i)} \frac{A_{ij}}{A_i} Y_{l,m}(\mathbf{r}_{ij}) \right|^2}, \quad (3)$$

where $A_{ij}$ denotes the area of the face of Voronoi cell connecting particle i and j and $A_i = \sum_j A_{ij}$. A similar idea had been suggested long ago in Ref. [74], to cure the artifacts mentioned above. In the following, we compute the weighted BOO employing the neighbors network determined via a radical Voronoi tessellation [76] obtained using Voro++ [77].

A standard indicator of structural ordering is the average $Q_6$, which takes large values for ordered, close-packed structures ($Q_6 = Q_6^w = 0.575$ for the fcc cell and $Q_6 = Q_6^w = 0.663$ for the icosahedron). In Figs. 4(a, b), we show the weighted BOO parameter $Q_6^w$ for both equilibrium fluid (open symbols) and jammed (closed symbols) states. The data are shown as a function of the corresponding volume fraction of jammed states, $\phi_J$, and as a function of the volume fraction of the parent equilibrium fluid, ($\phi_{fluid}$), in panels (a) and (b), respectively. In both fluid and jammed states, $Q_6^w$ increases monotonically with increasing packing fraction, indicating a systematic, progressive local ordering along the J-line. In agreement with the results of Ref. [19], we found that the values of the standard $Q_6$, obtained using the Voronoi neighbors network are smaller than for its weighted counterpart, $Q_6^w$ (by 10-15% in our case,
result not shown).

The different representations of the data in Fig. 4 convey two distinct messages: Panel (a) shows that fluid and jammed states can have a similar volume fraction and yet clearly different local structure, as illustrated by the vertical dashed line around $\phi_J = 0.655$ [44, 78]. Panel (b) demonstrates instead that above the onset volume fraction, $\phi_{onset} \approx 0.56$, around which glassy dynamics starts to manifest in the equilibrium fluid [79], the equilibrium fluid and the corresponding jammed states are characterized by very similar bond orientational order; at the local scale, jammed configurations retain essentially the local structure of the parent fluid. This result is confirmed by analysis of the other local structure metrics detailed below.

The evolution of $Q_6^w$ is generally not enough to detect the presence of local crystal structures in dense packings. In particular, it is difficult to disentangle local fcc structures from distorted icosahedra [19], which both give similar $Q_6$ values. To resolve fine details of the local structure, one has to resort to scatter plots of rotational invariants of different orders [74], and possibly introduce an additional averaging procedure over the neighbors [80]. Mickel et al. [19] have suggested to measure $Q_2^w$ invariants, which are zero for most crystal structures as a simpler approach to quantify the degree of local crystalline order. We found that $Q_2^w$ decreases monotonically by increasing volume fraction for both fluid and jammed configurations, and that its variation is smooth and continuous (not shown). Typical values of the average $Q_2^w$ range from about 0.1 at $\phi = 0.4$ to 0.05 at the largest packing fractions. These values are consistent with those observed in non-crystalline packings of monodisperse hard spheres [19]. This, together with the smooth evolution of the fluid equation of state and the regularity of the partial structure factors at small $k$ [79], confirms that our packings are fully amorphous along the whole J-line.

## 5.2 Locally favored structure

A slight but systematic change of bond orientational order along the J-line has been noted before [34, 36], albeit over a smaller range of volume fractions. While this suggests some sort of ordering of the packings, uncovering the precise nature of local order remains difficult and requires knowledge of a dictionary of “known” reference structures. In this section, we employ a strategy based on the statistics of Voronoi cell shapes [81, 82] which does not rely on any a priori reference to specific motifs. This approach has been employed successfully in the context of glasses, see e.g. Refs. [46, 83] for recent reviews, but its use in the context of jammed packings remains very limited [20].

As in Sec. 5.1, we perform a radical Voronoi tessellation [76] of each configuration. In this construction, which duly accounts for the size polydispersity, each particle $i$ in the system is enclosed in a polyhedral cell such that all points whose tangent distance from the surface of particle $i$ is smaller than the tangent distance from the surface of particle $j$ with $j \neq i$ [76]. The shape of a Voronoi cell encodes detailed information about the local arrangements of the neighbors around the central particle, hence the local structure. We characterize the shape of a cell through its signature $(n_3, n_4, n_5, ...)$, where $n_q$ is the number of faces of the cell with a given number $q$ of vertices.<fnref n="1" /> By averaging over an ensemble of configurations, we detect the most frequent signatures and monitor their concentration as a function of the relevant control parameter (here the volume fraction). This approach provides a fairly robust assessment of the preferred local order, often termed “locally favored structure” in simple models of glasses [46].

We find that at low volume fractions the fluid exhibits a variety of different Voronoi signatures, typical of highly disordered fluids. Upon increasing $\phi_{fluid}$, however, fluid samples become increasingly rich in (0, 0, 12) signatures, which are associated to icosahedral local order. This kind of signature is by far the most frequent one at largest volume fractions. In

<footnote n="1">¹We disregard null values of $n_q$ for $q > q'$, where $q'$ is the largest number of vertices such that $n_{q'} > 0$.</footnote>
<table><caption>Table 1: Most frequent Voronoi signatures and the corresponding average percentages in jammed configurations. The corresponding parent fluid packing fractions are $\phi_{\text{fluid}} = 0.500$ ($N = 8000$), 0.618 ($N = 8000$), and 0.655 ($N = 1000$), respectively.</caption><thead><tr><th colspan="2">φ<sub>J</sub> = 0.656</th><th colspan="2">φ<sub>J</sub> = 0.678</th><th colspan="2">φ<sub>J</sub> = 0.697</th></tr></thead><tbody><tr><td>(0,2,8,1)</td><td>2.8%</td><td>(0,2,8,1)</td><td>5.1%</td><td>(0,0,12)</td><td>14.2%</td></tr><tr><td>(0,3,6,4)</td><td>2.5</td><td>(0,2,8,2)</td><td>4.4</td><td>(0,1,10,2)</td><td>7.0</td></tr><tr><td>(0,3,6,3)</td><td>2.4</td><td>(0,2,8)</td><td>4.1</td><td>(0,2,8)</td><td>6.3</td></tr><tr><td>(0,2,8,2)</td><td>2.3</td><td>(0,0,12)</td><td>3.9</td><td>(0,2,8,1)</td><td>5.7</td></tr><tr><td>(0,3,6,1)</td><td>2.0</td><td>(0,1,10,2)</td><td>3.7</td><td>(0,2,8,2)</td><td>4.6</td></tr><tr><td>(0,4,4,3)</td><td>1.9</td><td>(0,3,6,4)</td><td>3.2</td><td>(0,3,6)</td><td>3.7</td></tr><tr><td>(0,2,8)</td><td>1.8</td><td>(0,3,6,3)</td><td>2.8</td><td>(0,3,6,1)</td><td>3.0</td></tr><tr><td>(0,3,6,2)</td><td>1.7</td><td>(0,3,6,1)</td><td>2.8</td><td>(0,3,6,4)</td><td>2.4</td></tr><tr><td>(0,3,6)</td><td>1.6</td><td>(0,3,6)</td><td>2.3</td><td>(0,3,6,3)</td><td>2.1</td></tr><tr><td>(1,2,5,3)</td><td>1.4</td><td>(0,2,8,4)</td><td>2.1</td><td>(0,2,8,4)</td><td>2.0</td></tr></tbody></table>

Figs. 4(c, d), we show the fraction $f_{(0,0,12)}$ of icosahedral structures for both fluid and jammed states. As for panels (a) and (b), we plot the data as a function of $\phi_J$ and of $\phi_{\text{fluid}}$. We find that $f_{(0,0,12)}$ for both fluid and jammed states is essentially negligible below $\phi_{\text{fluid}} \approx 0.56$, but then it increases rapidly with increasing volume fraction to reach a value of about 14% for our densest packings. We note that if we also take into account all particles connected to centers of icosahedral structures, as done for instance in Ref. [84], this fraction reaches about 80% at the largest density. Superficially, this behavior resembles the one found in the binary Lennard-Jones mixture introduced in Ref. [85], for which $f_{(0,0,12)}$ shows a marked increase around the onset of slow dynamics [86].

Turning our attention to jammed packings, no distinguishable local order is detectable below $\phi_{\text{onset}}$. Figure 4(d) demonstrates that the sudden emergence of the icosahedral local order around $\phi_{\text{onset}}$ in the equilibrium fluid is, again, inherited by the jammed packings. It is thus tempting to attribute the increase of $\phi_J$ along the J-line to the appearance of distorted icosahedral arrangements, which help maintaining an overall amorphous organization of the packings [20]. In Table 1, we report the occurrences of the most frequent Voronoi signatures in jammed packings at selected volume fractions. We note that, on average, (0,0,12) cells have contact numbers Z slightly smaller than 6 ($Z = 5.4$). We found that some less abundant signatures, such as (0,3,6) and (0,2,8), are associated to even smaller values of Z. These trends suggest a correlation between Z and the number of nearest neighbors, which we tentatively attribute to the size dispersity of the packings.

The presence of icosahedral order and, more generally, polytetrahedral order in fluid and jammed hard sphere packings has been discussed before, see e.g. Refs [17, 20, 87, 88]. The results presented above differ from earlier findings on two important aspects. First, the amount of (0, 0, 12) signatures found in our polydisperse packings substantially exceeds the one reported in monodisperse jammed packings [20] and is also appreciably larger than the one reported in less polydisperse hard sphere at equilibrium [88]. However, because of the large size polydispersity, the local arrangements associated to (0, 0, 12) signatures are often distorted and irregular. Such issues are usually neglected in the analysis of icosahedral order in simple binary mixtures [86].

To quantify the degree of asphericity of Voronoi cells, we analyzed the distribution of distances $r_{ij} = |\mathbf{r}_i - \mathbf{r}_j|$ separating a central particle i to its neighbors j. As a simple measure of cell
<figure><img src="image_8.png" /><figcaption>Figure 5: (a): Distribution of asphericity parameter s for polydisperse hard spheres for two different volume fractions, and for a monodisperse hard sphere system. (b, c, d): Typical icosahedral structures with varying degree of asphericity s.</figcaption></figure>

asphericity, we compute the normalized standard deviation of the distances from particle i,

$$s_i = \frac{1}{\tilde{r}_i} \sqrt{\frac{1}{n_b(i)} \sum_{j=1}^{n_b(i)} (r_{ij} - \tilde{r}_i)^2},$$

where $\tilde{r}_i = (\sum_{j=1}^{n_b(i)} r_{ij})/n_b(i)$ is the average nearest neighbors distance of particle i. See Ref. [89] for a more systematic approach based on Minkowski tensors. Of course $s \ge 0$ and the equality holds for a perfectly regular structure. In Fig. 5(a) we show the probability distribution of the asphericity parameter, $P(s)$, of icosahedral structures for jammed packings for a low ($\phi_J = 0.653$) and a high ($\phi_J = 0.685$) packing fraction. We find that the distribution is not very sensitive to the value of $\phi_J$. For comparison, we also include results for $P(s)$ measured in monodisperse hard sphere jammed packings. We obtained the latter by compressing Poisson distributed configurations at $\phi_{fluid} = 0.3$ using the algorithm described in Sec. 2.3. On average, icosahedral structures detected in polydisperse packings are more aspherical by about a factor two than those found in the monodisperse hard spheres. We found that the asphericity of icosahedra in a glassy binary mixture with modest size ratio [85], for which icosahedral order is most pronounced [45], is intermediate between those of the two systems in Fig. 5(a). Figures 5(b-d) show representative icosahedral structures detected in our densest packings. They are representative of different degrees of asphericity, ranging from fairly regular (b) to intermediate (c) and highly irregular (d) structures.

These results indicate that, despite the abundance of (0, 0, 12) Voronoi cells, the nature of icosahedral order in polydisperse packings might be different from the one observed in
<figure><img src="image_9.png" /><figcaption>Figure 6: (a) Diameter distribution for all particles, $f(\sigma)$, and diameter distribution for rattlers only, $f_R(\sigma)$. (b) Cluster size distribution for the rattlers (closed symbols), compared to the same fraction of randomly selected particles in the same packings (open symbols).</figcaption></figure>

simple binary mixtures. In fact, icosahedral structures in highly polydisperse spheres are more distorted and inherently more irregular due to the increased compositional freedom. In binary mixtures with sufficiently small size ratio icosahedral structures can strongly influence the local mobility of the particles in the supercooled regime [45, 84, 90]. Further work is needed to assess the relevance of icosahedral local order in dense, highly polydisperse packings, in particular in the context of glass transition studies.

## 5.3 Spatial organization of rattlers

In Fig. 2(b) we have shown that the fraction of rattlers grows steadily with increasing $\phi_J$. Naively, one might expect rattlers to occupy a larger free space, thus a large number of rattlers seems at odds with the increased efficiency of the packings.

First, we assess the size distribution of rattlers, $f_R(\sigma)$, along the J-line, see Fig. 6(a). By comparing $f_R(\sigma)$ to the distribution of all particles, $f(\sigma)$, we find rattlers are mostly small particles, which is an intuitive result. Interestingly, the shape of $f_R(\sigma)$ hardly changes along the J-line, whereas the fraction of rattlers increases substantially.

Next, we inspect the spatial organization of rattlers. To get a qualitative idea of their real space structure, we show in Figs. 7 some typical snapshots of rattlers found in jammed packings at $\phi_J = 0.656$ and $\phi_J = 0.688$. Rattlers appear to be distributed in rather homogeneous way in this representation. To get a more quantitative picture, we follow previous work [91] and study the spatial organization of rattlers into clusters. We define a rattler cluster as a group of rattlers where each rattler is neighbor to at least another member of the cluster. As in the previous sections, neighbors are identified through a radical Voronoi tessellation. Note that our analysis differs from the one of Ref. [91], in which rattler clusters were classified according to the connectivity of their corresponding cages.

It is sometimes assumed that rattlers are randomly distributed in the system [92]. To test this idea, we randomly pick a fraction of particles with the same concentration as the rattlers and compute the corresponding cluster size using the same definition given above. We show the distribution of rattler cluster sizes for jammed states with a low ($\phi_J = 0.656$) and a high ($\phi_J = 0.688$) volume fraction in Fig. 6(b). We find that the rattler clusters are smaller than the one formed by the randomly chosen particles, for both low and high $\phi_J$ values. This implies that overall the positions of the rattlers are only weakly correlated, displaying a slight tendency to form small clusters [91]. This is confirmed by visual inspection of snapshots of jammed
<figure><img src="image_10.png" /><figcaption>Figure 7: Typical snapshots of rattlers (red spheres) in jammed packings at $\phi_J = 0.656$ (left) and $\phi_J = 0.688$ (right). The white bonds connect neighboring rattlers that form a cluster.</figcaption></figure>

packings along the J-line, see Fig. 7 for representative examples. These relatively compact clusters might lead to small volume fraction fluctuations and may also affect hyperuniformity, as discussed in Sec. 7.

# 6 Structure at the large scale

Here we consider the structure of the packings at the large scale. The corresponding wave number regime is $k \to 0$, which means that we quantify fluctuations at very large length scales.

## 6.1 Hyperuniformity

Equilibrium fluids are characterized by a finite value of the static structure factor $S(k)$ at $k \to 0$, which is associated to a finite isothermal compressibility. However, it is known that hard sphere packings at the jamming transition show “unexpected” (that is, more complex) density fluctuations [22, 93]. It was reported that $S(k)$ in monodisperse jammed packings obeys a surprising linear behavior at low $k$, $S(k) \propto k$, which characterizes hyperuniform materials [22, 93]. It is sometimes assumed that hyperuniformity is related to the criticality of the jammed states [47], but the nature of this criticality and the connection to other signatures of jamming remain unclear [25, 26].

When the system is a mixture of several components or is continuously polydisperse, the structure factor $S(k)$ of jammed states does not show a linearly vanishing behavior. Instead, $S(k)$ takes a finite value at $k \to 0$ just as in fluid states [94, 95]. We show $S(k)$ of our polydisperse system for both fluid and jammed states in Figs. 8(a, b). These results confirm that the shape of $S(k)$ of jammed states is similar to the fluid ones, taking in particular a finite value, $S(k \to 0) \approx 0.5$, for all jammed states along the J-line, with no sign of a vanishing signal at low $k$.

In Refs. [23, 92], the concept of hyperuniformity was generalized from density to volume fraction fluctuations to discuss the hyperuniformity of multi-components system and of
<figure><img src="image_11.png" /><figcaption>Figure 8: (a, b): Static structure factor $S(k)$ reflecting number density fluctuations for fluid (a) and jammed (b) states. (c, d): Spectral density $\chi(k)$ reflecting volume fraction fluctuations for fluid (c) and jammed (d) states. Nearly hyperuniform behavior expected for $\chi(k)$ is only observed for the lowest $\phi_J$ value in (d).</figcaption></figure>

particles with non-spherical shapes. Thus, an appropriate observable to characterize the hyperuniformity of polydisperse systems is the spectral density $\chi(k)$, which quantifies the volume fraction fluctuations [23, 92]. It is defined as

$$ \chi(\mathbf{k}) = \frac{1}{V} \langle I_{\mathbf{k}}I_{-\mathbf{k}} \rangle, \qquad (4) $$

where $I_k$ is the Fourier transform of the indicator function $I(\mathbf{r})$. For spherical particle systems, $I(\mathbf{r})$ is given by

$$ I(\mathbf{r}) = \sum_{i=1}^{N} \theta(R_i - |\mathbf{r} - \mathbf{r}_i|), \qquad (5) $$

where $R_i$ is the radius of the i-th particle, $R_i = \sigma_i/2$. For homogeneous isotropic systems in $d=3$, $I_k$ becomes

$$ I_k = \sum_{i=1}^{N} \frac{4\pi}{k^3} [\sin(kR_i) - (kR_i) \cos(kR_i)] e^{-i\mathbf{k} \cdot \mathbf{r}_i}. \qquad (6) $$

Note that alternative definitions of the volume fraction fluctuations in the reciprocal space can be applied to study hyperuniformity in jammed packings [24, 96].

In Figs. 8(c, d), we show our numerical results for $\chi(k)$ for both fluid and jammed states. Before discussing $\chi(k)$ of jammed states at small $k$, where signs of hyperuniformity may be
<figure><img src="image_12.png" /><figcaption>Figure 9: Zoom in the small $k$ region of the spectral density $\chi(k)$ for jammed states shown in Fig. 8(d). The error bars represent the standard deviation obtained from sample to sample fluctuations. The gray dashed line is a linear fit for $\phi_J = 0.656$, which extrapolates to $\chi(k \to 0) = 7.9 \times 10^{-4}$ for the lowest $\phi_J$. Instead $\chi(k \to 0) = 4.2 \times 10^{-3}$ for the largest $\phi_J$.</figcaption></figure>

found, it is instructive to study $\chi(k)$ for the fluid states and its volume fraction dependence. We find that $\chi(k)$ for fluid states has a broad peak near $k \sim 2\pi/\sigma$ and a flat plateau at smaller $k$ [97]. In addition, the overall amplitude of $\chi(k)$ decreases as the volume fraction $\phi_{\text{fluid}}$ increases at $k \lesssim 2\pi/\sigma$, reflecting that volume fraction fluctuation are suppressed in order to achieve denser particle packings in equilibrium. This observation contrasts with the well-known increase of the first peak of $S(k)$ with increasing the density in dense fluids [98] which instead reflects the slight increase in local scale structural correlations as the density is increased.

In jammed states, the peak height of $\chi(k)$ is much smaller than in corresponding fluid states, and the peak amplitude decreases as $\phi_J$ increases, reflecting again that further optimization of local configurations has been realized. Note that $\chi(k)$ displays a rather strong volume fraction dependence in both fluid and jammed states, which contrasts with the relatively weak density dependence of $S(k)$ over the same range. This insensitivity of $S(k)$ in supercooled liquids is traditionally taken as a hallmark of glass physics [98]. However, our results demonstrate that $\chi(k)$ is a more sensitive probe of structural changes at the level of two-point correlations. This suggests that $\chi(k)$ might be a good observable to characterize the equilibrium structure of dense hard sphere systems.

Having described the main features of $\chi(k)$, we now concentrate on the low-$k$ behavior in jammed states. A zoom in this region is shown in Fig. 9, to investigate more carefully the possible existence of a hyperuniform behavior in this regime. At the lowest $\phi_J$ along the J-line for $N = 8000$, $\phi_J = 0.656$, $\chi(k)$ linearly decreases with decreasing $k$, as observed in previous numerical works [23, 24, 47, 92]. By fitting this linear behavior (shown with the dashed line), we can extrapolate the limit $k \to 0$ and we obtain a small finite value, $\chi(k \to 0) = 7.9 \times 10^{-4}$ for this volume fraction. Therefore, one might conclude that the system at $\phi_J = 0.656$ is very close to being hyperuniform [47]. However, as $\phi_J$ is increased from its lowest value, $\chi(k)$ quickly deviates from this linear behavior, and it even becomes flat, within our error bars, at low $k$. Simultaneously, $\chi(k \to 0)$ obtained by linear extrapolation increases systematically from $\chi(k \to 0) \approx 7.9 \times 10^{-4}$ to $\chi(k \to 0) \approx 4.2 \times 10^{-3}$ with increasing $\phi_J$, suggesting that deviations from hyperuniformity become stronger when $\phi_J$ increases.

Recently, deviations from hyperuniformity were reported in both two [24] and three [25, 26] spatial dimensions. Both works employed rapid quenches from fully random configura-
tions using energy minimization protocols. These studies showed that $\chi(k)$ (and thus $S(k)$ in monodisperse systems) linearly decreases with decreasing $k$, but it always saturates at a certain $k^*$ ($k^* \sim 0.2$ in Ref. [24] and $k^* \sim 0.4$ in Refs. [25,26]), indicating that hyperuniformity is eventually avoided at sufficiently small $k$. Ref. [47] attributes the breakdown of hyperuniformity at very small $k$ to a lack of numerical accuracy in conventional jamming algorithms. It is argued that producing truly jammed, hyperuniform states is a highly difficult numerical task and thus the observed saturation at $k^*$ using energy minimization could stem from numerical inaccuracy of a given algorithm.

However, the systematic deviation of $\chi(k)$ at denser $\phi_J$ of our system shown in Fig. 9 presumably has a different origin, because the system size we use, $N = 8000$, is not large enough to detect the saturation at $k^*$ reported in previous work. Also, the systematic increase of $\chi(k \to 0)$ is observed for a given system size and within a given numerical algorithm. Therefore, the strong deviations from hyperuniformity reported here must have a physical, rather than a computational, origin.

Our jammed packings at various $\phi_J$ are obtained with a similar degree of precision, in particular isostaticity is equally well satisfied for all configurations along the J-line. However, they display an increasing tendency to depart from hyperuniformity as $\phi_J$ is increased. Therefore, we conclude that the jamming criticality, which is well-obeyed all along the J-line, and hyperuniformity, which is increasingly violated, are distinct concepts. This conclusion is consistent with the fact that the linear dependence of $\chi(k)$ (and $S(k)$), which might be interpreted as a “precursor” of hyperuniformity [99], is robustly observed at densities well below [99] and above [24,25] $\phi_J$ and even in the presence of thermal fluctuations [25], whereas the jamming criticality originating from isostaticity is very quickly erased in similar conditions [100].

Furthermore, the nearly hyperuniform behavior observed at the lowest $\phi_J$ in this work implies that hyperuniformity, if it were to be observed, should characterize the lowest end of the J-line [32, 40]. Although the “maximally random jammed state” [2] and the lowest end point of the J-line are conceptually distinct, our observations lead us to speculate that these two concepts may be identical [101], and may both display the strongest (although presumably still imperfect) signature of hyperuniformity. This conclusion, based on our numerical findings, also directly contradicts an opposite theoretical prediction recently made in Ref. [102].

A hypothesis that could explain the present observations is that the increasing deviation from hyperuniformity observed on increasing $\phi_J$ correlates with the increasing number of rattlers in the corresponding packings. Rattlers can be seen as “defects” in the volume fraction field, and a finite fraction of rattlers could induce a finite amount of volume fraction fluctuations at large scale. This hypothesis is hard to test directly, as we cannot independently vary the fraction of rattlers and the degree of hyperuniformity in hard sphere packings. We notice that the fraction of rattlers increases by only a factor of 2 while the limit $\chi(k \to 0)$ increases by about an order of magnitude, which suggests that rattlers may not be the central explanation for this observation. This conclusion was indirectly tested in Ref. [24], which showed that upon compression of jammed packings the number of rattlers decreased but the behavior of $\chi(k)$ was essentially unchanged. We also confirmed this behavior in our simulations (not shown). We will discuss this issue further in Section 7.

## 6.2 Finite-size fluctuations of the critical density of jamming

Second order phase transitions are usually associated with diverging length scales [103]. Finite-size scaling represents a powerful tool to understand the nature of such transitions and of the corresponding diverging length scales. Although several distinct important length scales have been identified for the jamming transition [7, 25, 99, 100, 104–110], the particular length scale responsible for the finite-size effects of $\phi_J$ is not well understood [5, 105, 111].
<figure><img src="image_13.png" /><figcaption>Figure 10: (a): Probability distribution function of $\phi_J$ for several $N$ obtained by compressing dilute ($\phi_{fluid} = 0.4$) and dense ($\phi_{fluid} = 0.635$) equilibrium fluid configurations. (b): Finite-size scaling plot according to Eq. (7). Finite-size effects on the quantity $\phi_J(N)$ are considerably suppressed for $\phi_{fluid} = 0.635$.</figcaption></figure>

A qualitatively different aspect of the jamming transition compared to usual critical phenomena is the fact that the transition point is not unique and is protocol-dependent [11]. Thus, we examine the finite-size effect of $\phi_J$ for two different protocols which produce different averaged $\phi_J$ values. In practice, we compare results of two protocols, the non-equilibrium compression from dilute ($\phi_{fluid} = 0.4$) and dense ($\phi_{fluid} = 0.635$) fluids, which produce low ($\phi_J \approx 0.657$) and high ($\phi_J \approx 0.685$) jamming volume fractions. We expect a smooth evolution of the behavior with $\phi_J$, but these measurements are numerically demanding, so we limit ourselves to only two different critical points on the extremes of the J-line.

We show the probability distribution function $P(\phi_J)$ for several system sizes $N$ in Fig. 10(a). For the compression of the dilute fluid, $\phi_{fluid} = 0.4$, a broad Gaussian shape of $P(\phi_J)$ is obtained for the smallest system size. As $N$ increases, the distribution becomes narrower and the position of the peak shifts to higher volume fractions, as reported previously [5, 60, 112].

For the compression of the dense fluid, $\phi_{fluid} = 0.635$, the width of the distribution for a given $N$ is slightly narrower than the one for $\phi_{fluid} = 0.4$ [14, 34], but the distribution again becomes more sharply peaked as $N$ increases. Remarkably, the peak position hardly shifts with
$N$, in strong contrast with $\phi_{\text{fluid}} = 0.4$. Therefore, for $\phi_{\text{fluid}} = 0.635$, the finite size effect in terms of the mean value of $\phi_J$ is significantly suppressed and $\phi_J(N)$ quickly approaches the thermodynamic value $\phi_J^\infty$. For both protocols, the latter value can be extracted using standard finite-size scaling,

$$ \phi_J(N) = \phi_J^\infty - \delta N^{-1/\nu d}, \qquad (7) $$

where $\delta$ and $\nu$ are fitting parameters, and $d$ is the number of spatial dimensions. Note that we use the mean value for $\phi_J(N)$ instead of the peak value [5], since we do not have enough statistics to determine the peak position with high accuracy. The results of this fit are shown in Fig. 10(b). We obtain $\phi_J^\infty = 0.65698 \pm 0.00009$ and $\nu = 0.64 \pm 0.01$ for $\phi_{\text{fluid}} = 0.4$, and $\phi_J^\infty = 0.68548 \pm 0.00002$ and $\nu = 0.41 \pm 0.04$ for $\phi_{\text{fluid}} = 0.635$, respectively. The obtained $\nu$'s are slightly different from previous reports [5, 113]. Note that $\nu$ for $\phi_{\text{fluid}} = 0.4$ is compatible with very recent numerical work [102]. However, we do not wish to discuss these values quantitatively because $\nu$ might be protocol or algorithm dependent, and corrections to scaling should be included before drawing any strong conclusions [111]. We do not treat these corrections to scaling because they require huge statistics. Our main result here is more qualitative; finite size effects on $\phi_J$ are significantly suppressed for the compression from the dense fluid. The prefactor $\delta$ in Eq. (7) is more than 10 times smaller for the compression of the dense fluid, as can be directly seen in the data shown in Fig. 10(b). The smallness of the finite-size effect in fact makes a quantitative determination of the critical exponent very difficult for the largest $\phi_J$ values [36]. Our data conclusively demonstrate that the J-line is not due to a finite size effect and that it remains instead well-defined in the thermodynamic limit.

To get more physical insights into the finite-size effects for $\phi_J$, we monitor the probability distribution function $P(\Delta r)$ of single particle displacements during the compression from the fluid to the jammed states, $\Delta r = |\mathbf{r}_i^{\text{jam}} - \mathbf{r}_i^{\text{fluid}}|$, where $\mathbf{r}_i^{\text{fluid}}$ and $\mathbf{r}_i^{\text{jam}}$ are the positions of particle $i$ in the equilibrium parent fluid configuration and in the corresponding jammed configuration, respectively. The distribution $P(\Delta r)$ is computed for all the particles (including rattlers). $P(\Delta r)$ quantifies how much the particles need to move or rearrange during the compression until the system is jammed. In Figs. 11(a, b), we show $P(\Delta r)$ as a function of either $\Delta r$ or $(\Delta r)^2$, for three different system sizes, $N = 150, 1000$, and $8000$. For $\phi_{\text{fluid}} = 0.4$, $P(\Delta r)$ has a rather long tail, indicating that the particles may perform large displacements before finding the final jammed configuration. This tail can be fitted by a Gaussian form, $P(\Delta r) \propto \exp[-b(\Delta r)^2]$, where $b$ is a constant [see Fig 11(b)], which suggests that there is no strong positional correlation between before and after compression. Importantly, a noticeable finite-size effect is found in these tails, in the sense that smaller systems are characterized by smaller particle displacements. This observation provides a possible explanation for the observed finite-size effect in $\phi_J$: particles in smaller systems do not explore space as they do in large systems and thus get jammed in less optimized (i.e. less dense) configurations, leading to smaller $\phi_J$, as observed in Fig. 10.

Interestingly, for $\phi_{\text{fluid}} = 0.635$, the tail in $P(\Delta r)$ is significantly suppressed compared to the one for $\phi_{\text{fluid}} = 0.4$, indicating that particles get jammed in a position that is actually very close to their original position in the parent fluid. This observation is of course consistent with the fact that the parent fluid and its corresponding jammed state display very similar geometric structure when $\phi_J$ is large, as shown above in Figs. 4(b, d). Simultaneously, the finite-size effect on the tail in $P(\Delta r)$ almost disappears, which parallels the suppression of the finite-size effect in $\phi_J$ values. Furthermore, $P(\Delta r)$ is no longer Gaussian, but exhibits an exponential tail, $P(\Delta r) \propto \exp[-a\Delta r]$, where $a$ is a constant. This functional form has been reported in supercooled liquids undergoing a transition between nearby inherent structures [114–116]. This may suggest that during the compression of the denser fluid, the system may perform transitions among nearby locally stable configurations, while remaining firmly localised within a well-defined metabasin [64, 117–119].
<figure><img src="image_14.png" /><figcaption>Figure 11: (a): Probability distribution function $P(\Delta r)$ of particle displacements during a compression, $\Delta r = |\mathbf{r}_i^{\text{jam}} - \mathbf{r}_i^{\text{fluid}}|$. The straight line corresponds to $P(\Delta r) \propto \exp[-a\Delta r]$, where $a$ is a constant. (b): The same distributions $P(\Delta r)$ as a function of the quantity $(\Delta r)^2$. The straight line corresponds to $P(\Delta r) \propto \exp[-b(\Delta r)^2]$, where $b$ is a constant.</figcaption></figure>

# 7 Discussion and conclusions

Thanks to an efficient thermalisation algorithm which allowed us to equilibrate polydisperse hard spheres up to unprecedented packing fractions [44], we have significantly extended the study of the line of critical jamming transitions, or J-line. Our results demonstrate that isostaticity and the associated critical behavior of the pair correlation function remain unchanged along the entire J-line, while some other structural properties at larger length scales evolve qualitatively. Therefore, our results disentangle the structural properties which originate from isostaticity at the contact scale from other geometric properties at larger scale that seem unrelated or insensitive to the jamming criticality.

We found that the fraction of rattlers in jammed packings increases markedly with increasing $\phi_J$, which is a counterintuitive result. We confirm that a proper treatment of these rattlers is essential to reveal the critical behavior at contact scale. Rattlers also affect the organization of the packings at the neighbor scale and tend to form relatively compact clusters, which might induce subtle volume fraction fluctuations in the packing. In Ref. [47], it has been argued that
the presence of rattlers may interfere with hyperuniformity. We found that deviations from hyperuniform behavior in jammed polydisperse hard spheres become more pronounced as the volume fraction increases. However, hyperuniformity cannot be simply achieved by reducing the number of rattlers through additional compressions [24]. Moreover, hyperuniform behavior disappears completely when rattlers are removed from the computation of the spectral density, or volume fraction fluctuation [22]. Our results confirm this trend, which is in stark contrast with the jamming critical behavior seen in the pair correlation function. This would be another supporting evidence that hyperuniformity and jamming criticality are unrelated. These numerical results does not exclude the possibility that an ideal, rattler-free packing obtained through a more complex optimization process provides true hyperuniformity [47]. Along this line of thought, it would be interesting to study in more detail the relation between presence of rattlers and the violation of the hyperuniformity reported here.

We found that, above the onset volume fraction, the structure of jammed packings at the local scale closely tracks the one of the parent equilibrium fluid. This contrasts with the behavior at the contact scale, which is qualitatively very different in jammed and fluid states. In particular, as the volume fraction increases, both equilibrium fluid and jammed packings display a smooth but progressive local ordering, which we attributed to the emergence of distorted icosahedral structures. We attributed their irregular shape to the significant size dispersity of the system. We speculate that only the most regular structures, identified through the asphericity parameter used in this work or more sophisticated metrics [19], may actually provide locally stable arrangements in the dense equilibrium fluid. These ideas might find useful applications in the context of glass structure studies [46] and of investigations of the structure-dynamics relationship [90].

In a recent numerical study of polydisperse hard spheres [79], we demonstrated the growth of non-trivial static “point-to-set” correlations, which are at the core of thermodynamic pictures of the glass transition [120]. It has been suggested that such “amorphous order” could be identified with hyperuniformity [99, 102], although this link has been questioned [25]. Our results cast further doubts on this argument: point-to-set correlations grow as the volume fraction increases, while hyperuniform behavior is suppressed. The growth of amorphous order in the equilibrium fluid appears thus to compete with, rather than to enhance, hyperuniformity.

The significant suppression of finite size effects at higher $\phi_J$ suggests that the length scale causing this finite-size effect might not be related to the isostatic nature of the jamming transition, since all studied configurations are similarly isostatic. This conclusion seems to contrast with other diverging length scales near jamming [25, 100, 104, 107–109]. We notice that the finite-size effect for inherent structure energies in Lennard-Jones supercooled liquids is also suppressed when considering lower temperatures for the parent fluid [121], which is reminiscent of our observations for denser $\phi_{\text{fluid}}$ in hard spheres. Finite size effects for jammed and inherent structures may have a common physical origin, namely the different topography of the energy landscape probed in dilute or dense fluids [28, 122]. Dense fluids jam in a configuration that belongs to well-defined metabasins, whereas dilute fluids may explore a larger portion of the potential energy landscape before jamming, larger systems exploring a larger region of phase space.

We have used a system with continuous polydispersity to enhance thermalisation. However, we expect the conclusions drawn from the present model to hold more generally. Systems with discrete polydispersity, i.e. mixtures, are more easily amenable to theoretical investigations [101, 123–125]. Recently, it has been proposed that thermal glassy systems with continuous polydispersity can be mapped into an effective multi-component system to characterize their thermodynamic behavior [126]. Whether such an effective description is also applicable to the structural characterization of jammed states is an interesting open issue, which will likely require a careful treatment of the rattlers [124, 125].
Acknowledgments

We thank D. Durian, P. Charbonneau, A. Ikeda, H. Ikeda, K. Kim, A. Liu, A. Ninarello, M. Pica Ciamarra, P. Sollich, and D. Vågberg for helpful discussions. The numerical simulations were partially performed at Research Center of Computational Science (RCCS), Okazaki, Japan. We thank K. Kim for providing us with CPU time in RCCS. The research leading to these results has received funding from the European Research Council under the European Unions Seventh Framework Programme (No. FP7/2007-2013)/ERC Grant Agreement No.306845. This work was supported by a grant from the Simons Foundation (No. 454933, Ludovic Berthier).

References

[1] A. J. Liu and S. R. Nagel, *Nonlinear dynamics: Jamming is not just cool any more*, Nature **396**, 21 (1998), doi:10.1038/23819.

[2] S. Torquato and F. H. Stillinger, *Jammed hard-particle packings: From Kepler to Bernal and beyond*, Rev. Mod. Phys. **82**, 2633 (2010), doi:10.1103/RevModPhys.82.2633.

[3] M. Van Hecke, *Jamming of soft particles: geometry, mechanics, scaling and isostaticity*, J. Phys.: Condens. Matter **22**, 033101 (2009), doi:10.1088/0953-8984/22/3/033101.

[4] A. J. Liu and S. R. Nagel, *The jamming transition and the marginally jammed solid*, Annu. Rev. Condens. Matter Phys. **1**, 347 (2010), doi:10.1146/annurev-conmatphys-070909-104045.

[5] C. S. O'Hern, L. E. Silbert, A. J. Liu and S. R. Nagel, *Jamming at zero temperature and zero applied stress: The epitome of disorder*, Phys. Rev. E **68**, 011306 (2003), doi:10.1103/PhysRevE.68.011306.

[6] M. Wyart, L. E. Silbert, S. R. Nagel and T. A. Witten, *Effects of compression on the vibrational modes of marginally jammed solids*, Phys. Rev. E **72**, 051306 (2005), doi:10.1103/PhysRevE.72.051306.

[7] P. Olsson and S. Teitel, *Critical scaling of shear viscosity at the jamming transition*, Phys. Rev. Lett. **99**, 178001 (2007), doi:10.1103/PhysRevLett.99.178001.

[8] E. Lerner, G. Düring and M. Wyart, *A unified framework for non-Brownian suspension flows and soft amorphous solids*, Proc. Natl. Acad. Sci. **109**, 4798 (2012), doi:10.1073/pnas.1120215109.

[9] A. Donev, S. Torquato and F. H. Stillinger, *Pair correlation function characteristics of nearly jammed disordered and ordered hard-sphere packings*, Phys. Rev. E **71**, 011105 (2005), doi:10.1103/PhysRevE.71.011105.

[10] P. Charbonneau, E. I. Corwin, G. Parisi and F. Zamponi, *Universal microstructure and mechanical stability of jammed packings*, Phys. Rev. Lett. **109**, 205501 (2012), doi:10.1103/PhysRevLett.109.205501.

[11] C. P. Goodrich, A. J. Liu and J. P. Sethna, *Scaling ansatz for the jamming transition*, Proc. Natl. Acad. Sci. **113**, 9745 (2016), doi:10.1073/pnas.1601858113.

[12] M. Clusel, E. I. Corwin, A. O. Siemens and J. Brujić, *A 'granocentric' model for random packing of jammed emulsions*, Nature **460**, 611 (2009), doi:10.1038/nature08158.
[13] V. Ogarko, N. Rivas and S. Luding, Communication: Structure characterization of hard sphere packings in amorphous and crystalline states, J. Chem. Phys. **140**, 211102 (2014), doi:10.1063/1.4880236.

[14] J. M. Rieser, C. P. Goodrich, A. J. Liu and D. J. Durian, Divergence of voronoi cell anisotropy vector: A threshold-free characterization of local structure in amorphous materials, Phys. Rev. Lett. **116**, 088001 (2016), doi:10.1103/PhysRevLett.116.088001.

[15] S. Torquato, T. M. Truskett and P. G. Debenedetti, *Is random close packing of spheres well defined?*, Phys. Rev. Lett. **84**, 2064 (2000), doi:10.1103/PhysRevLett.84.2064.

[16] P. K. Morse and E. I. Corwin, *Geometric signatures of jamming in the mechanical vacuum*, Phys. Rev. Lett. **112**, 115701 (2014), doi:10.1103/PhysRevLett.112.115701.

[17] A. V. Anikeenko and N. N. Medvedev, *Polytetrahedral nature of the dense disordered packings of hard spheres*, Phys. Rev. Lett. **98**, 235504 (2007), doi:10.1103/PhysRevLett.98.235504.

[18] S. C. Kapfer, W. Mickel, K. Mecke and G. E. Schröder-Turk, *Jammed spheres: Minkowski tensors reveal onset of local crystallinity*, Phys. Rev. E **85**, 030301 (2012), doi:10.1103/PhysRevE.85.030301.

[19] W. Mickel, S. C. Kapfer, G. E. Schröder-Turk and K. Mecke, *Shortcomings of the bond orientational order parameters for the analysis of disordered particulate matter*, J. Chem. Phys. **138**, 044501 (2013), doi:10.1063/1.4774084.

[20] M. A. Klatt and S. Torquato, *Characterization of maximally random jammed sphere packings: Voronoi correlation functions*, Phys. Rev. E **90**, 052120 (2014), doi:10.1103/PhysRevE.90.052120.

[21] B. A. Klumov, Y. Jin and H. A. Makse, *Structural properties of dense hard sphere packings*, J. Phys. Chem. B **118**, 10761 (2014), doi:10.1021/jp504537n.

[22] A. Donev, F. H. Stillinger and S. Torquato, *Unexpected density fluctuations in jammed disordered sphere packings*, Phys. Rev. Lett. **95**, 090604 (2005), doi:10.1103/PhysRevLett.116.088001.

[23] C. E. Zachary, Y. Jiao and S. Torquato, *Hyperuniformity, quasi-long-range correlations, and void-space constraints in maximally random jammed particle packings*. *I. Polydisperse spheres*, Phys. Rev. E **83**, 051308 (2011), doi:10.1103/PhysRevE.83.051308.

[24] Y. Wu, P. Olsson and S. Teitel, *Search for hyperuniformity in mechanically stable packings of frictionless disks above jamming*, Phys. Rev. E **92**, 052206 (2015), doi:10.1103/PhysRevE.92.052206.

[25] A. Ikeda and L. Berthier, *Thermal fluctuations, mechanical response, and hyperuniformity in jammed solids*, Phys. Rev. E **92**, 012309 (2015), doi:10.1103/PhysRevE.92.012309.

[26] A. Ikeda, L. Berthier and G. Parisi, *The large-scale structure of randomly jammed particles*, Phys. Rev. E **95**, 052125 (2017), doi:10.1103/PhysRevE.95.052125.

[27] R. J. Speedy, *Random jammed packings of hard discs and spheres*, J. Phys.: Condens. Matter **10**, 4185 (1998), doi:10.1088/0953-8984/10/19/006.

[28] R. Mari, F. Krzakala and J. Kurchan, *Jamming versus glass transitions*, Phys. Rev. Lett. **103**, 025701 (2009), doi:10.1103/PhysRevLett.103.025701.
[29] I. Biazzo, F. Caltagirone, G. Parisi and F. Zamponi, *Theory of amorphous packings of binary mixtures of hard spheres*, Phys. Rev. Lett. **102**, 195701 (2009), doi:**10.1103/PhysRevLett.102.195701**.

[30] M. Hermes and M. Dijkstra, *Jamming of polydisperse hard spheres: The effect of kinetic arrest*, EPL **89**, 38005 (2010), doi:**10.1209/0295-5075/89/38005**.

[31] P. Chaudhuri, L. Berthier and S. Sastry, *Jamming transitions in amorphous packings of frictionless spheres occur over a continuous range of volume fractions*, Phys. Rev. Lett. **104**, 165701 (2010), doi:**10.1103/PhysRevLett.104.165701**.

[32] M. P. Ciamarra, M. Nicodemi and A. Coniglio, *Recent results on the jamming phase diagram*, Soft Matter **6**, 2871 (2010), doi:**10.1039/B926810C**.

[33] C. F. Schreck, C. S. O'Hern and L. E. Silbert, *Tuning jammed frictionless disk packings from isostatic to hyperstatic*, Phys. Rev. E **84**, 011305 (2011), doi:**10.1103/PhysRevE.84.011305**.

[34] D. Vågberg, P. Olsson and S. Teitel, *Glassiness, rigidity, and jamming of frictionless soft core disks*, Phys. Rev. E **83**, 031307 (2011), doi:**10.1103/PhysRevE.83.031307**.

[35] M. Otsuki and H. Hayakawa, *Critical scaling of a jammed system after a quench of temperature*, Phys. Rev. E **86**, 031505 (2012), doi:**10.1103/PhysRevE.86.031505**.

[36] M. Ozawa, T. Kuroiwa, A. Ikeda and K. Miyazaki, *Jamming transition and inherent structures of hard spheres and disks*, Phys. Rev. Lett. **109**, 205701 (2012), doi:**10.1103/PhysRevLett.109.205701**.

[37] V. Baranau and U. Tallarek, *Random-close packing limits for monodisperse and polydisperse hard spheres*, Soft Matter **10**, 3826 (2014), doi:**10.1039/c3sm52959b**.

[38] V. Baranau and U. Tallarek, *On the jamming phase diagram for frictionless hard-sphere packings*, Soft Matter **10**, 7838 (2014), doi:**10.1039/C4SM01439A**.

[39] T. Bertrand, R. P. Behringer, B. Chakraborty, C. S. O'Hern and M. D. Shattuck, *Protocol dependence of the jamming transition*, Phys. Rev. E **93**, 012901 (2016), doi:**10.1103/PhysRevE.93.012901**.

[40] N. Kumar and S. Luding, *Memory of jamming–multiscale models for soft and granular matter*, Granular Matter **18**, 1 (2016), doi:**10.1007/s10035-016-0624-2**.

[41] S. Luding, *So much for the jamming point*, Nat. Phys. **12**, 531 (2016), doi:**10.1038/nphys3680**.

[42] L. Berthier and T. A. Witten, *The glass transition of dense fluids of hard and compressible spheres*, Phys. Rev. E **80**, 021502 (2009), doi:**10.1103/PhysRevE.80.021502**.

[43] M. Maiti and S. Sastry, *Free volume distribution of nearly jammed hard sphere packings*, J. Chem. Phys. **141**, 044510 (2014), doi:**10.1063/1.4891358**.

[44] L. Berthier, D. Coslovich, A. Ninarello and M. Ozawa, *Equilibrium sampling of hard spheres up to the jamming density and beyond*, Phys. Rev. Lett. **116**, 238002 (2016), doi:**10.1103/PhysRevLett.116.238002**.

[45] D. Coslovich and G. Pastore, *Understanding fragility in supercooled Lennard-Jones mixtures. I. Locally preferred structures*, J. Chem. Phys. **127**, 124504 (2007), doi:**10.1063/1.2773716**.
[46] C. P. Royall and S. R. Williams, *The role of local structure in dynamical arrest*, Phys. Rep. **560**, 1 (2015), doi:10.1016/j.physrep.2014.11.004.

[47] S. Atkinson, G. Zhang, A. B. Hopkins and S. Torquato, *Critical slowing down and hyperuniformity on approach to jamming*, Phys. Rev. E **94**, 012902 (2016), doi:10.1103/PhysRevE.94.012902.

[48] M. P. Allen and D. J. Tildesley, *Computer simulation of liquids*, Oxford University Press, Oxford (1989).

[49] A. Santos, S. B. Yuste and M. López de Haro, *Contact values of the particle-particle and wall-particle correlation functions in a hard-sphere polydisperse fluid*, J. Chem. Phys. **123**, 234512 (2005), doi:10.1063/1.2136883.

[50] A. Ninarello, L. Berthier and D. Coslovich, *Models and algorithms for the next generation of glass transition studies*, Phys. Rev. X **7**, 021039 (2017), doi:10.1103/PhysRevX.7.021039.

[51] N. B. Wilding and P. Sollich, *Phase equilibria and fractionation in a polydisperse fluid*, EPL **67**, 219 (2004), doi:10.1209/epl/i2004-10064-2.

[52] N. B. Wilding and P. Sollich, *Phase behavior of polydisperse spheres: Simulation strategies and an application to the freezing transition*, J. Chem. Phys. **133**, 224102 (2010), doi:10.1063/1.3510534.

[53] P. Sindzingre, C. Massobrio, G. Ciccotti and D. Frenkel, *Calculation of partial enthalpies of an argon-krypton mixture by NPT molecular dynamics*, Chem. Phys. **129**, 213 (1989), doi:10.1016/0301-0104(89)80007-2.

[54] D. Gazzillo and G. Pastore, *Equation of state for symmetric non-additive hard-sphere fluids: An approximate analytic expression and new Monte Carlo results*, Chem. Phys. Lett. **159**, 388 (1989), doi:10.1016/0009-2614(89)87505-0.

[55] L. Santen and W. Krauth, *Liquid, glass and crystal in two-dimensional hard disks*, arXiv:cond-mat/0107459.

[56] T. S. Grigera and G. Parisi, *Fast Monte Carlo algorithm for supercooled soft spheres*, Phys. Rev. E **63**, 045102 (2001), doi:10.1103/PhysRevE.63.045102.

[57] Y. Brumer and D. R. Reichman, *Numerical investigation of the entropy crisis in model glass formers*, J. Phys. Chem. B **108**, 6832 (2004), doi:10.1021/jp037617y.

[58] L. Fernández, V. Martin-Mayor and P. Verrocchio, *Phase diagram of a polydisperse soft-spheres model for liquids and colloids*, Phys. Rev. Lett. **98**, 085702 (2007), doi:10.1103/PhysRevLett.98.085702.

[59] R. Gutiérrez, S. Karmakar, Y. G. Pollack and I. Procaccia, *The static lengthscale characterizing the glass transition at lower temperatures*, EPL **111**, 56009 (2015), doi:10.1209/0295-5075/111/56009.

[60] N. Xu, J. Blawzdziewicz and C. S. O'Hern, *Random close packing revisited: Ways to pack frictionless disks*, Phys. Rev. E **71**, 061306 (2005), doi:10.1103/PhysRevE.71.061306.

[61] K. W. Desmond and E. R. Weeks, *Random close packing of disks and spheres in confined geometries*, Phys. Rev. E **80**, 051305 (2009), doi:10.1103/PhysRevE.80.051305.
[62] J. Nocedal and S. Wright, *Numerical Optimization*, Springer Verlag, Berlin (1999), doi:10.1007/978-0-387-40065-5.

[63] S. Dagois-Bohy, B. P. Tighe, J. Simon, S. Henkes and M. van Hecke, *Soft-sphere packings at finite pressure but unstable to shear*, Phys. Rev. Lett. **109**, 095703 (2012), doi:10.1103/PhysRevLett.109.095703.

[64] L. Berthier, P. Charbonneau, Y. Jin, G. Parisi, B. Seoane and F. Zamponi, *Growing timescales and lengthscales characterizing vibrations of amorphous solids*, Proc. Natl. Acad. Sci USA **113**, 8397 (2016), doi:10.1073/pnas.1607730113.

[65] K. W. Desmond and E. R. Weeks, *Influence of particle size distribution on random close packing of spheres*, Phys. Rev. E **90**, 022204 (2014), doi:10.1103/PhysRevE.90.022204.

[66] S. Sastry, P. G. Debenedetti and F. H. Stillinger, *Signatures of distinct dynamical regimes in the energy landscape of a glass-forming liquid*, Nature **393**, 554 (1998), doi:10.1038/31189.

[67] Y. Brumer and D. R. Reichman, *Mean-field theory, mode-coupling theory, and the onset temperature in supercooled liquids*, Phys. Rev. E **69**, 041202 (2004), doi:10.1103/PhysRevE.69.041202.

[68] M. Müller and M. Wyart, *Marginal stability in structural, spin, and electron glasses*, Annu. Rev. Cond. Matt. Phys. **6**, 177 (2015), doi:10.1146/annurev-conmatphys-031214-014614.

[69] C. Zhang, C. B. O'Donovan, E. I. Corwin, F. Cardinaux, T. G. Mason, M. E. Möbius and F. Scheffold, *Structure of marginally jammed polydisperse packings of frictionless spheres*, Phys. Rev. E **91**, 032302 (2015), doi:10.1103/PhysRevE.91.032302.

[70] P. Charbonneau, J. Kurchan, G. Parisi, P. Urbani and F. Zamponi, *Fractal free energy landscapes in structural glasses*, Nat. Comm. **5** (2014), doi:10.1038/ncomms4725.

[71] L. E. Silbert, A. J. Liu and S. R. Nagel, *Structural signatures of the unjamming transition at zero temperature*, Phys. Rev. E **73**, 041304 (2006), doi:10.1103/PhysRevE.73.041304.

[72] E. Lerner, G. Düring and M. Wyart, *Low-energy non-linear excitations in sphere packings*, Soft Matter **9**, 8252 (2013), doi:10.1039/C3SM50515D.

[73] P. Charbonneau, E. I. Corwin, G. Parisi and F. Zamponi, *Jamming criticality revealed by removing localized buckling excitations*, Phys. Rev. Lett. **114**, 125504 (2015), doi:10.1103/physrevlett.114.125504.

[74] P. J. Steinhardt, D. R. Nelson and M. Ronchetti, *Bond-orientational order in liquids and glasses*, Phys. Rev. B **28**, 784 (1983), doi:10.1103/PhysRevB.28.784.

[75] D. Moroni, P. R. ten Wolde and P. G. Bolhuis, *Interplay between structure and size in a critical crystal nucleus*, Phys. Rev. Lett. **94**, 235703 (2005), doi:10.1103/PhysRevLett.94.235703.

[76] B. J. Gellatly and J. L. Finney, *Characterisation of models of multicomponent amorphous metals: The radical alternative to the Voronoi polyhedron*, J. Non-Cryst. Solids **50**, 313 (1982), doi:10.1016/0022-3093(82)90093-X.

[77] C. H. Rycroft, *Voro++: A three-dimensional Voronoi cell library in C++*, Chaos **19**, 041111 (2009), doi:10.1063/1.3215722.
[78] M. Z. Yamchi, S. Ashwin and R. K. Bowles, *Inherent structures, fragility, and jamming: Insights from quasi-one-dimensional hard disks*, Phys. Rev. E **91**, 022301 (2015), doi:10.1103/PhysRevE.91.022301.

[79] L. Berthier, P. Charbonneau, D. Coslovich, A. Ninarello, M. Ozawa and S. Yaida, *Configurational entropy measurements in extremely supercooled liquids that break the glass ceiling*, Proc. Natl. Acad. Sci. (to be published 2017), doi:10.1073/pnas.1706860114.

[80] W. Lechner and C. Dellago, *Accurate determination of crystal structures based on averaged local bond order parameters*, J. Chem. Phys. **129**, 114707 (2008), doi:10.1063/1.2977970.

[81] J. L. Finney, *Fine structure in randomly packed, dense clusters of hard spheres*, Mat. Sci. Eng. **23**, 199 (1976), doi:10.1016/0025-5416(76)90194-4.

[82] M. Tanemura, Y. Hiwatari, H. Matsuda, T. Ogawa, N. Ogita and A. Ueda, *Geometrical analysis of crystallization pf the soft-core model*, Prog. Theor. Phys. **58**, 1079 (1977), doi:10.1143/PTP.58.1079.

[83] E. Ma, *Tuning order in disorder*, Nat. Mater. **14**, 547 (2015), doi:10.1038/nmat4300.

[84] A. Malins, J. Eggers, C. P. Royall, S. R. Williams and H. Tanaka, *Identification of long-lived clusters and their link to slow dynamics in a model glass former*, J. Chem. Phys. **138**, 12A535 (2013), doi:10.1063/1.4790515.

[85] G. Wahnstrom, *Molecular-dynamics study of a supercooled two-component Lennard-Jones system*, Phys. Rev. A **44**, 3752 (1991), doi:10.1103/PhysRevA.44.3752.

[86] D. Coslovich, *Locally preferred structures and many-body static correlations in viscous liquids*, Phys. Rev. E **83**, 051505 (2011), doi:10.1103/PhysRevE.83.051505.

[87] M. Leocmach and H. Tanaka, *Roles of icosahedral and crystal-like order in the hard spheres glass transition*, Nat. Comm. **3**, 974 (2012), doi:10.1038/ncomms1974.

[88] A. J. Dunleavy, K. Wiesner, R. Yamamoto and C. P. Royall, *Mutual information reveals multiple structural relaxation mechanisms in a model glass former*, Nat. Comm. **6**, 6089 (2015), doi:10.1038/ncomms7089.

[89] G. E. Schröder-Turk, W. Mickel, M. Schröter, G. W. Delaney, M. Saadatfar, T. J. Senden, K. Mecke and T. Aste, *Disordered spherical bead packs are anisotropic*, EPL **90**, 34001 (2010), doi:10.1209/0295-5075/90/34001.

[90] G. M. Hocky, D. Coslovich, A. Ikeda and D. R. Reichman, *Correlation of local order with particle mobility in supercooled liquids is highly system dependent*, Phys. Rev. Lett. **113**, 157801 (2014), doi:10.1103/PhysRevLett.113.157801.

[91] S. Atkinson, F. H. Stillinger and S. Torquato, *Detailed characterization of rattlers in exactly isostatic, strictly jammed sphere packings*, Phys. Rev. E **88**, 062208 (2013), doi:10.1103/PhysRevE.88.062208.

[92] C. E. Zachary, Y. Jiao and S. Torquato, *Hyperuniform long-range correlations are a signature of disordered jammed hard-particle packings*, Phys. Rev. Lett. **106**, 178001 (2011), doi:10.1103/PhysRevLett.106.178001.

[93] S. Torquato and F. H. Stillinger, *Local density fluctuations, hyperuniformity, and order metrics*, Phys. Rev. E **68**, 041113 (2003), doi:10.1103/PhysRevE.68.041113.
[94] N. Xu and E. S. Ching, *Effects of particle-size ratio on jamming of binary mixtures at zero temperature*, Soft Matter **6**, 2944 (2010), doi:10.1039/b926696h.

[95] R. Kurita and E. R. Weeks, *Experimental study of random-close-packed colloidal particles*, Phys. Rev. E **82**, 011403 (2010), doi:10.1103/PhysRevE.82.011403.

[96] L. Berthier, P. Chaudhuri, C. Coulais, O. Dauchot and P. Sollich, *Suppressed compressibility at large scale in jammed packings of size-disperse spheres*, Phys. Rev. Lett. **106**, 120601 (2011), doi:10.1103/PhysRevLett.106.120601.

[97] M. A. Klatt and S. Torquato, *Characterization of maximally random jammed sphere packings. II. Correlation functions and density fluctuations*, Phys. Rev. E **94**, 022152 (2016), doi:10.1103/PhysRevE.94.022152.

[98] K. Binder and W. Kob, *Glassy materials and disordered solids: An introduction to their statistical mechanics*, World Scientific (2011), doi:10.1142/7300.

[99] A. B. Hopkins, F. H. Stillinger and S. Torquato, *Nonequilibrium static diverging length scales on approaching a prototypical model glassy state*, Phys. Rev. E **86**, 021505 (2012), doi:10.1103/PhysRevE.86.021505.

[100] A. Ikeda, L. Berthier and G. Biroli, *Dynamic criticality at the jamming transition*, J. Chem. Phys. **138**, 12A507 (2013), doi:10.1063/1.4769251.

[101] G. Parisi and F. Zamponi, *Mean-field theory of hard sphere glasses and jamming*, Rev. Mod. Phys. **82**, 789 (2010), doi:10.1103/RevModPhys.82.789.

[102] A. Coniglio, M. P. Ciamarra and T. Aste, *Universal behaviour of the glass and the jamming transitions in finite dimensions*, arXiv:1704.01231.

[103] H. Nishimori and G. Ortiz, *Elements of phase transitions and critical phenomena*, Oxford University Press, Oxford (2010), doi:10.1093/acprof:oso/9780199577224.001.0001.

[104] L. E. Silbert, A. J. Liu and S. R. Nagel, *Vibrations and diverging length scales near the unjamming transition*, Phys. Rev. Lett. **95**, 098301 (2005), doi:10.1103/PhysRevLett.95.098301.

[105] J. Drocco, M. Hastings, C. O. Reichhardt and C. Reichhardt, *Multiscaling at point J: Jamming is a critical phenomenon*, Phys. Rev. Lett. **95**, 088001 (2005), doi:10.1103/PhysRevLett.95.088001.

[106] C. Reichhardt and C. O. Reichhardt, *Aspects of jamming in two-dimensional athermal frictionless systems*, Soft Matter **10**, 2932 (2014), doi:10.1039/c3sm53154f.

[107] C. P. Goodrich, W. G. Ellenbroek and A. J. Liu, *Stability of jammed packings I: the rigidity length scale*, Soft Matter **9**, 10993 (2013), doi:10.1039/C3SM51095F.

[108] S. S. Schoenholz, C. P. Goodrich, O. Kogan, A. J. Liu and S. R. Nagel, *Stability of jammed packings II: the transverse length scale*, Soft Matter **9**, 11000 (2013), doi:10.1039/C3SM51096D.

[109] W. G. Ellenbroek, M. van Hecke and W. van Saarloos, *Jammed frictionless disks: Connecting local and global response*, Phys. Rev. E **80**, 061307 (2009), doi:10.1103/PhysRevE.80.061307.

[110] H. Ikeda and A. Ikeda, *One-dimensional KAC model of dense amorphous hard spheres*, EPL **111**, 40007 (2015), doi:10.1209/0295-5075/111/40007.
[111] D. Vågberg, D. Valdez-Balderas, M. Moore, P. Olsson and S. Teitel, *Finite-size scaling at the jamming transition: Corrections to scaling and the correlation-length critical exponent*, Phys. Rev. E **83**, 030303 (2011), doi:10.1103/PhysRevE.83.030303.

[112] G.-J. Gao, J. Bławzdziewicz and C. S. O'Hern, *Frequency distribution of mechanically stable disk packings*, Phys. Rev. E **74**, 061304 (2006), doi:10.1103/PhysRevE.74.061304.

[113] M. P. Ciamarra, A. Coniglio and A. de Candia, *Disordered jammed packings of frictionless spheres*, Soft Matter **6**, 2975 (2010), doi:10.1039/C001904F.

[114] T. B. Schrøder, S. Sastry, J. C. Dyre and S. C. Glotzer, *Crossover to potential energy landscape dominated dynamics in a model glass-forming liquid*, J. Chem. Phys. **112**, 9834 (2000), doi:10.1063/1.481621.

[115] M. Vogel, B. Doliwa, A. Heuer and S. C. Glotzer, *Particle rearrangements during transitions between local minima of the potential energy landscape of a binary Lennard-Jones liquid*, J. Chem. Phys. **120**, 4404 (2004), doi:10.1063/1.1644538.

[116] N. P. Bailey, T. B. Schrøder and J. C. Dyre, *Exponential distributions of collective flow-event properties in viscous liquid dynamics*, Phys. Rev. Lett. **102**, 055701 (2009), doi:10.1103/PhysRevLett.102.055701.

[117] P. Charbonneau, Y. Jin, G. Parisi, C. Rainone, B. Seoane and F. Zamponi, *Numerical detection of the Gardner transition in a mean-field glass former*, Phys. Rev. E **92**, 012316 (2015), doi:10.1103/PhysRevE.92.012316.

[118] A. Seguin and O. Dauchot, *Experimental evidences of the Gardner phase in a granular glass*, Phys. Rev. Lett. **117**, 228001 (2016), doi:10.1103/PhysRevLett.117.228001.

[119] Y. Jin and H. Yoshino, *Exploring the complex free energy landscape of the simplest glass by rheology*, Nat. Comm. **8**, 14935 (2016), doi:10.1038/ncomms14935.

[120] J.-P. Bouchaud and G. Biroli, *On the Adam-Gibbs-Kirkpatrick-Thirumalai-Wolynes scenario for the viscosity increase in glasses*, J. Chem. Phys. **121**, 7347 (2004), doi:doi:10.1063/1.1796231.

[121] B. Doliwa and A. Heuer, *Finite-size effects in a supercooled liquid*, J. Phys.: Condens. Matter **15**, S849 (2003), doi:10.1088/0953-8984/15/11/309.

[122] K. Kim and T. Munakata, *Glass transition of hard sphere systems: Molecular dynamics and density functional theory*, Phys. Rev. E **68**(2), 021502 (2003), doi:10.1103/PhysRevE.68.021502.

[123] A. Baule, F. Morone, H. J. Herrmann and H. A. Makse, *Edwards statistical mechanics for jammed granular matter*, arXiv:1602.04369.

[124] V. Ogarko and S. Luding, *Equation of state and jamming density for equivalent bi-and polydisperse, smooth, hard sphere systems*, J. Chem. Phys. **136**, 124508 (2012), doi:10.1063/1.3694030.

[125] V. Ogarko and S. Luding, *Prediction of polydisperse hard-sphere mixture behavior using tridisperse systems*, Soft Matter **9**, 9530 (2013), doi:10.1039/C3SM50964H.

[126] M. Ozawa and L. Berthier, *Does the configurational entropy of polydisperse particles exist?*, J. Chem. Phys. **146**, 014502 (2017), doi:10.1063/1.4972525.