Mean field theory of jamming of nonspherical particles

TL;DR

This paper extends mean field theory to analyze the jamming transition of nonspherical particles, revealing different scaling behaviors and providing detailed insights into the full scaling behavior near the jamming point.

Contribution

It develops a mean field model incorporating rotational degrees of freedom to better understand the jamming of nonspherical particles, expanding on previous work.

Findings

Different contact number and shear modulus scaling near jamming for nonspherical particles.

Detailed analysis of the model's scaling behavior in convex and non-convex phases.

Extension of perceptron-based models to include rotational degrees of freedom.

Abstract

Recent computer simulations have uncovered the striking difference between the jamming transition of spherical and non-spherical particles. While systems of spherical particles are isostatic at the jamming point, systems of nonspherical particles are not: the contact number and shear modulus of the former exhibit a square root singularity near jamming, while those of the latter are linearly proportional to the distance from jamming. Furthermore, while our theoretical understanding of jamming of spherical particles is well developed, the same is not true for nonspherical particles. To understand jamming of non-spherical particles, in the previous work [Brito, PNAS 115(46), 111736], we extended the perceptron model, whose SAT/UNSAT transition belongs to the same universality class of jamming of spherical particles, to include additional variables accounting for the rotational degrees of freedom of nonspherical particles. In this paper, we give more detailed investigations of the full scaling behavior of the model near the jamming transition point in both convex and non-convex phases.

Equations

$$h_{\mu }=\frac{\mathbf{x}\cdot {\mathbf{\xi }}^{\mu }}{\sqrt{N}}-\sigma .$$

$$h_{\mu }\ge 0\forall \mu =1,\dots ,M.$$

$$\mathcal{A}(\Delta )=\frac{S_d(\Delta )}{V_d(\Delta {)}^{\frac{d-1}{d}}}\frac{V_d(0{)}^{\frac{d-1}{d}}}{S_d(0)}-1,$$

$$h(r,\theta )\sim h(r)+\Delta f(\theta ),$$

$$\sigma \to \sigma +\Delta R_{\mu },$$

$$h_{\mu }=\frac{\mathbf{x}\cdot {\mathbf{\xi }}^{\mu }}{\sqrt{N}}-\sigma -\Delta R_{\mu }.$$

$$|\mathbf{R}{|}^2\equiv \sum _{\mu =1}^M{R}_{\mu }^2=M.$$

$$U_N$$

