Universality of jamming of non-spherical particles

Breathing particles model –

The BP model (43) was originally introduced to understand the physics of the Swap Monte Carlo algorithm (44), but here we will focus on its relation with the jamming of ellipsoids. The model consists of $N$ spherical particles with positions $\bm{x}_{i}$ in $d$-dimensions and radius $R_{i}\geq 0$, interacting via the potential energy:

where, defining $\theta(x)$ as the Heaviside theta function,

is the standard harmonic repulsive interaction potential of spherical particles such as bubbles and colloids (5), and the distribution of $R_{i}$, which can fluctuate around a reference value $R_{i}^{0}$, is controlled by the chemical potential term:

Here, $k_{R}$ is determined by imposing that the dimensionless standard deviation $\Delta\propto\sqrt{\sum_{i}(R_{i}-R_{i}^{0})^{2}/(NR_{0}^{2})}$ is constant, with $R_{0}=N^{-1}\sum_{i}R_{i}^{0}$. Note that $\Delta=0$ (corresponding to $k_{R}=\infty$) gives back the usual spherical particles (5), and that the full distribution of radii, $P(R)$, can generically change even if $\Delta$ is kept fixed. Upon approaching jamming, where the adimensional pressure $p$ (in units of $kR_{0}^{2-d}$) vanishes, it is found that $k_{R}=p/\Delta$ and $P(R)$ remains constant (43).

Because the BP model has $Nd$ translational degrees of freedom and $N$ radial degrees of freedom, the naive Maxwell stability condition requires $z\geq 2(d+1)$ in the thermodynamic limit (19, 45). However, a marginal stability argument and numerical simulations prove that the contact number at the jamming point $z_{J}$ increases continuously as $z_{J}-2d\propto\Delta^{1/2}$ (43) and the system is hypostatic for sufficiently small $\Delta$, i.e., the number of constraints is smaller than that required by the Maxwell’s stability condition. This is very similar to ellipsoids and motivates us to conjecture that the two models could belong to the same universality class. In the following, we show that this expectation is indeed true: hypostatic packings of the BP and ellipsoids are stabilized by a common mechanism and have the same critical exponents.

Mapping from ellipsoids to BP –

We now construct a mapping from a system of ellipsoids to the spherical BP model introduced above. Ellipsoids are described by their position $\bm{x}_{i}$ and by unit vectors $\hat{\bm{u}}_{i}$ along their principal axis, and for concreteness, we model them by the Gay-Berne potential (46, 31):

Here, $\hat{\bm{r}}_{ij}=(\bm{x}_{i}-\bm{x}_{j})/\left|\bm{x}_{i}-\bm{x}_{j}\right|$ is the unit vector connecting the $i$-th and $j$-th particles, $\varepsilon\sigma_{0}$ is the length of the principal axis, and $\chi=(\varepsilon^{2}-1)/(\varepsilon^{2}+1)$, where $\varepsilon$ denotes the aspect ratio. Because we are interested in the nearly spherical case, we expand the pair potential in small $\Delta=\varepsilon-1$ as

where $h_{ij}^{(0)}=r_{ij}/\sigma_{0}-1$ and $\Delta^{2}w_{ij}$ denotes the $O(\Delta^{2})$ term that we do not need to write explicitly. Substituting this in Eq. (4) and keeping terms up to $\Delta^{2}$, we obtain $V_{N}\approx U_{N}+\mu_{N}$, where

The stiffness matrix is $k_{i}^{ab}=-\Delta^{-1}\sum_{j(\neq i)}v^{\prime}(h_{ij}^{(0)})\hat{\bm{r}}_{ij}^{a}\hat{\bm{r}}_{ij}^{b}$, where $a,b=1,\cdots,d$. Note that near the jamming point, $k_{i}$ behaves as $k_{i}\sim v^{\prime}(h)/\Delta\sim p/\Delta$, which is the same scaling of the stiffness $k_{R}$ of the BP model, Eq. (3). Hence, if we identify $\Delta\hat{\bm{u}}_{i}$ with $R_{i}$, in the vicinity of jamming the potential for ellipsoids can be analyzed essentially in the same way as the BP model (43), as we discuss next.

Marginal stability –

