An exploratory study of the glassy landscape near jamming

TL;DR

This paper investigates the complex energy landscape of athermal soft spheres near the jamming transition, revealing a hierarchical structure and its implications for thermodynamic properties.

Contribution

It provides the first detailed analysis of the landscape structure of soft spheres near jamming, highlighting the hierarchical organization and its effect on system behavior.

Findings

Number of configurations increases steeply near jamming

Hierarchical landscape structure emerges close to transition

Predicts thermodynamic observables at the transition edge

Abstract

We present the study of the landscape structure of athermal soft spheres both as a function of the packing fraction and of the energy. We find that, on approaching the jamming transition, the number of different configurations available to the system has a steep increase and that a hierarchical organization of the landscape emerges. We use the knowledge of the structure of the landscape to predict the values of thermodynamic observables on the edge of the transition.

Figures

Equations

$$U(\{r\})=\frac{\epsilon}{2}\sum_{i<j}\Bigl{(}1-\frac{r_{ij}}{\sigma_{ij}}\Bigl{)}^{2}\theta{(\sigma_{ij}-r_{ij})}\,,$$

$${\delta }_n={\phi }_{n+1}-{\phi }_n\ge 0$$

$$0\le r_n=\frac{\text{min}\{{\delta }_n,{\delta }_{n-1}\}}{\text{max}\{{\delta }_n,{\delta }_{n-1}\}}\le 1$$

$$\Delta \equiv {\Delta }_{\alpha \beta }^2=\frac{1}{M}\sum _{i=1}^M({\mathbf{x}}_i^{\alpha }-{\mathbf{x}}_i^{\beta }-\mathbf{\delta }{)}^2$$

$$Q\equiv Q_{\alpha\beta}=\frac{1}{M}\sum_{i,j}^{1,M}\Theta\biggl{(}a-\mid\bm{x}_{i}^{\alpha}-\bm{x}_{j}^{\beta}-\bm{\delta}\mid\biggr{)}$$

$$H_{ij}^{\alpha \beta }=\frac{{\partial }^{2}U}{\partial r_i^{\alpha }\partial r_j^{\beta }}=$$

$$={\delta }_{ij}\sum _{k\in \partial i}\left[n_{ik}^{\alpha }{n}_{ik}^{\beta }+\frac{{\epsilon }_{ik}}{r_{ik}}\left(n_{ik}^{\alpha }{n}_{ik}^{\beta }-{\delta }^{\alpha \beta }\right)\right]+$$

$$-{\delta }_{⟨ij⟩}\left[n_{ij}^{\alpha }{n}_{ij}^{\beta }+\frac{{\epsilon }_{ij}}{r_{ij}}\left(n_{ij}^{\alpha }{n}_{ij}^{\beta }-{\delta }^{\alpha \beta }\right)\right],$$

$$\Gamma =\frac{1}{N_{dof}}\sum _{i=1}^{N_{dof}}\log ({\lambda }_i),$$

$$P(\Delta )=\frac{1}{Z^2}\int D{r}^{(1)}D{r}^{(2)}{e}^{-\beta H(r^{(1)})}{e}^{-\beta H(r^{(2)})}$$

$$\times \delta (\Delta -{\Delta }_{12}),$$

$$Z=\sum _{\alpha }{\int }_{r\in \alpha }Dr{e}^{-\beta H(r)}=\sum _{\alpha }{Z}_{\alpha },$$

$$Z_{\alpha }={\int }_{r\in \alpha }Dr{e}^{-\beta H(r)}.$$

$$P(\Delta )=\sum _{\alpha \beta }{w}_{\alpha }{w}_{\beta }\delta (\Delta -{\Delta }_{\alpha \beta }),$$