$$v(h)\{\begin{array}{ll}=0& h\ge 0\\ >0& h<0.\end{array}$$

$$V_N=U_N+{\mu }_N,U_N=\sum _{\mu =1}^{M}v(h_{\mu }),{\mu }_N=\frac{k_R}{2}\sum _{\mu =1}^M{R}_{\mu }^2$$

$$Z={\int }_{|\mathbf{x}{|}^2=N}\mathrm{d}\mathbf{x}{\int }_{|\mathbf{R}{|}^2=M}\mathrm{d}\mathbf{R}{e}^{-\beta V_N}$$

$$-\beta f=\underset{n\to 0}{\lim }\frac{\log \overline{Z^n}}{nN},$$

$$-\beta nf\equiv \frac{\log \overline{Z^n}}{N}\approx \frac{1}{2}\log \det Q-n\frac{\alpha }{2}\log \tilde{k}+\alpha \log \left(e^{\frac{1}{2}{\sum }_{ab}{Q}_{ab}\frac{{\partial }^2}{\partial h_a\partial h_b}}{\prod _{a=1}^n{e}^{-\beta {\varphi }_{\mathrm{eff}}(h_a)}|}_{h_a=0}\right),$$

$$e^{-\beta {\varphi }_{\mathrm{eff}}(r)}={\gamma }_{{\Delta }^2/\tilde{k}}\ast e^{-\beta v(r)}.$$

$$\frac{q_m}{\lambda (q_m{)}^2}+{\int }_{q_m}^q\mathrm{d}p\frac{1}{\lambda (p{)}^2}$$

$$\lambda (q)$$

$$\dot{f}(q,h)$$

$$\dot{P}(q,h)$$

$$f(q_M,h)$$

$$P(q_m,h)$$

$$\frac{1}{\lambda (q{)}^2}$$

$$x(q)$$

$$1$$

$$\rho (h)=\frac{1}{M}⟨\sum _{\mu =1}^M\delta (h-h_{\mu })⟩,$$

$$\rho (h)=\frac{1}{\alpha }\frac{\delta f}{\delta v(h)}=e^{-\beta v(h)}\int \mathrm{d}zP(q_M,z)e^{-f(q_M,z)}{\gamma }_{1-q_M+{\Delta }^2/k_R}(z-h).$$

$$z$$

$$p$$

$$g(h)$$

$$P(f)$$

$$\frac{q}{(1-q{)}^2}=\alpha \int \mathrm{d}h{P}_{\mathrm{RS}}(q,h)f_{\mathrm{RS}}^{\prime }(q,h{)}^2,$$

$$f_{\mathrm{RS}}(q,h)$$

$$P_{\mathrm{RS}}(q,h)$$

Full text

• March 2019

1 Introduction

A system consisting of macroscopic particles that are so large that the thermal fluctuations are negligible, at small packing fraction $\varphi$ easily deforms and flows under the action of external driving forces, such as shear and shake. Upon increasing $\varphi$, the constituent particles begin to be in contact and the system has a finite rigidity at a certain critical packing fraction $\varphi=\varphi_{J}$. This non-equilibrium fluid-solid transition is the so-called jamming transition and $\varphi_{J}$ is referred to as the jamming point [1, 2]. The jamming transition has been the focus of an intense research activity, as it is ubiquitously observed in a wide range of engineering and biological systems such as metallic balls [3], foams [4, 5], colloids [6], polymers [7], candies [8], dices [9] and tissues [10].

One of the most popular and simple models to study the jamming transition is a system consisting of frictionless spherical particles with purely repulsive interaction. A large number of numerical works have been performed on this model which have proven that (i) this system is isostatic at $\varphi_{J}$, meaning that the contact number between constituent particles is identical to the number the degrees of freedom, so that the system barely achieves mechanical stability [3], (ii) the distributions of inter particle distances and forces exhibit power law behaviors at $\varphi_{J}$ [11, 12, 13], (iii) various physical quantities, such as the shear modulus and contact number, exhibit power law scaling near $\varphi_{J}$ [5], and (iv) a set of lengthscales diverging at $\varphi_{J}$ [14, 15, 16]. These results suggest that the jamming transition point can be considered as a sort of out-of-equilibrium critical point.

The theoretical understanding of jamming of spherical particles is also quite advanced. On the one hand, M. Wyart et al. developed a variational argument showing that isostaticity plays a crucial role in determining the critical properties of the jamming transition [17, 18]. On the other hand, a system consisting of spherical particles has been solved exactly in the large dimension limit by using the replica liquid theory (RLT) [19, 20]. The RLT is a theory that combines the density functional theory of classical amorphous solids [21, 22, 23] with the replica method [24], which is a standard method to treat disordered systems originally developed to study spin glasses. The RLT predicts that a system near $\varphi_{J}$ is located in the (full) replica symmetric breaking (RSB) phase [25]. In the RSB phase, the system is marginally stable and has a gapless spectrum of excitations [26, 27]. Interestingly, this marginal stability is directly related to the isostaticity of the jamming transition point [28]. The exact solution of hard spheres provides the exact values of the critical exponents of the gap and force distribution function at $\varphi_{J}$ in infinite dimensions, whose values agree well with the numerical results in two and three dimensions within the numerical precision [28, 29, 30, 31].

Furthermore in recent years, Franz and Parisi have proposed to look at the jamming transition point as the satisfiability transition point for continuous constraint satisfaction problems [32]. In a constraint satisfaction problem one has to find a configuration of a set of dynamical variables that satisfy a set of constraints. As the number of constraints increases, one can observe a sharp phase transition from a satisfiable (SAT) phase, where there are configurations of the variables satisfying all the constraints, to an unsatisfiable (UNSAT) phase, where no such configurations exists and a finite fraction of constraints are violated [33, 34, 35]. The jamming transition of spherical particles can be considered as a SAT-UNSAT transition within a glass: in the unjammed phase, one can find disordered configurations satisfying the constraints that spheres do not overlap, meaning that the system is in the SAT phase, while in the jammed phase, one can not find such configurations and the number of contacts between particles is finite, meaning that the system is in the UNSAT phase 111Note that the jamming and glass transitions are different phenomena: the glass transition is defined as the point at which the system lose the ergodicity, while the jamming transition is the point at which the particle can not avoid overlap. . In [32] Franz and Parisi considered the simplest continuous constraint satisfaction problem, the perceptron, and showed that the gap and force distributions close to the non-convex SAT/UNSAT (jamming) threshold have the same critical exponents of those of hard spheres in the large dimension limit, implying that the two models belong to the same universality class. The simplicity of the perceptron also allows one to calculate the contact number and the density of states in the jammed (UNSAT) phase, and successfully reproduce the numerical results on jamming of spheres in finite dimensions (even if the origin of this universality is still unclear).

However, spherical particles are an idealized model. In general, asphericity is inevitably introduced in realistic situations (where also friction plays a role, although it will be always neglected here). It has been shown in the last years that the particle shape affects the properties of jamming. Numerical simulations on jamming of ellipsoids shows that the contact number $z_{J}$ at jamming smoothly increases as the asphericity increases, meaning that the system is not isostatic [8, 36]. The breakdown of isostaticity is also observed for spherocylinders [37], superballs [38], superellipsoids [39], other convex shaped particles [40] and even deformable polygons [41]. Furthermore, detailed numerical simulations showed that the critical exponents of the contact number and shear modulus near the jamming transition point differ quite significantly from those of spherical particles [36, 42]. This difference suggests that the universality class of the jamming transition of nonspherical particles is different from the one of spherical objects.

In recent years, a series of phenomenological and theoretical approaches have emerged to rationalize jamming of non-spherical particles. For instance, by analyzing the stability matrix, Donev et al. [43] argue that the breakdown of the isostaticity is caused by the positiveness of the pre-stress term. By using the Edwards ensemble approach, Baule et al. [44] calculated the contact number and the jamming transition point of nonspherical particles, which can qualitatively reproduce the numerical results. However, to the best of our knowledge there is no unified theoretical understanding that explains the origin of the different critical exponents and universality for nonspherical particles. In this work, we study a model that we have recently proposed in a shorter report [45] as a simple exactly solvable mean-field model for the jamming transition of nonspherical particles and we detail its analytical solution. We consider the perceptron model, introduced in [32] to study the jamming transition of spherical particles, and extend it to include additional dynamical variables, playing the same role as the rotational degrees of freedom of nonspherical particles [45]. The model and its scaling properties can be investigated analytically by using the replica method, as in the case of the original perceptron model [46]. We find that the gap and force distributions do not exhibit power law behavior and remain regular even at the jamming transition point. This is in marked constrast with the original perceptron model or spherical particles, where the gap and force distributions exhibit power laws [46]. We find that the absence of criticality of the gap and force distributions are originated from the breakdown of isostaticity at the jamming transition point, as predicted in [17, 18]. Furthermore, the regularity of these distributions leads to a trivial linear distance dependence of the physical quantities near the jamming transition point. We find that the critical exponents change discontinuously at zero asphericity, meaning that an infinitesimal asphericity is enough to alter the universality of the jamming transition.

The paper is organized as follows. In Sec. 2, we introduce the model and in Sec. 3, we derive the free energy and saddle point equations with the replica method. In Sec. 4, we calculate the phase diagram using the replica symmetric ansatz. In Sec. 5, we argue the scaling behavior of the model above and at the jamming transition point and in Sec. 6, we calculate the density of states in the UNSAT(jammed) phase and compare it with previous numerical results on ellipsoids. Finally in Sec. 7, we summarize our results and give some perspectives. We discuss the technical details in the appendix.

2 The polydisperse perceptron model

Jamming of particles can be regarded as a special case of a more general class of problems where one needs to find an assignement of a large set of continuous variables to satisfy a set of constraints [32]. These problems are therefore called Continuous Constraint Satisfaction Problems (CCSP). In this general setting, there is a precise dictionary between jamming of spheres and CCSP that allows to identify several physical quantities, from the pressure to gap variables and number of contacts, in both class of problems. Here we will not review this dictionary extensively but we will limit ourselves to discuss it online while presenting the application of this line of reasoning to study jamming of non-spherical particles. The interested reader can find more details in [26, 46].

The perceptron is a well known and studied linear classifier in machine learning [47, 48, 49, 50]. Here, following [32], we turn it into a continuous constraint satisfaction/optimization problem to study the jamming transition. The model studied in [32] is defined in terms of a state vector $\bm{x}=\{x_{1},\ldots,x_{N}\}$ which is $N$-dimensional and lives on the $(N-1)$-dimensional hypersphere defined by $|\bm{x}|^{2}=N$. Furthermore one defines a set of $M=\alpha N$ quenched $N$-dimensional random vectors $\{\bm{\xi}^{\mu}=\{\xi^{\mu}_{1},\ldots,\xi^{\mu}_{N}\}\}_{\mu=1,\ldots,M}$ whose components are Gaussian random variables with zero mean and unit variance. Given each random vector (often referred to as “pattern”) and the state vector of the system, one can define a gap variable as

[TABLE]

where $\sigma$ is a control parameter. The constraint satisfaction problem is defined by asking for a vector $\bm{x}$ that satisfies all the constraints

[TABLE]

On general grounds one can expect that if $\alpha$ is small enough it is possible to find a satisfiable (SAT) configuration of $\bm{x}$ while if $\alpha$ is large enough this is not possible and the problem becomes unsatisfiable (UNSAT). Indeed, in the thermodynamic limit $N\to\infty$, there exists a sharp SAT/UNSAT phase transition that marks the boundary between the two phases. The parameter $\sigma$ controls the convexity of the problem: if $\sigma>0$ the constraint satisfaction problem is convex while if $\sigma<0$ it is non-convex. In [32] it has been shown that for $\sigma<0$ the problem has a replica symmetry breaking phase near the jamming point as hard spheres in infinite dimensions [25], and the SAT/UNSAT transition is in the same universality class of jamming of spherical particles.

In this work, the perceptron model is extended to describe the jamming of nonspherical particles. We consider the polydisperse perceptron model that has been introduced in [45] and here we first recall its definition and then we construct its full analytical solution. For non-spherical particles like ellipsoids, the relative distance between constituent particles depends on their relative angles. Therefore the actual gap between particles is a function of both the distance of their centers and relative orientation. More precisely, consider non-spherical particles, with $\Delta$ being the linear deviation from spherical shape (for example, in ellipsoids $\Delta$ is the deviation of the aspect ratio from 1) and with the asphericity parameter being

[TABLE]

where $S_{d}(\Delta)$ and $V_{d}(\Delta)$ represent the surface and volume of a $d$-dimensional particle, respectively. Note that because the function $\mathcal{A}(\Delta)$ has a minimum at $\Delta=0$, it can be expanded as $\mathcal{A}(\Delta)=\mathcal{A}^{\prime\prime}(0)\Delta^{2}/2+\cdots$, which leads to $\Delta\sim\mathcal{A}^{1/2}$. Note also that for small $\Delta$, the gap between two particles can be written as [45]

[TABLE]

where $r$ is the distance between the centers, $h(r)=r-\sigma$ is the gap for spherical particles of diameter $\sigma$, and $f(\theta)$ is some function of the non-spherical degrees of freedom.

One can think about the orientation of a particle as a dynamical internal degree of freedom that each particle can change in order to satisfy better the constraints. In this way one can construct more general models in which beyond the standard degrees of freedom that are the position of the particles, one can consider internal degrees of freedom that must have the same role as orientations for non-spherical particles. In practice, any internal degree of freedom such that Eq. (4) holds should lead to similar results. For example, in [45] we have shown that in the small asphericity limit, the interaction potential of ellipsoids maps to the interaction potential of a model of breathing spheres, namely spheres that can change their diameter [51]. This suggests that one can extend the perceptron model studied in Ref. [32, 46], by replacing the effective diameter $\sigma$ by a fluctuating one [45]

[TABLE]

so that the new gap variables become

[TABLE]

The variables $\{R_{\mu}\}_{\mu=1,\ldots M}$ are additional internal dynamical variables that can be used to find solutions to the constraint satisfaction problem. We bound their amplitude by enforcing them to live on a $(M-1)$-dimensional hypersphere

[TABLE]

Because Eq. (6) has the same structure of Eq. (4), we expect this modified perceptron to fall into the same universality class of non-spherical particles with $\mathcal{A}\sim\Delta^{1/2}$. The control parameter $\Delta$ thus tunes the asphericity of the problem: in the limit $\Delta\to 0$, the model reduces to the standard perceptron model, which corresponds to the system consisting of spherical particles, except at the jamming transition point. However at the jamming transition point, the model leads to the completely different scaling from the standard perceptron even in the $\Delta\to 0$ limit, as we discuss in the rest of the manuscript.

In order to solve the constraint satisfaction problem defined by the polydisperse perceptron model we can define a cost function

[TABLE]

The local potential $v(h_{\mu})$ has to satisfy the property that

[TABLE]

The choice that we adopt here is the one describing harmonic particles $v(h)=h^{2}\theta(-h)/2$. The particular choice of the cost function does not affect the SAT phase as well as the SAT/UNSAT threshold, while it deeply affects the properties of the jammed phase [46, 52]. In order to impose the spherical constraint of Eq. (7) we add a chemical potential $\mu_{N}$ to the cost $U_{N}$ so that the final interaction potential reads

[TABLE]

where $k_{R}$ is a Lagrange parameter used to enforce Eq. (7). In the next section we study the partition function of the model at inverse temperature $\beta=1/T$ defined as

[TABLE]

and we study the corresponding average free energy using the replica method [24].

3 Free energy and thermodynamic quantities

The free energy of the model can be calculated by using the replica trick:

[TABLE]

where the overline denotes the average over the disorder $\{\bm{\xi}^{\mu}\}_{\mu=1,\ldots,M}$. As usual for integer $n$, one interprets the $n$th power of the partition function as the partition function of a system of $n$ copies or replicas of the original system, all subjected to the same realization of disorder. In this way, using standard manipulations, the free energy of the replicated system can be obtained as

[TABLE]

where $\tilde{k}=\beta k_{R}$, and $Q$ denotes the $n\times n$ overlap matrix $Q_{ab}=\left\langle\bm{x}^{a}\cdot\bm{x}^{b}\right\rangle/N$, $a,b=1,\cdots,n$ that has diagonal elements equal to one due to the spherical constraint on $\bm{x}$ while its off diagonal elements are obtained as a saddle point of Eq. (13). In Eq. (13) we have also introduced the effective potential

[TABLE]

Here $\gamma_{A}(x)$ denotes a Gaussian function of zero mean and variance $A$, and the star product denotes the convolution, $f*g(x)=\int_{-\infty}^{\infty}\mathrm{d}yf(y)g(x-y)$. The parameter $k_{R}$ should also be determined by a saddle point of Eq. (13), in order to enforce the constraint in Eq. (7). The saddle point equations for the overlap matrix $Q$ can be solved only assuming some ansatz that allows one to take the analytic continuation $n\to 0$. Following the standard strategy of replica theory, we assume a hierarchical ansatz for $Q$ [53], in which $Q$ is encoded by a continuous function $q(x)$ defined for $x\in[0,1]$ with the boundary conditions $q(x)=q_{m}$ for $x\in[0,x_{m})$ and $q(x)=q_{M}$ for $x\in(x_{M},1]$. It is practical to define $x(q)$ as the inverse function of $q(x)$ for $x\in[x_{m},x_{M}]$. After a standard but slightly lengthy procedure (see A and [46]), one obtains the saddle point condition for $x(q)$ as the solution of:

[TABLE]

where

[TABLE]

In Eq. (15), $f(q,h)$ and $P(q,h)$ are the solutions of the Parisi equations

[TABLE]

where $\dot{g}(q,h)=\partial g(q,h)/\partial q$ and $g^{\prime}(q,h)=\partial g(q,h)/\partial h$, and the boundary conditions are

[TABLE]

In the continuous RSB phase, $x(q)$ becomes a monotonic function of $q$, suggesting that Eq. (15) can be differentiated w.r.t. $q$ [46, 54]. The first and second derivatives lead to

[TABLE]

Finally, the spherical constraint for $R_{\mu}$, Eq. (7), reduces to

[TABLE]

Eqs. (15-21) represent the set of saddle point equations through which we can obtain the free energy of the model.

For later convenience, we introduce several thermodynamic quantities that we shall discuss in this manuscript. First we define the gap distribution as the following Boltzmann average

[TABLE]

which can be calculated from the first derivative of the free energy w.r.t. the interaction potential [46]:

[TABLE]

Using $\rho(h)$, one can calculate the isostaticity index as

[TABLE]

The value of $z$ equals one when the system is isostatic, namely, when the number of UNSAT gaps is equal to the dimension of $\bm{x}$. Note that the normalization of Eq. (24) is $N$, which is not the total number of degrees of freedom given by $N+M$ 222 The reason for this normalization is the following. One can easily show [17, 26, 45] that the matrix of second derivatives of $U_{N}$ has precisely $N\max(1-z,0)$ zero modes. For $\Delta=0$, i.e. in absence of $\mu_{N}$, stability then requires $z\geq 1$. Jamming has the minimal number of constraints and is then isostatic, $z=1$, corresponding to a number of constraints precisely equal to the number of degrees of freedom [17]. As soon as $\Delta>0$, one has $M+N$ degrees of freedom, and the matrix of second derivatives of $U_{N}$ has then $N\max(\alpha+1-z,0)$ zero modes. Isostaticity would then correspond to $z=1+\alpha$, but the term $\mu_{N}$ can stabilize some of the zero modes of $U_{N}$. Hence, the system is not constrained to be isostatic at jamming. We will show analytically that when $\Delta$ is small $0<z-1\sim\Delta^{1/2}<\alpha$ at the non-convex jamming point, and therefore the system is hypostatic (even if in our notation $z>1$).. The pressure $p$ is also calculated from $\rho(h)$:

[TABLE]

For comparison with numerical results, we introduce the positive gap distribution

[TABLE]

and the force distribution

[TABLE]

where $f=-h/p$ for negative $h$.

4 Phase diagram

4.1 Replica symmetric jamming transition point

We first consider the simplest ansatz for $Q$, namely the replica symmetric (RS) form $Q_{ab}=\delta_{ab}+(1-\delta_{ab})q$. In this case the corresponding saddle point equations reduce to

[TABLE]

where

[TABLE]

We want to consider what happens in the zero temperature limit where the partition function gives either the Gardner volume of solutions of the constraint satisfaction problem [48, 55] in the SAT phase or the ground state energy in the UNSAT phase. In this second case the overlap $q$ parametrizes the typical overlap between two configurations in the ground state basin and therefore for $T\to 0$ we have that $q\to 1$. Therefore we can expand $q$ as [46]

[TABLE]

where $\chi$ is a constant. At the jamming transition point, $\chi$ diverges to infinity since in the SAT phase $q<1$ [46]. For $\beta\gg 1$, we obtain (see Appendix C in Ref. [46])

[TABLE]

Substituting this into Eq. (28), we obtain

[TABLE]

where we have introduced an auxiliary function:

[TABLE]

With a similar calculation, one can show that Eq. (21) reduces to

[TABLE]

Eqs. (32) and (34) can be solved for $k_{R}$,

[TABLE]

which implies that $k_{R}$ vanishes on approaching the jamming transition point as $k_{R}\sim\chi^{-1}$. Substituting this into Eq. (32) and taking the $\chi\to\infty$ limit, we obtain the jamming transition point $\alpha_{J}$:

[TABLE]

The same equation is obtained by investigating the saddle point equations in the unjammed phase, see B. In Fig. 1, we show the typical behavior of $\alpha_{J}$ calculated by Eq. (36) for $\Delta=0.1$.

$\alpha_{J}$ increases with decreasing $\sigma$, which is a reasonable result considering that $\sigma$ corresponds to the effective diameter. For fixed $\sigma$, $\alpha_{J}$ increases as $\Delta$ increases, which is consistent with the numerical results of nonspherical particles of small asphericity [37, 8, 38, 40]. On the contrary, for larger asphericity, the jamming density of a particle system exhibits a non-monotonic behavior [37, 8, 38, 40], which is not captured by this extended perceptron model, which only represents nonspherical particles with small asphericity up to $O(\Delta^{2})=O(\mathcal{A})$ [45].

4.2 RS-RSB phase boundary

The RS ansatz breaks down in the (full) replica symmetric breaking (RSB) phase where $q(x)$ becomes a continuous function of $x$. In this case, the inverse function $x(q)$ is also a continuous function of $q$. Thus, at the RS-RSB phase boundary, $q$ calculated with the RS ansatz should satisfy Eq. (19), which leads to

[TABLE]

Using Eqs. (32) and (37), it can be shown that the RS-RSB boundary in the jammed phase is the vertical line defined by $\sigma=0$, see the blue vertical line in Fig. 1. Similarly, one can calculate the RS-RSB phase boundary by substituting the RS result in the unjammed phase into Eq. (19), see the blue line in the unjammed phase of Fig. 1. The resulting figure shows that the jamming transition line lies in the RSB region when $\sigma<0$, as in the case of the standard spherical perceptron [32].

5 Jamming scaling in the RSB phase

In this section, we discuss the scaling behavior of the jamming transition for $\sigma<0$, where jamming is described by the RSB equations and belongs to the same universality class of particle systems. For the sake of brevity, here we only discuss the scaling in the jammed phase, but the scaling solutions in the unjammed phase can be constructed in a very similar way, see C for the details.

5.1 Scaling solutions in the zero temperature limit

First we derive the asymptotic forms of the relevant equations in the zero temperature limit, $T\to 0$. As in the case of the RS analysis, we expand as $q_{M}=1-T\chi$, where $\chi$ has a finite limit for $T\to 0$ and diverges at the jamming transition point. We introduce the following scaling variables for $\chi\gg 1$ [46]:

[TABLE]

Then, the gap distribution $\rho(h)$, Eq. (23), is

[TABLE]

Using Eq. (24) and (19), we obtain

[TABLE]

Away from the jamming point, it has been shown that $y(q)$ is described by the following scaling solution near $q\sim 1$ [46]:

[TABLE]

which we refer to as the “regular” solution. For $m(q,h)$, we assume

[TABLE]

Using Eq. (20), one can calculate $\mathcal{Y}$ as

[TABLE]

5.2 Jamming for $\Delta=0$

For self-completeness, we first review the scaling behavior of the standard perceptron model ($\Delta=0$), which corresponds to the jamming of spherical particles and has been already well investigated in previous work [46]. In this case, the isostaticity index defined in Eq. (40) reduces to [46]

[TABLE]

At the jamming transition point, $\chi\to\infty$, one obtains

[TABLE]

meaning that the system becomes isostatic: the number of contacts is the same of that of the degree of freedom. Away from the jamming transition point, $z$ increases as

[TABLE]

For $\Delta=0$, Eq. (43) reduces to

[TABLE]

and for $\chi\gg 1$, we get

[TABLE]

implying that the scaling solution Eq. (41) breaks down when the system becomes isostatic and $\delta z=0$. Indeed in the isostatic limit, the scaling ansatz for $y(q)$ becomes:

[TABLE]

where the value of the critical exponent $\kappa=1.41574\ldots$ is obtained by solving the Parisi equations in the scaling regime [46]. Hereafter, we shall refer to Eq. (49) as the “critical” scaling solution. The matching argument between the regular and critical scaling solutions determines the asymptotic behavior of $P(1,h)$ as [46]

[TABLE]

where

[TABLE]

The critical exponents are

[TABLE]

Using Eqs. (25), (39) and (50), the pressure is calculated as

[TABLE]

leading to

[TABLE]

This is consistent with the numerical results of the jamming transition of spherical particles interacting with the harmonic potential [5]. Using Eqs. (26), (27), (39), and (50), we arrive at the asymptotic forms of the positive gap and force distributions:

[TABLE]

where

[TABLE]

Eqs. (55) show that the gap and force distributions exhibit a power law behavior if the system is isostatic, $\delta z=0$, while they remain finite and regular if the system is not isostatic, $\delta z\neq 0$.

5.3 Jamming for $\Delta>0$

We now discuss the scaling behavior for $\Delta>0$, corresponding to the jamming of nonspherical particles. We first study the equations in the jamming limit, $\chi\to\infty$. The isostaticity index at the jamming transition point is

[TABLE]

implying that

[TABLE]

Eq. (43) at the jamming point becomes

[TABLE]

Note that, unlike for jamming of spherical particles, $\mathcal{Y}\sim KP(1,0)\sim\delta z^{-1}P(1,0)$ remains finite even at jamming, implying that the regular scaling solution, described by Eq. (41) persists even at the jamming transition. The regular solution connects to the critical solution described by Eq. (49), in the spherical limit $\Delta\to 0$. It is worth noting that Eqs. (44) and (47) can be identified with Eqs. (57) and (59), if one replaces $\chi$ with $K$, implying that the scaling form of $P(1,h)$ for $K\gg 1$ is obtained by simply replacing $\chi$ in Eq. (50) with $K$ :

[TABLE]

Using this scaling form Eq. (21) reduces to

[TABLE]

which, for $\Delta\ll 1$, leads to

[TABLE]

From this equation and Eq. (58), we have

[TABLE]

which agrees with numerical results for jamming of nonspherical particles [36, 40] if one identifies $\Delta$ with the square root of the asphericity [45].

Next, we discuss the scaling behavior above the jamming transition point, where $\chi$ is large but finite. The pressure can be represented using $\chi$ and $K$ as

[TABLE]

implying

[TABLE]

Then, the isostaticity index, Eq. (40), can be expanded as

[TABLE]

where $c_{1}$ and $c_{2}$ are constants. Note that $z$ linearly depends on $p$. This is consistent with numerical results for the jamming transition in ellipsoids [56] and in marked contrast with the standard perceptron, corresponding to spherical particles, where $\delta z\sim p^{1/2}$ according to Eq. (54). Note that $\delta z$ should converge to Eq. (54) in the $\Delta\to 0$ limit, which requires the following scaling form:

[TABLE]

where $\mathcal{Z}(x)\to\rm{const}$ for $x\ll 1$ and $\mathcal{Z}(x)\to x^{1/2}$ for $x\gg 1$ [45]. Eq. (67) implies that, upon approaching the jamming transition point, the scaling behavior of $z-z_{J}$ is changed qualitatively at $p\sim\Delta$ from the scaling of spherical particles $z-z_{J}\propto p^{1/2}$ to that of nonspherical particles $z-z_{J}\propto p$, as shown in Fig. 2.

From Eqs. (38), (62), and (65), we have

[TABLE]

Finally, by substituting Eq. (65) into Eq. (60), we can calculate $g(h)$ and $P(f)$ at the jamming transition point. Interestingly, the resulting equations are identical to Eqs. (55), meaning that the criticality of $g(h)$ and $P(f)$ is controlled by $\delta z$ only. Because $\delta z>0$ for nonspherical particles, $g(h)$ and $P(f)$ remain finite and regular functions even at the jamming transition point.

6 Vibrational density of states

In this section, we discuss the scaling of the vibrational density of states $D(\omega)$ of the polydisperse perceptron. We focus on the region $\sigma<0$ where the model falls into the same universality class of nonspherical particles, as shown in the previous sections. In the UNSAT phase the ground state is a configuration of the state vectors $\bm{x}$ and $\bm{R}$ that minimizes the potential $V_{N}$. By defining $\rho(\lambda)$ as the eigenvalue distribution of the Hessian of $V_{N}$ computed in the ground state, the density of states is given by

[TABLE]

To calculate $\rho(\lambda)$, we consider a slightly modified interaction potential:

[TABLE]

where $\zeta$ denotes the Lagrange multiplier that is needed to impose the spherical constraint $|\bm{x}|^{2}=N$. In the UNSAT phase, the first term of Eq. (70) is a sum over the UNSAT constraints, which are the $Nz$ contacts defined as the gaps that satisfy $\theta(-h_{\mu})=1$. For convenience, we hereafter reassign the index $\mu$ only to these gaps $\mu=1,\cdots Nz$. Note that in the following, we neglect the $(\alpha-z)N$ degrees of freedom $R_{\mu}$ associated to the SAT constraints, such that $\theta(-h_{\mu})=0$. These degrees of freedom are trivially decoupled from the system and, if included, would give rise to a delta function $\delta(\lambda-k_{R})$ in the density of states. Using this convention, we can write a reduced $N(1+z)\times N(1+z)$ Hessian of the potential $V_{N}$ as

[TABLE]

We define the Hessian matrix as

[TABLE]

One can calculate the eigenvalue distribution $\rho(\lambda)$ of $\mathcal{H}$ by using the Edwards-Jones formula [57]:

[TABLE]

where the overline denotes the average over the quenched disorder $\xi_{i}^{\mu}$, and we have introduced a partition function as

[TABLE]

Here $I_{d}$ denotes the $d$-dimensional unit matrix. Due to the mean field nature of the model, we can derive the analytical expression of $\rho(\lambda)$ (see D):

[TABLE]

where the explicit expressions of $\lambda_{1},\cdots,\lambda_{4}$ are given in D. To calculate $\rho(\lambda)$, one needs to calculate $\zeta$, $z$, and $k_{R}$ for given $p$ and $\Delta$. In the RSB phase, there is a useful relation to express $\zeta$ (see E):

[TABLE]

To calculate $z$ and $k_{R}$, one should solve the RSB equation numerically, which is a difficult task. Since we are mainly interested in the scaling property of the density of states, here we assume a suitable ansatz for $z$ and $k_{R}$ of the same form discussed in the previous sections. The scaling of $z$ given by Eq. (67), can be satisfied by assuming

[TABLE]

where $c_{1}$ and $c_{2}$ are constants. To satisfy Eq. (68), we assume

[TABLE]

where $c_{3}$ is another constant. Because these constants do not affect the scaling behavior, we set them to $c_{1}=c_{2}=c_{3}=1$. Substituting Eqs. (77) and (78) into Eq. (75), and using Eq. (69), one can calculate $D(\omega)$.

In Fig. 3, we show the typical behavior of $D(\omega)$, which consists of three regions separated by finite gaps:

The lowest band is quasi-gapless and ends at $\omega_{0}$. The existence of gapless excitations is a direct consequence of RSB as already discussed for the standard perceptron [26]. The weight of the lowest band is $f_{0}\equiv\int_{0}^{\omega_{0}^{2}}d\lambda\rho(\lambda)=\int_{0}^{\omega_{0}}d\omega D(\omega)=1/(1+z)$. 2. 2.

The delta function at $\omega_{1}$. The weight is $f_{1}=(z-1)/(1+z)$. 3. 3.

The highest band starting at $\omega_{2}$. The weight of this band is $f_{2}=1/(1+z)$.

Note that $D(\omega)$ is normalized so that $\int_{0}^{\infty}d\omega D(\omega)=f_{0}+f_{1}+f_{2}=1$. The behavior of $D(\omega)$ that we find in the polydisperse spherical perceptron resembles very closely to that obtained in numerical simulations of ellipsoids [42], where again $D(\omega)$ consists of three well-separated regions. To uncover further similarity, we discuss the scaling behavior of the characteristic frequencies, $\omega_{0}$, $\omega_{1}$, and $\omega_{2}$.

In Fig. 4, we show the $p$ and $\Delta$ dependence of the characteristic frequencies, from which one can deduce the following scaling behavior:

[TABLE]

The same scaling can also be obtained directly by the asymptotic analysis of $\lambda_{1},\cdots,\lambda_{4}$. For $p\ll\Delta$ we find the same scaling reported in Ref. [42] for the case of ellipsoids. Instead it seems that the aspect ratios used in Ref. [42] are too large to observe the scaling for $p\gg\Delta$. It would be interesting to redo the numerical simulation to further investigate this case. For the lowest frequency regime, we get that $\rho(\lambda)\sim\sqrt{\lambda}\theta(\lambda)$, which leads to $D(\omega)\sim\omega^{2}$, see Fig. 3 as in the case of the spherical perceptron [26].

7 Summary and discussions

In this manuscript, we constructed a mean field theory for the jamming transition of nonspherical particles based on the analytical solution of the polydisperse spherical perceptron model. This model is an extension of the spherical perceptron model introduced in [32] to study jamming of spherical particles. In order to take into account the particles shape, which results in internal degrees of freedom such as the orientations, we introduced additional internal degrees of freedom in the perceptron model [45]. We can parametrize the asphericity of particles using an additional control parameter $\Delta$. We solved the model through the replica method and we determined the phase diagram and the asymptotic behavior of the gap (and force) distributions near and at the jamming transition point. The resulting generic picture is showed in the left panel of Fig. 5, where we summarize the phase diagram predicted by our theoretical calculation. For $\Delta=0$, which corresponds to spherical particles, the system is isostatic at the jamming point, and the gap and force distributions show a power law behavior [32]. On the contrary, for $\Delta>0$, which corresponds to nonspherical particles of finite asphericity, the system is not isostatic at the jamming transition point, and the gap and force distributions remain regular and finite. In the right panel of Fig. 5, we show the critical exponents of the shear modulus and contact number. Due to the regularity of the distribution function, the critical exponent takes a rather simple value for $\Delta>0$: the shear modulus and contact number are linearly proportional to the proximity to the jamming transition point. On the contrary, they are proportional to the square root of the proximity to the jamming transition point for spherical particles. Thus, the exponent jumps from one-half to one at $\Delta=0$. We also calculated the density of states $D(\omega)$ and found that our model reproduces the scaling behavior of $D(\omega)$ of ellipsoids near the jamming transition point, previously reported in numerical simulations. Our results highlight the specificity of the universality of jamming of spherical particles, which holds only when the asphericity is precisely zero, whereas the universality of nonspherical particles is far more general as it holds no matter how small the asphericity is.

There are still several important points that deserve further investigation. Here we give a tentative list:

• •

In previous research, we confirmed our theoretical prediction for the gap distribution function by computer simulation of breathing particles, which belong to the same universality class of nonspherical particles [45]. It would be important to repeat the same analysis for various shapes of nonspherical particles to confirm the absence of criticality also in those cases.

• •

$D(\omega)$ of nonspherical particles differs dramatically from that of spherical particles. An interesting question is how this difference affects the steady-state rheology near the jamming transition point. For spherical particles, it has been shown that the minimal eigenfrequency of $D(\omega)$ controls the viscosity near the jamming transition point [58]. It would be interesting to extend the theory in Ref. [58] to the present case to study how the rheology of nonspherical particles differs from that of spherical ones.

• •

For the low frequency regime of the density of states, our mean field model predicts a quadratic scaling $D(\omega)\sim\omega^{2}$. However, in finite dimensions, it has been shown theoretically that spatial heterogeneity modifies the mean-field result below some characteristic frequency, and the quadratic scaling is replaced by a quartic scaling $D(\omega)\sim\omega^{4}$ [59]. This is indeed consistent with numerical observations in breathing particles, which can be considered as nonspherical particles of small asphericity [60]. It would be interesting future work to test if the same quartic scaling holds for ellipsoids and other nonspherical models.

• •

In this manuscript, we investigated the fully connected mean-field model, which does not have a spatial structure. A natural next step is to include spatial fluctuations and study the divergence of the correlation length near the jamming transition point. This is possible by considering a Kac version of the model, as previously done for the $p$-spin spherical model [61] and hard spheres [62].

• •

Although the scaling behavior of our model near the jamming transition point is entirely different from that of the original perceptron model, the two models have qualitatively the same phase diagram, in particular, the jamming line is surrounded by the replica symmetric breaking (RSB) line for $\sigma<0$, where the model can be mapped into a particle system. This means that RSB occurs before the jamming transition point is reached. It would be interesting to see if signatures of such a transition are present for models of nonspherical particles in finite dimension.

• •

Our approach is justified only for small enough asphericity. For larger asphericity, completely new effects may appear such as (partial) orientational ordering, orientational jamming vs. positional jamming, additional many-body correlations, higher order corrections, those of which can not be captured by our model. To discuss those effects, it is tempting to extend the present calculation for the perceptron to more realistic models of nonspherical particles. A possibility in this direction is to solve nonspherical particle models in the large dimensional limit where one can write down the analytical expression of the free energy as a functional of the density of particle positions and angles, as shown in Ref. [63].

This paper is dedicated to Giorgio Parisi in his 70th birthday; Giorgio is a great source of inspiration for all of us and discussions with him, as well as previous joint work, have been crucial for the development of this research. We particularly thank C. Brito and M. Wyart for previous joint research [45] on which this work is based. We also thank B. Chakraborty, A. Ikeda, J. Kurchan, S. Nagel, S. Franz, and H. Yoshino for interesting discussions related to this work. 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 “Investissements d’Avenir” LabEx PALM (ANR-10-LABX-0039-PALM)(P. Urbani).

Appendix A Full RSB free energy and saddle point equations

Here, we briefly explain the derivation of the replicated free-energy and saddle point equations. For a more complete discussion see Ref. [46].

Substituting a hierarchical replica ansatz [24] into the free energy of Eq. (13), we obtain

[TABLE]

where the dot denotes derivation with respect to $x$, and

[TABLE]

The function $f(x,h)$ is obtained as the solution of the Parisi equation [64]:

[TABLE]

with the boundary condition:

[TABLE]

We assume that for $x\not\in[x_{m},x_{M}]$, $\dot{q}(x)=0$ while $q(x)$ takes a constant value, $q(x)=q_{m}$ for $x\in[0,x_{m})$ and $q(x)=q_{M}$ for $x\in(x_{M},1]$. For $x\in[x_{m},x_{M}]$, $q(x)$ is a monotonic function, and thus one can define the inverse function $x(q)$. Using $x(q)$, Eq. (82) reduces to

[TABLE]

with

[TABLE]

The replicated free energy is

[TABLE]

where

[TABLE]

In order to compute the saddle point for $q(x)$ or equivalently for $x(q)$, one should extremize Eq. (86) w.r.t $x(q)$, with the constraint that $f(q,h)$ satisfies Eq. (84) and (85). To do this one can introduce the Lagrange multiplier $P(q,h)$ [54] as,

[TABLE]

Note that the saddle point conditions for $P(q,h)$ and $P(q_{M},h)$ correctly reproduce Eq. (85) and (83), respectively. The equations for the Lagrange multiplier $P(q,h)$ is obtained taking the functional derivatives w.r.t. $f(q,h)$ and $f(q_{m},h)$:

[TABLE]

The saddle point condition for $x(q)$ is

[TABLE]

In the continuous RSB phase, $x(q)$ is a continuous function for $x\in[q_{m},q_{M}]$, which allows us to calculate derivation of Eq. (90) w.r.t $q$. The first and second derivatives lead to

[TABLE]

Also, the spherical constraint $\sum_{\mu=1}^{M}R_{\mu}^{2}=M$, gives the last saddle point equation for $\tilde{k}$

[TABLE]

Appendix B Replica symmetric analysis in the unjammed phase

Here we investigate the replica symmetric (RS) saddle point equations in the unjammed phase. In this case one has that $e^{-\beta v(h)}=\theta(h)$ so that

[TABLE]

where $\Theta(x)=({\rm erf}(x)+1)/2$. Approaching the jamming point one has that $q\to 1$ and $\tilde{k}\to\infty$. In this limit, using the asymptotic expansion of the error function, we can show that

[TABLE]

Substituting this expression into Eqs. (28) and Eq. (93), we obtain

[TABLE]

and

[TABLE]

which coincides with Eq. (36) obtained from the analysis in the jammed phase.

Appendix C Scaling in the unjammed phase

Here we discuss the scaling behavior when approaching to the jamming transition point from the unjammed phase.

C.1 Failure of the critical solution

We first show that the critical jamming solution for the spherical perceptron does not work for the nonspherical perceptron model when $\Delta>0$. For this purpose, we introduce the following scaling variables:

[TABLE]

where $\varepsilon$ is the linear distance from the jamming transition point. In the critical solution [46], one assumes

[TABLE]

[TABLE]

and

[TABLE]

where

[TABLE]

Substituting these relations into Eq. (93), we get

[TABLE]

which leads to

[TABLE]

This equation can be satisfied only when $\Delta=0$, i.e., in the case of the standard perceptron. A similar equation can be also obtained on approaching to the jamming transition from the jammed phase. We deduce that the critical jamming solution is inconsistent when $\Delta>0$.

C.2 Scaling solution for $\Delta>0$

The jamming scaling for $\Delta>0$ is describe by the regular full RSB solution. We define the scaling solutions as

[TABLE]

For $\varepsilon\ll 1$, we have

[TABLE]

Also we assume that $P(q,h)$ is a regular and finite function. Then, Eq. (93) reduces to

[TABLE]

Thus, the regular scaling solution gives a well-defined expression in the $\varepsilon\to 0$ limit, unlike the critical scaling. Similarly, the pressure $p$ can be obtained from

[TABLE]

which implies that $p\sim\varepsilon^{-1}$. Combining this with Eqs. (105), we can determine the pressure dependence of $q_{M}$ and $\tilde{k}$ for $p\ll 1$:

[TABLE]

Appendix D Derivation of the eigenvalue distribution

We want to compute the eigenvalue distribution $\rho(\lambda)$ of the $(N+zN)\times(N+zN)$ dimensional Hessian defined by:

[TABLE]

where

[TABLE]

for $i,j=1,\cdots,N$ and $\mu,\nu=1,\cdots,zN$. We perform the computation using the Edwards-Jones formula [57]:

[TABLE]

where the overline denotes the average over the quenched disorder $\xi_{i}^{\mu}$, and we have introduced the partition function as

[TABLE]

Here $I_{d}$ denotes the $d$-dimensional identity matrix. Performing the Gaussian integration we get

[TABLE]

where we have introduced a $N\times N$ matrix as

[TABLE]

Replacing the quenched average by the annealed one and using the saddle point method, we can show that

[TABLE]

where $q=N^{-1}\sum_{i=1}^{N}u_{i}^{2}$ is to be determined by the saddle point condition:

[TABLE]

This can be solved for $q$:

[TABLE]

where the sign is to be chosen so that $\rho(\lambda)$ is positive and normalized. Substituting Eqs. (114) and (116) into Eq. (73), we have

[TABLE]

After a straightforward calculation, we finally get

[TABLE]

where

[TABLE]

Since $z>1$ in the RSB jammed phase, Eq. (120) can be slightly simplified as

[TABLE]

Appendix E Calculation of $\zeta$

We here determine the Lagrange multiplier $\zeta$ introduced in Sec. 3. For this purpose, we rewrite the free energy Eq. (88) as

[TABLE]

where $q_{d}=N^{-1}\sum_{i}x_{i}^{2}$. The saddle point condition for $q_{d}$ leads to

[TABLE]

After some manipulations, this can be rewritten as

[TABLE]

Substituting $q_{M}\sim q_{d}-T\chi$, $f(q_{M},h)\sim-\beta h^{2}\theta(-h)/(1+\chi+\Delta^{2}/k_{R})$, and taking the zero temperature limit $T\to 0$, we get

[TABLE]