# Jamming is a first-order transition with quenched disorder in amorphous
materials sheared by cyclic quasistatic deformations

Yue Deng  CAS Key Laboratory of Theoretical Physics, Institute of Theoretical
Physics, Chinese Academy of Sciences, Beijing 100190, China  School of
Physical Sciences, University of Chinese Academy of Sciences, Beijing 100049,
China  Deng Pan  CAS Key Laboratory of Theoretical Physics, Institute of
Theoretical Physics, Chinese Academy of Sciences, Beijing 100190, China
Yuliang Jin<EMAIL_ADDRESS>CAS Key Laboratory of Theoretical
Physics, Institute of Theoretical Physics, Chinese Academy of Sciences,
Beijing 100190, China  School of Physical Sciences, University of Chinese
Academy of Sciences, Beijing 100049, China  Center for Theoretical
Interdisciplinary Sciences, Wenzhou Institute, University of Chinese Academy
of Sciences, Wenzhou, Zhejiang 325001, China

###### Abstract

Jamming is an athermal transition between flowing and rigid states in
amorphous systems such as granular matter, colloidal suspensions, complex
fluids and cells. The jamming transition seems to display mixed aspects of a
first-order transition, evidenced by a discontinuity in the coordination
number, and a second-order transition, indicated by power-law scalings and
diverging lengths. Here we demonstrate that jamming is a first-order
transition with quenched disorder in cyclically sheared systems with
quasistatic deformations. Particle models are simulated in two and three
dimensions, which undergo reversible-irreversible transitions under cyclic
shear. The ensemble of configurations, with partially crystallized and fragile
states excluded, is generated in the irreversible phase, where the system is
stationary and particles are diffusive. Detailed scaling analyses are
performed on the distribution of coordination numbers. The fluctuation of the
jamming density in finite-sized systems has important consequences on the
finite-size effects of various quantities, resulting in a square relationship
between disconnected and connected susceptibilities, a key signature of first-
order transitions with quenched disorder. This study puts the jamming
transition into the category of a broad class of transitions in disordered
systems, including the first-order transition in the random-field Ising model,
brittle yielding of amorphous materials, and melting of ultra-stable glasses.

Introduction.

Jamming in athermal particles is a paradigm of transitions between fluids and
amorphous solids Makse  _et al._ (2000); O’hern  _et al._ (2003); Liu and
Nagel (2010); van Hecke (2009); Pan  _et al._ (2023a), with a deep connection
to the glass transition in thermal systems Parisi and Zamponi (2010);
Charbonneau  _et al._ (2017); Parisi  _et al._ (2020). Recent studies have
revealed extremely rich features in the jamming phenomenon, but the nature of
the jamming transition remains inconclusive.

(i) Mechanical marginality, related diverging length scales and power-law
scalings. At the jamming transition density (volume fraction) $\varphi_{\rm
J}$, the isostatic condition needs to be satisfied for the coordination number
$Z$ (average number of contacts per particle): $Z=Z_{\rm iso}=2d$ in
frictionless, infinite systems, where $d$ is the dimensionality Maxwell
(1870); Alexander (1998). Isostaticity implies that at $\varphi_{\rm J}$ the
system is marginally stable, inspiring the search for diverging length scales.
According to the “cutting argument”, the isostaticity gives rise to a
diverging isostatic length scale at jamming, $l^{*}\sim\Delta
Z^{-1}\sim(\varphi-\varphi_{\rm J})^{-1/2}$, below which mechanical stability
of the bulk system is affected by the cutting boundaries Wyart (2005); wya
(2005); Silbert  _et al._ (2005), where $\Delta Z=Z-Z_{\rm iso}$ is the excess
coordination number. The effective medium theory gives a scattering length
scale diverging at $\varphi_{\rm J}$, $l_{c}\sim\Delta
Z^{-1/2}\sim(\varphi-\varphi_{\rm J})^{-1/4}$, below which continuum
elasticity breaks down Wyart (2010); DeGiuli  _et al._ (2014a); Lerner  _et
al._ (2014). Other related diverging length scales include the transverse
wavelength, $\xi_{T}\sim(\varphi-\varphi_{\rm J})^{-\nu_{T}}$ with
$\nu_{T}\approx 0.24$ Silbert  _et al._ (2005), and the longitudinal
wavelength, $\xi_{L}\sim(\varphi-\varphi_{\rm J})^{-\nu_{L}}$ with
$\nu_{L}\approx 0.48$ Silbert  _et al._ (2005).

At jamming, the marginal stability analysis provides relationships between the
exponents $\theta$ appearing in the power-law distribution of weak inter-
particle forces $P(f)\sim f^{\theta}$ and $\alpha$ in the distribution of
small inter-particle gaps $P(h)\sim h^{-\alpha}$: $\theta=1/\alpha-2$ for
extensive modes Wyart (2012) and $\theta=1-2\alpha$ for localized buckling
modes Lerner  _et al._ (2013); DeGiuli  _et al._ (2014b). Above jamming, the
marginal boundary between unstable and stable phases is defined by a scaling
relation, $\Delta Z\sim(\varphi-\varphi_{\rm J})^{1/2}$ Wyart (2005); Müller
and Wyart (2015). Other scalings have been established for over-jammed systems
near $\varphi_{\rm J}$ Makse  _et al._ (2000); O’hern  _et al._ (2003);
Goodrich  _et al._ (2016); Liu and Nagel (2010); Baule  _et al._ (2018); Pan
_et al._ (2023a). For example, a relationship between the shear modulus $G$
and $\Delta Z$ can be derived by microscopic elastic theories, $G\sim\Delta
Z\sim(\varphi-\varphi_{\rm J})^{1/2}$ Wyart (2005); Zaccone and Scossa-Romano
(2011); Pan  _et al._ (2023b).

(ii) Hyperuniformity and associated diverging length scales. Recent studies
reveal the spatial distribution of the single-particle contact number,
$Z_{i}$, is hyperuniform at jamming, implying a diverging hyperuniform
correlation length, $\xi_{\rm H}\sim\Delta Z^{-\nu}$ Hexner  _et al._ (2018,
2019). The value of exponent $\nu$ appears to depend on the way to extract the
correlation length and dimensionality Hexner  _et al._ (2018).

(iii) Gardner glass phase, landscape marginality and associated criticality.
The marginality at jamming has been established by an independent approach
within the framework of replica symmetry breaking Charbonneau  _et al._
(2014); Parisi  _et al._ (2020); Charbonneau  _et al._ (2015a); Berthier  _et
al._ (2016, 2019); Li  _et al._ (2021); Urbani  _et al._ (2023). Unjammed hard
sphere glasses undergo a Gardner transition where the free energy landscape
becomes fractal and marginal, and the caging susceptibility diverges Berthier
_et al._ (2016). The entire Gardner phase, including the jamming limit, is
critical. In other words, the caging correlation length $\xi_{\rm
G}\sim(\varphi_{\rm J}-\varphi)^{0}$ remains infinite near jamming. The mean-
field replica theory predicts the values of exponents in the weak force and
small gap distributions, $\alpha=0.41269$ and $\theta=0.42311$, coinciding
with the relationship $\theta=1/\alpha-2$ given by the mechanical marginal
stability analysis Charbonneau  _et al._ (2014, 2015b); Wang  _et al._ (2022).
The theory provides an additional scaling relationship between the cage size
$\Delta$ and the entropic pressure $p$, $\Delta\sim p^{-\kappa}$ with
$\kappa=1.41574$. The exponents appear to be independent of the dimensionality
$d$ for $d\geq d_{\rm u}$ Charbonneau  _et al._ (2014, 2015b), where $d_{\rm
u}=2$ is the conjectured upper critical dimension Goodrich  _et al._ (2012).

(iv) Criticality in shear rheology. The criticality of the jamming transition
is suggested by scaling analysis of the rheological data obtained in finite-
rate shear simulations of flowing states near $\varphi_{\rm J}$ Olsson and
Teitel (2007, 2011); Kawasaki  _et al._ (2015). Combing the power-law
divergence of the viscosity, $\eta=(\varphi_{\rm J}-\varphi)^{-\beta}$ and the
vanishing of yield stress, $\sigma_{\rm Y}\sim(\varphi-\varphi_{\rm
J})^{\Delta}$, a critical scaling function is proposed. A diverging
rheological correlation length can be extracted from velocity correlations
Olsson and Teitel (2007) or non-affine displacements Heussinger and Barrat
(2009), $\xi_{\rm R}\sim(\varphi_{\rm J}-\varphi)^{-\nu}$, with $\nu=0.6-1.1$
Olsson and Teitel (2007, 2011); Heussinger and Barrat (2009); Vagberg  _et
al._ (2011). The criticality seems to be at odds with the hyperuniformity
discussed in (ii) – in the thermodynamical limit, the fluctuations diverge in
the former and vanish in the latter.