$$\{\begin{array}{l}{P}_p(\Delta )={\mathcal{N}}_p{\sum }_{\alpha \beta }{e}^{-Np/{\phi }_J^{\alpha }}{e}^{-Np/{\phi }_J^{\beta }}\delta (\Delta -{\Delta }_{\alpha \beta })\\ P_T(\Delta )={\mathcal{N}}_T{\sum }_{\alpha \beta }{e}^{-\beta E_{\alpha }}{e}^{-\beta E_{\beta }}\delta (\Delta -{\Delta }_{\alpha \beta })\end{array}$$

$$P_{T,vib}(\Delta )={\mathcal{N}}_{T,vib}\sum _{\alpha \beta }{e}^{-\beta F_{\alpha }^{(h)}}{e}^{-\beta F_{\beta }^{(h)}}\delta (\Delta -{\Delta }_{\alpha \beta }),$$

$$F^{(h)}_{\alpha}=E_{\alpha}+T\bigl{(}3\Gamma_{\alpha}-3N\log(T)\bigr{)}\,.$$

$$({\mathbf{x}}_i-{\mathbf{x}}_j{)}^2+2({\mathbf{x}}_i-{\mathbf{x}}_j)({\mathbf{\Delta }}_i-{\mathbf{\Delta }}_j)+$$

$$+({\mathbf{\Delta }}_i-{\mathbf{\Delta }}_j{)}^2-\alpha {\sigma }_{ij}^2\le 0$$

$$F(\mathbf{h},ϵ)=-\sum _{i=1}^{N^{\prime }}\{\begin{array}{ll}\ln (h_i)& \text{if hi>ϵ}\\ [\frac{h_i-ϵ}{ϵ}+\ln (ϵ)]& \text{if hi<ϵ}\end{array}$$

$$\Delta =\frac{1}{M}\sum _{i=1}^M({\mathbf{x}}_i^{\alpha }-{\mathbf{x}}_i^{\beta }-\mathbf{\delta }{)}^2$$

$$Q=\frac{1}{M}\sum_{i,j}^{1,M}\Theta\biggl{(}a-\mid\bm{x}_{i}^{\alpha}-\bm{x}_{j}^{\beta}-\bm{\delta}\mid\biggr{)}$$

Full text

I introduction

In this exploratory study, we investigate the properties of the landscape near jamming starting from independent configurations in the same cage. We follow the line of research of Charbonneau et al. (2015a, b); Berthier et al. (2016). We look at the landscape local minima both at jamming and in the overcompressed region. The overcompressed phase near the jamming point is expected to present a similar physics to hard spheres at jamming Charbonneau et al. (2012); Franz et al. (2017); Berthier et al. (2011): in particular, we expect to find the same values for most of the critical exponents. Therefore, in both cases, we expect to find a Gardner-like landscape Kurchan et al. (2013); Charbonneau et al. (2014a, b, 2017); Biroli and Urbani (2018); Scalliet et al. (2019a), in which each cage is broken into a fractal hierarchy of subcages. Beyond the Gardner transition, the replica symmetry becomes continuously broken (full-RSB) Mézard et al. (1987) and the subcages are organized in an ultrametric structure.

The presence of the Gardner transition in finite-dimensional glasses is still a matter of debate: renormalization group studies yield controversial results Moore and Bray (2011); Urbani and Biroli (2015); Charbonneau and Yaida (2017), while numerical simulations suggest that the existence of the Gardner transition may be model dependent DeGiuli et al. (2014); Berthier et al. (2016); Scalliet et al. (2017); Seoane and Zamponi (2018); Liao and Berthier (2019); Scalliet et al. (2019b). However, in hard spheres systems the jamming transition presents universal critical properties which seems to be independent both on the preparation protocol and on the dimension of the system 111Anyway, in finite dimensions new effects might come into play, due, for instance, to the rattlers and the bucklers Charbonneau et al. (2015a)., showing the same features from $d=\infty$ to $d=2$ Charbonneau et al. (2017).

We present a direct inspection of the existence of the Gardner phase in spheres systems in three dimensions near jamming, looking at the properties of the landscape local minima. We study the packing fraction and the energy distribution of the local minima packings, and their relative distance in the configuration space, both at jamming and in the overcompressed region. We look for the ultrametric structure of the landscape at jamming. Furthermore, we predict how the shape of key thermodynamical observables modifies when approaching the transition point both in temperature and in pressure.

Our analysis shows that the system undergoes a roughness transition, which brings about isostaticity on approaching jamming Charbonneau et al. (2014b). The transition is characterized by a steep increase in the number of local minima that, at the jamming point, are organized in an ultrametric way. The deepest minima are close both in terms of the packing fraction and of the distance in the configuration space. They also have large basins of attraction. Moreover, in the overcompressed region, the cumulative of the number of local minima at small energies behaves as a power-law with a packing fraction dependent exponent.

The plan of the paper is as follows. In Sec. II, we describe the general framework and the model. The results of the numerical simulations on the landscape structure at jamming are presented in Sec. III. In Sec. IV, we describe the results in the overcompressed region. Sec. V is devoted to the extrapolation of the landscape properties at finite temperatures and pressures and, in Sec. VI, we summarize our findings and we discuss the perspectives future studies.

II General Framework

This Section is devoted to summarizing the main features of the hard spheres phase diagram in infinite dimensions Parisi and Zamponi (2010); Charbonneau et al. (2017).

The control parameters of a hard spheres system are the pressure and the packing fraction $\varphi$ 222The packing fraction $\varphi$ is defined as the fraction of the volume occupied by the particles. In a monodisperse system, $\varphi=\frac{4}{3}\pi r^{3}\rho$, where $\rho=\frac{N}{V}$ is the density and $r$ is the radius of the spheres.. Different regions of the phase diagram are identified by the behavior of the mean-square-displacement (MSD) ${\Delta(t)=\frac{1}{N}\sum_{i=1}^{N}|\textbf{r}_{i}(0)-\textbf{r}_{i}(t)|^{2}}$. When the pressure is small, the system is ergodic and $\underset{t\to\infty}{\lim}\Delta(t)=\infty$. An equilibrium compression of the liquid can be carried out only up to $\varphi_{d}$, the dynamical transition point. At $\varphi_{d}$, the phase space becomes clustered into an exponential number of glassy states. These clusters are mutually inaccessible and trap the dynamics at infinite times: $\underset{t\to\infty}{\lim}\Delta(t)=\Delta_{\text{liq}}<\infty$. For $\varphi_{g}>\varphi_{d}$, it is possible to define a restricted equilibrium state Franz and Parisi (1995), known as the stable glass phase: the system can completely relax inside the metastable state but long-time diffusion is forbidden. The particles of a stable glass are caged by their neighbors and vibrate around their “amorphous equilibrium positions” in cages whose typical size is $\Delta_{\text{liq}}$.

Compressing further the stable glass up to $\varphi_{G}(\varphi_{g})$, one encounters the Gardner transition where even the restricted equilibrium is lost: the stable glass state breaks into a hierarchical structure of marginal states (landscape marginal stability, LMS). This implies the existence of delocalized soft modes, diverging susceptibilities Biroli and Urbani (2016) and anomalous rheological properties Yoshino and Zamponi (2014).

Finally, compressing the system up to the point where the pressure diverges, the system jams. At jamming, the packings are mechanically rigid and isostatic, meaning that the number of mechanical constraints is equal to the number of degrees of freedom O’Hern et al. (2003). Isostaticity implies that the system is mechanically marginally stable (MMS). Hence, at the jamming point, the number of soft modes is enhanced.

Validating the mean-field picture for finite-dimensional systems would greatly increase the global understanding of the glass transition Debenedetti and Stillinger (2001); Cavagna (2009); Berthier and Biroli (2011). In finite dimensions, the dynamical transition reduces to a crossover because the energy barriers between metastable states are finite. For this reason, it is possible to numerically generate equilibrated glassy configurations even at $\varphi_{g}>\varphi_{d}$ via improved Monte Carlo simulations, known as Swap Monte Carlo Gazzillo and Pastore (1989); Grigera and Parisi (2001); Ninarello et al. (2017). Nevertheless, in conventional dynamics, the time spent by the system in a metastable state can be, to a good extent, considered infinite.

II.1 Methods

Beyond the Gardner transition the number of minima of each cage is expected to diverge in a way that depends exponentially on the number of components of the system Stillinger (1999). Therefore, we choose to restrict our study to 100 3$d$ spheres. In order to enhance the equilibration process, the spheres diameters are drawn from the continuous probability distribution $P(\sigma)\propto\sigma^{-3}$ with $\sigma_{min}/\sigma_{max}=0.45$ Ninarello et al. (2017). The equilibration has been achieved via the constant-pressure Swap Monte Carlo algorithm 333The equilibration beyond the dynamical transition point via the Swap Monte Carlo dynamics has been carried out by the Montpellier research group and L. Berthier, that we warmly thank. to produce 5 glassy configurations at ${\varphi_{g}=0.647}$ ($\varphi_{d}=0.594$), which are in different cages and whose position in the phase diagram is represented by the green square of Fig. 1. We study the landscape near jamming starting from independent glassy configurations in the same cage.

The absence of spatial order and long-time diffusion are key features of glasses. Hence, we firstly study the structure of the 5 initial glassy configurations, computing the radial distribution function Hansen and McDonald (1990), and we investigate their dynamical properties, measuring the MSD as a function of time. Our analysis does not show any sign of crystalline order, phase separation, and long-time diffusion.

II.1.1 Samples Generation and Compression Protocols

Using an NVT Monte Carlo dynamics 444In a sweep of the NVT Monte Carlo dynamics, we propose the displacement of each particle of the system. If the displaced particle overlaps another particle, the movement is refused. We measure the time in units of Monte Carlo sweep., we evolve in time each of the starting configurations. At $t=\tau_{cage}$, the system enters the caging regime, signaled by a plateau in the MSD. Saving the time-evolved configurations every $3\cdot\tau_{cage}$ sweeps, we create a set of independent configurations belonging to the same cage, called clones. The set of clones generated from the same initial configuration is called sample. To prevent long-time diffusion and the breaking of the cage, we periodically restart the dynamics from the starting configuration. Indeed, due to the stochasticity of the Monte Carlo dynamics, we are guaranteed to sample different trajectories at each restart. We end up with a set of $2\times 10^{5}$ independent glassy configurations in the same cage.

Furthermore, via a fast compression 555We move the particles via the NVT Monte Carlo dynamics and, after a fixed number of steps, we increase all the diameters of a factor $\gamma=10^{-3}$. If in the new packing some spheres do overlap, we further move the particles. The procedure stops when the packing fraction of the configuration reaches $\varphi=0.68$. We call this as fast compression because it carries the system out-of-equilibrium. of the starting configuration, we generate a new glassy configuration at higher packing fraction, $\varphi=0.68$. The fast compression does not distort the jamming landscape, but pushes the system into one of the subcages, if, at that packing fraction value, the cage is already broken into subcages. From the new configuration, we produce another set of $2\times 10^{5}$ packings in the same (sub)cage. Therefore, studying the jamming landscape from a sufficiently compressed configuration, we expect to retrieve a subset of the local minima found starting from the equilibrated glassy configuration 666This method can also be employed to study the heights of the landscape barriers. Given a $\varphi>\varphi_{g}$, one can generate many compressed configurations at $\varphi$ via independent fast compressions and, from each of the compressed configurations, generate a sample. Each sample can then be brought to the jamming transition via LP. Repeating the procedure at increasing values of $\varphi$ and studying which landscape local minima survive in different samples, a full map of the landscape can be constructed, including the heights of landscape barriers..

Thus, we follow three compression protocols: (i) taking each sample at $\varphi=0.647$, we bring each clone to the jamming point, via an instantaneous LP compression (App. A); (ii) we repeat the same procedure for each sample at $\varphi=0.68$; (iii) we take $10^{5}$ clones from each sample at $\varphi=0.647$ and we bring them at several packing fraction values beyond the jamming point. Then, we locally relax them via the FIRE algorithm Bitzek et al. (2006).

II.1.2 Studying the Landscape at Jamming

To study the landscape at jamming, we reach the jamming point via the Linear Programming (LP) algorithm Donev et al. (2004); Torquato and Jiao (2010) (App. A). LP works in the uncompressed region 777The uncompressed region is known as the SAT phase in the context of constraint satisfaction problems, while the overcompressed region is called the UNSAT phase Franz et al. (2017)., meaning that overlaps among particles are not allowed during the compression protocol.

In the uncompressed region, the probability of having hard spheres packing at $\varphi$ is $\propto e^{-\frac{Np}{\varphi}}$, where $p$ is a proxy for the pressure and $N$ is the system size. Hence, the jamming transition is reached at $p=\infty$ and the jammed packings maximize $\varphi$.

In this exploratory study, we present the analyses of the jamming landscape of 4 cages, reaching the jamming point both from $\varphi=0.647$ and $\varphi=0.68$.

II.1.3 Studying the Landscape in the Overcompressed Region

To investigate the overcompressed region it is necessary to introduce a soft sphere potential. We choose to employ a harmonic repulsive potential O’Hern et al. (2003)

[TABLE]

where $r_{ij}=|\bm{r}_{i}-\bm{r}_{j}|$ is the distance between the centers of two particles, $\sigma_{ij}=\frac{\sigma_{i}+\sigma_{j}}{2}$ is the sum of their radii and $\epsilon$ is a constant that fixes the energy unit.

We can study how the system behaves in the overcompressed phase when a finite temperature $\beta^{-1}$ is introduced, weighting each energy minimum with its Gibbs measure $\propto e^{-\beta E}$. Temperature is measured in units of $\epsilon/k_{B}$.

III Landscape at Jamming

In this Section, we present our results on the study of the landscape local minima at jamming. Notice that larger values of the jamming packing fraction $\varphi_{J}$ correspond to minima which are deeper in the landscape; while smaller values of $\varphi_{J}$ coincide with minima higher in the landscape.

Most of the jammed packings obtained using the LP algorithm are isostatic (App. B), i.e. they verify $N_{c}=N_{c}^{\text{iso}}\equiv(N-1)d+1$ 888 $N_{c}=(N-1)d+1$ is the isostaticity condition for a finite system under periodic boundary conditions (PBC) Charbonneau et al. (2015a). Here, $N_{c}$ is the number of contacts of the jammed packing, $N$ is the number of particles (excluding the rattlers (App. A)) and $d$ is the dimension of the system., consistently with the mean field solution. Since several works Lerner et al. (2013); Wyart (2012); Goodrich et al. (2012); Hopkins et al. (2013); Charbonneau et al. (2015a) have highlighted the importance of having $N_{c}=N_{c}^{\text{iso}}$ to observe key aspects of jamming criticality, in this study we restrict our analysis to the packings with $N_{c}=N_{c}^{\text{iso}}$.

III.1 Local Minima

Compressing the clones of a sample up to jamming, we find many jammed packings with the same $\varphi_{J}$. This is especially true for high values of $\varphi_{J}$. Furthermore, the packings with the same $\varphi_{J}$ present also the same arrangement of the particles but for a rigid translation 999The rigid translation is due to the translational invariance of the system. The effect has been taken into account in all the results shown in the present study. In particular, in computing the distance $\Delta$ (Eq. 4) and the overlap Q (Eq. 5), we subtracted the displacement of the center of mass among the two configurations., i.e. they represent the same local minimum at jamming.

In each different cage, we find that the jamming local minima coming from $\varphi=0.68$ are a subset of those coming from $\varphi=0.647$, meaning that at $\varphi=0.68$ the starting cage is already broken into subcages.

For each sample, we compute the local minima distribution at jamming. The local minima distributions, shown in Fig. 2, have different supports, depending on the depth of the basin of that cage. The $\varphi_{J}$’s distributions are not self-averaging quantities. The average $\varphi_{J}$ values, merging the data from different cages, are $\overline{\varphi}_{J,\,0.647}=0.6955$ and $\overline{\varphi}_{J,\,0.68}=0.6957$.

Two important features of the landscape are summarized by Fig. 2: many deep local minima are found with high probability and are characterized by similar $\varphi_{J}$ values. We argue that the first feature means that their basins of attraction are large.

We find that the $\varphi_{J}$ distributions may be wider or narrower. The average differences between the highest and the lowest jamming packing fractions are $\overline{\Delta\varphi}_{J,\,0.647}^{\,\text{max-min}}=0.0035$ and $\overline{\Delta\varphi}_{J,\,0.68}^{\,\text{max-min}}=0.0016$. The average ratios between the number of distinct local minima found and the total number of clones in the same sample are, respectively, $f_{0.647}=0.32$ and $f_{0.68}=0.20$.

Ordering in a progressive way the jamming packing fractions, we can define

[TABLE]

If the $\varphi_{J}$’s are uniformly distributed, one expects the level statistics to be Poissonian Oganesyan and Huse (2007) and, thus, $P_{\text{Poisson}}(r)=2/(1+r)^{2}$ with $\langle r\rangle_{\,\text{Poisson}}\simeq 0.386$.

Computing $\langle r\rangle$ in each sample, the Poissonian statistics turns out to be almost valid. Indeed, the average on all the samples is $\langle r\rangle_{0.647}=0.384\pm 0.001$ and $\langle r\rangle_{0.68}=0.382\pm 0.002$.

III.2 Structure in the Configuration Space

Since at jamming, in each sample, we find a huge number of distinct local minima ($\simeq 10^{4}-10^{5}$), here we study in detail the arrangement in the landscape only of a subset of them, the 1000 deepest minima.

First of all, we want to understand how the deepest minimum (i.e., the one with the highest $\varphi_{J}$) of each sample is located in the configuration space. To do so, we compute its square distance

[TABLE]

from the other 200 deepest and the 200 highest minima of the same sample. Here, $\alpha$ and $\beta$ are the minima index and $i$ is the particle index; $M$ is the total number of particles excluding the rattlers of both $\alpha$ and $\beta$ and $\bm{\delta}$ is the distance of the centers of mass of the two minima (Apps. B, C). The $\Delta$ values have to be compared with the average values of the MSD in the caging regime, which are $\overline{\Delta_{\text{liq}}}\equiv\overline{\text{MSD}_{\text{plateau}}^{0.647}}=7.3\times 10^{-3}$ and $\overline{\Delta_{\text{liq}}^{0.68}}\equiv\overline{\text{MSD}_{\text{plateau}}^{0.68}}=2.6\times 10^{-3}$. In Fig. LABEL:sub@fig-delta-conf1-647-deep, is shown the result for cage 1 generating the clones at $\varphi=0.647$. The deepest minimum of this sample appears to be located in a deep well where there are many other very deep minima of the cage (blue line). Moreover, the deepest minimum is faraway from the highest minima (red line). We find almost the same behavior in all the samples.

In Fig. LABEL:sub@fig-delta-1000-200-reciproche are shown the $\Delta$ distributions among the 200 and 1000 deepest minima, having joined the results of all the cages. It is clear that the deepest minima of the landscape in each sample are usually close one to the other.

In a sample, the configuration correspondent to the deepest minimum is usually found many times; but, it is not the one with the greatest number of occurrences. We denote the latter as the most probable minimum.

Performing the same analyses of Fig. 3 on the most probable minimum, we find that it is close to many other highly probable minima; hence, in the landscape at jamming, there is a large basin of attraction made up by large subbasins of attraction.

Furthermore, studying in the same way the properties of the highest minima, we find that they are far apart from all the other minima and do not form a well. They have small basins of attraction. We can imagine them as narrow ponds on the walls of the landscape basins.

Another observable useful to understand the relative position of two local minima at jamming is the overlap Ninarello et al. (2017). It is defined as

[TABLE]

As in Eq. 4, the sum is restricted to the $M$ particles which are non-rattlers in none of the two configurations $\alpha$ and $\beta$; $\bm{\delta}$ is the distance of the centers of mass of $\alpha$ and $\beta$ (Apps. B, C). $\Theta$ is the Heaviside step function. We choose $a=0.03$. $Q\rightarrow 1$ when the number of displaced particles decreases.

$Q$ gives complementary information to $\Delta$. For instance, a high value of $\Delta$ can be due to the presence of only one particle which has very different positions in the two packings $\alpha$ and $\beta$; if it is so, $Q$ has a large value. Otherwise, if many particles are displaced by a small amount, the same $\Delta$ can correspond to a small value of $Q$. In Fig. 4 is shown the behavior of the overlap $Q$ for the same minima used in Fig. LABEL:sub@fig-delta-1000-200-reciproche. In our results, given a couple of minima $\alpha$ and $\beta$ there is not a direct correspondence between a high value of $Q$ and a small value of $\Delta$. However, Fig. LABEL:sub@fig-delta-1000-200-reciproche and Fig. 4 show that the deepest minima of each sample on average are close in terms of both $Q$ (few particles are displaced) and $\Delta$ (small distance). This confirms our description of the landscape structure in terms of deep wells in which the deepest minima are contained.

III.3 Hierarchical Structure

According to the mean-field picture, we expect to find an ultrametric structure of the basins. To verify this hypothesis, we construct the heatmaps 101010The heatmaps are made with R R Development Core Team (2008), using the average method. of the 1000 deepest minima of each sample, using the distance $\Delta$ as the dissimilarity measure.

In Fig. 5 are shown the heatmaps of cage 1 and cage 3 for the samples at $\varphi=0.647$ and $\varphi=0.68$, respectively. To provide a quantitative characterization of the clustering properties of the selected minima, we compute the agglomerative coefficient (AC) for each heatmap 111111$\langle\text{AC}=1-m(i)\rangle$. For each observation $i$, $m(i)$ is the dissimilarity to the first cluster it is merged with divided by the dissimilarity of the merger in the final step of the algorithm.. The average on all the samples is $\overline{\text{AC}_{0.647}}=0.994$ and $\overline{\text{AC}_{0.68}}=0.995$. These high values of AC tell us that the minima have a good clustering structure and, so, present a nearly ultrametric structure.

III.4 Interim Discussion: The Main Features of the Landscape at Jamming

Summarizing, the main features of the landscape at jamming are:

• •

All the samples show a huge number of distinct local minima at jamming. The landscape appears very complex.

• •

The deepest minima of the landscape are found more often than the highest, and so we argue that their basins of attraction are larger.

• •

The dendrograms constructed with the 1000 deepest minima of each sample have good clustering properties (high values of AC). This is evidence in favor of an ultrametric structure of the landscape. In particular, the landscape seems to be more ultrametric when we reach the jamming point coming from $\varphi=0.68$. Notice that a rigorous ultrametric structure can be seen only considering all the minima of a basin. Indeed, when we reach jamming from $\varphi=0.68$, we have pushed the initial equilibrated configuration into a subcage: we are sampling a smaller basin and we expect to detect better the ultrametricity.

• •

The arrangement of the local minima is compatible with several studies on disordered systems Parisi (2003): the deepest minima are close in the configuration space and form a large deep well; instead, the highest minima do not form a well and have small probabilities of being found.

• •

A small square distance $\Delta$ often means a high $Q$ value: there are few displaced particles between the two configurations. Anyway, there are cases in which a lot of particles are slightly moved.

• •

The level statistics at jamming seems to be determined by Poisson statistics, meaning that there is no level-repulsion.

In the next Section, we will study the landscape in the overcompressed region where, near the jamming point, we expect to find a similar physics to hard spheres at jamming Charbonneau et al. (2012); Franz et al. (2017); Berthier et al. (2011).

IV Landscape in the Overcompressed Region

Since the point at which a sample jams depends on the preparation protocol Chaudhuri et al. (2010), we independently determine the jamming packing fraction $\varphi_{J}$ for each cage through compression-decompression cycles Charbonneau et al. (2015a), extrapolating the point at which the energy reaches zero. We then instantaneously inflate each clone to an excess packing fraction $\varphi_{e}=\varphi-\varphi_{J}$, ranging from $\varphi_{e}=5\times 10^{-3}\,\varphi_{J}$ to $\varphi_{e}=100\times 10^{-3}\,\varphi_{J}$, and we locally relax it to the closest energy minimum, i.e., the inherent structure Stillinger and Weber (1982), via a FIRE routine Bitzek et al. (2006).

IV.1 Local Minima

In Fig. 6 we depict the number of different inherent structures, $N_{\text{IS}}$, as a function of $\delta\varphi=\frac{\varphi_{e}}{\varphi_{J}}$. Fig. 6 clearly shows that, as jamming is approached, the number of distinct inherent structures steeply increases, passing from a few units at the highest packing fractions to several thousand, consistently with a roughening of the landscape at jamming.

In Fig. 7(a) and Fig. 8, we show the cumulative distributions of the number of inherent structures as a function of the inherent structure energy 121212Note that $\log{N(E)}$ is related to the configurational entropy of the system., $N(E-E_{0})$ ($E_{0}$ denotes the energy of the deepest minimum of the cage). Fig. 7(a) presents the results averaged over all the samples for several packing fraction values, ranging from $\delta\varphi=5\times 10^{-3}$ to $\delta\varphi=40\times 10^{-3}$, while in Fig. 8 we show the behavior in different samples at $\delta\varphi=5\times 10^{-3}$. Three regimes are present: (i) seemingly exponential at small energies, probably in connection with finite-size effects Kallus (2016), (ii) power-law with a packing fraction dependent exponent $\alpha(\varphi)$ at intermediate energies, (iii) plateau at high energies.

To better understand the origin of the plateau, we look at the probability distribution of the inherent structure energies, in Fig. 7(b). It appears highly concentrated around values close to the deepest minimum of the landscape. It follows that the highest energy minima have very narrow basins of attraction Sciortino (2005); Martiniani et al. (2016) that are visited with very small probability by the equilibrium dynamics, even though they contribute significantly to the configurational entropy of the system Parisi (2003). Therefore, in Figs. 7(a)-8 the saturation of $N(E-E_{0})$ at high energies is likely due to the inability of the algorithm to sample all the minima in the cage. As a consequence, Figs. 7(a)-8 provide a good representation of $N(E-E_{0})$ only at low and intermediate energies.

Hence, Fig. 2 and Fig. 7 represent landscapes with the same features: the deepest minima of the cage are more likely to be found, while the highest have narrow basins of attraction and are hardly sampled by gradient-descent-like algorithms.

Moreover, Fig. 8 gives evidence of the presence of the Gardner physics: even though the cumulative distributions $N(E-E_{0})$ of the 5 starting cages have the same qualitative behavior, the total number of detected states has very large sample-to-sample fluctuations. This leads to the loss of the self-averaging property when the jamming transition is approached.

IV.2 Harmonic Properties of the Local Minima

The analysis of the harmonic properties of the minima is carried out by studying the Hessian matrix:

[TABLE]

with $\alpha,\beta=x,y,z$ are position vector components, $\varepsilon_{ij}=\sigma_{ij}-r_{ij}$ is the overlap between two spheres, $\vec{n}_{ij}=(\vec{r}_{j}-\vec{r}_{i})/r_{ij}$ is a unit vector, both $\delta_{ij}$ and $\delta^{\alpha\beta}$ are Kronecker deltas, $\delta_{\langle ij\rangle}$ indicates a contact between a pair of particles, and $\partial i$ denotes the set of neighbors of $i$.

Matrix IV.2 has 3 zero modes corresponding to the global translational invariance of the system plus a set of clone-dependent zero modes connected to the presence of rattlers, which are particles that are not part of the force network Charbonneau et al. (2015a), and therefore play no role in the thermodynamics of the system 131313The rattlers are removed before computing the Hessian..

Given the eigenvalues of the Hessian $\{\lambda\}$, we define the logarithm of the curvature as:

[TABLE]

where $N_{dof}=3(N-1-n_{ratt})$ is the total number of degrees of freedom, $n_{ratt}$ being the number of rattlers.

Fig. 9 tracks the evolution of the curvature as a function of the packing from $\delta\varphi=35\times 10^{-3}$, where a Gaussian profile first appears, to $\delta\varphi=5\times 10^{-3}$. At larger packing fractions, the Gaussian profile is not present.

Gaussianity implies that large finite-size effects $O(1/\sqrt{N})$ take place approaching the jamming transition. We can also notice that the distribution shifts towards smaller values of $\Gamma$ as the packing fraction is decreased, compatibly with the onset of mechanical marginal stability at jamming. Indeed, smaller values of $\Gamma$ are associated with larger basins of attraction of the inherent structures. Error bars are computed with the jackknife method Amit and Martin-Mayor (2005) on the 5 samples and their height is a consequence of the large sample-to-sample fluctuations, as already highlighted in Sec. IV.1.

V Thermodynamical Properties of the Landscape

Building on the knowledge of the landscape structure, we try to predict how physical observables would change on approaching the jamming transition at finite temperature in the overcompressed region and at finite pressure in the undercompressed one. In the first case, we start our analysis from the potential energy landscape of Sec. IV, while in the second case we start from the landscape at jamming studied in Sec. III.

For two systems with the same disorder (in this context, two clones) we define the probability distribution of their distance (see Eq. 4 and App. C) as

[TABLE]

where $Dr$ represents the integration over all the configurational degrees of freedom and $Z=\int Dre^{-\beta H({r})}$ is the partition function.

In presence of ergodicity breaking, the Gibbs measure can be decomposed over pure states Mézard et al. (1987); Parisi (1988); Castellani and Cavagna (2005), yielding

[TABLE]

where

[TABLE]

and $\alpha$ is the pure state index. Whenever Eq. 9 is valid, we can rewrite Eq. 8 by taking explicitly into account the contribution of different pure states:

[TABLE]

where $w_{\alpha}=Z_{\alpha}/Z$ is the statistical weight of the state $\alpha$. The sum runs over all the minima and includes the case $\alpha=\beta$. We define its average over different realizations of the disorder as $\overline{P(\Delta)}=\overline{\sum_{\alpha\beta}w_{\alpha}w_{\beta}\delta(\Delta-\Delta_{\alpha\beta})}$.

In the low temperature limit, we can expand the integral in Eq. 10 around the energy $E_{\alpha}$ of the correspondent inherent structure getting, to the first order, $Z_{\alpha}\sim e^{-\beta E_{\alpha}}$. Analogously, by repeating the same analysis for the uncompressed region, where the control parameter is the pressure, we get $Z_{\alpha}\sim e^{-Np/\varphi^{\alpha}_{J}}$ where $p$ is a proxy for the pressure, $N$ is the system size, and $\varphi^{\alpha}_{J}$ is the jamming packing fraction of state $\alpha$.

Substituting these approximations into Eq. 11, we obtain that the distance distribution is given by

[TABLE]

where $\mathcal{N}_{p}$ and $\mathcal{N}_{T}$ are normalization constants, in both cases equal to the square of the partition function.

Qualitatively, $P(\Delta)$ displays a delta peak at $\Delta=0$, coming from the self-part of the summation ($\alpha=\beta$), and a continuous band that originates from the exchange part ($\alpha\neq\beta$). The former can be used to predict a crossover at which the weight of the entropic term becomes more important than the energetic one. The presence of the continuous band from the exchange part is in agreement with the $d=\infty$ picture which predicts the presence of a huge number of states at jamming due to the Gardner transition. The continuous band contains the details of the organization of the states.

We define the coefficient of the self-part as $W_{T,p}^{0}=\sum_{\alpha}w^{2}_{\alpha}$. From Eq. 12, we see that $W^{0}_{T,p}$ becomes more and more important when $T$ decreases or $p$ increases. Indeed, if $T=0$ or $p=\infty$ the only contribution comes from the deepest minimum of the basin with itself and $W^{0}_{T,p}=1$. On the contrary, at high $T$ or low $p$, we expect to have also the continuous part of the distribution $P(\Delta\neq 0)$. When $T=\infty$ or $p=0$ all the minima contribute with the same weight and $W^{0}_{T,p}=N_{minima}^{-1}$.

In the overcompressed phase, we can easily take into account the harmonic corrections to get

[TABLE]

where $F^{(h)}_{\alpha}$ denotes the free energy in the harmonic approximation

[TABLE]

V.1 Finite Temperature

In Fig. 10(a), we show the behavior of $W^{0}_{T}$, the coefficient of the self-part, as a function of the temperature. We observe the presence of a crossover between the low and the high temperature behavior, where we go from a situation in which the thermodynamic properties are determined by a finite number of states, whose energy is close to that of the deepest minimum, to a regime in which all the minima are equally important.

We define the crossover temperature $T^{*}$ via the condition $W^{0}_{T^{*}}=1/3$. The result is shown in Fig. 10(b): we see that, as the packing fraction decreases, and consequently the number of states increases, smaller temperatures are needed to overcome the entropic contribution. From Fig. 10(a) we also see that, decreasing $\delta\varphi$, the crossover becomes steeper.

We notice that, since this crossover involves just a finite number of states of low energy, it is not affected by the plateau at high energies found in the cumulative distributions of the number of inherent structures (Fig. 7); in other words, the presence in the cage of a number of inherent structures larger than what can be sampled by our algorithm, affects only Fig. 10(a) at very small values of $\beta$, but leaves the qualitative picture and the determination of $\beta^{*}=1/T^{*}$ intact.

Now, let us focus on the continuous branch of the displacement distribution $P_{T}(\Delta)$. In Fig. 11(a) we show the distributions at $\beta=0$ and $\delta\varphi=5\times 10^{-3}$ for all the samples. We find that the distribution has support over a continuous range of $\Delta$ values, covering nearly 5 orders of magnitude, compatibly with the existence of many states as predicted by the mean-field theory. $P_{T}(\Delta)$ exhibits a nearly flat behavior. We notice that, as already highlighted in Fig. 8, there are large sample-to-sample fluctuations, but the overall shape of the probability distribution is preserved.

In Figs.  11(b), 12(a) we consider how the distribution $\overline{P_{T}(\Delta)}$ modifies when the packing fraction or the temperature is changed. Here, $\overline{\bullet}$ indicates the average over all the samples. We see that the probability of finding two local minima of the landscape at smaller $\Delta$ increases when the temperature is lowered or the packing fraction is increased. The former can be easily explained by observing that the Boltzmann weight $e^{-\beta E}$ strongly suppresses the pairs of configurations with the highest energy, that in our system also correspond to those at a larger distance. The dependence on the packing fraction can be understood by noticing that, given any minimum, the number of neighboring minima increases with the distance. Since, decreasing $\delta\varphi$, the landscape roughens and the number of local minima increases, there will be more configurations at large distances giving rise to a peak at high $\Delta$ values.

We also notice that in Fig. 12(a) the probability distribution is rather insensible to temperature changes up to the crossover temperature, where the states close to the ground state start to dominate the partition function.

Finally, in Fig. 12(b) we present $\overline{P_{T,vib}(\Delta)}$ for several temperatures, where we have included the harmonic terms as prescribed by Eqs. 13-14. Their effect is to make more uniform the behavior of the distribution at different $\beta$’s, by enhancing small $\Delta$ probabilities. However, the qualitative picture of $\overline{P_{T}(\Delta)}$ is not altered.

V.2 Finite Pressure

In analogy with Secs. III.2-III.3, in the uncompressed region we compute $P_{p}(\Delta)$ considering only the 1000 deepest minima of each sample.

In Fig. 13 we present the behavior of $W^{0}_{p}$ as a function of the pressure for two samples reaching the jamming point from $\varphi=0.68$. As in Sec. V.1, this plot allows us to determine the crossover from the entropy ruled and the ground state ruled landscape.

In Fig. 14 is shown the $P_{p}(\Delta)$ for the same samples of Fig. 13. In both cages, the probability $P_{p}(\Delta)$ at low $p$ is almost flat. Furthermore, as expected, the distribution presents isolated peaks when $p$ increases.

How the distribution changes with pressure is different in each cage because of the high heterogeneity of the landscape basins. For instance, in cage 1 at $p=10^{6}$ $P(\Delta\neq 0)\neq 0$ because the sum in Eq. 12 has a residual contribution from the second deepest minimum, which has $\varphi_{J}$ very close to the one of the ground state. At low pressure, all the minima contribute.

VI Conclusions and Outlook

This exploratory study aims to investigate the structure of the landscape near the jamming transition in three-dimensional soft spheres systems, in the light of the existence of the Gardner transition in the exact solution of hard spheres in $d=\infty$. We followed two different but complementary approaches.

In the first one, we brought the system to the jamming point from the uncompressed phase, starting from two different packing fraction values. This approach pointed out the presence of subcages and of a huge number of local minima at jamming, in agreement with the fractal nature of the landscape in $d=\infty$. We found that, in many cases, the deepest minima are close in the landscape and located in a deep well. They usually have large basins of attraction. Instead, the minima corresponding to smaller values of $\varphi_{J}$ are scattered in all the jamming landscape and they usually have small basins of attraction. The ultrametric structure of the landscape and the clustering properties of the jammed configurations seem to be verified by the high value of the AC coefficients, even though tested on a small bunch of minima. Furthermore, our study shows the possibility of using in an efficient way a Linear Programming (LP) algorithm to reach the jamming point in finite-dimensional hard spheres system.

On approaching the jamming transition from the overcompressed region, we found a new kind of criticality: the number of inherent structures as a function of the packing fraction seems to diverge in a power-law fashion when $\varphi\rightarrow\varphi_{J}$, and the cumulative of the number of inherent structures at low energies behaves as a power-law in a wide range of packing fraction values, with a $\varphi$-dependent exponent. Moreover, the onset of the critical behavior is accompanied by the emergence of large sample-to-sample fluctuations, which can be seen both in the cumulative of the inherent structures and in the distribution of the curvatures from the Hessian matrix. Furthermore, the latter displays a Gaussian behavior.

The study of the distribution of the distances both in the uncompressed and overcompressed region, as a function of $p$ and $T$ respectively, revealed the existence of a crossover between ground state ruled and entropy ruled landscapes. We found the presence of an isolated peak at zero distance and of a continuous part of the distribution at larger distances. The latter displays an almost flat behavior across several values of pressure, temperature and packing fraction. The existence of a continuous part extended over five orders of magnitude shows that the landscape is made up of a huge number of states which have a continuous distribution of distances, in agreement with the $d=\infty$ picture. The inclusion of harmonic corrections does not change in a significant way the behavior of the distribution that is, therefore, correctly captured already at a purely energetical level.

The pictures of the landscape that emerge from the study at jamming and in the overcompressed region are in agreement and consistent with the mean-field predictions. However, we observe that quantities like $P(\Delta)$ are very atypical and, therefore, a conclusive analysis should be done over a larger number of local minima and averaging over several different cages. Indeed, increasing the number of starting metabasins, would allow to perform better disorder averages and have deeper control over sample-to-sample fluctuations.

Moreover, repeating the same analysis with systems of different sizes and dimensions would allow us to determine the entity of finite-size effects and to properly extrapolate the thermodynamic limit of the observables.

Furthermore, since it is known Berthier et al. (2016); Charbonneau et al. (2015b) that the starting point of the compression plays an important role in determining the strength of the Gardner transition, it would also be important to repeat the analysis by starting from different points along the equation of state and identify, in a clear way, the universal features of the behavior at jamming

VII Acknowledgements

We thank L. Berthier for providing the initial configurations and useful exchanges. We thank P. Charbonneau, G. Folena, C. Scalliet, A. Scardicchio, G. Sicuro and G. Tsekenis for many useful discussions related to this work. The research has been supported by the Simons Foundation (grant No. 454949, G. Parisi).

Appendix A The Linear Programming Algorithm

A Linear Programming algorithm is a method to solve linear optimization problems. An LP problem is the problem of minimization or maximization of a linear function subjected to linear constraints Dantzig and Thapa (2006a, b). Our LP problem has the objective to maximize the particles’ radii. The inflation rate has to be equal for all the particles to leave unchanged the degree of polydispersity. The variables are the inflation rate $\alpha$ and the particles’ displacements $\bm{\Delta}$ ($3(N-1)$ scalar variables) 141414We fix the position of a sphere to avoid the rigid translation of all the system.. The constraints of the problem are

[TABLE]

where $\sigma_{ij}=\frac{\sigma_{i}+\sigma_{j}}{2}$ is the sum of the radii, $\mid\Delta^{x,y,z}_{i}\mid\leq c$ and $0\leq\alpha\leq c^{\prime}$ for $c$ and $c^{\prime}$ reasonable values with respect to the linear dimensions of the system.

To use LP we need a linear problem. We can neglect the term $(\bm{\Delta}_{i}-\bm{\Delta}_{j})^{2}$ supposing that the displacements are small. This assumption is asymptotically justified because, after some iterations, the increase of the particles’ radii will reduce the magnitude of the possible displacements. In the end, the problem can be written as a maximization problem with objective function $\alpha$.

For each couple of particles, the constraint in Eq. 15 is satisfied and the inflation $\alpha$ is maximized when the spheres go as far as possible along the direction orthogonal to the constraint and they are inflated until the constraint is saturated (Fig. 15). Since each sphere has many neighbors and all the constraints have to be satisfied simultaneously, the displacements would not be only along this orthogonal direction and the constraint will not be saturated after only one iteration. We need some iterations ($\sim\mathcal{O}(10)$) to reach the jamming point.

The LP algorithm we used is in the GLPK library, version 4.55 151515GLPK uses the symplectic method. The Reference Manual can be found here glp ..

Appendix B Rattlers

In jammed hard spheres packings in finite dimensions some particles are not part of the force network. They are called rattlers and they have $Z<d+1$ contacts. Typically, their neighbors are part of the force network and form a cage in which the rattlers can freely move. In this study, we compute the observables excluding the rattlers. The distance between two packings is computed as the distance of their backbones intersection (i.e. we excluded the rattlers in common and not in common between the two packings).

In a polydisperse system, the number of rattlers is bigger than in a monodisperse one. Indeed, it is easier for a particle with a small radius to be trapped by the others. At jamming, we find that the average number of rattlers per jammed configuration is $\overline{N}_{rattlers}=15$.

With the LP algorithm of the GLPK library, version 4.55, after each iteration, one finds that most of the rattlers computed on the LP constraints have $3$ contacts. This means that they are leaned against the walls of the virtual box of the constraints of the LP problem. This is not what one would expect because rattlers should be free to move and, so, they typically have 0 contacts. Furthermore, the final jammed packings are not isostatic on the physical contacts. We believe this result is due to a computer rounding on the forces acting on the rattlers, which moves them arbitrary.

To solve this problem, we define a smooth potential which, minimized at each LP iteration, moves the rattlers in the middle of their virtual box. This potential is defined for each group of nearest rattlers. In the case of an isolated rattler, the potential is defined as

[TABLE]

where the sum is over all the $N^{\prime}$ non-rattlers of the packing, $h_{i}$ is the gap between the considered rattler and the particle $i$ and $\epsilon=10^{-12}$. The generalization to the case of a group of rattlers is straightforward. For the minimization of the potential, we use the nmsimplex2 algorithm in the GSL library gsl .

In this way, the resulting jammed configurations are isostatic in the physical gaps, except for a small percentage of non-isostatic packings 161616The percentage of non-isostatic packings is $9\%$. The percentage is much smaller excluding from the average the second cage (both starting from $\varphi=0.647$ and $\varphi=0.68$) and becomes $2\%$. The jammed configurations of cage 2 are more difficult to reach with our LP algorithm. This might be due to the presence of more rattlers than in the other cages.. Although the potential expedient, we find that the rattlers’ positions in the jammed packings are random. Indeed, looking at the rattlers in different clones that have reached the same local minimum, they have different positions. Hence, to avoid spurious contributions from the rattlers we compute the distances between two packings as the distance of the backbones intersection (see App. C).

Appendix C Distance among Jammed Packings

The distance of two configurations is defined as (Eq. 4)

[TABLE]

while the overlap is (Eq. 5)

[TABLE]

$\alpha$ and $\beta$ are the minima indeces, $i$ is the particle index, $M$ is the number of spheres in the backbones intersection, i.e. $M$ is the total number of particles which are non-rattlers neither in packing $\alpha$ nor in $\beta$, and $\bm{\delta}$ is the distance of the centers of mass of $\alpha$ and $\beta$.

Indeed, as highlighted in Sec. III.1, when we compute observables such as $\Delta$ and $Q$ we need to avoid spurious effects coming from the system translational invariance. Hence, to compute $\Delta$ and $Q$ we perform a rigid translation in order to have $\alpha$ and $\beta$ with coincident centers of mass: we compute the average displacement of the particles among the two configurations as $\bm{\delta}=\frac{1}{M}\sum_{i=1}^{M}(\bm{x}_{i}^{\alpha}-\bm{x}_{i}^{\beta})$ and we subtract $\bm{\delta}$ in each term of the sum $\bm{x}^{\alpha}_{i}-\bm{x}^{\beta}_{i}$. From Eqs. 16, 17 it is clear that, in this way, packings representing the same local minima have $\Delta=0$ and $Q=1$.

In the LP protocols, before doing any analysis, we computed the distance distribution among the 1000 deepest minima in each sample. Studying them in log scale, we found that at very small distances the distributions present some isolated peaks. Moreover, the same feature is present in the log distribution of the packing fraction differences $\Delta\varphi_{J}$. After a careful analysis, we argued that configurations with both a very small $\Delta$ and a very small $\Delta\varphi_{J}$ should be considered degenerate because they represent the same local minimum of the landscape. Looking at those log distributions, we fixed two threshold values $\Delta^{*}\leq 10^{-9}$ and $\Delta\varphi_{J}^{*}\leq 10^{-9}$. Whenever two configurations have $\Delta\leq\Delta^{*}$ and $\Delta\varphi_{J}\leq\Delta\varphi_{J}^{*}$, we considered them as degenerate.