Vibrational density of states of jammed packing at high dimensions: mean-field theory

Harukuni Ikeda [email protected] Department of Physics, Gakushuin University, 1-5-1 Mejiro, Toshima-ku, Tokyo 171-8588, Japan Masanari Shimada Department of Physics, Toronto Metropolitan University, M5B 2K3, Toronto, Canada

Abstract

Several mean-field theories predict that the Hessian matrix of amorphous solids converges the Wishart matrix in the limit of the large spatial dimensions $d\to\infty$ . Motivated by these results, we here calculate the density of states of random packing of harmonic spheres by mapping the Hessian of the original system to the Wishart matrix. We compare our result with that of previous numerical simulations of harmonic spheres in several spatial dimensions $d=3$ , $5$ , and $9$ . For small pressure $p\ll 1$ (near jamming), we find a good agreement even in $d=3$ , and obtain better agreements in larger $d$ , suggesting that the approximation becomes exact in the limit $d\to\infty$ .

^†^†preprint: AIP/123-QED

I Introduction

The vibrational density of states $D(\omega)$ plays a central role to characterize the low-temperature properties of solids. For both crystals and amorphous solids, $D(\omega)$ for small $\omega$ eventually follows the prediction of the Debye model $D(\omega)\sim\omega^{d-1}$ , suggesting that the vibrational excitation is dominated by phonon modes Kittel and McEuen (1976); Mizuno et al. (2017). However, for amorphous solids, in addition to the phonon modes, there arise excess non-phonon modes for small $\omega$ . This phenomenon, often referred to as the boson peak, is considered as a universal feature of amorphous solids Phillips and Anderson (1981).

From the theoretical point of view, a first step to tackle the problem is to consider mean-field models/theories. Several mean-field models, such as the $p$ -spin spherical model Biroli and Bouchaud (2012) and perceptron Franz et al. (2015), and theories, such as the effective medium theory DeGiuli et al. (2014); Shimada and De Giuli (2020), cavity method Parisi (2014), etc. Beltukov (2015); Cicuta et al. (2018); Baggioli et al. (2019); Baggioli and Zaccone (2019), suggest that Hessian matrices of amorphous solids are approximated by the Wishart matrix. However, somewhat surprisingly, the functional form of $D(\omega)$ of particle systems has not been calculated yet, even in the large dimensional limit $d\to\infty$ , where the mean-field theory becomes exact. As a consequence, one should introduce fitting parameters to compare the theory and numerical results Manning and Liu (2015); Beltukov (2015); Baggioli et al. (2019), even in large $d$ Charbonneau et al. (2016).

In this work, we focus on frictionless spherical particles interacting with the harmonic potential O’Hern et al. (2003). Since the harmonic potential is a purely repulsive potential, the system gets unstable in the zero pressure limit $p\to 0$ , which is known as the (un)jamming transition O’Hern et al. (2003). Near the jamming transition point ( $p\ll 1$ ), several physical quantities, such as the contact number $z$ , exhibit the power-law behavior O’Hern et al. (2003). The critical exponents near the jamming transition are calculated by several mean-field theories Wyart et al. (2005); DeGiuli et al. (2014); Franz et al. (2017); Parisi et al. (2020). However, again, the detailed functional form of $z$ is still undetermined, even in $d\to\infty$ .

Recently, one of the present authors performed an extensive numerical simulation of harmonic spheres and calculated $D(\omega)$ and $z$ in spatial dimensions from $d=3$ to $d=9$ Shimada et al. (2020). Therefore, it is now desirable to directly compare the numerical results in large $d$ with the predictions of the mean-field theory.

Here, we theoretically calculate $D(\omega)$ and $z$ of harmonic spheres in large $d$ , and compare them with the previous numerical results. For this purpose, inspired by the previous mean-field calculations, we assume that the Hessian of harmonic spheres converges to the (shifted) Wishart matrix in the mean-field limit $d\to\infty$ . We determine the pre-factors of the Wishart matrix so that its trace is consistent with that of the Hessian of the original model. For small pressure, our results well agree with the previous numerical results Shimada et al. (2020) even in $d=3$ , and obtain better agreements in larger $d$ , suggesting that our theory becomes exact in the limit of $d\to\infty$ .

