JJamming below upper critical dimension
Harukuni Ikeda ∗ Graduate School of Arts and Sciences, The University of Tokyo 153-8902, Japan (Dated: July 14, 2020)Extensive numerical simulations in the past decades proved that the critical exponents of thejamming of frictionless spherical particles are the same in two and three dimensions. This impliesthat the upper critical dimension is d u = 2 or lower. In this work, we study the jamming transitionbelow the upper critical dimension. We investigate a quasi-one-dimensional system: disks confinedin a narrow channel. We show that the system is isostatic at the jamming transition point as in thecase of standard jamming transition of the bulk systems in two and three dimensions. Nevertheless,the scaling of the excess contact number shows the linear scaling. Furthermore, the gap distributionremains finite even at the jamming transition point. These results are qualitatively different fromthose of the bulk systems in two and three dimensions. Introduction. –
When compressed, particles interact-ing with finite ranged potential undergo the jammingtransition at the critical packing fraction ϕ = ϕ J atwhich particles start to touch, and the system acquiresrigidity without showing apparent structural changes [1].One of the most popular models of the jamming transi-tion is a system consisting of frictionless spherical par-ticles [2]. The nature of the jamming transition of themodel is now well understood due to experimental andnumerical investigations in the past decades [1]. A fewremarkable properties are the following: (i) the system isnearly isostatic at ϕ J ; namely, the number of constraintsis just one greater than the number of degrees of free-dom [3, 4], (ii) the excess contact number δz from theisostatic value exhibits the power-law scaling δz ∼ δϕ a where δϕ = ϕ − ϕ J denotes the excess packing fraction [2],(iii) the distribution of the gap between particles g ( h ) ex-hibits the power-law divergence g ( h ) ∼ h − γ at ϕ J [5], and(iv) the critical exponents, a = 1 / γ = 0 .
41, do notdepend on the spatial dimensions d for d ≥ a and γ agree with theresults of the mean-field theories, such as the replicamethod [6–8], variational argument [9, 10], and effec-tive medium theory [11]. This implies that the uppercritical dimension d u , above which the mean-field theoryprovides correct results, is d u ≤
2. An Imry-Ma-typeargument [12] and recent finite-size scaling analysis [13]also suggest d u ≤ d <
2. How-ever, the jammed configuration of a true d = 1 systemis trivial: for ϕ ≥ ϕ J , the number of contacts per parti-cle is just z = 2, unless next nearest neighbor particlesbegin to interact at very high ϕ . To obtain non-trivialresults, we consider a quasi-one-dimensional system asshown in Fig. 1, where particles are confined betweenthe walls at y = 0 and y = L y . In the thermodynamiclimit with fixed L y , the model can be considered as aone-dimensional system, but the jammed configurationis still far from trivial. FIG. 1. A configuration at ϕ J for N = 32 and L y = 2 σ max .Gray circles represent particles, and the solid lines denote thecontacts. In the previous works, quasi-one-dimensional systemshave been studied to elucidate the effect of confinementon the jamming transition [14, 15]. These studies uncoverhow the confinement changes the transition point ϕ J [15]and the distribution of the stress near the walls [14].However, the investigation of the critical properties islimited for the systems with very small L y where thejammed configuration is similar to that of the true d = 1system: each particle contact with at most two particles,and therefore one can not discuss the scaling of δz [16–18]. To our knowledge, the scaling of δz for an interme-diate value of L y has not been studied before.In this work, by means of extensive numerical simu-lations, we show that the system is always isostatic atthe jamming transition point for all values of L y , as inthe case of the jamming in d ≥
2. Nevertheless, the criti-cal behavior of the jamming of the quasi-one-dimensionalsystem is dramatically different from the jamming tran-sition in d ≥
2. We find that the excess contact number δz , and the excess constraints δc , which plays a similarrole as δz , exhibit the linear scaling δz ∼ δc ∼ δϕ . Fur-thermore, we find that g ( h ) remains finite even at ϕ J .These results prove that the jamming transition of thequasi-one-dimensional system indeed shows the distinctscaling behaviors from those in d ≥ Model. –
Here we describe the details of our model.We consider two dimensional disks in a L x × L y box. Forthe y -direction, particles are confined between the wallsat y = 0 and y = L y . For the x -direction, we impose theperiodic boundary condition. The interaction potential a r X i v : . [ c ond - m a t . s o f t ] J u l of the model is given by V N = ,N (cid:88) i We perform numerical simulations for N = 1024 disks. We find ϕ J by combining slow compres-sion and decompression as follows [2]. We first generatea random initial configuration at a small packing fraction ϕ = 0 . y = 0 and y = L y . Then, weslowly compress the system by performing an affine trans-formation along the x -direction. For each compressionstep, we increase the packing fraction with a small incre-ment δϕ = 10 − , and successively minimize the energywith the FIRE algorithm [19] until the squared force act-ing on each particle becomes smaller than 10 − . Afterarriving at a jammed configuration with V N /N > − ,we change the sign and amplitude of the increment as δϕ → − δϕ/ 2. Then, we decompress the system until weobtain an unjammed configuration with V N /N < − .We repeat this process by changing the sign and ampli-tude of the increment as δϕ → − δϕ/ V N /N ∈ (10 − , × − ). Wedefine ϕ J as a packing fraction at the end of the abovealgorithm.After obtained a configuration at ϕ J , we re-compressthe system to obtain configurations above ϕ J . As re-ported in Ref. [20], some fraction of samples become un-stable during the compression (compression unjamming).We neglect these samples. We remove the rattlers thathave less than three contacts before calculating physicalquantities. Hereafter, we refer the number of the non-rattler particles as N nr . To improve the statistics, weaverage over 50 independent samples. ϕ J and z J . – First, we discuss the L y dependence ofthe jamming transition point ϕ J and the contact num-ber par particle at that point z J . In Fig. 2 (a), we show ϕ J as a function of σ max /L y . For intermediate values of σ max /L y , ϕ J shows a non-monotonic behavior. A simi-lar non-monotonic behavior has been reported in a pre-vious numerical simulation for a binary mixture [15]. In σ max / L y φ J ( a ) σ max / L y z J ( b ) FIG. 2. L y dependence of (a) the jamming transition point ϕ J and (b) the contact number per particle at the jammingtransition point z J . Markers denote numerical results, andsolid lines denote the guide to the eye. The dashed linesdenote the linear fits ϕ J = 0 . − . σ max /L y and z J =4 − . σ max /L y . the limit σ max /L y → ϕ J converges to its bulk value ϕ bulk J = 0 . 84 as ϕ bulk J − ϕ J ∝ /L y , see the dashed line inFig. 2 (a). The same scaling has been observed in the pre-vious simulation for the binary mixture [15]. The scalingimplies the growing length scale ξ ∼ ( ϕ bulk J − ϕ ) − ν with ν = 1. It is worth mentioning that this is the same expo-nent observed by a correction to scaling analysis [21] andalso our replica calculation for a confined system [22].In Fig. 2 (b), we show z J as a function of σ max /L y . It iswell known that z J = z bulk J = 4 for bulk two dimensionaldisks [2]. In the L y → ∞ limit, z J converges to the bulkvalue as z bulk J − z J ∼ /L y , see the dashed line in Fig. 2(b). Isostaticity. – Next we discuss the isostaticity of ourmodel at ϕ J . The number of degrees of freedom of thenon-rattler particles is N f = 2 N nr − N nr denotesthe number of non-rattler particles, and we neglect theglobal translation along the x -axis. The number of con-strains is N c = N nr z − N w N w = N nr z N w , (2)where z denotes the number of contacts per particle, N w denotes the number of contacts between particles andwalls, and ( N nr z − N w ) / c = N c /N nr . When the system is isostatic N c = N f , we get c = c iso = 2 in the thermodynamic limit. In Fig. 3,we show our numerical result of c at ϕ J as a functionof σ max /L y . This plot proves that the system is alwaysisostatic, irrespective of the value of L y .Now we shall discuss the behavior above ϕ J . As men-tioned in the introduction, we will investigate the modelmainly for L y > σ min so that some fraction of disks canpass through, and thus the contact network undergoes anon-trivial rearrangement on the change of ϕ . σ max / L y c J FIG. 3. L y dependence of the number of constraints per par-ticle at the jamming transition point c J . Markers denote thenumerical results, and the dashed line denotes the isostaticnumber c iso = 2. Ly = σ max Ly = σ max Ly = σ max Ly = σ max Ly = σ max - - - - - - δφ V N / N ( a ) - - - - - - - δφ p ( b ) FIG. 4. (a) δϕ dependence of the energy per particle V N /N .Maker denote numerical results, and the solid line denotes δϕ . (b) δϕ dependence of the pressure p . Maker denotenumerical results, and the solid line denotes δϕ . Energy and pressure. – For ϕ > ϕ J , the particlesoverlap each other. As a consequence, the energy V N and pressure p have finite values. Since we only considerthe compression along the x -axis, we define the pressureas p = − V ∂V N ( { x (cid:48) i } ) ∂ε (cid:12)(cid:12)(cid:12)(cid:12) ε =0 = − V (cid:88) i Next we ob-serve the density dependence of the number of con-straints. For this purpose, we introduce the excess con-straints as δc = N c − ( N f + 1) N nr . (4)where N f + 1 denotes the minimal number of constraintsto stabilize a system consisting of frictionless spherical Ly = σ max Ly = σ max Ly = σ max Ly = σ max Ly = σ max - - - - - - δφ δ c ( a ) - - - - - - ( L y / σ max ) δφ ( L y / σ m a x ) δ c ( b ) FIG. 5. (a) δc as a function of δϕ . Markers denote numericalresults. The solid and dashed lines denote δc ∼ δϕ / and δc ∼ δϕ , respectively. (b) Scaling plot for the same data. particles [4, 12]. For the bulk limit L y ∼ σ max √ N , δc can be identified with the excess contact number δz . Inthis case, the extensive finite size scaling analysis provedthe following scaling form [4]: δc = N − C (cid:0) N δϕ (cid:1) , (5)where the scaling function C ( x ) behaves as C ( x ) ∼ (cid:40) x / x (cid:29) x x (cid:28) . (6)This implies that the square root behavior δc ∼ δϕ / is truncated at δϕ ∼ N − for a finite N system.For δϕ (cid:28) N − , one observes a linear scaling behavior δc ∼ N δϕ [21].To investigate how the behavior changes for L y (cid:28) σ max √ N , in Fig. 5 (a), we show the δϕ dependence of δc for several L y . For large L y and intermediate δϕ , weobserve the square root scaling δc ∼ δϕ / . On the con-trary, for small L y and δϕ , δc shows the linear behavior δc ∼ δϕ . To discuss the scaling behavior more closely,we assume the following scaling form: δc = l αy C (cid:48) (cid:0) l βy δϕ (cid:1) , (7)where l y = L y /σ max , and C (cid:48) ( x ) shows the same scalingbehavior as C ( x ), Eq. (6). When l y ∼ √ N , the scalingshould converge to that of the bulk d = 2 system, Eq. (5).This requires α = − β = 4. In Fig. 5, we test thisprediction. A good scaling collapse verifies the scalingfunction Eq. (7).Note that for a bulk system in d ≥ 2, the system ex-hibits the linear scaling only for δϕ (cid:28) N − : the lin-ear regime vanishes in the thermodynamic limit. Con-trary, Eq. (7) implies that the linear scaling regime per-sists even in the thermodynamic limit for the quasi-one-dimensional system as long as L y is finite. Therefore,the quasi-one-dimensional system indeed has a distinctcritical exponent from that of the bulk systems in d ≥ Ly = σ max Ly = σ max Ly = σ max Ly = σ max Ly = σ max - - - δφ z ( a ) - - - - - - - δφ δ z ( b ) - - - - - ( L y / σ max ) δφ ( L y / σ m a x ) δ z ( c ) FIG. 6. (a) z as a function of δϕ . Markers denote numericalresults. (b) δz = z − z J as a function of δϕ . The solid anddashed lines denote δz ∼ δϕ / and δz ∼ δϕ , respectively. (c)Scaling plot for the same data. Ly = σ max Ly = σ max Ly = σ max Ly = σ max Ly = σ max - - - - - h CD F ( a ) - - - - ( L y / σ max ) μ h ( L y / σ m a x ) CD F ( b ) FIG. 7. (a) CDF of the gap function h . Markers denotenumerical results. The solid and dashed lines denote h − γ and h , respectively. (b) Scaling plot for the same data. In Figs.(a)–(c), we also show the behaviors of the con-tact number per particle z , excess contacts δz = z − z J ,and its scaling plot. The data for δz are more noisy than δc , presumably due to the fluctuation of z J , but still wefind a reasonable scaling collapse by using the same scal-ing form as δc . Gap distribution. – Another important quantity tocharacterize the critical property of the jamming transi-tion is the gap distribution g ( h ). For the bulk systems in d ≥ g ( h ) exhibits the power-law divergence at ϕ J : g ( h ) ∼ h − γ (8)with γ = 0 . 41 [6]. In order to improve the statistics,we observe the cumulative distribution function (CDF)of the gap functions ( h ij and h t , b i ), instead of g ( h ) itself.In this case, the power-law divergence Eq. (8) appears asCDF ∼ h − γ . In Fig. 7 (a), we show our numerical resultsof CDF for several L y . We find that for small L y and h ,CDF ∼ h meaning that g ( h ) remains finite g ( h ) ∼ h even at ϕ J . On the contrary, for large L y , there appearsthe intermediate regime where CDF ∼ h − γ , as in d ≥ ∼ h to CDF ∼ h − γ ,we assume the following scaling form:CDF = l ζy F (cid:48) (cid:0) l ηy h (cid:1) , (9) Ly = σ max Ly = σ max bulk - - - - - - - - δφ δ c ( a ) - - - - - - - h CD F ( b ) FIG. 8. Results for the quasi-two-dimensional system andbulk three dimensional system. (a) δc as a function of δϕ .Markers denote numerical results. The solid and dashed linesdenote δc ∼ δϕ / and δc ∼ δϕ , respectively. (b) CDF ofthe gap function h . Markers denote numerical results. Thesolid and dashed lines denote CDF ∼ h − γ and CDF ∼ h ,respectively. where the scaling function F (cid:48) ( x ) behaves as F (cid:48) ( x ) ∼ (cid:40) x − γ x (cid:29) x x (cid:28) . (10)When l y ∼ √ N , this should converge to the scal-ing form for finite N , CDF( h ) = N − F ( N µ h ), where µ = 1 / (1 − γ ), and F ( x ) shows the same scaling asEq. (10) [23]. This requires ζ = − η = 2 µ . InFig. 7 (b), we check this prediction. The excellent col-lapse of the data for h (cid:28) Quasi-two-dimensional system. – One may suspectthat the distinct scaling of the quasi-one-dimensionalsystem is due to the effect of the boundary condition,not the spatial dimensions. To investigate this possibil-ity, we conduct a numerical simulation for a quasi- two -dimensional system. We consider the same interactionpotential as Eq. (1) with the same system size N = 1024and polydispersity σ i ∈ [1 . , . L x × L y × L z box. As before, particlesare confined between the walls at y = 0 and y = L y ,and the periodic boundary conditions are imposed alongthe x and z directions. We fix L y and change L x = L z to control ϕ . For comparison, we also perform numericalsimulations for the bulk three dimensional system, where L x = L y = L z and the periodic boundary conditions areimposed for all directions. In Fig. 8, we summarize ourresults for δc and CDF of the gaps. One can see thatthe scaling of the quasi-two-dimensional system is thesame as that of the bulk three dimensional system. Thisresult implies that the different scaling of the quasi-one-dimensional system is indeed a consequence of the factthat one dimension is lower than the upper critical di-mension. Conclusions. – In this work, we showed that thejamming transition in a quasi-one-dimensional system isqualitatively different from that in d ≥ δc ∼ δz ∼ δϕ , instead of the square root scaling δz ∼ δϕ / , and the gap distribution g ( h ) remains fi-nite even at ϕ J , instead of the power-law divergence g ( h ) ∼ h − γ .Important future work is to test the robustness of ourresults for other shapes of the quasi-one-dimensional ge-ometries such as a d -dimensional box with an infinitelength in only one direction and fixed lengthes in theother d − Acknowledgements. – We warmly thank M. Ozawa,A. Ikeda, K. Hukushima, Y. Nishikawa, F. Zamponi,P. Urbani, and M. Moore for discussions related to thiswork. We would like in particular to thank the anony-mous referee and M. Ozawa for suggesting the numeri-cal simulation of the quasi-two-dimensional system. Thisproject has received funding from the European Re-search Council (ERC) under the European Union’s Hori-zon 2020 research and innovation program (grant agree-ment n. 723955-GlassUniversality) and JSPS KAKENHIGrant Number JP20J00289. ∗ [email protected][1] A. J. Liu and S. R. Nagel, Annu. Rev. Condens. MatterPhys. , 347 (2010).[2] C. S. O’Hern, L. E. Silbert, A. J. Liu, and S. R. Nagel,Phys. Rev. E , 011306 (2003).[3] J. Bernal and J. Mason, Nature , 910 (1960).[4] C. P. Goodrich, A. J. Liu, and S. R. Nagel, Phys. Rev.Lett. , 095704 (2012).[5] A. Donev, S. Torquato, and F. H. Stillinger, Phys. Rev.E , 011105 (2005).[6] P. Charbonneau, J. Kurchan, G. Parisi, P. Urbani, andF. Zamponi, Nat. Commun. , 3725 (2014).[7] S. Franz and G. Parisi, J. Phys. A , 145001 (2016).[8] S. Franz, G. Parisi, M. Sevelev, P. Urbani, and F. Zam-poni, SciPost Phys. , 019 (2017).[9] M. Wyart, L. E. Silbert, S. R. Nagel, and T. A. Witten,Phys. Rev. E , 051306 (2005).[10] L. Yan, E. DeGiuli, and M. Wyart, EPL , 26003(2016).[11] E. DeGiuli, A. Laversanne-Finot, G. D¨uring, E. Lerner,and M. Wyart, Soft Matter , 5628 (2014).[12] M. Wyart, arXiv preprint cond-mat/0512155 (2005).[13] D. Hexner, P. Urbani, and F. Zamponi, Phys. Rev. Lett. , 068003 (2019).[14] J. W. Landry, G. S. Grest, L. E. Silbert, and S. J. Plimp-ton, Phys. Rev. E , 041303 (2003).[15] K. W. Desmond and E. R. Weeks, Phys. Rev. E ,051305 (2009).[16] S. S. Ashwin and R. K. Bowles, Phys. Rev. Lett. ,235701 (2009).[17] S. S. Ashwin, M. Zaeifi Yamchi, and R. K. Bowles, Phys. Rev. Lett. , 145701 (2013).[18] M. J. Godfrey and M. A. Moore, Phys. Rev. E , 032111(2014).[19] E. Bitzek, P. Koskinen, F. G¨ahler, M. Moseler, andP. Gumbsch, Phys. Rev. Lett. , 170201 (2006).[20] K. VanderWerf, A. Boromand, M. D. Shattuck, and C. S.O’Hern, Phys. Rev. Lett. , 038004 (2020).[21] D. V˚agberg, D. Valdez-Balderas, M. A. Moore, P. Olsson,and S. Teitel, Phys. Rev. E , 030303(R) (2011).[22] H. Ikeda and A. Ikeda, EPL , 40007 (2015).[23] H. Ikeda, C. Brito, and M. Wyart, J. Stat. Mech.: The-ory Exp. , 033302 (2020).[24] Note that our scaling prediction does not work for h ∼109