(v) Discontinuity in the coordination number. The coordination number $Z$,
which is considered as an order parameter of the jamming transition, jumps
discontinuously from $Z=0$ below jamming, to $Z\geq Z_{\rm iso}$ above, in
quasi-static compression or decompression Makse  _et al._ (2000); O’hern  _et
al._ (2003). Apparently, as the signature of a first-order transition, this
discontinuity is inconsistent with the viewpoints of a continuous transition
described above.

It is clear that various diverging length scales have been suggested
throughout the literature. However, none of the lengths in (i-iv) can explain
the finite-size scaling behavior of the jamming fraction $F_{\rm
J}(\varphi,N)$ O’hern  _et al._ (2003); Bertrand  _et al._ (2016); Baity-Jesi
_et al._ (2017); Jin and Yoshino (2021): the data can be reasonably collapsed
by a master curve, $F_{\rm J}(\varphi,N)=\mathcal{F}[(\varphi-\varphi_{\rm
J})N^{1/2}]$, valid for both compression and shear jamming in two (2D) and
three (3D) dimensions. In this study, we show that this scaling can be fully
explained by a first-order transition scenario of the jamming transition with
quenched disorder. The form $(\varphi-\varphi_{\rm J})N^{1/2}$ originates from
the disorder-induced fluctuation of the jamming density itself in finite-sized
systems of $N$ particles, which follows the standard central limit theorem.
Thus this finite-size scaling is independent of isostaticity, marginality,
criticality and hyperuniformity.

Three important differences between previous approaches (i-v) and ours shall
be denoted. First, in (i-iii), the jamming limit is approached from one side
only, i.e., the over-jammed side ($\varphi>\varphi_{\rm J}$) in (i) and (ii),
and the unjammed side ($\varphi<\varphi_{\rm J}$) in (iii). Here we consider a
well-defined ensemble including both over-jammed and unjammed states, whose
ratio is essential in the scaling analysis. Note that once an ensemble average
is taken, the (v) discontinuity in $Z$ turns into a smooth function, and thus
a finite-size analysis becomes essential to see the asymptotic behavior in the
thermodynamic limit.

Second, in conventional compression quenching protocols, the generated
ensemble depends on the initial conditions Chaudhuri  _et al._ (2010); Ozawa
_et al._ (2012); Jin and Yoshino (2021) and the basins of attraction Bertrand
_et al._ (2016); Ashwin  _et al._ (2012); Martiniani  _et al._ (2017). Here we
instead consider an ensemble prepared by cyclic shear, where the states are
sampled by well-controlled dynamics similar to that in thermal systems.
Recently, the response of amorphous assemblies of particles to cyclic shear
has attracted great interest, due to the presence of a nonequilibrium phase
transition, called the reversible-irreversible (RI) transition Pine  _et al._
(2005); Corte  _et al._ (2008); Schreck  _et al._ (2013); Nagasawa  _et al._
(2019); Das  _et al._ (2020). Interestingly, the jamming transition lies in
the irreversible regime where particle trajectories are asymptotically
diffusive Nagasawa  _et al._ (2019); Das  _et al._ (2020). In this study, we
restrict our ensemble within the irreversible phase.

Third, in our ensemble, we carefully exclude partially crystallized and
fragile states with $Z<Z_{\rm iso}$ that are sensitive to mechanical
perturbations or protocol parameters, and regard them as unjammed states with
$Z=0$. If such states are included, the discontinuity in $Z$ diminishes and
the jamming transition looks continuous, similar to the results reported in
Ref. Heussinger and Barrat (2009) obtained by uniform shear. Thus we expect
the impactibility between the criticality viewpoint in (iv) and our first-
order picture originating from finite-rate effects. For finite-time scales,
the existence of transient states with $Z<Z_{\rm iso}$ can lead to a
continuous jamming transition Olsson and Teitel (2007, 2011); Heussinger and
Barrat (2009). However, our results suggest that, after a sufficiently long
time, any configurations with $Z<Z_{\rm iso}$ would eventually relax to
unjammed states with zero energy and inter-particle contacts.

An ensemble generated by cyclic athermal quasi-static shear.

We apply cyclic athermal quasistatic shear (CAQS) to standard models of soft,
frictionless particles in 2D and 3D (see Models and Shear protocol in
Methods). We present 2D data in the main text and 3D data in the Extended
Data. For the 2D model, previous studies report a jamming density (J-point
density) $\varphi_{\rm J}\approx 0.842-0.843$ O’Hern  _et al._ (2002); Vagberg
_et al._ (2011). The phase diagram of RI transitions near $\varphi_{\rm J}$ is
plotted in Fig. 1a. The RI transition is characterized by the one-cycle
displacement averaged over particles Nagasawa  _et al._ (2019),

$\delta
r(t)=\frac{1}{N}\sum_{i=1}^{N}\left|\boldsymbol{\mathrm{r}}_{i}(t+1)-\boldsymbol{\mathrm{r}}_{i}(t)\right|,$
(1)

and the mean-squared displacement (MSD),

$\langle\Delta
r^{2}(t)\rangle=\frac{1}{N}\left\langle\sum_{i=1}^{N}\left|\boldsymbol{\mathrm{r}}_{i}(t)-\boldsymbol{\mathrm{r}}_{i}(0)\right|^{2}\right\rangle,$
(2)

where $t$ is the number of cycles playing a similar role as the time in
thermal systems, $\boldsymbol{\mathrm{r}}_{i}(t)$ the position of particle $i$
at time $t$ and zero strain $\gamma=0$, and $\langle x\rangle$ the average
over $N_{\rm s}$ samples.

The RI dynamics near $\varphi_{\rm J}$ are systematically studied in Ref.
Nagasawa  _et al._ (2019), according to which $\delta r(t)$ displays two-step
relaxation behavior typically appearing in glassy systems. For $\tau_{\rm
R}<t<\tau_{\rm L}$, $\delta r(t)$ develops a plateau at $\delta r_{\rm s}$. In
the interested irreversible phase, $\delta r_{\rm s}>0$, $\tau_{\rm R}\sim 0$
and $\tau_{\rm L}>t_{\rm max}$ with $t_{\rm max}$ the maximum simulation time,
suggesting that the system reaches a stationary state (see Fig. 1b). In
addition, the MSD in the irreversible phase shows typical diffusive behavior
$\langle\Delta r^{2}(t)\rangle\sim t$ (see Fig. 1c). In contrast, in the
reversible phase, $\delta r_{\rm s}=0$, which means that the system is
“absorbed” into an invariant state. In practice, one defines $\delta
r_{\infty}=\delta r(t_{\rm max})$ and distinguish between the reversible and
irreversible phases by comparing the value of $\delta r_{\infty}$ to a
threshold $\delta r_{\rm th}$. In this study, we set $t_{\rm max}=4000$ and
$\delta r_{\rm th}=0.1$, giving the boundary between irreversible (yellow) and
reversible (blue) regions in Fig. 1a.

The above results imply that the configurational space is effectively explored
by the dynamical trajectory of the system in the irreversible phase,
encouraging us to consider a statistical ensemble generated by shear dynamics.
In the main text below, we fix $\gamma_{\rm max}={0.7}$, and vary $\varphi$
systematically in the window of $[0.833,0.849]$. At each $\varphi$, in total
$t_{\rm max}\times N_{\rm s}$ independent configurations are collected to
construct the ensemble. In the Extended Data, we present additional results
obtained for $\gamma_{\rm max}={1}$. All configurations are typically
isotropic since they are collected at $\gamma=0$.