The organization of the paper is as follows. In Sec. II, we introduce the model and several physical quantities. In Sec. III, we calculate $D(\omega)$ in the limit of $d\to\infty$ . In Sec. IV, we summarize the results.

II Settings

Here we introduce the model and several physical quantities. We consider a system consisting of frictionless spherical particles interacting with the harmonic potential O’Hern et al. (2003):

\displaystyle V=\sum_{i<j}^{1,N}k\frac{h_{ij}^{2}}{2}\theta(-h_{ij}),\ h_{ij}=\left|\bm{r}_{i}-\bm{r}_{j}\right|-R_{i}-R_{j}

(1)

where $N$ denotes the number of particles, $k$ denotes the spring constant, and $\theta(x)$ denotes the Heaviside step function. $\bm{r}_{i}=\{x_{i1},\dots,x_{id}\}$ and $R_{i}$ denote the position and radius of the $i$ -th particle, respectively. To simplify the notation, hereafter, we set $k=1$ .

The Hessian of the potential is

	$\displaystyle\mathcal{H}_{ia,jb}=\frac{\partial^{2}V}{\partial x_{ia}\partial x_{jb}}={\mathcal{H}}_{ia,jb}^{(1)}+{\mathcal{H}}_{ia,jb}^{(2)},$
	$\displaystyle{\mathcal{H}}_{ia,jb}^{(1)}=\sum_{\mu=1}^{Nz/2}\frac{\partial h_{\mu}}{\partial x_{ia}}\frac{\partial h_{\mu}}{\partial x_{jb}},\ {\mathcal{H}}_{ia,jb}^{(2)}=\sum_{\mu=1}^{Nz/2}h_{\mu}\frac{\partial^{2}h_{\mu}}{\partial x_{ia}\partial x_{jb}},$		(2)

where

\displaystyle z=\frac{1}{N}\sum_{i<j}\theta(-h_{ij}).

(3)

denotes the number of contacts per particle, and $\sum_{\mu=1}^{Nz/2}$ denotes the sum of all pairs $\mu=(ij)$ for which $h_{ij}<0$ . Once we have the eigenvalue distribution of ${\mathcal{H}}$ , $\rho(\lambda)$ , the vibrational density of states $D(\omega)$ is calculated as

\displaystyle D(\omega)=2\omega\rho(\lambda=\omega^{2}).

(4)

For the control parameter, we use the pre-stress defined as Shimada et al. (2020)

\displaystyle e=-\frac{2(d-1)}{Nz}\sum_{\mu=1}^{Nz/2}\frac{h_{\mu}}{r_{\mu}}=(d-1)\left\langle\frac{R_{i}+R_{j}}{r_{ij}}-1\right\rangle,

(5)

where $\left\langle\bullet\right\rangle$ denotes the average for the all contacts $\left\langle\bullet\right\rangle=(Nz/2)^{-1}\sum_{i<j}^{Nz/2}\theta(-h_{ij})\bullet$ . The right most expression in Eq. (5) clearly shows that $e$ proportional to the average overlap of particles. The proportional constant $(d-1)$ has been chosen so that $e$ remains finite in the limit $d\to\infty$ Shimada et al. (2020). Near the jamming transition point, $e$ is proportional to the pressure, $e\sim p$ and vanishes at the jamming transition point. In a previous numerical study Charbonneau et al. (2016), the packing fraction was used as a control parameter. However, it has been pointed out that $e$ is a more natural control parameter Bi et al. (2015); Shimada et al. (2020). Below, we calculate $z$ and $D(\omega)$ as functions of $e$ .

III Theory

III.1 Summary of previous works

Here we briefly review the previous works. The seminal work has been done by G. Parisi Parisi (2014). He showed that the eigenvalue distribution of the Hessian of harmonic spheres converges to the Marcenko-Pastur distribution in the limit $d\to\infty$ , meaning that the Hessian is identified with the Wishart matrix ${\mathcal{H}}\sim{\mathcal{W}}$ Parisi (2014), but the effects of the pre-stress $e$ were neglected at that time. More recently, G. Parisi with coauthors performed a more complete analysis for the perceptron, which is a mean-field model of random sphere packing of harmonic spheres Franz and Parisi (2016); Franz et al. (2017). The analysis of the perceptron suggests that the Hessian of harmonic spheres in large $d$ is written as

