Infinitesimal asphericity changes the universality of the jamming transition

TL;DR

This paper investigates how infinitesimal asphericity in particles alters the universality class of the jamming transition, affecting critical exponents and the nature of singularities in force and gap distributions.

Contribution

It provides a comprehensive review and new numerical evidence on how internal degrees of freedom influence the critical behavior of jamming in non-spherical particles.

Findings

Infinitesimal asphericity smooths out power-law singularities in force and gap distributions.

Critical exponents for contact number and shear modulus are derived using variational methods.

Numerical data support the theoretical predictions about the universality class of jamming with internal degrees of freedom.

Abstract

The jamming transition of non-spherical particles is fundamentally different from the spherical case. Non-spherical particles are hypostatic at their jamming points, while isostaticity is ensured in the case of the jamming of spherical particles. This structural difference implies that the presence of asphericity affects the critical exponents related to the contact number and the vibrational density of states. Moreover, while the force and gap distributions of isostatic jamming present power-law behaviors, even an infinitesimal asphericity is enough to smooth out these singularities. In a recent work [PNAS 115(46), 11736], we have used a combination of marginal stability arguments and the replica method to explain these observations. We argued that systems with internal degrees of freedom, like the rotations in ellipsoids, or the variation of the radii in the case of the \textit{breathing} particles fall in the same universality class. In this paper, we review comprehensively the results about the jamming with internal degrees of freedom in addition to the translational degrees of freedom. We use a variational argument to derive the critical exponents of the contact number, shear modulus, and the characteristic frequencies of the density of states. Moreover, we present additional numerical data supporting the theoretical results, which were not shown in the previous work.

Equations

$$V_N$$

$${\sigma }_i(\hat{r},{\mathbf{u}}_i)={\sigma }^{(0)}+f(\hat{r},{\mathbf{u}}_i)\Delta +O({\Delta }^2),$$

$$h_{ij}(\Delta )-h_{ij}(0)$$

$$h_{ij}(0)=r_{ij}-{\sigma }_i^{(0)}-{\sigma }_j^{(0)}.$$

$$V_N$$

$$U_N$$

$$Q_N$$

$$g_i({\mathbf{u}}_i)$$

$$Q_N\sim g_i\sim p\Delta ,$$

$$\Delta \sim (\mathcal{A}-1{)}^{1/2}.$$

$$\delta V_N=\sum _{ij=1}^N\delta {\mathbf{x}}_i\frac{\partial V_N}{\partial {\mathbf{x}}_i\partial {\mathbf{x}}_j}\delta {\mathbf{x}}_j=\sum _{ij=1}^N\delta {\mathbf{x}}_i\left[\frac{{\partial }^{2}v(h_{ij}(0))}{\partial h_{ij}(0{)}^2}\frac{\partial h_{ij}(0)}{\partial {\mathbf{x}}_i}\frac{\partial h_{ij}(0)}{\partial {\mathbf{x}}_j}+\frac{\partial v(h_{ij}(0))}{\partial h_{ij}(0)}\frac{{\partial }^2{h}_{ij}(0)}{\partial {\mathbf{x}}_i\partial {\mathbf{x}}_j}\right]\delta {\mathbf{x}}_j.$$

$$\delta V_N=\frac{k}{2}\sum _{\alpha }^{Nz/2}\delta r_{\alpha }^2,$$

$${\lambda }_{\min }\sim k\underset{\delta {\mathbf{x}}_i}{\min }\frac{{\sum }_{\alpha =1}^{Nd}\delta r_{\alpha }^2+{\sum }_{\alpha =Nd+1}^{N\delta z/2}{\left(\delta r_{\alpha }-y_{\alpha }\right)}^2}{{\sum }_{i=1}^N\delta {\mathbf{x}}_i\cdot \delta {\mathbf{x}}_i},$$

$${\lambda }_{\min }\sim k\delta z^2.$$

$${\lambda }_{\min }\sim k\delta z^2-cp,$$

$$\delta z\sim p^{1/2}.$$

$$N_V^c$$

$$N_0=N(d+d_{\mathrm{rot}})-N_V^c=N{d}_{\mathrm{rot}}-\frac{N}{2}\delta z,$$

$$N_Q^c-N_0=\frac{N}{2}\delta z>0.$$