Figure 1d shows the time evolution of the coordination number $Z(t)$ in the
irreversible phase, obtained from a typical simulation near $\varphi_{\rm J}$.
Rattlers (particles with fewer than $d+1$ contacts) are removed in the
computation of $Z$. At first glance, one sees the coexistence of jammed
($Z\approx 2d=4$) and unjammed states ($Z=0$), similar to the coexistence of
two ferromagnetic states (positive $m$ and negative $m$) in the time evolution
of the magnetization $m(t)$ in an Ising model near a first-order phase
transition Binder (1987a). However, a more careful examination reveals four
distinct states, which we discuss below.

Figure 1: Cyclic athermal quasi-static shear simulations. (a) Reversible
($\delta r_{\infty}<0.1$, blue region) and irreversible ($\delta
r_{\infty}>0.1$, yellow region) phases ($N={1000}$). The solid line represents
$\varphi_{\rm J}^{\infty}={0.8432}$, and the dashed line represents
$\gamma_{\rm max}=0.7$ used to generate the ensemble in the main text. (b)
One-cycle displacement $\delta r(t)$ and (d) average coordination number
$Z(t)$ for a typical sample. (c) The MSD data show diffusive behavior
$\langle\Delta r^{2}(t)\rangle\sim t$ (line). Data in (b-c) are obtained at a
fixed $\varphi=0.841$.

Four states: unjammed, jammed, partially crystallized and fragile.

Figure 2: Four states. (a) Semi-log plot of the probability distribution
$p(Z)$. The four peaks of unjammed (U), partially crystallized (C), fragile
(F) and jammed (J) states, are indicated. The inset shows a linear plot of
$p(Z)$ without the unjammed peak. Fractions are plotted as a function of
$\varphi$ in (b), and of $t$ in (c). In (c), each data point is obtained for a
small window $[t-\delta t,t+\delta t]$ with $\delta t={250}$. (d) Peak
coordination numbers $Z_{\rm J}^{*}$ of jammed states (upper branch, open
symbols) and $Z_{\rm F}^{*}$ of fragile states (lower branch, filled symbols)
as functions of $\varphi$, for $N=256$ (diamonds), 512 (triangles), 1000
(squares), and 2000 (circles). In (e-h), we show typical configurations of (e)
partially crystallized, (f) fragile percolating in one direction, (g) fragile
percolating in both directions, and (h) jammed states. Contact forces are
represented by bonds, whose width is proportional to the magnitude of force.
The red and blue disks are non-rattler and rattler particles respectively. We
set $\varphi=0.841$ for (a) and (c), $N=256$ for (b) and (c), and $N=1000$ for
(e-h).

In Fig. 2a, we plot the probability distribution $p(Z)$ of the states in the
considered ensemble, at a fixed $\varphi=0.841$ in the irreversible phase. The
following four states, represented by peaks in $p(Z)$, can be identified (see
Figs. 2e-h).

(i) Unjammed states. The left-most peak is a delta-function $p_{\rm
U}(Z)\sim\delta(Z)$, corresponding to unjammed states. All unjammed states
have strictly zero contacts, $Z=0$, once rattlers are removed.

(ii) Jammed states. The right-most peak $p_{\rm J}(Z)$ at $Z\geq 4$
corresponds to jammed states. Their average coordination numbers satisfy
minimally the isostatic condition, $Z\geq 2d=4$.

(iii) Partially crystallized states. The delta-peak at $p_{\rm
C}(Z)\sim\delta(Z-24/7)$ represents the states with a single unit cell of the
hexagonal crystal (see Fig. 2e). Here $24/7\approx 3.4286$ is the average
number of contacts of the seven particles forming the unit cell. Occasionally,
states with two or three unit crystal cells can be also found but they are
rare.

(iv) Fragile states. We define the states in the broad peak $3<Z<4$, excluding
the crystalline peak $p_{\rm C}(Z)$ as the fragile states.

The fractions, $F_{\rm U}$, $F_{\rm J}$, $F_{\rm C}$ and $F_{\rm F}$, of the
above four states (i-iv), are obtained by integrating corresponding peaks in
$p(Z)$. In Fig. 2b, the fractions are plotted as functions of $\varphi$,
showing that $F_{\rm J}(\varphi)$ increases from zero to one across
$\varphi_{\rm J}$. This behavior is quantitatively similar to $F_{\rm
J}(\varphi)$ obtained in previous rapid quench simulations, where finite-size
analyses have been carried out to precisely determine the asymptotic jamming
transition density $\varphi_{\rm J}^{\infty}$ in the thermodynamic limit
$N\rightarrow\infty$ O’hern  _et al._ (2003); Vagberg  _et al._ (2011). We
will perform such finite-size analyses later, after discussing the nature of
fragile states.

In Fig. 2c, we report various fractions as functions of $t$ at
$\varphi={0.841}$. The fractions are independent of $t$, confirming that the
system is stationary. We further divide $F_{\rm C}(t)$ into two parts, $F_{\rm
C}(t)=F_{\rm C}^{1}(t)+F_{\rm C}^{>}(t)$, where $F_{\rm C}^{1}(t)$ and $F_{\rm
C}^{>}(t)$ are respectively the fractions of partially crystallized states
with one and more than one crystal unit cells. Both $F_{\rm C}^{1}(t)$ and
$F_{\rm C}^{>}(t)$ are independent of $t$, indicating that the growth of seed
crystals is not observed. Because $F_{\rm C}$ is generally several orders of
magnitude lower than the fractions of other types, partially crystallized
states will be ignored in the following analyses.

In Fig. 2d, we plot the maximum coordination numbers, $Z_{\rm J}^{*}$ and
$Z_{\rm U}^{*}$, of jammed and fragile peaks. Around $\varphi_{\rm J}$, there
is a small gap $\delta Z_{\rm gap}\approx 0.2$ between $Z_{\rm J}^{*}$ (upper
branch) and $Z_{\rm U}^{*}$ (lower branch). The results are very similar to
those obtained by quasistatic uniform shear in Ref. Heussinger and Barrat
(2009) (see Fig. 4 therein). In Ref. Heussinger and Barrat (2009), it is
suggested that $\delta Z_{\rm gap}$ vanishes in the thermodynamic limit, and
thus the jamming transition is continuous. However, as we show below, the
fragile states are generated due to incomplete energy minimization and are
mechanically unstable. Once such fragile states are excluded, the lower branch
only contains unjammed states with $Z_{\rm U}^{*}=0$, which is separated by a
large gap $\delta Z_{\rm gap}\approx 4$ from the upper branch $Z_{\rm
J}^{*}\geq 4$.

Instability of fragile states.

Previously, the fragile states were obtained by uniform shear at non-zero
strains $\gamma>0$ below $\varphi_{\rm J}$ in experiments Bi  _et al._ (2011)
and simulations Kumar and Luding (2016). It was proposed that fragile and
jammed states differ in the percolation of the strong force network: in
fragile states the percolation occurs anisotropically in one direction only,
while in jammed states it occurs isotropically in both directions. In this
study, all states are collected during cyclic shear at $\gamma=0$, without
systematically introducing anisotropy. Indeed, the fragile states with $3<Z<4$
can have percolated force networks of non-rattlers in one or two directions
(see Fig. 2f and g). Thus in our case, anisotropy cannot effectively
distinguish fragile from jammed states.

We demonstrate that the essential difference between fragile and jammed states
is their mechanical stability. In fragile states, the potential energy per
particle is negligibly low (see Stopping criterion for energy minimization in
Methods), $\overline{e}_{i}<e_{\rm th}=10^{-20}$ (see Fig. 3a), but the
fraction of non-rattler particles is non-zero. These non-rattlers experience
forces, which can form a transient, percolated network. Such networks are
highly heterogeneous (see Fig. 2f and g), compared to those in jammed states
(see Fig. 2h). More importantly, the force networks in fragile states are
unstable, revealed by the non-zero net force per particle,
$\overline{f}_{i}>f_{\rm th}=10^{-14}$ (see Fig. 3b). Thus fragile states are
not strictly equilibrated; they turn into unjammed states by sufficiently long
relaxation (accurate energy minimization) or mechanical perturbations.

To demonstrate the instability of fragile states, two tests are carried out.
First, we perform a compression-decompression perturbation,
$\varphi\rightarrow\varphi+\delta\varphi\rightarrow\varphi$, where
$\delta\varphi=10^{-8}$, and the energy is minimized after each step. All
fragile states ($3<Z<4$) become unjammed ($Z=0$) after this perturbation,
while jammed states remain. In Fig. 3c, we plot the distribution
$p^{\prime}(Z)$ after the compression-decompression perturbation.
Independently, without any perturbation we simply regard all fragile states
($3<Z<4$) as unjammed ($Z=0$) and recalculate the distribution $\tilde{p}(Z)$.
Figure 3c shows that the two distributions $p^{\prime}(Z)$ and $\tilde{p}(Z)$
perfectly coincide.