\displaystyle{\mathcal{H}}_{\rm MF}=a{\mathcal{W}}+be{\mathcal{I}},

(6)

where $a$ and $b$ denote constants, ${\mathcal{I}}_{ia,jb}=\delta_{ia,jb}$ denotes the identity matrix, and

\displaystyle{\mathcal{W}}_{ia,jb}=\frac{2}{Nz}\sum_{\mu=1}^{Nz/2}\xi_{ia}^{\mu}\xi_{jb}^{\mu}

(7)

denotes the Wishart matrix. $\xi_{ia}^{\mu}$ denotes the i.i.d gaussisan random variable with zero mean and unit variance Franz et al. (2015). The replica calculation of the perceptron also proves the marginal stability Müller and Wyart (2015): the minimal eigenvalue $\lambda_{\rm min}$ of the Hessian ${\mathcal{H}}$ vanishes $\lambda_{\rm min}=0$ near the jamming transition point Franz et al. (2015). Interestingly, several other mean-field theories also suggest that the Hessian of harmonic spheres is written as Eq. (6), though the precise values of $a$ and $b$ are still unknown DeGiuli et al. (2014); Beltukov (2015); Cicuta et al. (2018); Baggioli et al. (2019); Ikeda (2020).

In this work, we use Eq. (6) as an Ansatz. We determine $a$ and $b$ by requiring that the trace of ${\mathcal{H}}_{\rm MF}$ is consistent with that of the Hessian of the original model, Eq. (2).

III.2 Calculation of $a$ and $b$

For simplicity, we first discuss the case without pre-stress. By setting $e=0$ , we get

\displaystyle{\mathcal{H}}\to{\mathcal{H}}^{(1)},\ {\mathcal{H}}_{\rm MF}\to a{\mathcal{W}}.

(8)

Then, we determine $a$ from the following condition:

\displaystyle{\rm Tr}{\mathcal{H}}^{(1)}=a{\rm Tr}{\mathcal{W}}.

(9)

The LHS in Eq. (9) can be calculated as

\displaystyle{\rm Tr}{\mathcal{H}}^{(1)}=\sum_{i=1}^{N}\sum_{a=1}^{d}\sum_{\mu=1}^{Nz/2}\left(\frac{\partial h_{\mu}}{\partial x_{ia}}\right)^{2}=Nz,

(10)

where we used

\displaystyle\frac{\partial h_{ij}}{\partial x_{ka}}=\left(\delta_{ik}-\delta_{jk}\right)\frac{x_{ia}-x_{ja}}{\left|\bm{r}_{i}-\bm{r}_{j}\right|},

\displaystyle\sum_{k=1}^{N}\sum_{a=1}^{d}\left(\frac{\partial h_{ij}}{\partial x_{ka}}\right)^{2}=2.

(11)

The RHS is

\displaystyle a{\rm Tr}{\mathcal{W}}=\frac{2a}{Nz}\sum_{i=1}^{N}\sum_{a=1}^{d}\sum_{\mu=1}^{Nz/2}\left(\xi_{ia}^{\mu}\right)^{2}=aNd.

(12)

By using Eqs. (9), (10) and (12), we get

\displaystyle a=\frac{z}{d}.

(13)

Next we consider the full matrix including the term proportional to $e$ . Assuming that ${\rm Tr}{\mathcal{H}}={\rm Tr}{\mathcal{H}}_{\rm MF}$ , we get

	$\displaystyle{\rm Tr}\left({\mathcal{H}}^{(1)}+{\mathcal{H}}^{(2)}\right)={\rm Tr}\left(a{\mathcal{W}}+be{\mathcal{I}}\right)$
	$\displaystyle\to{\rm Tr}{\mathcal{H}}^{(2)}=be{\rm Tr}{\mathcal{I}}.$		(14)

The LHS can be calculated as