$${\lambda }_{\min }\sim k_R\delta z^2\sim c_{1}p\Delta \delta z^2,$$

$${\lambda }_{\min }\sim p\left[c_1\Delta \delta z^2+c_2{\Delta }^2+O({\Delta }^3)\right],$$

$$\delta z\sim {\Delta }^{1/2}\sim (\mathcal{A}-1{)}^{1/4}.$$

$$\delta z={\Delta }^{1/2}\mathcal{Z}({\Delta }^{-1}p),$$

$$\mathcal{Z}(x)$$

$$z-z_J\sim \frac{p}{{\Delta }^{1/2}},$$

$$\delta V_N\sim \delta Q_N\sim \underset{{\mathbf{y}}_i}{\min }\sum _{i=1}^N\sum _{\alpha =1}^{d_{\mathrm{rot}}}{k}_i^{\alpha }{\left(\delta {\tilde{u}}_i^{\alpha }+y_i^{\alpha }\right)}^2,$$

$$\delta V_N\sim k_R(N{d}_{\mathrm{rot}}-N_0){\left(\frac{\delta \gamma }{\Delta }\right)}^2\sim \frac{Np\delta z\delta {\gamma }^2}{\Delta }.$$

$$G\sim \frac{\delta V_N}{N\delta {\gamma }^2}\sim \frac{p}{{\Delta }^{1/2}}.$$

$$G\sim p^{1/2}.$$

$$G={\Delta }^{1/2}\mathcal{G}({\Delta }^{-1}p),$$

$$\mathcal{G}(x)$$

Full text

• March 2019

1 Introduction

The jamming transition of non-spherical particles is qualitatively different from that of spherical particles [1, 2]. Several experimental and numerical investigations uncover that (i) systems consisting of non-spherical particles are not isostatic at their jamming transition point [3], while systems of spherical particles are isostatic [4], (ii) the pair correlation of non-spherical particles does not exhibit the power law singularity at the jamming transition point [5], while that of spherical particles does [6], and (iii) the critical exponents of non-spherical particles are different from those of spherical particles [7, 8].

The theoretical understanding of the jamming of non-spherical particles is challenging because particles do not hold rotational symmetry. In the previous work [5, 9], we proposed a way to bypass this difficulty by considering the mapping from non-spherical particles to the breathing particles (BP), defined as a system of spherical particles for which their radii are allowed to fluctuate [10]. An advantage of the BP particles is that the model holds the rotational symmetry, and thus, one can apply the same technique developed for the spherical particle without any difficulty. Using the BP, we theoretically and numerically confirmed that the gap and force distributions of non-spherical particles are regular and finite even at the jamming transition point, while those quantities exhibit the power law in the case of spherical particles. Furthermore, we showed that the critical exponents of several physical quantities, such as the contact number, shear modulus, and characteristic frequencies of the density of states, have different values from those of spherical particles. This confirms that the jamming of non-spherical particles belongs to a different universality class from that of spherical particles.

This paper is a longer version of our previous work [5]. We shall give a more straight forward derivation of the scaling functions without mapping to the BP particles, and additional numerical data supporting the theoretical results. The organization of the remaining paper is as follows. In Sec. 2, we develop a variational argument for non-spherical particles. In Sec. 3, we discuss the connection between nonspherial particles and BP. In Sec. 4, we discuss the universal form of the gap and force distributions near the isostatic point. In Sec. 5, we discuss the scaling behavior of the density of states of the BP and show that the characteristic frequencies exhibit the same scaling as non-spherical particles. In Sec. 6, we summarize and conclude the work.

2 Variational argument

Here we derive the scaling functions of non-spherical particles for a small asphericity by using the variational argument [11, 12]. In the previous work [5], we performed this calculation by mapping the Gay-Berne potential, which is a model for ellipsoids, to the breathing particles (BP), which is the model consisting of spherical particles where the radii of particles can vary continuously [10]. In this paper, instead, we present a more direct derivation of the scaling functions of non-spherical particles without using the mapping to the BP model.

2.1 Interaction potential

For concreteness, we consider the following interaction potential:

[TABLE]

where $h_{ij}$ denotes the minimal distance between the $i$-th and $j$-th particles, and $v(h)$ denotes a purely repulsive and finite ranged potential, such as the harmonic potential $v(h)=h^{2}\theta(-h)/2$, where $\theta(x)$ denotes the Heaviside function.