In the second test, we repeat CAQS simulations by systematically varying the
threshold $f_{\rm th}$ in the criterion of the energy minimization algorithm.
The fraction $F_{\rm F}$ of fragile states grows and $F_{\rm U}$ of unjammed
states decays with increasing $f_{\rm th}$ (Fig. 3d), suggesting that many
unjammed states become fragile under a looser force-balance condition
$\overline{f}_{i}<f_{\rm th}$ with a larger $f_{\rm th}$. In contrast, the
constant $F_{\rm J}$ shows that the definition of jammed states is insensitive
to the algorithm parameter. Note that in granular experiments Bi  _et al._
(2011); Zhao  _et al._ (2019), the friction may play the role of $f_{\rm th}$,
i.e., the net inter-particle force could be balanced by the frictional force
between particles and the supporting plate. According to this assumption, the
probability of observing fragile states in experiments would depend on the
particle-plate friction that can be changed by replacing the materials. It
provides a protocol to examine our scenario of fragile states in future
experimental studies.

Because fragile states are unstable, from now on we count them as unjammed
states. More specifically, we replace the original distribution $p(Z)$ with
the modified distribution $\tilde{p}(Z)$. For simplicity, we omit the tilde
and denote $\tilde{p}(Z)$ by $p(Z)$ below.

Figure 3: Stability tests. Average single-particle (a) potential energy
$\overline{e}_{i}$ and (b) net force $\overline{f}_{i}$ as functions of $Z$.
Lines in (a) and (b) are guides to the eye. (c) Comparison between
$p^{\prime}(Z)$ obtained after the compression-decompression perturbation and
$\tilde{p}(Z)$ with fragile states counted as unjammed states, for $N=1000$ at
$\varphi=0.844$. (d) Fractions of four states as functions of $f_{\rm th}$.
Data are obtained for $N=256$ systems, except for (c). Figure 4: Distribution
of the coordination number. (a) Probability distribution $p(Z)$ at a few
different $\varphi$ for $N=1000$. (b) $p(Z)$ at $\varphi=0.843$ for a few
different $N$. For better visualization, we do not show the unjammed delta-
peak at $Z=0$. (c) The fraction of jammed states $F_{\rm J}$ as a function of
$\varphi$ for a few different $N$. The intersection of curves gives
$\varphi_{\rm J}^{\infty}=0.8432$. (d) The data points of $F_{\rm J}$ with
different $N$ collapse as a function of $(\varphi-\varphi_{\rm
J}^{\infty})N^{1/2}$. The solid line represents fitting to Eq. (4), with two
fitting parameters $u=0.041$ and $\sigma_{\varphi}=0.043$. The average
coordination number $\langle Z\rangle$ is plotted as a function of $\varphi$
in (e) and of $(\varphi-\varphi_{\rm J}^{\infty})N^{1/2}$ in (f). The solid
line in (f) represents $\langle Z\rangle=F_{\rm J}Z_{\rm J}^{*}$ using
$u=0.041$, $\sigma_{\varphi}=0.043$ and $Z_{\rm J}^{*}=4.1$. Symbols in (b-f)
have the same meanings.

Scaling analysis of $p(Z)$ near the jamming transition.

The new distribution $p(Z)$ contains two peaks of unjammed and jammed states,
whose fractions satisfy $F_{\rm U}(\varphi)+F_{\rm J}(\varphi)=1$ at any
$\varphi$. In Figs 4a and b, we plot $p(Z)$ for a few different $\varphi$ and
$N$, showing that $p(Z)$ is always double-peaked across the jamming
transition, which is typical behavior of a first-order rather than a second-
order transition.

To analyze the scaling behavior of $p(Z)$, we consider a general first-order
form of first-order phase transitions Binder (1987a, b),

$p(Z)=(1-F_{\rm J})\delta(Z)+F_{\rm J}\,p_{\rm J}\left[(Z-Z_{\rm
J}^{*})N^{\eta}\right],$ (3)

where $\delta(Z)$ and $p_{\rm J}(Z)$ correspond to the unjammed and jammed
peaks respectively. The asymptotic jamming density $\varphi_{\rm
J}^{\rm\infty}=0.8432$ for $N\rightarrow\infty$ is determined by the
intersection of $F_{\rm J}(\varphi)$ curves with different $N$ (see Fig. 4c),
close to the J-point density $\varphi_{\rm J}\approx 0.842-0.843$ reported in
previous studies O’hern  _et al._ (2003); Vagberg  _et al._ (2011). The
difference between $\varphi_{\rm J}^{\rm\infty}=0.8432$ by $\gamma_{\rm
max}=0.7$ and $\varphi_{\rm J}^{\rm\infty}=0.8435$ by $\gamma_{\rm max}=1$
(see Extended Data Fig. E1) is too small to conclude any systematic dependence
on $\gamma_{\rm max}$. In Ref. Das  _et al._ (2020), an unjamming pocket is
reported in the phase diagram for $\gamma_{\rm max}<0.17$, while for larger
$\gamma_{\rm max}$, $\varphi_{\rm J}$ is independent of $\gamma_{\rm max}$,
consistent with our observations. These results also suggest that the jamming
density obtained by cyclic shear is the lowest jamming density on the J-line
of all possible jammed states Parisi and Zamponi (2010); Chaudhuri  _et al._
(2010); Babu  _et al._ (2021).

As shown in Fig. 2d, $Z_{\rm J}^{*}$ is independent of $N$ near $\varphi_{\rm
J}$, and weakly depends on $\varphi$. In the following scaling analysis,
$Z_{\rm J}^{*}$ is approximated by a constant value $Z_{\rm J}^{*}\approx 4.1$
at $\varphi_{\rm J}^{\rm\infty}=0.8432$. To keep the expressions concise, we
ignore the corrections to the scaling functions from the $\varphi$-dependence
of $Z_{\rm J}^{*}$.

We assume the fraction of jammed states $F_{\rm J}$ having the following
scaling form, $F_{\rm J}(\varphi,N)=\mathcal{F}\left[(\varphi-\varphi_{\rm
J}^{\infty})N^{\lambda}\right]$. The value of the exponent $\lambda$ is
important for understanding the nature of the transition. If the system were
thermal, the fraction would be determined by the Boltzmann distribution,
$F_{\rm J}\sim\exp\left(\frac{N\delta f}{k_{\rm B}T}\right)$, where $\delta f$
is the single-particle free energy difference between two phases Binder
(1987b, a). Because the free energy is non-singular around a first-order phase
transition, it can be expanded, giving $\delta f\sim(T-T_{\rm c})$ to the
lowest order. Thus, if the jamming transition were a standard first-order
phase transition, then $F_{\rm J}(\varphi,N)=\mathcal{F}_{\rm
J}\left[(\varphi-\varphi_{\rm J}^{\infty})N\right]$, i.e., $\lambda=1$ (note
that $\varphi$ in the granular system is analogous to the temperature $T$ in
thermal systems). However, our numerical data can not be collapsed using
$\lambda=1$; in contrast, they can be perfectly collapsed using $\lambda=1/2$
(see Fig. 4d).

To understand the origin of $\lambda=1/2$ scaling, recall that the athermal
jamming transition is not driven by the free energy difference between the two
phases. Below we explain the scaling by the scenario of a first-order
transition with quenched disorder. For a finite $N$, due to the presence of
disorder in the packing structure, the jamming density $\varphi_{\rm J}^{N}$
should fluctuate around the asymptotic value $\varphi_{\rm J}^{\infty}$. Let
us assume a simple Gaussian form of the distribution, $\rho(\varphi_{\rm
J}^{N})\sim\exp\left[-\frac{(\delta\hat{\varphi}_{\rm
J}+u)^{2}}{2\sigma_{\varphi}^{2}}\right],$ where $\delta\hat{\varphi}_{\rm
J}=(\varphi_{\rm J}^{N}-\varphi_{\rm J}^{\infty})N^{1/2}$ follows the standard
central limit theorem. Note that $\rho(\varphi_{\rm J}^{N})$ has been
explicitly measured in the compression protocol O’Hern  _et al._ (2002);
O’hern  _et al._ (2003): the width of $\rho(\varphi_{\rm J}^{N})$ scales as
$w\sim N^{-0.55}$ in both 2D and 3D, supporting our assumption. The above
assumption also predicts, $\varphi_{\rm J}^{\infty}-\varphi_{\rm J}^{N}\sim
N^{-1/2}$, independent of dimensionality. The jammed states are defined by
$\varphi>\varphi_{\rm J}^{N}$, and thus