$\displaystyle{\rm Tr}{\mathcal{H}}^{(2)}$	$\displaystyle=\sum_{i=1}^{N}\sum_{a=1}^{d}\sum_{\mu=1}^{Nz/2}h_{\mu}\frac{\partial^{2}h_{\mu}}{\partial x_{ia}^{2}}$
	$\displaystyle=2(d-1)\sum_{\mu=1}^{Nz/2}\frac{h_{\mu}}{r_{\mu}}$
	$\displaystyle=-Nze,$	(15)

where we used

	$\displaystyle\frac{\partial^{2}h_{ij}}{\partial x_{ka}^{2}}=\left(\delta_{ik}+\delta_{jk}\right)\frac{r_{ij}^{2}-(x_{ia}-x_{jb})^{2}}{r_{ij}^{3}},$
	$\displaystyle\sum_{k=1}^{N}\sum_{a=1}^{d}\frac{\partial^{2}h_{ij}}{\partial x_{ka}^{2}}=2\frac{d-1}{r_{ij}}.$		(16)

From Eqs. (14), (15) and ${\rm Tr}{\mathcal{I}}=Nd$ , we get

\displaystyle b=-\frac{z}{d}.

(17)

In summary, we get

\displaystyle{\mathcal{H}}_{\rm MF}=\frac{z}{d}{\mathcal{W}}-\frac{z}{d}e{\mathcal{I}}.

(18)

III.3 Eigenvalue distribution

It is well-known that the eigenvalue distribution of the Wishart matrix ${\mathcal{W}}$ follows the Marchencko-Pastur law Livan et al. (2018)

\displaystyle\rho_{\rm MP}(\lambda)=\frac{z}{2d}\frac{\sqrt{(\lambda_{+}-\lambda)(\lambda-\lambda_{-})}}{2\pi\lambda},\ \lambda_{\pm}=\left(1\pm\sqrt{\frac{2d}{z}}\right)^{2}.

(19)

Let $\bm{e}_{n}$ be an eigenvector of ${\mathcal{W}}$ , and $\lambda_{n}^{\rm MP}$ be the corresponding eigenvalue. Then, we have

\displaystyle{\mathcal{H}}_{\rm MF}\cdot\bm{e}_{n}=\left(\frac{z}{d}\lambda_{n}^{\rm MP}-\frac{z}{d}e\right)\bm{e}_{n},

(20)

meaning that $\bm{e}_{n}$ is also an eigenvector of ${\mathcal{H}}_{\rm MF}$ and the corresponding eigenvalue is $\lambda_{n}=\frac{z}{d}\lambda_{n}^{\rm MP}-\frac{z}{d}e$ . Therefore, the eigenvalue distribution of ${\mathcal{H}}_{\rm MF}$ is calculated as

\displaystyle\rho(\lambda)=\rho_{\rm MP}(\lambda_{\rm MP})\frac{d\lambda_{\rm MP}}{d\lambda}=\frac{d}{z}\rho_{\rm MP}(d\lambda/z+e).

(21)

In particular, the minimal eigenvalue is

\displaystyle\lambda_{\rm min}=\frac{z}{d}\left(1-\sqrt{\frac{2d}{z}}\right)^{2}-\frac{z}{d}e.

(22)

III.4 Marginal stability and contact number

The replica calculation in the limit $d\to\infty$ predicts that the system is marginally stable near the jamming transition point, $\lambda_{\rm min}=0$ ¹¹1 Here we take the limit $d\to\infty$ so that $z/d$ remains finite. This is consistent with the numerical result the jamming transition point where $z\approx 2d$ O’Hern et al. (2003). Franz et al. (2015); Kurchan et al. (2013); Parisi et al. (2020). By using the marginal stability and Eq. (22), we get

\displaystyle\frac{z(e)}{2d}=\frac{1}{(1-e^{1/2})^{2}}.

(23)

Here we assumed that $z/d$ takes a finite value in the limit $d\to\infty$ , because the numerical results suggest that $z\to 2d$ at the jamming transition point O’Hern et al. (2003). For $e\ll 1$ , we reproduce the well-known scaling observed by numerical simulations O’Hern et al. (2003)

\displaystyle z/2d-1\sim 2e^{1/2}\sim p^{1/2}.

(24)