2.2 Perturbation around spherical particles

We first derive a convenient expression of interaction potential for a small asphericity. Non-spherical particle have rotational degrees of freedom in addition to the translational degrees of freedom. We assign the vectors $\bm{x}_{i}$ and a unit vector $\bm{u}_{i}$ to express the position and direction of the $i$-th particle, respectively. The radius $\sigma_{i}$ of a non-spherical particle along the direction $\hat{r}$ varies depending on both $\bm{u}_{i}$ and $\hat{r}$. We shall assume that there is a small parameter $\Delta$ representing the deviation from a sphere. We expand the radius using $\Delta$ as

[TABLE]

where $\sigma^{(0)}$ represents the radius of the reference sphere, and $f(\hat{r},\bm{u})$ represents the coefficient of the first order term. Following the similar procedure, one can expand the gap function, which is the minimal distance between the $i$-th and $j$-th particles. The first order correction of $\Delta$ comes from the change of radii of the $i$-th and $j$-th particles along the direction $\hat{r}_{ij}=(\bm{x}_{i}-\bm{x}_{j})/\left|\bm{x}_{i}-\bm{x}_{j}\right|$, namely,

[TABLE]

where $h_{ij}(0)$ is the gap function of the reference spherical particles:

[TABLE]

Substituting this into Eq. (1) and expanding by $\Delta$, we have

[TABLE]

where

[TABLE]

Here $w_{ij}\Delta^{2}$ denotes the $O(\Delta^{2})$ term of the interaction potential, and we have introduced the auxiliary function as

[TABLE]

From Eqs. (6) and (7), one can show that

[TABLE]

where $p\sim-v^{\prime}(h)$ denotes the pressure.

To understand the physical meaning of $\Delta$, it is convenient to clarify the relation between $\Delta$ and the sphericity $\mathcal{A}$, which represents how far is the shape of a particle from a perfect sphere; $\mathcal{A}=1$ for a perfect sphere and $\mathcal{A}>1$ otherwise. By definition, $\mathcal{A}$ takes a minimal value at $\Delta=0$, which allows us to expand it as $\mathcal{A}(\Delta)=1+\frac{1}{2}\mathcal{A}^{\prime\prime}(0)\Delta^{2}+O(\Delta^{3})$, leading to

[TABLE]

This is a useful relation to compare with numerical and experimental results.

2.3 Variational argument for spherical particles

In this work, we derive the scaling behaviors of non-spherical particles by using the variational argument and assumption of the marginal stability. The variational argument gives a typical amplitude of the minimal eigenvalue near the jamming transition point, while the marginal stability requires that the minimal eigenvalue vanishes. For spherical particles, this approach can reproduce the correct scaling of the excess contact number and shear modulus [11]. Since our argument is very similar to that of spherical particles, we first give a summary of the variational argument of spherical particles.

For spherical particles, the interaction potential is given by the $\Delta\to 0$ limit of Eqs. (5). We consider a quadratic expansion of the potential around an equilibrium position:

[TABLE]

The second term in the square bracket is the so-called pre-stress, which is proportional to the pressure and reduces the eigenvalues [13]. If we neglect the pre-stress, the quadratic expansion of the potential can be rewritten as

[TABLE]

where $k$ denotes the characteristic stiffness of the potential, $\alpha=(i,j)$ denotes the contact pair, $z$ denotes the contact number per particle, and $\delta r_{ij}=\hat{r}_{ij}\cdot(\delta\bm{x}_{i}-\delta\bm{x}_{j})$. We shall start from an isostatic configuration $z=2d$, where the system has a zero mode and thus $\lambda_{\rm min}=0$. To obtain a configuration above the jamming transition point, we add $N\delta z/2$ extra contacts. With this setting, one can express the minimal eigenvalue in a variational form:

[TABLE]

where $y_{\alpha}=\varepsilon$ denotes the rest length of the additional contacts. Combining an appropriate linear transformation and finite $N$ scaling, one can estimate the typical amplitude of $\lambda_{\rm min}$ for $\delta z\ll 1$ as [12]

[TABLE]

So far, we have neglected the effects of the pre-stress. For spherical particles, one can show that the pre-stress always gives a negative contribution to the eigenvalues [11]. As the pre-stress is proportional to the pressure, $\partial v(h_{ij}(0))/\partial h_{ij}(0)\sim-p$, we get