$\begin{split}F_{\rm J}(\varphi,N)=\int_{0}^{\varphi}\rho(\varphi_{\rm
J}^{N})\,d\varphi_{\rm
J}^{N}\approx\frac{1}{2}+\frac{1}{2}\text{erf}\left[\frac{\delta\hat{\varphi}+u}{\sqrt{2}\sigma_{\varphi}}\right],\end{split}$
(4)

where $\delta\hat{\varphi}=(\varphi-\varphi_{\rm J}^{\infty})N^{1/2}$.
Equation (4) agrees well with the data in Fig. 4d.

Figure 5: Susceptibilities. We plot (a) the disconnected susceptibility
$\chi_{\rm dis}$ rescaled by $N$, (b) the connected susceptibility $\chi_{\rm
con}$ rescaled by $N^{1/2}$, and (c) the ratio between $\chi_{\rm dis}$ and
$\chi_{\rm con}^{2}$, as functions of $(\varphi-\varphi_{\rm
J}^{\infty})N^{1/2}$. The solid black lines in (a-c) are Eqs. (5), (6) and (7)
respectively. The dashed red line in (c) is the constant term in Eq. (7),
i.e., $\chi_{\rm dis}/\chi_{\rm con}^{2}=\pi\sigma_{\varphi}^{2}/2$. To draw
these lines, we use $u=0.036$ and $\sigma_{\varphi}=0.043$ obtained from the
fitting in Fig. 4d, and $Z_{\rm J}^{*}=4.1$ determined at $\varphi_{\rm
J}^{\infty}=0.8432$ in Fig. 2d. No fitting is performed here. Symbols in (a-c)
have the same meanings as indicated in the legend of (a). In (d) and (e), the
susceptibilities of jammed states $\chi_{\rm J}$ are plotted as functions of
$N$ for a few different $\varphi$, and are fitted to Eq. (8). The open squares
in (d) are $\chi_{\rm J}$ data obtained in a small pressure window
$0.0009<P<0.0011$, corresponding to $\varphi\approx 0.843$. The fitting
parameter $\mu$ is plotted in (f) as a function of $\varphi$, where the error
bar represents the fitting error.

Once the scaling behavior of $F_{\rm J}(\varphi,N)$ is obtained, it is easy to
derive scaling forms of susceptibilities. In the random field Ising model
(RFIM) Gofman  _et al._ (1993); Nattermann (1998), two susceptibilities have
been introduced: a connected susceptibility $\chi_{\rm con}$ associated to
thermal fluctuations, and a disconnected susceptibility $\chi_{\rm dis}$
associated to fluctuations caused by disorder. The two susceptibilities are
related by $\chi_{\rm dis}=\chi_{\rm con}^{2}$ around the critical point and
first-order transitions in RFIM. We similarly define two susceptibilities in
the jamming transition. The disconnected susceptibility is defined as
$\chi_{\rm dis}\equiv N\sigma_{Z}^{2}=N\langle(Z-\langle
Z\rangle)^{2}\rangle$, where $\sigma_{Z}^{2}$ is the variance of $p(Z)$ and
$\langle\ldots\rangle$ is the average over all states. From Eq. (3), we
obtain,

$\frac{\chi_{\rm dis}(\varphi,N)}{N}\approx(Z_{\rm J}^{*})^{2}[1-F_{\rm
J}(\varphi,N)]F_{\rm J}(\varphi,N),$ (5)

where $F_{\rm J}(\varphi,N)$ is given by Eq. (4). Equation (5) suggests a
scaling form $\chi_{\rm dis}(\varphi,N)/N=\mathcal{X}_{\rm
dis}\left[(\varphi-\varphi_{\rm J}^{\infty})N^{1/2}\right]$, in agreement with
the data in Fig. 5a.

The alternative connected susceptibility is defined as $\chi_{\rm
con}=d\langle Z\rangle/d\varphi$. From Eq. (3), one finds that $\langle
Z\rangle\approx F_{\rm J}Z_{\rm J}^{*}$ and thus $\langle Z\rangle$ has the
same scaling form as $F_{\rm J}$, i.e., $\langle
Z\rangle(\varphi,N)=\mathcal{Z}[(\varphi-\varphi_{\rm J}^{\infty})N^{1/2}]$,
consistent with the data in Fig. 4e and f. Using $\chi_{\rm con}=Z_{\rm
J}^{*}\frac{dF_{\rm J}}{d\varphi}$ and Eq. (4), we derive,

$\frac{\chi_{\rm con}(\varphi,N)}{N^{1/2}}\approx\frac{Z_{\rm
J}^{*}}{\sigma_{\varphi}\sqrt{2\pi}}\exp\left[-\left(\frac{\delta\hat{\varphi}+u}{\sqrt{2}\sigma_{\varphi}}\right)^{2}\right],$
(6)

which is verified by the simulation data in Fig. 5b.

Now we can look at the relationship between $\chi_{\rm dis}$ and $\chi_{\rm
con}$. Expanding Eqs. (5) and (6) around the maxima, where
$x\equiv\frac{\delta\hat{\varphi}+u}{\sqrt{2}\sigma_{\varphi}}=0$, we obtain,
up to the quadratic order,

$\frac{\chi_{\rm dis}}{\chi_{\rm
con}^{2}}\approx\frac{\pi\sigma^{2}_{\varphi}}{2}\left[1+\left(2-\frac{4}{\pi}\right)\left(\frac{\delta\hat{\varphi}+u}{\sqrt{2}\sigma_{\varphi}}\right)^{2}\right].$
(7)

To the lowest order, Eq. (7) gives a scaling, $\chi_{\rm dis}\sim\chi_{\rm
con}^{2},$ which is a key signature of the presence of quenched disorder. This
square relationship appears generally at first-order transitions in disordered
systems, such as the first-order phase transition in the RFIM Nattermann
(1998), the brittle yielding of amorphous solids Ozawa  _et al._ (2018) and
the melting of ultrastable hard-sphere glasses Zhang  _et al._ (2022).
Comparing to the simulation data, Eq. (7) works well around the extreme point
(see Fig. 5c).

It seems that the non-Gaussian effect in the distribution $\rho(\varphi_{\rm
J}^{N})$, which is neglected in the present analysis, is amplified in the data
of $\chi_{\rm dis}/\chi_{\rm con}^{2}$, resulting in slight asymmetry.
Extended Data Figs. E1-E4 show that the scaling functions Eqs. (4)-(7) work in
3D, and for a different max strain $\gamma_{\max}=1$ in 2D.

Finally let us discuss the finite-size effects of the jammed peak $p_{\rm
J}(Z)$ in Eq. (3). For standard first-order phase transitions, $\eta=1/2$
Binder (1987b, a). Thus in that scenario the susceptibility of the jamming
peak would scale as, $\chi_{\rm J}\equiv N\sigma_{\rm J}^{2}=N\langle(Z-Z_{\rm
J}^{*})^{2}\rangle_{\rm J}\sim N^{0}$, where $\sigma_{\rm J}^{2}$ is the
variance of $p_{\rm J}(Z)$ and $\langle\ldots\rangle_{\rm J}$ represents the
average restricted to the jammed states only. However, our $\chi_{\rm J}$ data
disagree with this scaling (see Figs. 5d-f). At different $\varphi$,
$\chi_{\rm J}$ can be fitted to a pow-law form,

$\chi_{\rm J}\sim N^{-\mu/d},$ (8)

where $\mu=d(2\eta-1)$ is a fitting parameter. At large or small $\varphi$
away from $\varphi_{\rm J}^{\rm\infty}$ (see Fig. 5f), the exponent is close
to zero ($\mu\approx-0.2$), suggesting that the local coordination numbers
$Z_{i}$ are uniformly distributed. The deviation from the uncorrelated
behavior ($\mu=0$) is significant around $\varphi_{\rm J}^{\rm\infty}$, where
$\mu\approx-0.8$ reaches the minimum.