The critical exponent $1/2$ was previously derived by using the variational argument Wyart et al. (2005), effective medium theory DeGiuli et al. (2014), and replica theory Franz et al. (2017), but our result Eq. (23) also allows us to access the pre-factor and non-linear terms. Somewhat surprisingly, Eq. (23) suggests that $z(e)$ depends only on $e$ , and does not depend on the preparation protocols. It is interesting future work to see if this property survives in finite $d$ .

In Fig. 1 (a), we compare Eq. (23) with numerical results in several spatial dimensions $d$ obtained by rapid quench from high temperature random configurations. See Ref. Shimada et al. (2020) for the details of the numerical simulations. The theory well agrees with the numerical results for small $e$ . For more quantitative discussion, in Fig. 1 (b), we show the difference between the results of the theory $z_{\rm the}$ and simulation $z_{\rm sim}$ :

\displaystyle\epsilon=\frac{z_{\rm the}/2d-1-\left(z_{\rm sim}/2d-1\right)}{z_{\rm the}/2d-1}=\frac{z_{\rm the}-z_{\rm sim}}{z_{\rm the}-2d}.

(25)

The data collapse onto a single curve if we rescale the vertical axis by $d$ (Fig. 1 (c)), meaning that the deviation scales as $\epsilon\sim 1/d$ .

Refer to caption — Figure 1: (a) $e$ dependence of $z$ . Markers denote numerical results taken from Ref. Shimada *et al.* (2020), while the solid line denotes the theoretical prediction. (b) $\epsilon$ for the same data. (c) $d\epsilon$ for the same data.

III.5 Vibrational density of states

By using Eq. (21), the vibrational density of states $D(\omega)$ is calculated as ${D(\omega)=2\omega\rho(\lambda=\omega^{2})}$ . Although $D(\omega)$ depends on both $z$ and $e$ , Eq. (23) allows us to eliminate the dependency on $z$ . After some manipulations, we get

\displaystyle D(\omega)=\frac{\omega^{2}\sqrt{(1-e^{1/2})^{3}\left\{8-(1-e^{1/2})\omega^{2}\right\}}}{2\pi\left\{2e+(1-e^{1/2})^{2}\omega^{2}\right\}}.

(26)

In Fig. 2, we compare the theoretical prediction Eq. (26) and numerical results. The results are consistent near jamming $e=0.01$ even in $d=3$ , while there is a visible deviation for small $\omega$ far from jamming $e=0.2$ even in $d=9$ . It is an interesting future work to see if a better agreement is obtained in higher $d$ .

For $e\ll 1$ and $\omega\ll 1$ , we get the following scaling:

\displaystyle D(\omega)\sim\frac{\omega^{2}\sqrt{\omega_{\rm max}^{2}-\omega^{2}}}{2\pi(\omega^{2}+\omega_{*}^{2})}\sim\begin{cases}{\rm const}&\omega\gg\omega_{*}\\ \delta z^{-2}\omega^{2}&\omega\ll\omega_{*},\end{cases}

(27)

where $\omega_{\rm max}=2\sqrt{2}$ , and $\omega_{*}=\sqrt{2e}\propto z/2d-1$ . In particular, $D(\omega)={\rm const}$ at the jamming transition point $e=0$ . The similar results have been previously derived by applying the effective medium theory to the disordered lattices DeGiuli et al. (2014), and the replica method to the mean-field models Franz et al. (2015); Ikeda (2020).

IV Summary

In this work, we calculated the contact number $z$ and vibrational density of states $D(\omega)$ for harmonic spheres in the large spatial dimensions $d\to\infty$ . Our theoretical results well agree with the results of the previous numerical simulation in large $d$ .

Our theoretical results are relied on the Ansatz Eq. (6), that is, the Hessian of harmonic spheres has the form of the shifted Wishart matrix. The consistency between our theoretical results and previous numerical results suggests that the Ansatz becomes exact in the limit $d\to\infty$ . This result motivates us to develop more rigorous calculation in $d\to\infty$ without using the Ansatz, as done in the previous exact calculations for hard spheres Kurchan et al. (2012, 2013); Charbonneau et al. (2014). We left it as future work.

Acknowledgements.

We thank F. Zamponi for useful comments. This project has received JSPS KAKENHI Grant Numbers 21K20355 and 19J20036.

References