[TABLE]

where $c$ denotes a positive constant. At the jamming transition point, $\delta z=0$ and $p=0$, meaning that the system is marginal stable $\lambda_{\rm min}=0$. We assume that this marginal stability persists even above the jamming transition point; we require $\lambda_{\rm min}\sim 0$ or at most $\lambda_{\rm min}\sim p$. From this condition and Eq. (14), we have [11]

[TABLE]

This reproduces the results of numerical simulations in two and three dimensions [14].

2.4 Variational argument for non-spherical particles

Now we apply the variational argument to non-spherical particles. We first discuss the stability of a system consisting of non-spherical particles. At the jamming transition point $p=0$, $Q_{N}$ vanishes, and the number of constraints given by $V_{N}$ is

[TABLE]

where $N$ denotes the number of particles, and $z$ denotes the number of contacts per particle. For non-spherical particles, $N_{V}^{c}$ is smaller than the number of degrees of freedom, meaning that there are unconstrained modes, which we hereafter refer to as the zero modes. The number of the zero modes is

[TABLE]

where $\delta z=z-2d$, and $d_{\rm rot}$ denotes the number of rotational degrees of freedom per particle. For $p>0$, $Q_{N}$ has a finite value that would stabilize some of the zero modes. $Q_{N}$ gives $N_{Q}^{c}=Nd_{\rm rot}$ number of constraints, the typical stiffness of which is $k_{R}=\partial_{\bm{u}_{i}}\partial_{\bm{u}_{j}}Q_{N}\sim p\Delta$. This indeed suffices to stabilize the zero modes:

[TABLE]

When $\delta z\ll 1$, the system is considered to be nearly isostatic if one only takes into account the zero modes. In this case, we can apply the variational argument in Ref. [12] to the zero modes, as in the case of spherical particles near the jamming transition point. By repeating a similar argument used in the previous sub-section, we obtain

[TABLE]

where $c_{1}$ is a constant. Here we assume that $c_{1}>0$, which can be validated by the explicit calculations for ellipsoids interacting with the Gay-Berne potential, and breathing particles [10]. However, this assumption is not validated for some shapes of particles. We shall discuss this point in Sec. 6. If we take into account the second order term of $\Delta$, we have

[TABLE]

where all terms should be proportional to $p$, because the zero mode vanishes as $p$. Although the first order term of $\Delta$ can be positive, the pre-stress, in general, destabilizes the system due to the structural buckling [13]. Therefore, it is natural to assume that the second order term gives a negative contribution $c_{2}<0$. The marginal stability requires $\lambda_{\rm min}\sim 0$, meaning that the first and second terms in the RHS of Eq. (20) should be canceled each other, which leads to the scaling of $\delta z$:

[TABLE]

In Fig. 1, we compare the theoretical prediction, Eq. (21), with the numerical results of various shapes of non-spherical particles in Ref. [15] where the jammed configurations were generated by the binary search combining isotropic compression and decompression. We get an excellent agreement with the theory and numerical result, though there are visible deviations for $\mathcal{A}-1\ll 1$, which might be originated from the lack of statistics, numerical precision, or finite size effects.

In the $\Delta\to 0$ limit, Eq. (21) should be smoothly connected to the result of spherical particles Eq. (15). From this condition, one can decide the scaling form of $\delta z$ as

[TABLE]

where

[TABLE]

$\mathcal{Z}(x)$ is a finite and regular function at $x=0$, and thus one can expand it as $\mathcal{Z}(x)=\mathcal{Z}(0)+\mathcal{Z}^{\prime}(0)x+\cdots$, which leads to

[TABLE]

where $z_{J}=2d+\Delta^{1/2}\mathcal{Z}(0)$. This is again consistent with a numerical result of ellipsoids [16].

We now turn our attention to the scaling of the shear modulus $G$. In the standard numerical procedure to calculate $G$, one first imposes the small strain and then minimizes the energy. Comparing the resultant energy with that of the undeformed one, one can calculate $G$ [14]. Previous numerical and theoretical investigations prove that the square root singularity $G\sim p^{1/2}$ appears near the jamming transition of spherical particles [14, 11]. We here discuss how this scaling is altered for non-spherical particles. We assume that the imposed shear excites only the zero modes because the typical energy of those modes is much smaller than the other near the jamming transition point. When $\Delta\ll 1$, the zero modes mainly consists of the rotational degrees of freedom. The typical displacement $\delta u_{i}$ caused by the imposed shear strain $\delta\gamma$ is roughly $\delta u_{i}\sim\delta\gamma/\Delta$. After the minimization, the energy difference caused by the shear can be expressed as