The result of $\mu$ shall be interpreted with care. In general, (i) $\mu=0$
corresponds to uniformity of $Z_{i}$, (ii) $\mu>0$ to hyperuniformity with a
vanishing $\chi_{\rm J}$ in the thermodynamic limit, and (iii) $\mu<0$ to
hyperfluctuations with a diverging $\chi_{\rm J}$ in the thermodynamic limit
that typically appears at the critical point in a second-order phase
transition. However, the negative $\mu$ in Fig. 5f is not due to the
criticality of the jamming transition. Here it is essential to consider the
volume fluctuations. The negative $\mu$ is obtained from a
$\varphi$-controlled setup in our CAQS simulations. If we instead select
configurations around a constant pressure $P$, then the fluctuations become
significantly smaller, and $\mu\approx 0$ (Fig. 5d). This observation is
consistent with a previous study Ikeda (2023), suggesting that the large
fluctuation of $Z$ near $\varphi_{\rm J}^{\rm\infty}$ in the
$\varphi$-controlled protocol might be originated from the fluctuation of
$\varphi_{\rm J}^{N}$. In fact, careful measurement gives $\mu\simeq 1$ in a
$P$-controlled protocol at jamming, confirming hyperuniformity Hexner  _et
al._ (2018, 2019). We emphasize that the finite-size effects of the jammed
peak $p_{\rm J}(Z)$ contribute negligibly to the scaling of $F(\varphi,N)$,
$\chi_{\rm dis}(\varphi,N)$ and $\chi_{\rm con}(\varphi,N)$. In other words,
the first-order nature of the jamming transition is independent of how $Z_{i}$
is spatially distributed in jammed packings.

Conclusion.

In this study, we understand the nature of jamming transition through an
ensemble approach analogous to the statistical mechanics in equilibrium
systems Baule  _et al._ (2018); Bi  _et al._ (2015). Within such a framework,
the order of the jamming transition is independent of the complex properties
of jammed packings, including isostaticity, marginality and hyperuniformity
(see (i-iii) in Introduction). As a first-order jamming transition, the
finite-size scaling exponents measured in this study are naturally universal
in 2D and 3D, without the need to introduce the upper critical dimension
$d_{\rm u}=2$ Parisi  _et al._ (2018), a concept relevant to continuous
transitions only. It is of particular interest to reconcile the first-order
transition established here under quasi-static shearing and the second-order
transition observed in previous finite-rate rheology (see (iv) in
Introduction). The jamming transition appears to be a rare example with the
simultaneous existence of first-order and second-order transitions between
liquid and solid states. Conventionally, the two kinds of transitions only
coexist in gas-liquid systems, but not in liquid-crystalline solid systems.

Methods.

Models. We study models of frictionless, bidisperse particles in two and three
dimensions. The number ratio between large and small particles is 1:1, and the
diameter ratio is 1.4:1. The potential energy between two particles is:

$U(r_{ij})=\frac{\epsilon}{2}\left(1-\frac{r_{ij}}{\sigma_{ij}}\right)^{2}\Theta\left(1-\frac{r_{ij}}{\sigma_{ij}}\right),$
(9)

where $\epsilon=1$ is the energy unit, $r_{ij}$ the distance between particles
$i$ and $j$, $\sigma_{ij}=\frac{\sigma_{i}+\sigma_{j}}{2}$ the mean diameter,
and $\Theta(x)$ the Heaviside step function. We set unit particle mass, and
the diameter of small particles as unit length.

Shear protocol. Particles are randomly distributed at an initial volume
fraction $\varphi_{0}=0.02$. The system is then rapidly quench compressed to
the target $\varphi$. The CAQS is performed under the Lees-Edwards boundary
conditions, with a fixed $\varphi$. We use a strain step $\delta\gamma=0.1$ to
generate the phase diagram in Fig. 1a, and $\delta\gamma=0.01$ for other
results. At each step, particle positions are shifted according to
$x_{i}\rightarrow x_{i}+\delta\gamma y_{i}$, and then the system’s energy is
minimized using the FIRE algorithm Bitzek  _et al._ (2006). During a cycle,
the strain is varied as $\gamma=0\rightarrow\gamma_{\rm max}\rightarrow 0$. We
set $\gamma_{\rm max}=0.7$ in the main text and $\gamma_{\rm max}=1$ in the
Extended Data for the 2D model, and $\gamma_{\rm max}=0.5$ for the 3D model
(Extended Data). The number of cycles is represented by $t$ with unit
oscillation period ($t=1$). In 2D, the maximum number of cycles is $t_{\rm
max}=250$ for $0.839\leq\varphi\leq{0.849}$ and $t_{\rm max}=4000$ for
${0.833}\leq\varphi<0.839$, while in 3D, $t_{\rm max}=144$ for $2000$, and
$t_{\rm max}=50$ for $N=512$ and $1000$. We generate $N_{\rm s}=4000-12000$
independent samples at each $\varphi$. At each $\varphi$, in total $t_{\rm
max}\times N_{\rm s}$ configurations are collected to compute statistical
quantities.

Stopping criterion for energy minimization. We terminate the energy
minimization when the potential energy per particle
$\overline{e}_{i}=\frac{1}{N}\sum_{i=1}^{N}e_{i}$ falls below a threshold
$e_{\rm th}$ or the average single-particle net force
$\overline{f}_{i}=\frac{1}{N}\sum_{i=1}^{N}f_{i}$ falls below $f_{\rm th}$,
whichever is satisfied earlier. We set $e_{\rm th}=10^{-20}$ and $f_{\rm
th}=10^{-14}$ unless otherwise specified.

###### Acknowledgements.

We warmly thank Hajime Yoshino for inspiring discussions. We acknowledge
financial support from NSFC (Grants 12161141007, 11935002, 11974361 and
12047503), from Chinese Academy of Sciences (Grants ZDBS-LY-7017 and
KGFZD-145-22-13), and from Wenzhou Institute (Grant WIUCASQD2023009). In this
work access was granted to the High-Performance Computing Cluster of Institute
of Theoretical Physics - the Chinese Academy of Sciences.