The distinctive feature of both BP and ellipsoids is that the total potential, and thus the Hessian matrix, can be split in two parts: one having finite stiffness, and the second having vanishing stiffness $p/\Delta$ by dimensional arguments. The zero modes of the first term are stabilized by the second, as recognized in Refs. (29, 32). We now provide additional insight on this structure by generalizing a marginal stability argument discussed for the BP in Ref. (43). At jamming, $p=0$ and $V_{N}=U_{N}$ because $\mu_{N}\propto p$. The ${\cal N}_{3}\equiv Nz/2$ constraints coming from $U_{N}$, one per mechanical contact, stabilize the same number of vibrational modes. Because the system is hypostatic, there remain ${\cal N}_{0}\equiv N(d+d_{\rm ex})-Nz/2=N(d_{\rm ex}-\delta z/2)$ zero-frequency modes, where $\delta z=z-2d$ and $d_{\rm ex}$ is the number of extra degree of freedom per particle, i.e., $d_{\rm ex}=1$ for the BP and $d_{\rm ex}=d-1$ for ellipsoids. Above jamming, where $p>0$, the ${\cal N}_{0}$ zero modes are stabilized by the “soft” constraint coming from $\mu_{N}$ whose characteristic stiffness is $k_{R}\sim k_{i}\sim k(p/\Delta)\ll k$, where $k$ is the stiffness associated to $U_{N}$. Hence, the energy scale of these modes remains well separated from that of the ${\cal N}_{3}$ other modes, and we can restrict to the ${\cal N}_{0}$-dimensional subspace of the soft modes. In this space, we have ${\cal N}_{0}=N(d_{\rm ex}-\delta z/2)$ degrees of freedom, and $\mu_{N}$ provides $Nd_{\rm ex}$ constraints, hence the number of degrees of freedom is $N\delta z/2$ less than the number of constraints. When $\delta z\ll 1$, a variational argument was developed in (17, 47) to describe the low-frequency spectrum. It shows that the soft modes are shifted above a characteristic frequency $\omega_{*}^{2}\sim k_{i}\delta z^{2}\sim k_{R}\delta z^{2}\sim\Delta^{-1}p\,\delta z^{2}$, which is reduced by $\sim-p$ by the so-called pre-stress terms, resulting in $\omega_{*}(p)^{2}=c_{1}\Delta^{-1}p\delta z^{2}-c_{2}p$, where $c_{1}$ and $c_{2}$ are unknown constants. Assuming that the system is marginally stable, $\omega_{*}(p)=0$, results in (43)

This explains the universal square root singularity of the contact number $z_{J}$ observed in ellipsoids, BP and several other models (9, 29, 43), as illustrated in Fig. 1.

Eq. (8) holds when $p\ll\Delta$, because in the argument we assumed to be close to jamming ($p\sim 0$) at fixed $\Delta$. On the contrary, when $\Delta\ll p$, the contact number should have the same scaling of spherical particles:

Eqs. (8) and (9) imply that $p$ and $\Delta$ have the same scaling dimension and the following scaling holds:

In the $\Delta\rightarrow 0$ limit, Eq. (10) reduces to Eq. (9), which requires $\gamma=1/2$ and $f(x)\rightarrow x^{1/2}$ for $x\gg 1$. In the $p\rightarrow 0$ limit, we should recover Eq. (8), which requires $f(x)\rightarrow const$ for $x\ll 1$. For the BP, Eq. (10) is confirmed by numerical simulations (43). Assuming that $f(x)$ is a regular function around $x\sim 0$, one can expand it as $f(x)=c_{0}+c_{1}x+\cdots$ and obtains

where $z_{J}=2d+c_{0}\Delta^{1/2}$. This is compatible with previous numerical results of ellipsoids, where $z-z_{J}\sim\Delta^{-0.35\pm 0.1}p$ (48). We can also study the response to shear deformation, which mainly excites the zero-modes (30). Applying the argument in Ref. (18) to the zero-modes and using Eq. (8), the shear modulus $G$ behaves as $G\sim\delta zk_{R}\sim\delta zk_{i}\sim p/\sqrt{\Delta}$, in perfectly agreement with the numerical result (30).

Vibrational spectrum –

The marginal stability argument suggests that ${\cal N}_{0}$ soft vibrational modes can be found in the frequency range $\omega^{*}\lesssim\omega\lesssim\sqrt{k_{R}}$, with $\omega^{*}\sim 0$ due to marginal stability and $k_{R}\sim p/\Delta$, while the remaining ${\cal N}_{3}$ modes have finite frequency at jamming. We now refine the argument to discuss in more details the vibrational density of states $D(\omega)$. It is convenient to define the ${\cal N}\times{\cal N}$ Hessian matrix of the BP model, with ${\cal N}=N(d+d_{\rm ex})$, as the second derivative of the interaction potential $V_{N}$ w.r.t. $\bm{x}_{i}$ and $R_{i}/\Delta$, in such a way that it has a similar scaling of the one of ellipsoids, where $R_{i}/\Delta$ is mapped onto the angular degrees of freedom $\hat{\bm{u}}$.