[TABLE]

where $k_{i}$ and $\delta\tilde{\bm{u}}_{i}$ denotes the eigenvalue and eigenvector of $\partial_{\bm{u}_{i}^{\alpha}}\partial_{\bm{u}_{i}^{\beta}}g_{i}(\bm{u}_{i})$, respectively. $\bm{y}_{i}=\{y_{i}^{1},\cdots,y_{i}^{d_{\rm rot}}\}$ denotes the vector spanned by the $N_{0}$ zero modes to be chosen to minimize the energy. Using the standard technique of the linear algebra, one can eliminate $N_{0}$ terms among $Nd_{\rm rot}$ terms in Eq. (25) (see. Ch.7 in Ref. [17]). Thus, the typical amplitude of $\delta V_{N}$ is

[TABLE]

The shear modulus $G$ is then calculated as

[TABLE]

The result is consistent with the numerical results of ellipsoids [7]. In the $\Delta\to 0$ limit, Eq. (27) smoothly connects to the result of spherical particles [14, 11],

[TABLE]

This requires the following scaling form:

[TABLE]

where the scaling function $\mathcal{G}(x)$ satisfies

[TABLE]

In Fig. 2, we confirm our scaling prediction for ellipsoids interacting with harmonic potential, where $\Delta$ can be identified with the aspect ratio, and $p\propto\delta\varphi\equiv\varphi_{J}-\varphi$. One can see that the data of different aspect ratios are collapsed on a single curve, proving the validity of the scaling prediction Eq. (29).

3 Connection between non-spherical particles and breathing particles

The above argument can be generally applied for models having extra degrees of freedom, in addition to the translational degrees of freedom. Besides non-spherical particles, another interesting model that belongs to the same universality class is the so-called breathing particle (BP) model [5]. The model consists of $N$ polydisperse particles, the radii of which can change continuously. The interaction potential of the model is given by

[TABLE]

where

[TABLE]

and

[TABLE]

We chose the stiffness $k$ so that the standard deviation of the radii is proportional to $\Delta$:

[TABLE]

From the saddle point condition, $\partial_{R_{i}}V_{N}=0$, one can infer that

[TABLE]

By introducing the new variable $u_{i}\equiv(R_{i}-R_{i}^{0})/\Delta$, Eqs. (32) and (33) can be rewritten as

[TABLE]

where

[TABLE]

Note that $k_{R}$ now has the same order as that of non-spherical particles. Thus, one can repeat the same arguments in the previous sections for non-spherical particles, which leads to the same critical exponents [5]. The numerical implementation of the BP model is rather simpler, and the calculation time is shorter than those of non-spherical particles, as the extra degrees of freedom are simple scalar variables. For this reason, we shall use the BP in the numerical experiments in the following sections, instead of non-spherical particles. To obtain jammed configurations of the BP system, we use the FIRE algorithm to find the inherent structures of the potential Eq.(36). We use the Barendsen barostat to find them at a fixed pressure. All the details are explained in Ref. [5].

4 Universal scaling of the gap and force distributions near isostatic point

Here we show that the gap and force distributions exhibit the universal scaling behavior near the isostatic point.

4.1 Definition of the distribution functions

Here we investigate the gap distribution

[TABLE]

At the zero temperature $T=0$, $\rho(h)$ has a gap at $h=0$ [14]. For this reason, it is convenient to define distributions for the positive and negative $h$, separately. We define a positive gap distribution

[TABLE]

and a force (normalized negative gap) distribution

[TABLE]

where $\theta(x)$ is the Heviside function, and $f=-h/p$. It is well known that at the jamming transition point of spherical particles, $g(h)$ and $P(f)$ exhibit the power laws for small $h$ and $f$ [6, 18]:

[TABLE]