Kittel and McEuen (1976) C. Kittel and P. McEuen, Introduction to solid state physics, Vol. 8 (Wiley New York, 1976).
Mizuno et al. (2017) H. Mizuno, H. Shiba, and A. Ikeda, Proceedings of the National Academy of Sciences 114, E9767 (2017).
Phillips and Anderson (1981) W. A. Phillips and A. Anderson, Amorphous solids: low-temperature properties, Vol. 24 (Springer, 1981).
Biroli and Bouchaud (2012) G. Biroli and J.-P. Bouchaud, Structural Glasses and Supercooled Liquids: Theory, Experiment, and Applications , 31 (2012).
Franz et al. (2015) S. Franz, G. Parisi, P. Urbani, and F. Zamponi, Proceedings of the National Academy of Sciences 112, 14539 (2015).
DeGiuli et al. (2014) E. DeGiuli, A. Laversanne-Finot, G. Düring, E. Lerner, and M. Wyart, Soft Matter 10, 5628 (2014).
Shimada and De Giuli (2020) M. Shimada and E. De Giuli, arXiv preprint arXiv:2008.11896 (2020).
Parisi (2014) G. Parisi, arXiv preprint arXiv:1401.4413 (2014).
Beltukov (2015) Y. Beltukov, JETP Letters 101, 345 (2015).
Cicuta et al. (2018) G. M. Cicuta, J. Krausser, R. Milkus, and A. Zaccone, Phys. Rev. E 97, 032113 (2018).
Baggioli et al. (2019) M. Baggioli, R. Milkus, and A. Zaccone, Phys. Rev. E 100, 062131 (2019).
Baggioli and Zaccone (2019) M. Baggioli and A. Zaccone, Phys. Rev. Research 1, 012010 (2019).
Manning and Liu (2015) M. L. Manning and A. J. Liu, EPL (Europhysics Letters) 109, 36002 (2015).
Charbonneau et al. (2016) P. Charbonneau, E. I. Corwin, G. Parisi, A. Poncet, and F. Zamponi, Phys. Rev. Lett. 117, 045503 (2016).
O’Hern et al. (2003) C. S. O’Hern, L. E. Silbert, A. J. Liu, and S. R. Nagel, Phys. Rev. E 68, 011306 (2003).
Wyart et al. (2005) M. Wyart, L. E. Silbert, S. R. Nagel, and T. A. Witten, Phys. Rev. E 72, 051306 (2005).
Franz et al. (2017) S. Franz, G. Parisi, M. Sevelev, P. Urbani, and F. Zamponi, (2017).
Parisi et al. (2020) G. Parisi, P. Urbani, and F. Zamponi, Theory of simple glasses: exact solutions in infinite dimensions (Cambridge University Press, 2020).
Shimada et al. (2020) M. Shimada, H. Mizuno, L. Berthier, and A. Ikeda, Physical Review E 101, 052906 (2020).
Bi et al. (2015) D. Bi, S. Henkes, K. E. Daniels, and B. Chakraborty, Annu. Rev. Condens. Matter Phys. 6, 63 (2015).
Franz and Parisi (2016) S. Franz and G. Parisi, Journal of Physics A: Mathematical and Theoretical 49, 145001 (2016).
Müller and Wyart (2015) M. Müller and M. Wyart, (2015).
Ikeda (2020) H. Ikeda, Phys. Rev. Research 2, 033220 (2020).
Livan et al. (2018) G. Livan, M. Novaes, and P. Vivo, Introduction to random matrices: theory and practice, Vol. 26 (Springer, 2018).
Note (1) Here we take the limit $d\to\infty$ so that $z/d$ remains finite. This is consistent with the numerical result the jamming transition point where $z\approx 2d$ O’Hern et al. (2003).
Kurchan et al. (2013) J. Kurchan, G. Parisi, P. Urbani, and F. Zamponi, The Journal of Physical Chemistry B 117, 12979 (2013).
Kurchan et al. (2012) J. Kurchan, G. Parisi, and F. Zamponi, Journal of Statistical Mechanics: Theory and Experiment 2012, P10012 (2012).
Charbonneau et al. (2014) P. Charbonneau, J. Kurchan, G. Parisi, P. Urbani, and F. Zamponi, Journal of Statistical Mechanics: Theory and Experiment 2014, P10009 (2014).