Vibrational density of states of jammed packing: mean-field theory and beyond
VVibrational density of states of jammed packing: mean-field theory and beyond
Harukuni Ikeda ∗ and Masanari Shimada Graduate School of Arts and Sciences, The University of Tokyo 153-8902, Japan (Dated: November 13, 2020)Several mean-field theories predict that Hessian matrices of amorphous solids can be written byusing the random matrix in the limit of the large spatial dimensions d → ∞ . Motivated by theseresults, we here propose a way to map a Hessian of the amorphous solid to a random matrix. Thisis possible by determining the coefficients of a random matrix so that the trace of the randommatrix coincides with the Hessian of the original system. We compare our result with that ofprevious numerical simulations of harmonic spheres in several spatial dimensions d = 3, 5, and 9.For small pressure p (cid:28) d = 3, and obtainbetter agreements in larger d , suggesting that the approximation indeed becomes exact in the limitof large spatial dimensions. We also discuss the finite dimensional effects, which are not consideredin the mean-field theory and lead to the coexistence of the Debye density of state D ( ω ) ∼ ω d − andthe quartic mode D ( ω ) ∼ ω for small ω . I. INTRODUCTION
The vibrational density of states D ( ω ) plays a cen-tral role to characterize the low-temperature propertiesof solids. For both crystals and amorphous solids, D ( ω )eventually follows the prediction of the Debye model D ( ω ) ∼ ω d − in the low frequency limit ω →
0, sug-gesting that the vibrational excitation is dominated byphonon modes [1, 2]. However, for amorphous solids, inaddition to the phonon modes, there arise excess non-phonon modes for small ω . This phenomenon, often re-ferred to as the boson peak, is considered as a universalfeature of amorphous solids [3].From the theoretical point of view, a first step to tacklethe problem is to consider mean-field models/theories.Several mean-field models, such as the p -spin sphericalmodel [4] and perceptron [5], and theories, such as theeffective medium theory [6, 7], cavity method [8], etc. [9–11], suggest that Hessian matrices of amorphous solidsare approximated by random matrices. However, some-what surprisingly, the functional form of D ( ω ) of particlesystems has not been calculated yet, even in the largedimensional limit d → ∞ , where the mean-field theorybecomes exact. As a consequence, one should introducefitting parameters to compare the theory and numericalresults [9, 11, 12], even in large d [13].In this work, we focus on frictionless spherical particlesinteracting with the harmonic potential [14]. Since theharmonic potential is a purely repulsive potential, thesystem gets unstable in the zero pressure limit p → p (cid:28) z , exhibit thepower-law behavior [14]. The critical exponents near thejamming transition are calculated by several mean-fieldtheories [6, 15–17]. However, again, the detailed func-tional form of z is still undetermined, even in d → ∞ . ∗ [email protected] Recently, one of the present authors performed an ex-tensive numerical simulation of harmonic spheres and cal-culated D ( ω ) and z in spatial dimensions from d = 3 to d = 9 [18]. Therefore, it is now desirable to directly com-pare the numerical results in large d with the predictionsof the mean-field theory.Here, we theoretically calculate D ( ω ) and z of har-monic spheres in large d , and compare them with theprevious numerical results. For this purpose, inspired bythe previous mean-field calculations, we assume that theHessian of harmonic spheres converges to a random ma-trix in the mean-field limit d → ∞ , and then determinethe pre-factors by requiring that the trace of the randommatrix coincides with the Hessian of the original model.For small pressure, our results well agree with the pre-vious numerical results [18] even in d = 3, and obtainbetter agreements in larger d , suggesting that our theoryindeed becomes exact in the limit of d → ∞ .Although the main purpose of this work is to calcu-late the physical quantities in d → ∞ , we aware thatthe the finite d effects alter the mean-field predictionsin finite d . For instance, the mean-field theory near thejamming transition point predicts D ( ω ) ∼ ω in the lowfrequency limit ω → d = 3 reported the coexistence of the phonon mode D ( ω ) ∼ ω d − and the quartic mode D ( ω ) ∼ ω in thatlimit [2]. We briefly discuss that the discrepancy betweenthe mean-field theory in d → ∞ and numerical results infinite d can be reconciled by using the fluctuated effec-tive medium theory, which is an extended version of theeffective medium theory [6] to take into account the fluc-tuation of the pre-stress [19].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 ( ω ) in the limit of d → ∞ .In Sec. IV, we discuss how the finite d effects alter themean-field result in d → ∞ . In Sec. V, we summarizethe results. a r X i v : . [ c ond - m a t . d i s - nn ] N ov II. SETTINGS
Here we introduce the model and several physicalquantities. We consider a system consisting of frictionlessspherical particles interacting with the harmonic poten-tial [14]: V = ,N (cid:88) i Nz/ (cid:88) µ =1 h µ r µ . (5) e is proportional to the pressure, e ∼ p and vanishesat the jamming transition point. In a previous numericalsimulation [13], the packing fraction was used as a controlparameter. However, it has been recently realized that e is a more natural control parameter [18]. Hereafter, wecalculate z and D ( ω ) for a given e . III. MEAN-FIELD THEORYA. Summary of previous works Here we briefly review the previous works. The sem-inal work has been done by G. Parisi [8]. He showed that the eigenvalue distribution of the Hessian of har-monic spheres converges to the Marcenko-Pastur distri-bution in the limit d → ∞ , meaning that the Hessianis identified with the Wishart matrix H ∼ W [8], butthe effects of the pre-stress e were neglected at thattime. More recently, G. Parisi with coauthors performeda more complete analysis for the perceptron, which is amean-field model of random sphere packing of harmonicspheres [16, 20]. The analysis of the perceptron suggeststhat the Hessian of harmonic spheres in large d is writtenas H MF = a W + be I , (6)where a and b denote constants, I ia,jb = δ ia,jb denotesthe identity matrix, and W ia,jb = 2 N z Nz/ (cid:88) µ =1 ξ µia ξ µjb (7)denotes the Wishart matrix. ξ µia denotes the i.i.d gaus-sisan random variable with zero mean and unit vari-ance [5]. The replica calculation of the perceptron alsoproves the marginal stability [21]: the minimal eigen-value λ min of the Hessian H vanishes λ min = 0 near thejamming transition point [5]. Interestingly, several othermean-field theories also suggest that the Hessian of har-monic spheres is written as Eq. (6), though the precisevalues of a and b are still unknown [6, 9–11, 22].In this section, we use Eq. (6) as an Ansatz that maybe exact in the limit d → ∞ . Below, we determine a and b by requiring that the trace of H MF is consistentwith that of the Hessian of the original model, Eq. (2).In Appendix. A, we provide another complemental andmore physically motivated derivation. B. Calculation of a and b For simplicity, we first discuss the case without pre-stress. By artificially setting e = 0, we get H → H (1) , and , H MF → a W . (8)Then, by requiring Tr H = Tr H MF , we get a = Tr H (1) Tr W = zd , (9)where we used ∂h ij ∂x ka = ( δ ik − δ jk ) x ia − x ja | r i − r j | , and (cid:80) ka (cid:16) ∂h ij ∂x ka (cid:17) = 2.Now we consider full matrix including the term pro-portional to e . Assuming that Tr H = Tr H MF , we getTr (cid:16) H (1) + H (2) (cid:17) = Tr ( a W + be I ) → b = Tr (cid:0) H (2) (cid:1) Tr ( e I ) = − zd , (10)where we used ∂ h ij ∂x ka = ( δ ik + δ jk ) r ij − ( x ia − x jb ) r ij , and (cid:80) ,dNka ∂ h ij ∂x ka = 2 d − r ij . C. Eigenvalue distribution In summary, we get H MF = zd W − zd e I . (11)It is well-known that the eigenvalue distribution ofthe Wishart matrix W follows the Marchencko-Pasturlaw [23] ρ MP ( λ ) = z d (cid:112) ( λ + − λ )( λ − λ − )2 πλ , λ ± = (cid:32) ± (cid:114) dz (cid:33) . (12)Let e n be an eigenvector of W , and λ MP n be the corre-sponding eigenvalue. Then, we have H · e n = (cid:16) zd λ MP n − zd e (cid:17) e n , (13)meaning that e n is also an eigenvector of H and the cor-responding eigenvalue is λ n = zd λ MP n − zd e . Therefore, theeigenvalue distribution of H is calculated as ρ ( λ ) = ρ MP ( λ MP ) dλ MP dλ = dz ρ MP ( dλ/z + e ) . (14)In particular, the minimal eigenvalue is λ min = zd (cid:32) − (cid:114) dz (cid:33) − zd e. (15) D. Marginal stability and contact number The replica calculation in the limit d → ∞ predictsthat the system is marginally stable near the jammingtransition point, λ min = 0 [5, 17, 24]. By using themarginal stability and Eq. (15), we get z ( e ) = 2 d (1 − e / ) . (16)For e (cid:28) 1, we reproduce the well-known scaling observedby numerical simulations [14] z − d ∼ de / ∼ p / . (17)The critical exponent 1 / z ( e ) depends only on e , and does not depend on the ●●●●●●●●●● ●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●■■■■■■■■■■ ■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■◆◆◆◆◆◆◆◆◆◆ ◆ ◆ ◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆● d = ■ d = ◆ d = - - - - - e z / - ( a ) ●●●●●●●●●● ●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●■■■■■■■■■■ ■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■◆◆◆◆◆◆◆◆◆◆ ◆ ◆ ◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆● d = ■ d = ◆ d = - - - - - e ϵ ( b ) ●●●●●●●●●● ●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●■■■■■■■■■■ ■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■◆◆◆◆◆◆◆◆◆◆ ◆ ◆ ◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆● d = ■ d = ◆ d = - - - - e d ϵ ( c ) FIG. 1. (a) e dependence of z . Markers denote numericalresults taken from Ref. [18], while the solid line denotes thetheoretical prediction. (b) (cid:15) for the same data. (c) d(cid:15) for thesame data. ●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆● d = ■ d = ◆ d = - - - - - ω D ( ω ) ( a ) e = ●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆● d = ■ d = ◆ d = - - - - - ω D ( ω ) ( b ) e = ●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆● d = ■ d = ◆ d = - - - - - ω D ( ω ) ( c ) e = ●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆◆● d = ■ d = ◆ d = - - - - - ω D ( ω ) ( d ) e = FIG. 2. Density of states D ( ω ). Markers denote numericalresults taken from Ref. [18]. The solid lines denote theoreticalpredictions. preparation protocols. It is interesting future work tosee if this property survives in finite d .In Fig. 1 (a), we compare Eq. (16) with numerical re-sults in several spatial dimensions d obtained by rapidquench from high temperature random configurations.See Ref. [18] for the details of the numerical simulations.The theory well agrees with the numerical results forsmall e . For more quantitative discussion, in Fig. 1 (b),we show the difference between the results of the theory z the and simulation z sim : (cid:15) = z the / d − − ( z sim / d − z the / d − z the − z sim z the − d . (18)The data collapse onto a single curve if we rescale thevertical axis by d (Fig. 1 (c)), meaning that (cid:15) ∼ /d . E. Vibrational density of states By using Eq. (14), the vibrational density of states D ( ω ) is calculated as D ( ω ) = 2 ωρ ( λ = ω ). Although D ( ω ) depends on both z and e , Eq. (16) allows us to elim-inate the dependency on z . After some manipulations,we get D ( ω ) = ω (cid:113) (1 − e / ) (cid:8) − (1 − e / ) ω (cid:9) π (cid:8) e + (1 − e / ) ω (cid:9) . (19)In Fig. 2, we compare the theoretical prediction Eq. (19)and numerical results. Results are consistent near jam-ming e = 0 . 01 even in d = 3, while there is a visibledeviation for small ω far from jamming e = 0 . d = 9. It is an interesting future work to see if a betteragreement is obtained in higher d .For e (cid:28) ω (cid:28) 1, we get the following scaling: D ( ω ) ∼ ω (cid:112) ω − ω π ( ω + ω ∗ ) ∼ (cid:40) const ω (cid:29) ω ∗ δz − ω ω (cid:28) ω ∗ , (20)where ω max = 2 √ 2, and ω ∗ = √ e ∝ z/ d − 1. Thesimilar results have been previously derived by applyingthe effective medium theory to the disordered lattices [6],and the replica method to the mean-field models [5, 22]. IV. BEYOND MEAN-FIELD THEORY Here we discuss how the finite dimensional effects alterthe mean-field result. In this section, we only focus on thescaling behavior for small ω , and neglect the proportionalconstants. A. Effective medium theory The mean-field theory does not consider the spatialstructure. As a consequence it fails to reproduce the well-known Debye density of states D ( ω ) ∼ ω d − . One cancure this problem by considering the effective mediumtheory (EMT), where a disordered system is approxi-mated as an effective medium with the effective stiff-ness k eff , namely, H → k eff δ ( p + q ) in the momentumspace [6, 25, 26]. In general, k eff may depend on λ andhave the imaginary part. In the continuum limit, theGreen function at the origin is written as G ( λ ) = Tr( H − λ ) − → G EMT ( λ ) = C (cid:90) q D dq q d − k eff ( λ ) q − λ , (21)where C and q D are constants. Below, we determine k eff ( λ ) by requiring for λ (cid:28) G EMT ( λ ) ∼ G MF , (22)where G MF ( λ ) denotes the Green function calculated bythe mean-field theory. For λ (cid:28) 1, we get G EMT ( λ ) ∼ ( k eff ) − C (cid:48) , (23) where C (cid:48) = C (cid:82) q D dqq d − . Hereafter, we only focus on d > C (cid:48) takes a finite value. On the contrary,the mean-field theory predicts for λ (cid:28) δz (cid:28) G MF ( λ ) = Tr (cid:104) zd W − zed I − λ I (cid:105) − ∼ (cid:104) c δz + ic (cid:112) λ − λ min + O ( λ, δz ) (cid:105) − ,λ min = c δz − c e, (24)where c , . . . , c denote constants, and λ min denotes theminimal eigenvalue of the Hessian. The mean-field theorypredicts that the system is marginally stable λ min = 0in d → ∞ [5], but, the numerical results suggest thatthe marginal stability does not hold λ min > d [6, 27]. Assuming that G EMT ( λ ) ∼ G MF ( λ ) for λ (cid:28) k eff ( λ ) ∼ c δz + ic (cid:112) λ − λ min . (25)The similar result was obtained by a seemingly differentmethod: the coherent potential approximation for a dis-ordered lattice [6]. This is, however, not very surprising,because Eq. (25) is a consequence of the jamming scaling,see Appendix. B for details. The remaining calculationis actually almost identical with a previous work [6]. So,here we just summarize the main steps, see Ref. [6] fordetails.The eigenvalue distribution ρ ( λ ) is calculated as ρ ( λ ) = lim ε → π Im G ( λ − iε ) ∼ lim ε → π Im (cid:90) q D dq q d − k eff q − λ + iε . (26)For λ (cid:28) 1, the asymptotic form of ρ ( λ ) is ρ ( λ ) ∼ (cid:40) Im [ k eff ( λ )] − λ > λ min [ k eff ( λ )] − d/ λ d/ − λ ≤ λ min , (27)where we used the fact that Im k eff = 0 and thuslim ε → Im( k eff q − λ + iε ) − = πδ ( k eff q − λ ) for λ ≤ λ min .The density of states D ( ω ) = 2 ωρ ( λ = ω ) for ω (cid:28) D ( ω ) ∼ (cid:40) ω √ ω − ω ω + ω ∗ ω min (cid:28) ω (cid:28) c δz + c ω min ) − d/ ω d − ω ≤ ω min , ∼ const ω ∗ (cid:28) ω (cid:28) δz − ω ω min (cid:28) ω (cid:28) ω ∗ δz − d/ ω d − ω (cid:28) ω min , (28)where ω min = λ / , and ω ∗ = (cid:112) c δz − c λ min . Herewe assumed ω min (cid:28) ω ∗ , which is consistent with theresults of numerical simulations. If ω min = 0, we recoverthe qualitatively same result as mean-field Eq. (19) for ω (cid:28) e (cid:28) 1. Contrarily, if ω min > 0, we recoverthe well-known Debye density of states D ( ω ) ∼ ω d − for ω ≤ ω min . B. Fluctuated effective medium theory Another important finite d effect is the spatial hetero-geneity: the pre-stress e fluctuate over the space. Toconsider this effect, we write e as e = e c (1 − χ ) , (29)where e c = c δz /c denotes the critical pre-stress de-fine by λ min ( e c ) = 0, and χ denotes a random positivevariable [19]. In the limit d → ∞ near the jamming tran-sition point, the system is marginal stable, meaning thatthe distribution of χ is the delta function P ( χ ) = δ ( χ ).For finite d , on the contrary, P ( χ ) may have a finite widthof distribution around χ = 0: P ( χ ) ∼ (cid:40) ∆ − χ (cid:28) ∆0 χ (cid:29) ∆ , (30)where ∆ denotes a small constant. Using Eq. (27), weget ρ ( λ ) = (cid:90) ∞ dχP ( χ ) ρ ( λ ) ∼ (cid:90) ∞ dχP ( χ ) θ ( λ − c δz χ )Im [ k eff ( λ )] − + (cid:90) ∞ dχP ( χ ) θ ( c δz χ − λ ) [ k eff ( λ )] − d/ λ d/ − ∼ λ − / λ ∗ (cid:28) λ (cid:28) δz − λ / c e ∆ (cid:28) λ (cid:28) λ ∗ δz − λ / + δz − d/ λ d/ − λ (cid:28) c e ∆ , (31)and the density of states is D ( ω ) ∼ const ω ∗ (cid:28) ω (cid:28) δz − ω ω (cid:28) ω (cid:28) ω ∗ ∆ − δz − ω + δz − d/ ω d − ω (cid:28) ω , (32)where ω = √ c e ∆. Note that for ω (cid:28) ω , there arisesthe quartic mode D ( ω ) ∼ ω in addition to the Debyemode D ( ω ) ∼ ω d − . For d < 5, the Debye mode over-whelms the quartic mode. However, detailed numericalinvestigations reported that the quartic modes indeedarise even in d = 3, if the Debye modes are carefullyremoved [2, 28, 29]. V. SUMMARY In this work, we first calculated the contact number z and vibrational density of states D ( ω ) for harmonicspheres in the large spatial dimensions d → ∞ . Our the-oretical results well agree with the results of the previ-ous numerical simulation in large d . Next, we discussed how the finite d effects alter the scaling of D ( ω ). Weshowed that the mean-field result D ( ω ) ∼ ω is replacedby D ( ω ) ∼ c ω + c ω d − in finite d . This is also consis-tent with the numerical results. ACKNOWLEDGMENTS We thank F. Zamponi for useful comments. Thisproject has received JSPS KAKENHI Grant NumbersJP20J00289 and 19J20036. Appendix A: Another derivation of H of harmonicspheres in d → ∞ Here we provide another derivation of Eq. (11). Westart from the second term H (2) . Since H (2) is a realand symmetric matrix, we can diagonalize it using theorthogonal matrix O :( OH (2) O ) ia,jb = λ ia δ ia,jb . (A1)We now use the first approximation. In the spirit of themean-field theory, we replace λ ia by its mean value. λ ia → N d ,dN (cid:88) jb λ jb = 1 N d ,dN (cid:88) jb H (2) jb,jb = zd e (A2)where we have defined the pre-stress [18]: e = − d − N z Nz/ (cid:88) µ =1 h µ r µ . (A3) e is proportional to the pressure, e ∼ p , and vanishesat the jamming transition point. In a previous numericalsimulation [13], the packing fraction was used as a controlparameter, however, it has been recently realized that e is a more natural control parameter [18].The first term H (1) can be calculated as (cid:16) OH (1) O (cid:17) ia,jb = Nz/ (cid:88) µ =1 ( O · n µ ) ia ( O · n µ ) jb , (A4)where n µ is a N d dimensional vector. For µ = ij , its ka component is defined as (cid:0) n ij (cid:1) ka ≡ ∂h ij ∂x ka = ( δ ik − δ jk ) x ia − x ja | r i − r j | . (A5)For amorphous solids both O and n µ are random. So,we now use the second approximation:( O · n µ ) ia → (cid:114) N d ξ µia , (A6)where ξ µia denotes an i.i.d random variable of zero meanand unit variance. The pre-factor in Eq. (A6) has beenchosen so that (cid:80) Ndia ( O · n µ ) ia = (cid:80) Ndia ( n µia ) = 2 onaverage.Using the approximations, Eqs. (A2) and (A6), we geta simple expression:( OHO ) ia,jb → zd W ia,jb − zd eδ ia,jb , (A7)where W ia,jb = 2 N z Nz/ (cid:88) µ =1 ξ µia ξ µjb (A8)denotes the so-called Wishart matrix. Appendix B: Heuristic derivation of k eff ( λ ) Here we calculate k eff by assuming the jamming scal-ing. For the effective medium theory or any other ap-proximation, k eff is calculated by solving a self-consistentequation: F ( k eff , δz, e, λ ) = 0 , (B1) where the specific form of F may depend on the details ofthe approximation. On the other hand, the scaling nearjamming is e ∼ δz , k eff ∼ δz, λ ∼ δz . (B2)If F is a regular function, we can expand F with respectto k eff , δz , e , and λ . The O ( δz ) order terms may bewritten as λ = ak + bδzk eff + cδz + de, (B3)where a , b , c , and d are unknown constants. Solving theabove equation, we get k eff ( λ ) = c (cid:104) c δz + ic (cid:112) λ − λ min (cid:105) ,λ min = c δz − c e, (B4)where c , . . . c are constants. Herewith we get essentiallythe same result as Eq. (25). [1] C. Kittel and P. McEuen, Introduction to solid statephysics , Vol. 8 (Wiley New York, 1976).[2] H. Mizuno, H. Shiba, and A. Ikeda, Proceedings of theNational Academy of Sciences , E9767 (2017).[3] W. A. Phillips and A. Anderson, Amorphous solids: low-temperature properties , Vol. 24 (Springer, 1981).[4] G. Biroli and J.-P. Bouchaud, Structural Glasses andSupercooled Liquids: Theory, Experiment, and Appli-cations , 31 (2012).[5] S. Franz, G. Parisi, P. Urbani, and F. Zamponi, Pro-ceedings of the National Academy of Sciences , 14539(2015).[6] E. DeGiuli, A. Laversanne-Finot, G. D¨uring, E. Lerner,and M. Wyart, Soft Matter , 5628 (2014).[7] M. Shimada and E. De Giuli, arXiv preprintarXiv:2008.11896 (2020).[8] G. Parisi, arXiv preprint arXiv:1401.4413 (2014).[9] Y. Beltukov, JETP Letters , 345 (2015).[10] G. M. Cicuta, J. Krausser, R. Milkus, and A. Zaccone,Phys. Rev. E , 032113 (2018).[11] M. Baggioli, R. Milkus, and A. Zaccone, Phys. Rev. E , 062131 (2019).[12] M. L. Manning and A. J. Liu, EPL (Europhysics Letters) , 36002 (2015).[13] P. Charbonneau, E. I. Corwin, G. Parisi, A. Poncet, andF. Zamponi, Phys. Rev. Lett. , 045503 (2016).[14] C. S. O’Hern, L. E. Silbert, A. J. Liu, and S. R. Nagel,Phys. Rev. E , 011306 (2003).[15] M. Wyart, L. E. Silbert, S. R. Nagel, and T. A. Witten,Phys. Rev. E , 051306 (2005). [16] S. Franz, G. Parisi, M. Sevelev, P. Urbani, and F. Zam-poni, (2017).[17] G. Parisi, P. Urbani, and F. Zamponi, Theory of sim-ple glasses: exact solutions in infinite dimensions (Cam-bridge University Press, 2020).[18] M. Shimada, H. Mizuno, L. Berthier, and A. Ikeda,Physical Review E , 052906 (2020).[19] H. Ikeda, Phys. Rev. E , 050901 (2019).[20] S. Franz and G. Parisi, Journal of Physics A: Mathemat-ical and Theoretical , 145001 (2016).[21] M. M¨uller and M. Wyart, (2015).[22] H. Ikeda, Phys. Rev. Research , 033220 (2020).[23] G. Livan, M. Novaes, and P. Vivo, Introduction to ran-dom matrices: theory and practice , Vol. 26 (Springer,2018).[24] J. Kurchan, G. Parisi, P. Urbani, and F. Zamponi, TheJournal of Physical Chemistry B , 12979 (2013).[25] W. Schirmacher, G. Ruocco, and T. Scopigno, Phys.Rev. Lett. , 025501 (2007).[26] M. Wyart, EPL (Europhysics Letters) , 64001 (2010).[27] C. Scalliet, L. Berthier, and F. Zamponi, Phys. Rev.Lett. , 205501 (2017).[28] E. Lerner, G. D¨uring, and E. Bouchbinder, Phys. Rev.Lett. , 035501 (2016).[29] L. Angelani, M. Paoluzzi, G. Parisi, and G. Ruocco,Proceedings of the National Academy of Sciences115