## References

  * Makse  _et al._ (2000) Hernán A Makse, David L Johnson, and Lawrence M Schwartz, “Packing of compressible granular materials,” Physical review letters 84, 4160 (2000).
  * O’hern  _et al._ (2003) Corey S O’hern, Leonardo E Silbert, Andrea J Liu, and Sidney R Nagel, “Jamming at zero temperature and zero applied stress: The epitome of disorder,” Physical Review E 68, 011306 (2003).
  * Liu and Nagel (2010) Andrea J Liu and Sidney R Nagel, “The jamming transition and the marginally jammed solid,” Annu. Rev. Condens. Matter Phys. 1, 347–369 (2010).
  * van Hecke (2009) Martin van Hecke, “Jamming of soft particles: geometry, mechanics, scaling and isostaticity,” Journal of Physics: Condensed Matter 22, 033101 (2009).
  * Pan  _et al._ (2023a) Deng Pan, Yinqiao Wang, Hajime Yoshino, Jie Zhang, and Yuliang Jin, “A review on shear jamming,” Physics Reports 1038, 1–18 (2023a), a review on shear jamming.
  * Parisi and Zamponi (2010) Giorgio Parisi and Francesco Zamponi, “Mean-field theory of hard sphere glasses and jamming,” Reviews of Modern Physics 82, 789 (2010).
  * Charbonneau  _et al._ (2017) Patrick Charbonneau, Jorge Kurchan, Giorgio Parisi, Pierfrancesco Urbani, and Francesco Zamponi, “Glass and jamming transitions: From exact results to finite-dimensional descriptions,” Annual Review of Condensed Matter Physics 8, 265–288 (2017).
  * Parisi  _et al._ (2020) Giorgio Parisi, Pierfrancesco Urbani, and Francesco Zamponi,  _Theory of simple glasses: exact solutions in infinite dimensions_ (Cambridge University Press, 2020).
  * Maxwell (1870) J Clerk Maxwell, “I.—on reciprocal figures, frames, and diagrams of forces,” Earth and Environmental Science Transactions of the Royal Society of Edinburgh 26, 1–40 (1870).
  * Alexander (1998) Shlomo Alexander, “Amorphous solids: their structure, lattice dynamics and elasticity,” Physics reports 296, 65–236 (1998).
  * Wyart (2005) M Wyart, “On the rigidity of amorphous solids,” in  _Annales de Physique_ , Vol. 30 (EDP Sciences, 2005) pp. 1–96.
  * wya (2005) “Geometric origin of excess low-frequency vibrational modes in weakly connected amorphous solids,” Europhysics Letters 72, 486 (2005).
  * Silbert  _et al._ (2005) Leonardo E Silbert, Andrea J Liu, and Sidney R Nagel, “Vibrations and diverging length scales near the unjamming transition,” Physical review letters 95, 098301 (2005).
  * Wyart (2010) Matthieu Wyart, “Scaling of phononic transport with connectivity in amorphous solids,” Europhysics Letters 89, 64001 (2010).
  * DeGiuli  _et al._ (2014a) Eric DeGiuli, Adrien Laversanne-Finot, Gustavo Düring, Edan Lerner, and Matthieu Wyart, “Effects of coordination and pressure on sound attenuation, boson peak and elasticity in amorphous solids,” Soft matter 10, 5628–5644 (2014a).
  * Lerner  _et al._ (2014) Edan Lerner, Eric DeGiuli, Gustavo Düring, and Matthieu Wyart, “Breakdown of continuum elasticity in amorphous solids,” Soft Matter 10, 5085–5092 (2014).
  * Wyart (2012) Matthieu Wyart, “Marginal stability constrains force and pair distributions at random close packing,” Physical review letters 109, 125502 (2012).
  * Lerner  _et al._ (2013) Edan Lerner, Gustavo Düring, and Matthieu Wyart, “Low-energy non-linear excitations in sphere packings,” Soft Matter 9, 8252–8263 (2013).
  * DeGiuli  _et al._ (2014b) Eric DeGiuli, Edan Lerner, Carolina Brito, and Matthieu Wyart, “Force distribution affects vibrational properties in hard-sphere glasses,” Proceedings of the National Academy of Sciences 111, 17054–17059 (2014b).
  * Müller and Wyart (2015) Markus Müller and Matthieu Wyart, “Marginal stability in structural, spin, and electron glasses,” Annu. Rev. Condens. Matter Phys. 6, 177–200 (2015).
  * Goodrich  _et al._ (2016) Carl P Goodrich, Andrea J Liu, and James P Sethna, “Scaling ansatz for the jamming transition,” Proceedings of the National Academy of Sciences 113, 9745–9750 (2016).
  * Baule  _et al._ (2018) Adrian Baule, Flaviano Morone, Hans J Herrmann, and Hernán A Makse, “Edwards statistical mechanics for jammed granular matter,” Reviews of modern physics 90, 015006 (2018).
  * Zaccone and Scossa-Romano (2011) Alessio Zaccone and Enzo Scossa-Romano, “Approximate analytical description of the nonaffine response of amorphous solids,” Physical Review B 83, 184205 (2011).
  * Pan  _et al._ (2023b) Deng Pan, Fanlong Meng, and Yuliang Jin, “Shear hardening in frictionless amorphous solids near the jamming transition,” PNAS nexus 2, pgad047 (2023b).
  * Hexner  _et al._ (2018) Daniel Hexner, Andrea J Liu, and Sidney R Nagel, “Two diverging length scales in the structure of jammed packings,” Physical review letters 121, 115501 (2018).
  * Hexner  _et al._ (2019) Daniel Hexner, Pierfrancesco Urbani, and Francesco Zamponi, “Can a large packing be assembled from smaller ones?” Physical review letters 123, 068003 (2019).
  * Charbonneau  _et al._ (2014) Patrick Charbonneau, Jorge Kurchan, Giorgio Parisi, Pierfrancesco Urbani, and Francesco Zamponi, “Fractal free energy landscapes in structural glasses,” Nature communications 5, 3725 (2014).
  * Charbonneau  _et al._ (2015a) Patrick Charbonneau, Yuliang Jin, Giorgio Parisi, Corrado Rainone, Beatriz Seoane, and Francesco Zamponi, “Numerical detection of the gardner transition in a mean-field glass former,” Physical Review E 92, 012316 (2015a).
  * Berthier  _et al._ (2016) Ludovic Berthier, Patrick Charbonneau, Yuliang Jin, Giorgio Parisi, Beatriz Seoane, and Francesco Zamponi, “Growing timescales and lengthscales characterizing vibrations of amorphous solids,” Proceedings of the National Academy of Sciences 113, 8397–8401 (2016).
  * Berthier  _et al._ (2019) Ludovic Berthier, Giulio Biroli, Patrick Charbonneau, Eric I Corwin, Silvio Franz, and Francesco Zamponi, “Gardner physics in amorphous solids and beyond,” The Journal of chemical physics 151 (2019).
  * Li  _et al._ (2021) Huaping Li, Yuliang Jin, Ying Jiang, and Jeff ZY Chen, “Determining the nonequilibrium criticality of a gardner transition via a hybrid study of molecular simulations and machine learning,” Proceedings of the National Academy of Sciences 118, e2017392118 (2021).
  * Urbani  _et al._ (2023) Pierfrancesco Urbani, Yuliang Jin, and Hajime Yoshino, “The gardner glass,” in  _Spin Glass Theory and Far Beyond: Replica Symmetry Breaking After 40 Years_ (World Scientific, 2023) pp. 219–238.
  * Charbonneau  _et al._ (2015b) Patrick Charbonneau, Eric I Corwin, Giorgio Parisi, and Francesco Zamponi, “Jamming criticality revealed by removing localized buckling excitations,” Physical review letters 114, 125504 (2015b).
  * Wang  _et al._ (2022) Yinqiao Wang, Jin Shang, Yuliang Jin, and Jie Zhang, “Experimental observations of marginal criticality in granular materials,” Proceedings of the National Academy of Sciences 119, e2204879119 (2022).
  * Goodrich  _et al._ (2012) Carl P Goodrich, Andrea J Liu, and Sidney R Nagel, “Finite-size scaling at the jamming transition,” Physical review letters 109, 095704 (2012).
  * Olsson and Teitel (2007) Peter Olsson and Stephen Teitel, “Critical scaling of shear viscosity at the jamming transition,” Physical review letters 99, 178001 (2007).
  * Olsson and Teitel (2011) Peter Olsson and S Teitel, “Critical scaling of shearing rheology at the jamming transition of soft-core frictionless disks,” Physical Review E 83, 030302 (2011).
  * Kawasaki  _et al._ (2015) Takeshi Kawasaki, Daniele Coslovich, Atsushi Ikeda, and Ludovic Berthier, “Diverging viscosity and soft granular rheology in non-brownian suspensions,” Physical Review E 91, 012203 (2015).
  * Heussinger and Barrat (2009) Claus Heussinger and Jean-Louis Barrat, “Jamming transition as probed by quasistatic shear flow,” Physical review letters 102, 218303 (2009).
  * Vagberg  _et al._ (2011) Daniel Vagberg, Daniel Valdez-Balderas, MA Moore, Peter Olsson, and Stephen Teitel, “Finite-size scaling at the jamming transition: Corrections to scaling and the correlation-length critical exponent,” Physical Review E 83, 030303 (2011).
  * Bertrand  _et al._ (2016) Thibault Bertrand, Robert P Behringer, Bulbul Chakraborty, Corey S O’Hern, and Mark D Shattuck, “Protocol dependence of the jamming transition,” Physical Review E 93, 012901 (2016).
  * Baity-Jesi  _et al._ (2017) Marco Baity-Jesi, Carl P Goodrich, Andrea J Liu, Sidney R Nagel, and James P Sethna, “Emergent so (3) symmetry of the frictionless shear jamming transition,” Journal of Statistical Physics 167, 735–748 (2017).
  * Jin and Yoshino (2021) Yuliang Jin and Hajime Yoshino, “A jamming plane of sphere packings,” Proceedings of the National Academy of Sciences 118, e2021794118 (2021).
  * Chaudhuri  _et al._ (2010) Pinaki Chaudhuri, Ludovic Berthier, and Srikanth Sastry, “Jamming transitions in amorphous packings of frictionless spheres occur over a continuous range of volume fractions,” Physical review letters 104, 165701 (2010).
  * Ozawa  _et al._ (2012) Misaki Ozawa, Takeshi Kuroiwa, Atsushi Ikeda, and Kunimasa Miyazaki, “Jamming transition and inherent structures of hard spheres and disks,” Physical review letters 109, 205701 (2012).
  * Ashwin  _et al._ (2012) SS Ashwin, Jerzy Blawzdziewicz, Corey S O’Hern, and Mark D Shattuck, “Calculations of the structure of basin volumes for mechanically stable packings,” Physical Review E 85, 061307 (2012).
  * Martiniani  _et al._ (2017) Stefano Martiniani, K Julian Schrenk, Kabir Ramola, Bulbul Chakraborty, and Daan Frenkel, “Numerical test of the edwards conjecture shows that all packings are equally probable at jamming,” Nature physics 13, 848–851 (2017).
  * Pine  _et al._ (2005) David J Pine, Jerry P Gollub, John F Brady, and Alexander M Leshansky, “Chaos and threshold for irreversibility in sheared suspensions,” Nature 438, 997–1000 (2005).
  * Corte  _et al._ (2008) Laurent Corte, Paul M Chaikin, Jerry P Gollub, and David J Pine, “Random organization in periodically driven systems,” Nature Physics 4, 420–424 (2008).
  * Schreck  _et al._ (2013) Carl F Schreck, Robert S Hoy, Mark D Shattuck, and Corey S O’Hern, “Particle-scale reversibility in athermal particulate media below jamming,” Physical Review E 88, 052205 (2013).
  * Nagasawa  _et al._ (2019) Kentaro Nagasawa, Kunimasa Miyazaki, and Takeshi Kawasaki, “Classification of the reversible–irreversible transitions in particle trajectories across the jamming transition point,” Soft matter 15, 7557–7566 (2019).
  * Das  _et al._ (2020) Pallabi Das, HA Vinutha, and Srikanth Sastry, “Unified phase diagram of reversible–irreversible, jamming, and yielding transitions in cyclically sheared soft-sphere packings,” Proceedings of the National Academy of Sciences 117, 10203–10209 (2020).
  * O’Hern  _et al._ (2002) Corey S O’Hern, Stephen A Langer, Andrea J Liu, and Sidney R Nagel, “Random packings of frictionless particles,” Physical Review Letters 88, 075507 (2002).
  * Binder (1987a) K Binder, “Finite size effects on phase transitions,” Ferroelectrics 73, 43–67 (1987a).
  * Bi  _et al._ (2011) Dapeng Bi, Jie Zhang, Bulbul Chakraborty, and Robert P Behringer, “Jamming by shear,” Nature 480, 355–358 (2011).
  * Kumar and Luding (2016) Nishant Kumar and Stefan Luding, “Memory of jamming–multiscale models for soft and granular matter,” Granular Matter 18, 58 (2016).
  * Zhao  _et al._ (2019) Yiqiu Zhao, Jonathan Barés, Hu Zheng, Joshua ES Socolar, Robert P Behringer, _et al._ , “Shear-jammed, fragile, and steady states in homogeneously strained granular materials,” Physical review letters 123, 158001 (2019).
  * Binder (1987b) Kurt Binder, “Theory of first-order phase transitions,” Reports on progress in physics 50, 783 (1987b).
  * Babu  _et al._ (2021) Varghese Babu, Deng Pan, Yuliang Jin, Bulbul Chakraborty, and Srikanth Sastry, “Dilatancy, shear jamming, and a generalized jamming phase diagram of frictionless sphere packings,” Soft Matter 17, 3121–3127 (2021).
  * Gofman  _et al._ (1993) Michael Gofman, Joan Adler, Amnon Aharony, A Brooks Harris, and Moshe Schwartz, “Evidence for two exponent scaling in the random field ising model,” Physical review letters 71, 1569 (1993).
  * Nattermann (1998) Thomas Nattermann, “Theory of the random field ising model,” in  _Spin glasses and random fields_ (World Scientific, 1998) pp. 277–298.
  * Ozawa  _et al._ (2018) Misaki Ozawa, Ludovic Berthier, Giulio Biroli, Alberto Rosso, and Gilles Tarjus, “Random critical point separates brittle and ductile yielding transitions in amorphous materials,” Proceedings of the National Academy of Sciences 115, 6656–6661 (2018).
  * Zhang  _et al._ (2022) Kaihua Zhang, Xinyang Li, Yuliang Jin, and Ying Jiang, “Machine learning glass caging order parameters with an artificial nested neural network,” Soft Matter 18, 6270–6277 (2022).
  * Ikeda (2023) Harukuni Ikeda, “Control parameter dependence of fluctuations near jamming,” The Journal of Chemical Physics 158 (2023).
  * Bi  _et al._ (2015) Dapeng Bi, Silke Henkes, Karen E Daniels, and Bulbul Chakraborty, “The statistical physics of athermal materials,” Annu. Rev. Condens. Matter Phys. 6, 63–83 (2015).
  * Parisi  _et al._ (2018) Giorgio Parisi, Yoav G Pollack, Itamar Procaccia, Corrado Rainone, and Murari Singh, “Robustness of mean field theory for hard sphere models,” Physical Review E 97, 063003 (2018).
  * Bitzek  _et al._ (2006) Erik Bitzek, Pekka Koskinen, Franz Gähler, Michael Moseler, and Peter Gumbsch, “Structural relaxation made simple,” Physical review letters 97, 170201 (2006).