Then, $D(\omega)$ near jamming can be separated into the following three regions. (i) The lowest band corresponds to the ${\cal N}_{0}=N(d_{\rm ex}-\delta z/2)$ zero modes stabilized by $\mu_{N}$. Their typical frequency is $\omega_{0}^{2}\sim\partial^{2}\mu_{N}/\partial(\Delta^{-1}R_{i})^{2}\sim k_{R}\Delta^{2}\sim\Delta p$. The remaining ${\cal N}_{3}={\cal N}-{\cal N}_{0}=Nz/2$ modes can be split into two bands: (ii) an intermediate band corresponding to the extra (rotational or radial) degrees of freedom ${\cal N}_{1}=N\delta z/2$, with typical frequency $\omega_{1}^{2}\sim\partial^{2}V_{N}/\partial(\Delta^{-1}R_{i})^{2}\sim\Delta^{2}$, and (iii) the highest band corresponding to the ${\cal N}_{2}=Nd$ translational degree of freedom. For $\Delta\ll 1$, the additional degrees of freedom do not strongly affect these modes, and one can apply the standard variational argument of spherical particles (17, 47), which predicts that their typical frequency is $\omega_{2}^{2}\sim\delta z^{2}\sim\Delta$. The resulting $D(\omega)$ differs significantly from that of isostatic packings of spherical particles, which displays a single translational band.

Numerical results for $D(\omega)$ of ellipsoids from (30), of the BP from (43), and analytical results for the perceptron model to be introduced below, are reported in Fig. 2. Details about the simulations of the BP are explained in (43); here we show data for $N=484$ particles, averaged over at least 1000 samples for each state point. As predicted by our theory, $D(\omega)$ consists of three separated bands with characteristic peak frequencies $\omega_{0,1,2}$. Their scaling with $\Delta$, also reported in Fig. 2 at fixed $p$, follows the theoretical predictions $\omega_{0}\propto\Delta^{1/2}$, $\omega_{1}\propto\Delta$ and $\omega_{2}\propto\Delta^{1/2}$. We also find that $\omega_{0}\propto p^{1/2}$ for small $p$, while $\omega_{1,2}$ do not change much with $p$, which is again consistent with the theory. Finally, in Fig. 3 we report the fraction $f_{i}={\cal N}_{i}/{\cal N}$ of modes in each band for the BP, which also follow the theoretical prediction as a function of $\Delta$ and $p$.

Mean field model –

The universality class of isostatic jamming is well understood: it can be described analytically by particles in $d\rightarrow\infty$ (15) or, equivalently, by the perceptron model (24, 25, 26): both models reproduce the critical exponents of isostatic jamming in all dimensions $d$, leading to the conjecture that its lower critical dimension is $d=2$ (49).

We now introduce a new mean field model which describes the universality class of hypostatic jamming in the BP, ellipsoids and many other models of non-spherical particles. The model, which is a generalization of the perceptron, can be solved analytically and, as we shall show, the solution reproduces all the critical exponents of hypostatic jamming. It consists of one tracer particle with coordinate $\bm{x}$ on the surface of the $N$ dimensional hypersphere of radius $\sqrt{N}$, and $M$ obstacles of coordinates $\bm{\xi}_{\mu}$ and “size” $\sigma+R_{\mu}$. The interaction potential between the tracer particle and the obstacles is

where $v(h)=h^{2}\theta(-h)/2$ and the gap variable $h_{\mu}$ is defined as

The $\bm{\xi}_{\mu}$ are frozen variables, and each of their components follows independently a normal distribution of zero mean and unit variance. The dynamical variables are $\bm{x}$ and the $R_{\mu}$, whose variance is controlled by the chemical potential $\mu_{N}$. We fix the value of $k_{R}$ so that $\sum_{\mu=1}^{M}R_{\mu}^{2}=M\Delta^{2}$. In the $\Delta\rightarrow 0$ limit, the system reduces to the standard perceptron model investigated in Ref. (26), while for $\Delta>0$ the $R_{\mu}$ play the same role of the particle sizes in the BP model.