$\gamma$ does not depend on the spatial dimensions $d$ for $d\geq 2$ and follows the mean-field prediction $\gamma=0.41$ [19], while $\theta$ exhibits the weak $d$ dependence due to the localized excitations [20]. Below, we discuss that the power law is generally truncated at finite $h$ and $f$ if the system is not isostatic.

4.2 Finite size scaling

Here we first describe the distribution functions of spherical particles when the number of particles $N$ is finite. Then, in the next subsection, we show that the scaling of finite $N$ can be generalized to the scaling of non-isostatic systems, including non-spherical particles at the jamming transition point.

The minimal gap $h_{\rm min}$ of the $N$ particle system is calculated by using the extreme statistics [21]

[TABLE]

When $h\ll h_{\rm min}$, $g(h)$ quickly decreases, implying that the power law divergence of $g(h)$, Eq. (41), is truncated at $h\sim h_{\rm min}$. Thus, the scaling form of $g(h)$ at finite $N$ would be [21]

[TABLE]

where $p_{0}^{+}(x)$ is a regular and finite function, and

[TABLE]

We perform the numerical simulation for two-dimensional harmonic spheres to test the above conjecture. Instead of $g(h)$, we observe the cumulative distribution $Z(h)=\int_{0}^{h^{\prime}}dh^{\prime}g(h^{\prime})$ to improve the statistics. From Eq. (43), $Z(h)$ should satisfy

[TABLE]

In Fig. 3, we plot $Z(h)$ and its scaling form for several $N$. The excellent scaling collapse justifies our scaling argument.

4.3 General scaling form of the distribution functions near isostatic point

We want to generalize the above argument for more general systems close to the isostatic point $\delta z\ll 1$. For this purpose, we shall consider some function $F(h)$ which has the following scaling form for $\delta z\ll 1$:

[TABLE]

where $\alpha$ and $\beta$ denote the critical exponents we want to determine from the finite size scaling. The extensive numerical simulations of spherical particles prove that the scaling like Eq. (46) persists up to $\delta z=1/N$ [22, 23], suggesting that for a finite size system at the jamming transition point, we have

[TABLE]

In other words, the scaling for $\delta z\ll 1$ can be obtained by substituting $N=\delta z^{-1}$ into the result of the finite size scaling. From Eqs. (43) and (45), we have

[TABLE]

and

[TABLE]

We propose that the above equations hold for any system sufficiently near the isostatic point, i.e., $\delta z\ll 1$. We shall test this conjecture for the BP at the jamming transition point, at which we have shown that $\delta z\sim\Delta^{1/2}$ [5]. In Fig. 4, we show $Z(h)$ and its scaling form of the BP at the jamming transition point. The excellent scaling collapse justifies the validity of Eqs. (48) and (49). Note that the same equation of Eq. (48) holds exactly in the case of the mean-field model of non-spherical particles [5, 9] and the spherical particles slightly above the jamming transition point, where $\delta z\sim p^{1/2}$ [24].

For the force distribution $P(f)$, one can apply a similar argument, leading to

[TABLE]

where $p_{0}^{-}(x)$ denotes a finite and regular function, and we have introduced the critical exponent by

[TABLE]

The numerical justification of $P(f)$ is rather tricky because one should carefully separate the localized and extended modes to compare them with the theoretical prediction [20], which we leave for future work.

5 Vibrational density of states

The vibrational density of states $D(\omega)$ characterizes the low temperature physics of solids [25]. In this section, we investigate $D(\omega)$ of the BP model near the jamming transition point, which exhibits the same scaling as non-spherical particles. The results shown in the following are made for a system with $N=484$ particles.

5.1 Characteristic frequencies

We first describe the qualitative shape of $D(\omega)$ and define the characteristic frequencies.

In Fig. 5 (a), we show the typical behavior of $D(\omega)$ of the BP model. $D(\omega)$ consists of the three separate parts: (i) the lowest band at $\omega_{0}$, (ii) the intermediate band at $\omega_{1}$, and (iii) the highest band starts from $\omega_{2}$. In Fig 5 (b), we show the $\Delta$ dependence of the characteristic frequencies, $\omega_{0}$, $\omega_{1}$, and $\omega_{2}$. The characteristic frequencies are well fitted by the following power laws (see solid lines):

[TABLE]

In Fig. 5 (c), we show the $p$ dependence of $\omega_{0}$. We found that

[TABLE]