Extended Data

Figure E1: Distribution of the coordination number in two dimensions for
$\gamma_{\rm max}=1$. (a) Fraction of jammed states $F_{\rm J}$ as a function
of $\varphi$ for a few different $N$. The intersection of curves gives
$\varphi_{\rm J}^{\infty}=0.8435$. (b) The data points of $F_{\rm J}$ with
different $N$ collapse as a function of $(\varphi-\varphi_{\rm
J}^{\infty})N^{1/2}$. The solid line represents the fitting according to Eq.
(4), with two fitting parameters $u=0.051$ and $\sigma_{\varphi}=0.043$.
Figure E2: Susceptibilities in two dimensions for $\gamma_{\rm max}=1$. (a)
The disconnected susceptibility $\chi_{\rm dis}$ rescaled by $N$, (b) the
connected susceptibility $\chi_{\rm con}$ rescaled by $N^{1/2}$, and (c) the
ratio between $\chi_{\rm dis}$ and $\chi_{\rm con}^{2}$ are plotted as
functions of $(\varphi-\varphi_{\rm J}^{\infty})N^{1/2}$. The solid black
lines in (a-c) are Eqs. (5), (6) and (7) respectively. The dashed red line in
(c) is the constant term in Eq. (7), i.e., $\chi_{\rm dis}/\chi_{\rm
con}^{2}=\pi\sigma_{\varphi}^{2}/2$. To draw these lines, we use $u=0.051$ and
$\sigma_{\varphi}=0.043$ obtained from the fitting in Fig. E1b. Symbols in
(a-c) have the same meanings as indicated in the legend of (a).  Figure E3:
Distribution of the coordination number in three dimensions. (a) Fraction of
jammed states $F_{\rm J}$ as a function of $\varphi$ for a few different $N$.
The intersection of curves gives $\varphi_{\rm J}^{\infty}=0.648$. (b) The
data points of $F_{\rm J}$ with different $N$ collapse as a function of
$(\varphi-\varphi_{\rm j}^{\infty})N^{1/2}$. The solid line represents the
fitting according to Eq. (4), with two fitting parameters $u=0.077$ and
$\sigma_{\varphi}=0.051$.  Figure E4: Susceptibilities in three dimensions.
(a) The disconnected susceptibility $\chi_{\rm dis}$ rescaled by $N$, (b) the
connected susceptibility $\chi_{\rm con}$ rescaled by $N^{1/2}$, and (c) the
ratio between $\chi_{\rm dis}$ and $\chi_{\rm con}^{2}$ are plotted as
functions of $(\varphi-\varphi_{\rm J}^{\infty})N^{1/2}$. The solid black
lines in (a-c) are Eqs. (5), (6) and (7) respectively. The dashed red line in
(c) is the constant term in Eq. (7), i.e., $\chi_{\rm dis}/\chi_{\rm
con}^{2}=\pi\sigma_{\varphi}^{2}/2$. To draw these lines, we use $u=0.077$ and
$\sigma_{\varphi}=0.051$ obtained from the fitting in Fig. E3b. Symbols in
(a-c) have the same meanings as indicated in the legend of (a).