Because the model can be solved by the same procedure of the standard perceptron model, here we just give a brief sketch of our calculation, which will be given in a longer publication. The free energy of the model at temperature $T=1/\beta$ can be calculated by the replica method, $-\beta f=\lim_{n\rightarrow 0}\frac{1}{nN}\log\overline{Z^{n}}$, where $Z=\int d^{N}\bm{x}d^{M}Re^{-\beta V_{N}}$ and the overline denotes the averaging over the quenched randomness $\bm{\xi}_{\mu}$. Here we are interested in the athermal limit $T\rightarrow 0$. Using the saddle point method, the free energy can be expressed as a function of the overlap $q_{ab}=\left\langle\bm{x}^{a}\cdot\bm{x}^{b}\right\rangle/N$, where $\bm{x}^{a}$ and $\bm{x}^{b}$ denote the positions of the tracer particles of the $a$-th and $b$-th replicas, respectively. In the $n\rightarrow 0$ limit, $q_{ab}$ is parametrized by a continuous variable $x\in$, $q_{ab}\rightarrow q(x)$. The function $q(x)$ plays the role of the order parameter and characterizes the hierarchical structure of the metastable states (50). We first calculate the phase diagram assuming a constant $q(x)=q$, which is the so-called replica symmetric (RS) ansatz that describes an energy landscape with a single minimum. The result for $\Delta=0.1$ is shown in Fig. 4. The control parameters are the obstacle density $\alpha=M/N$ and size $\sigma$. If $\alpha$ is small, the tracer particle can easily find islands of configurations $\bm{x}$ that satisfy all the constraints $h_{\mu}>0$: the total potential energy $U_{N}$ and the pressure vanish and the system is unjammed. The overlap $q<1$ measures the typical distance between two zero-energy configurations. Upon increasing $\alpha$, $q$ increases and eventually reaches $q=1$ at $\alpha_{J}$, which is the jamming transition point (Fig. 4). Naturally, due to the additional degrees of freedom when $\Delta>0$, we have $\alpha_{J}(\Delta)>\alpha_{J}(0)$ for equal $\sigma$. For $\sigma>0$, the RS ansatz is stable for all values of $\alpha$ and it describes the jamming transition. For $\sigma<0$ instead, the jamming line is surrounded by a replica symmetry broken (RSB) region where the RS ansatz is unstable. The jamming transition should thus be described by the RSB ansatz where $q(x)$ is not constant, corresponding to a rough energy landscape. The qualitative behavior of the phase diagram is independent of $\Delta$, in particular the jamming line for $\sigma<0$ is always surrounded by a RSB region.

An important observable to characterize jamming is the gap distribution $\rho(h)\equiv\frac{1}{N}\sum_{\mu=1}^{M}\left\langle\delta(h_{\mu}-h)\right\rangle$ that also gives the contact number $z=\int_{-\infty}^{0}dh\rho(h)$. At jamming, $z$ counts the gaps $h_{\mu}$ that are exactly equal to zero. For comparison with numerical results, we introduce the positive gap distribution $g(h)\equiv\theta(h)\rho(h)/\int_{0}^{\infty}dh\rho(h)$, and the force distribution $P(f)\equiv\theta(-h)\rho(h)\frac{\partial h}{\partial f}/\int_{-\infty}^{0}\rho(h)\frac{\partial h}{\partial f}df$, where $f=-h/p$ (corresponding to negative gaps), both normalized to 1. For the standard perceptron model with $\Delta=0$ and $\sigma<0$, jamming is isostatic with $z=1$ (26), and both $g(h)$ and $P(f)$ exhibit a power law behavior (24, 25, 26). In the jammed phase and $\alpha\gtrsim\alpha_{J}$, the system is described by a “regular” full RSB solution where $1-q(x)\sim y_{\chi}^{2}x^{-2}$ for $q(x)\sim 1$, and $g(h)$ and $P(f)$ are regular and finite functions. The prefactor $y_{\chi}$ is predominantly controlled by the contact number $z$, and diverges at isostaticity when $z=1$ (26) and the regular solution breaks down. At $\alpha_{J}$, the model is described by the “jamming” solution where $1-q(x)\sim x^{-\kappa}$, $g(h)\sim h^{-\gamma}$ and $P(f)\sim f^{\theta}$, with critical exponents $\kappa\simeq 1.42$, $\gamma=(2-\kappa)/\kappa$ and $\theta=(3\kappa-4)/(2-\kappa)$ (15, 24, 25, 26). Near $\alpha_{J}$, the regular solution should connect to the jamming solution. This matching argument leads to $z-1\sim p^{1/2}$, which is the same scaling behavior of spherical particles (6).

The situation is completely different if $\Delta>0$. One can show that the contact number at jamming is $z_{J}\geq 1$, meaning that the regular solution persists even at $\alpha_{J}$. Consequently, $g(h)$ and $P(f)$ are finite and regular functions at jamming, and the square-root behavior of the contact number is replaced by $z-z_{J}=c_{\Delta}p$. At $\alpha_{J}$, the regular solution should connect to the jamming solution in the limit of $\Delta\rightarrow 0$. Using the form of the scaling solution derived for $\Delta\rightarrow 0$ in (26) and $z-z_{J}\sim p$ this matching argument leads to the scaling behavior of $g(h)$ and $P(f)$ at $\alpha_{J}$:

with new critical exponents $\mu=\kappa/(4\kappa-4)=0.851$, $\nu=\mu-1/2$, and a universal scaling function $p_{0}(x)$. The scaling analysis also leads to $z_{J}-1\sim\Delta^{1/2}$ and $c_{\Delta}\sim\Delta^{-1/2}$, consistently with the marginal stability argument, Eqs. (8), (11).

The simplicity of the model allows us to derive the analytical form of the density of states $D(\omega)$. As before, we define the Hessian matrix as the second derivatives of the interaction potential $V_{N}$, Eq. (12) w.r.t $x_{i}$ and $R_{\mu}/\Delta$. Using the Edwards-Jones formula for the eigenvalue density $\rho(\lambda)$ (51, 52), the density of states $D(\omega)=2\omega\rho(\omega^{2})$ can be expressed analytically in closed form as a function of $z$, $k_{R}$ and $p$. These quantities should be obtained by solving numerically the full RSB equations but for simplicity, because here we are interested only in the scaling properties of $D(\omega)$, to obtain Fig. 2 we used arbitrary functions $z$, $k_{R}$ and $p$ which are compatible with the analytical scaling derived from the full RSB equation. We find that $D(\omega)$ displays three separate bands (Fig. 2). As in the standard perceptron (24), marginal stability in the full RSB phase implies that the lowest band starts from $\omega=0$ and for small $\omega$, $D(\omega)\sim\omega^{2}$. The lowest band terminates at $\omega_{0}\sim\Delta^{1/2}p^{1/2}$ near which $D(\omega)$ exhibits a sharp peak. At $\omega_{1}\sim\Delta$ a delta peak is found, while the highest band starts from $\omega_{2}\sim\Delta^{1/2}$. The qualitative behavior of $D(\omega)$, and the scaling of $\omega_{0}$, $\omega_{1}$ and $\omega_{2}$ are the same of all the models displaying hypostatic jamming, such as ellipsoids (31, 32) and BP (43). This confirms that the generalised perceptron can reproduce analytically all the critical properties of the hypostatic jamming transition.

As a final check of universality, we test the prediction for the $\Delta$ dependence of the gap distribution function $g(h)$ at the jamming point, (14). In Fig. 5, we show numerical results (obtained as in (43)) for $g(h)$ of the BP model at $p=10^{-5}$, a value small enough to observe the critical behavior. Here, as usual for particle systems, $g(h)$ is normalized by $g(h)\rightarrow 1$ for larger $h$. When $\Delta=0$, $g(h)$ exhibits a power law divergence, $g(h)\sim h^{-\gamma}$, where $\gamma=0.413$, consistently with previous numerical observation (6, 14, 15). For finite $\Delta$, on the contrary, the divergence of $g(h)$ is cutoff (Fig. 5), consistently with the theoretical prediction of (14).

Conclusions –

Using a marginal stability argument, we derived the scaling behavior of the contact number $z$ and the density of states $D(\omega)$ of ellipsoids and breathing particles. Our theory predicts that the scaling behaviors of the two models are identical, which we confirmed numerically. Many other models of non-spherical particles display the same jamming criticality (40), which defines a new universality class of hypostatic jamming. We introduced an analytically solvable model which allows us to derive analytically the critical exponents associated to the new universality class.

One of the most surprising output of our theory is the universality of the density of states $D(\omega)$ (Fig. 2). This might be relevant for some colloidal experiments where the constituents are non-spherical (53), in which the vibrational modes could be experimentally extracted from the fluctuations of positions (54, 55). Another relevant question is how non-spherical particles would flow under shear (30). The divergence of the viscosity at jamming is related to the low eigenvalues of $D(\omega)$ (56), which suggests that the shear flow of non-spherical particles should be quite different from that of spherical particles, in agreement with recent experiments (57).

We thank B. Chakraborty, A. Ikeda, J. Kurchan, S. Nagel and S. Franz for interesting discussions. This project has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement n.723955-GlassUniversality). This work was supported by a grant from the Simons Foundation (#454953, Matthieu Wyart and #454955, Francesco Zamponi) and by “Investissements d’Avenir” LabEx PALM (ANR-10-LABX-0039-PALM, P. Urbani). We thank the authors of Refs. (40) and (32) for sharing their data used in Figs. 1 and 2.

References