while $\omega_{1}$ and $\omega_{2}$ remain constant (not shown). The above scaling is the same as ellipsoids if we identify $\Delta$ with the aspect ratio [8], which is another evidence that the BP and ellipsoids belong to the same universality class.

Using the previous theoretical analysis [5, 9], we can understand the above scalings in the three regions. (i) The lowest band corresponds to the zero modes stabilized by the positive part of the pre-stress. As the pre-stress scales as $k_{R}\sim p\Delta$, Eq. (37), the characteristic frequency of the mode is $\omega_{0}\sim\sqrt{k_{R}}\sim p^{1/2}\Delta^{1/2}$. (ii) The intermediate band corresponds to the breathing motion of the BP or the rotation of ellipsoids. Therefore, the characteristic frequency is $\omega_{1}\sim\sqrt{\partial_{R}^{2}V_{N}}\sim\Delta$. (iii) The highest band corresponds to the translational degrees of freedom. As in the case of spherical particles, the characteristic frequency is proportional to $\delta z$ [17]. Using Eq. (21), we get $\omega_{2}\sim\Delta^{1/2}$. Herewith we recover the above numerical results. To give further evidence to support the above picture, we calculate the weights of each band by numerically integrating $D(\omega)$. If the above description is correct, one should have the following equations:

[TABLE]

where $N_{i}$ denotes the number of the modes included in the $i$-th band. Since the total number of modes is $3N$, the fraction $f_{i}$ of modes in each band is given by $f_{i}=N_{i}/3N$. In Fig. 6, we show the $\Delta$ and $p$ dependencies of $f_{i}$ and compare with the theoretical prediction Eqs. (54). We obtain quite good agreement.

6 Summary and discussions

In this paper, we investigated the jamming transition of non-spherical particles and breathing particles. Using both numerical and scaling arguments, we confirmed that the critical behavior of the jamming of non-spherical particles and breathing particles is qualitatively different from that of spherical particles. In the left panel of Fig. 7, we summarize our scaling prediction for the shear modulus $G$. Note that, for non-spherical particles ($\Delta>0$), $G$ always shows the linear pressure $p$ dependence sufficiently near the jamming transition point ($p\ll\Delta$), while it exhibits the square root dependence for spherical particles ($\Delta=0$). This means that the critical exponent of $G$ changes discontinuously at $\Delta=0$ from one to one half; in other words, the small asphericity is enough to change the universality class of the jamming transition. We also show that non-spherical particles and breathing particles are not critical at the jamming transition point in terms of the gap distribution $g(h)$, see the right panel of Fig. 7. The power law divergence of $g(h)$ is truncated at finite $h$, and thus the gap distribution is finite and analytic even at the jamming transition point. This is a sharp contrast to spherical particles, where the power law divergence of $g(h)$ persists up to $h=0$. Furthermore, we fully characterized the scaling of the characteristic frequencies of the density of states near the jamming transition point, which are again dramatically different from those of spherical particles.

There are still several open questions. A tentative list is the following:

• •

It is important to understand the rheological properties of the system near the jamming transition point. It has been shown that the divergence of the viscosity is strongly connected to the lowest excitation of the density of states [26]. As discussed in Sec. 5, the density of states of non-spherical particles is very different from that of spherical particles, which would change the critical exponent of the viscosity compared to that of spherical particles. A further study of this point is left as an open problem.

• •

In this work, we assumed that the pre-stress gives a positive contribution for the first order of $\Delta$, $c_{1}>0$. This assumption is violated for non-spherical particles consisting of spherical particles, such as dimers. In this case, the rotational motions of non-spherical particles can be identified with the translational motions of particles that consists non-spherical particles. Therefor, the system becomes isostatic and exhibits the same scaling of that of spherical particles [16, 27]. Furthermore, Platonic solids are also known to be isostatic at the jamming transition point [28]. Further studies are necessary to uncover when the assumption $c_{1}>0$ can be validated.

• •

In this work, we assumed that the two particles could have at most one contact. This assumption is correct for particles of convex shape. However, for particles of non-convex shape, the two particles can have more than two contacts. The extension of our work for such non-convex shape particles is an interesting future work. We believe that the study along this direction would be a promising way to introduce the effect of the friction, which is considered to be originated from the surface roughness of the constituent particles.

We thank F. Zamponi and P. Urbani for previous joint research [5] on which this work is based. 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).