Jamming transition and normal modes of polydispersed soft particle packing

Abstract

The jamming transition of soft particles characterized by narrow size distributions has been well studied by physicists. However, polydispersed systems are more relevant to engineering, and the influence of polydispersity on jamming phenomena is still unexplored. Here, we numerically investigate jamming transitions of polydispersed soft particles in two dimensions. We find that polydispersity strongly influences contact forces, local coordination, and the jamming transition density. In contrast, the critical scaling of pressure and elastic moduli is not affected by the particle size distribution. Consistent with this observation, we find that the vibrational density of states is also insensitive to the polydispersity. Our results suggest that, regardless of particle size distributions, both mechanical and vibrational properties of soft particle packings near jamming are governed by the distance to jamming.

---

Introduction

Soft particles such as foams, emulsions, and granular materials are ubiquitous in our daily lives. They are of great importance to technologies, including food, granular, and pharmaceutical products.^1,2 It is now well known that soft particles exhibit a rigidity transition, i.e., the so-called jamming transition, at critical packing fraction ϕ_c.^3–5 In recent years, critical behavior of their mechanical, geometrical, and rheological properties (e.g., pressure, elastic moduli, excess coordination number, the first peak of radial distribution function, and viscosity) near jamming has well been tested by numerous experiments and simulations.^6–11 Furthermore, disordered configurations of jammed soft particles contrast sharply with periodic structures of regular lattices such that their normal modes are distinct from those of usual solids.¹² For instance, the vibrational density of states (VDOS) of jammed soft particles exhibits a plateau extending down to zero frequency as the system approaches the unjamming transition.^13–18 In addition, quasi-localized modes coexist with low-frequency vibrations^19–30 and the non-Debye scaling of VDOS is observed in between the low-frequency and plateau regimes.^31–33 It is also theoretically and numerically confirmed that the elastic^34–37 and complex moduli^38–40 are directly linked to the VDOS, hence linear (visco)elasticity of soft particle packings can be predicted from knowledge of low frequency (long wavelength) vibrations of the particles.^36–40

Though the jamming transition and normal modes of soft particle packings have extensively been explored by the theories, experiments, and numerical simulations, most previous works assumed that the particles are either monodispersed, bidispersed, or weakly polydispersed. Since the seminal work by O’Hern et al.³ employed monodisperse systems in three dimensions and bidisperse mixtures of soft particles with a size ratio of 1.4 in two dimensions, these systems have become canonical reference points. Much less attention has been paid to systems with broadly distributed particle sizes. Nevertheless, polydisperse systems are intrinsically relevant to geophysics and civil engineering,⁴¹ because grain sizes in a seismic fault are power law-distributed.^42,43 Moreover, the particle size distribution of Apollonian packings is given by a power law.⁴⁴ There are also indications that polydispersity plays an important role in jamming and elasticity. Both the critical packing fraction ϕ_c and bulk modulus of bidisperse mixtures are sensitive to the size ratio.⁴⁵ In addition, recent study of droplets with a power law distribution showed that pressure and ϕ_c are controlled by distribution's exponent.⁴⁶ In contrast, the same study found that the distribution of coordination number is independent of the same exponent. Similarly, an experimental study of polydispersed granular particles revealed that the macroscopic friction coefficient in the critical state is not affected by polydispersity.⁴⁷ These contrasting results highlight the need for a systematic study of the interplay between polydispersity and its interplay with critical scaling near the jamming transition.

In this communication, we report results of simulations of polydispersed soft particles in d = 2 spatial dimensions generated with molecular dynamics (MD). We analyze packings’ statistics, geometry, and mechanics by systematically varying their packing fraction ϕ and their polydispersity, quantified by a ratio between the maximum size to the minimum size λ. We find that some features are directly controlled by λ, and we quantify the form of this dependence. These “sensitive” features include the distributions of local forces and coordination (which broaden dramatically) and the critical packing fraction (which increases). Other features are surprisingly insensitive to λ. These include the mean coordination, pressure, elastic moduli, and VDOS.

Numerical methods

We study two-dimensional polydispersed particles using MD simulations. We model a repulsive force between the particles, i and j, in contact by a linear spring as f_ij = kδ_ijn_ij with the stiffness k. Here, n_ij = (r_i − r_j)/r_ij with the center-to-center distance, r_ij ≡ |r_i − r_j|, is a unit vector parallel to the normal direction, where r_i(r_j) is the position of the i-th (j-th) particle. In addition, δ_ij ≡ d_ij − r_ij is introduced as the overlap between the particles, with d_ij ≡ R_i + R_j. We randomly sample each particle radius R_i from a power-law size distribution, P(R_i) ∝ R_i^−ν, with the exponent ν. The distribution function is limited to the range, R_min < R_i < R_max, so that we can control polydispersity of the system by changing a size ratio, λ ≡ R_max/R_min. In this study, we fix ν = 3, which is typical of the size distributions of grains in seismic faults,‡ and vary λ from 2 to 20. In ESI,† we show the size distributions P(R_i). We also examined the influence of the power-law exponent ν and confirmed that our results are qualitatively the same if the exponent is in the range, 3 ≤ ν ≤ 4.

To make a static packing of polydispersed particles, we randomly distribute the N = 2048 particles in a L × L square periodic box such that packing fraction of the particles is given by . We then minimize the elastic energy using the FIRE algorithm⁴⁸ with all particle masses set to unity.

Fig. 1 displays snapshots of static packings after minimization. The packing fraction is given by ϕ = 0.9 and the size ratio increases from (a) λ = 2 to (c) 20. See ESI† for a full image of (c). If the size ratio is small, the force network (solid lines) is homogeneous in space and the local coordination number for each particle varies little (Fig. 1(a)). However, as λ increases, one observes that the forces become heterogeneous and the coordination numbers for the largest particles are much larger than for smaller particles (Fig. 1(b) and (c)).

Fig. 1 Snapshots of polydispersed particles (circles), where the packing fraction is ϕ = 0.9 and the size ratio increases as λ = (a) 2, (b) 10, and (c) 20. The solid lines represent force-chains, where their width is proportional to the magnitude of repulsive force between the particles in contact, i.e. kδ_ij. See ESI† for a full image of (c).

Statistics of contact forces and coordination number

The force heterogeneity evident in Fig. 1 is reflected in the distribution function of contact forces P(f), which broadens with increasing polydispersity. Fig. 2(a) displays P(f) for size ratios from λ = 2 to 20 (see ESI† for the full data set). Here, f represents the magnitude of repulsive force f_ij, i.e. kδ_ij, and the horizontal axis is scaled by the average 〈f〉 for each λ. If λ is small, e.g. λ = 2, P(f) is well fitted to a Gaussian distribution (solid line). However, the tail broadens in highly polydispersed packings, such that P(f) is better described with an exponential at large f, P(f) ∼ exp(−f/〈f〉) (dashed line). A similar crossover from compressed-exponential to exponential tails has been observed as a result of other physical parameters, including increasing stress anisotropy,^49,50 increasing spatial dimension,^51,52 increasing particle asphericity,⁵³ and decreasing distance to the unjamming transition.⁵⁴ Our results add polydispersity to this list.

Fig. 2 (a) Semi-logarithmic plots of the distribution function of contact forces f and (b) double logarithmic plots of the distribution function of coordination number z. The packing fraction of the particles is given by ϕ = 0.90 and the size ratio λ increases as listed in the legend of (a) and indicated by the arrows. In (a), the horizontal axis is scaled by the average 〈f〉 for each λ and the solid line represents a Gaussian fit to the data of λ = 2. The dashed lines indicate (a) the exponential tail, P(f) ∼ exp(−f/〈f〉), and (b) power-law, P(z) ∼ z^−4.2, for the data of λ = 20. The inset to (b) shows that the cutoff scales as z* ∼ λ^0.74 (dashed line).

The distribution of coordination number, P(z), is also affected by the polydispersity. Fig. 2(b) shows that P(z) broadens with increasing λ (see ESI† for the full data set). For sufficiently large λ, the distribution approaches a power law, P(z) ∼ z^−4.2 (dashed line), with a large-z cutoff. One might expect the cutoff to be proportional to the perimeter of the largest disks, and therefore linear in λ; instead we find z* ∼ λ^0.74 (see inset), which grows more slowly but still diverges. These results complement a previous study,⁴⁶ which found power law decay for varying exponent ν. In addition, we note that the “granocentric” model^55,56 successfully reproduces the contact number distribution in narrowly polydisperse packings of emulsion droplets. It may be possible to extend the model to broadly polydisperse packings; however the calculation is challenging and beyond the scope of this paper.

In the ESI,† we also examine the radial distribution function, g(r), of polydispersed packings. We find that the first peak of g(r) gets higher and both the first and second peaks shift to shorter distances with the increase of λ.

Jamming transition and critical packing fraction

We have shown that distribution functions are sensitive to polydispersity, as quantified by λ. It is therefore natural to ask to what extent this polydispersity-dependence is inherited by macroscopic (averaged) quantities. In canonical bidisperse packings (λ = 1.4), the pressure p and excess coordination number, Δz ≡ 〈z〉 − z_c, scale as p/k ∼ (ϕ − ϕ_c) and Δz ∼ (ϕ − ϕ_c)^1/2, respectively. Here ϕ_c is the critical packing fraction,^3–5 〈z〉 is the mean coordination number, and z_c = 2d − 2d/N is the central force isostatic value (for N particles in d-dimensions).⁵⁷ We also calculate p and Δz of polydisperse packings to examine the effect of polydispersity on their scaling near jamming. Fig. 3 displays (a) the scaled pressure p/k and (b) Δz as functions of the packing fraction, where λ increases as indicated by the arrows. As can be seen, both p/k and Δz start to increase from zero at ϕ = ϕ_c, where ϕ_c shifts to higher values with the increase of λ. When calculating Δz, we first remove mechanically unstable particles (“rattlers”) from the system. In the ESI,† we show that the fraction of rattlers linearly increases with λ except for the case of λ = 1 (monodisperse packings).

Fig. 3 (a) The scaled pressure p/k and (b) excess coordination number, Δz ≡ 〈z〉 − z_c, as functions of the packing fraction ϕ. The size ratio λ increases as listed in the legend of (a) and indicated by the arrows. The solid lines in (a) represent fitting functions (see the text).

In the ESI,† we also show the dependence of elastic energy E and mean overlap 〈δ〉 on the packing fraction. As p/k and Δz, both E and 〈δ〉 start to increase from zero at ϕ_c, where ϕ_c increases with the increase of λ.

It is apparent from Fig. 3 that the jamming transition density is dependent on polydispersity. ϕ_c is higher in systems with higher polydispersity, because small particles can fill voids between larger ones. To quantify the λ-dependence of ϕ_c, we fit a power law to each p/k dataset and extrapolate the x-intercept. The results are shown in Fig. 4, and are well described by the power-law, (solid line). Here, λ* = 1 and indicate the size ratio and critical density for monodispersed particles, respectively. While we have been unable to find a theoretical explanation for the specific value of the exponent, we note that similar power-law shifts in the critical packing fraction also occur in packings of ellipsoidal particles⁵⁸ and sticky particles.^37,59

Fig. 4 The critical packing fraction ϕ_c(λ) extracted from the data of p/k (Fig. 3(a)). The inset is a double logarithmic plot of and λ − λ*. The solid lines represent the power-law, with λ* = 1 and .

Elastic moduli

Though the dependence of p/k and Δz on the packing fraction is influenced by the polydispersity (Fig. 3), the relation between p/k and Δz is independent of the size ratio λ. Fig. 5(a) displays scatter plots of p/k and Δz, where λ increases as listed in the legend. Strikingly, all the data are nicely collapsed onto the critical scaling, p/k ∼ Δz² (dashed line). In Fig. 5(b)–(d), we also show elastic moduli, i.e. shear modulus G and bulk modulus B, of polydispersed particles and their ratio, i.e. G/B, as functions of Δz. (A discussion of how to calculate the moduli is presented in the ESI.†) All the data in (b)–(d) are well collapsed. The scaled shear modulus exhibits the critical scaling, G/k ∼ Δz (dashed line in (b)). Moreover, the scaled bulk modulus B/k converges to a constant (dashed line in (c)) and the critical scaling, G/B ∼ Δz (dashed line in (d)), can be confirmed as the system approaches the unjamming transition, i.e. as Δz → 0. Therefore, we conclude that the scaling relations p/k ∼ Δz², G/k ∼ Δz, and G/B ∼ Δz, which are all hallmarks of jamming transition, are insensitive to polydispersity. Instead, linear elasticity near jamming is controlled only by the mean coordination number. This is especially surprising in light of the observation that the coordination distribution P(z) is sensitive to λ (Fig. 2).

Fig. 5 The (a) scaled pressure p/k, (b) scaled shear modulus G/k, (c) scaled bulk modulus B/k, and (d) the ratio G/B as functions of the excess coordination number Δz. The size ratio λ increases as listed in the legend of (a). The dashed lines in (a), (b), and (d) represent the critical scaling, i.e. (a) p/k ∼ Δz², (b) G/k ∼ Δz, and (d) G/B ∼ Δz, respectively. In (c), B/k converges to a constant (dashed line) as the system approaches the unjamming transition, Δz → 0.

Normal modes

A system's elastic moduli are determined by its vibrational properties.^38,60 We therefore examine whether jammed systems’ vibrational properties display the same insensitivity to polydispersity seen in G and B. We calculate eigenfrequencies ω of polydispersed packings by diagonalizing their dynamical matrix (see ESI†). Fig. 6(a) displays the vibrational density of states (VDOS, i.e. distribution function of ω) with the horizontal axis non-dimensionalized by a time unit, (m₀ is the particle mass). The size ratio is λ = 20, and the excess coordination number increases as indicated by the arrow. As in the case of bidispersed packings,^3–5 the VDOS exhibits a plateau (horizontal dashed line) above a characteristic frequency, ω_* < ω. As shown in Fig. 6(b), the characteristic frequency is linear in the excess coordination number as ω_* ∼ Δz. In contrast, the VDOS is unaffected by the polydispersity (Fig. 6(c)). The vibrational properties of polydispersed packings are therefore governed only by the coordination number – just like their elastic moduli. In the ESI,† we show our results of the participation ratio P_r(ω), where P_r(ω) in intermediate frequencies (the plateau regime in the VDOS) slightly decreases with the increase of λ.

Fig. 6 (a) The VDOS of polydispersed particle packings, where the size ratio is given by λ = 20. The excess coordination number increases from Δz = 10⁻¹ to 10^0.1 (symbols) as indicated by the arrow. (b) A scaling data collapse of the VDOS, where the horizontal axis is divided by the excess coordination number Δz. (c) The VDOS with Δz = 10^−0.5, where λ increases as listed in the legend.

Summary & outlook

In this study, we have numerically investigated polydispersed particle packings in two dimensions. We found that the distributions of contact force broadens with increasing polydispersity λ. The contact number distribution also broadens, developing a power law tail with a characteristic cutoff z* ∼ λ^0.74. The jamming density ϕ_c monotonically increases as . In contrast, the critical scaling of pressure and elastic moduli are unaffected by λ. Furthermore, the VDOS is independent of polydispersity and the critical scaling of characteristic frequency, ω_* ∼ Δz, is the same with the results of monodisperse and bidisperse systems. Therefore, the mechanical properties and normal modes of soft particle packings are not affected by the particle size distribution and governed only by the excess coordination number.

In our MD simulations, we employed a minimal model of polydispersed particles, where every contact force is proportional to the same stiffness k and every mass m₀ is unique. Because the variance of particle mass, i.e., m_i (i = 1,…,N), merely rescales each row of the dynamical matrix, we do not expect that distributions of particle mass significantly change the VDOS and elastic moduli. As demonstrated in random elastic networks,⁶¹ however, the scaling of VDOS and characteristic frequency ω_* is controlled by the stiffness distribution. Thus, the effect of stiffness distribution on our results have to be examined, which we leave as a future work. Moreover, further analysis in three dimensions is useful for practical applications of this work and linear viscoelastic properties of polydispersed particles are also an interesting topic for future works.

Author contributions

KS and BPT designed the research and wrote the article. KS performed the research.

Data availability

The data supporting this article have been included as part of the ESI.†

Conflicts of interest

There are no conflicts of interest to declare.

Acknowledgements

We thank K. Yoshii, H. Katsuragi, and H. Ikeda for fruitful discussions. This work was financially supported by JSPS KAKENHI Grant Number 22K03459 and the Information Center of Particle Technology.

References

1. R. G. Larson, The Structure and Rheology of Complex Fluids, Oxford University Press, New York, 1999.
2. D. Bonn, M. M. Denn, L. Berthier, T. Divoux and S. Manneville, Rev. Mod. Phys., 2017, 89, 035005.
3. C. S. O'Hern, L. E. Silbert, A. J. Liu and S. R. Nagel, Phys. Rev. E, 2003, 68, 011306.
4. M. van Hecke, J. Phys.: Condens. Matter, 2010, 22, 033101.
5. A. J. Liu and S. R. Nagel, Annu. Rev. Condens. Matter Phys., 2010, 1, 347.
6. T. S. Majmudar, M. Sperl, S. Luding and R. P. Behringer, Phys. Rev. Lett., 2007, 98, 058001.
7. W. G. Ellenbroek, E. Somfai, M. van Hecke and W. van Saarloos, Phys. Rev. Lett., 2006, 97, 258001.
8. L. E. Silbert, A. J. Liu and S. R. Nagel, Phys. Rev. E, 2006, 73, 041304.
9. P. Olsson and S. Teitel, Phys. Rev. Lett., 2007, 99, 178001.
10. G. Katgert, B. P. Tighe and M. van Hecke, Soft Matter, 2013, 9, 9739–9746.
11. K. N. Nordstrom, E. Verneuil, P. E. Arratia, A. Basu, Z. Zhang, A. G. Yodh, J. P. Gollub and D. J. Durian, Phys. Rev. Lett., 2010, 105, 175701.
12. N. W. Ashcroft and N. D. Mermin, Solid State Physics, Saunders College Publishing, US, 1976.
13. L. E. Silbert, A. J. Liu and S. R. Nagel, Phys. Rev. Lett., 2005, 95, 098301.
14. M. Wyart, L. E. Silbert, S. R. Nagel and T. A. Witten, Phys. Rev. E, 2005, 72, 051306.
15. M. Wyart, S. R. Nagel and T. A. Witten, Europhys. Lett., 2005, 72(3), 486.
16. N. Xu, M. Wyart, A. J. Liu and S. R. Nagel, Phys. Rev. Lett., 2007, 98, 175502.
17. L. E. Silbert, A. J. Liu and S. R. Nagel, Phys. Rev. E, 2009, 79, 021308.
18. N. Xu, V. Vitelli, A. J. Liu and S. R. Nagel, Europhys. Lett., 2010, 90, 56001.
19. Z. Zeravcic, W. van Saarloos and D. R. Nelson, Europhys. Lett., 2008, 83, 44001.
20. N. Xu, V. Vitelli, A. J. Liu and S. R. Nagel, Europhys. Lett., 2010, 90, 56001.
21. K. Chen, W. G. Ellenbroek, Z. Zhang, D. T. N. Chen, P. J. Yunker, S. Henkes, C. Brito, O. Dauchot, W. van Saarloos, A. J. Liu and A. G. Yodh, Phys. Rev. Lett., 2010, 105, 025501.
22. K. Chen, M. L. Manning, P. J. Yunker, W. G. Ellenbroek, Z. Zhang, A. J. Liu and A. G. Yodh, Phys. Rev. Lett., 2011, 107, 108301.
23. E. Lerner, G. Düring and E. Bouchbinder, Phys. Rev. Lett., 2016, 117, 035501.
24. H. Mizuno, H. Shiba and A. Ikeda, Proc. Natl. Acad. Sci. U. S. A., 2017, 114, E9767.
25. G. Kapteijns, E. Bouchbinder and E. Lerner, Phys. Rev. Lett., 2018, 121, 055501.
26. M. Shimada, H. Mizuno and A. Ikeda, Phys. Rev. E, 2018, 97, 022609.
27. C. Rainone, E. Bouchbinder and E. Lerner, Proc. Natl. Acad. Sci. U. S. A., 2020, 117, 5228.
28. S. Mahajan and M. P. Ciamarra, Phys. Rev. Lett., 2021, 127, 215504.
29. L. Wang, G. Szamel and E. Flenner, Phys. Rev. Lett., 2021, 127, 248001.
30. L. Wang, A. Ninarello, P. Guan, L. Berthier, G. Szamel and E. Flenner, Nat. Commun., 2019, 10, 26.
31. P. Charbonneau, E. I. Corwin, G. Parisi, A. Poncet and F. Zamponi, Phys. Rev. Lett., 2016, 117, 045503.
32. J. Lin, I. Jorjadze, L.-L. Pontani, M. Wyart and J. Brujic, Phys. Rev. Lett., 2016, 117, 208001.
33. F. Arceri and E. I. Corwin, Phys. Rev. Lett., 2020, 124, 238002.
34. C. Maloney and A. Lemaitre, Phys. Rev. Lett., 2004, 93, 016001.
35. C. E. Maloney and A. Lemaitre, Phys. Rev. E, 2006, 74, 016118.
36. H. Mizuno, K. Saitoh and L. E. Silbert, Phys. Rev. E, 2016, 93, 062905.
37. D. J. Koeze, L. Hong, A. Kumar and B. P. Tighe, Phys. Rev. Res., 2020, 2, 032047.
38. B. P. Tighe, Phys. Rev. Lett., 2011, 107, 158303.
39. K. Baumgarten and B. P. Tighe, Soft Matter, 2017, 13, 8368.
40. K. Saitoh, T. Hatano, A. Ikeda and B. P. Tighe, Phys. Rev. Lett., 2020, 124, 118001.
41. T. Ueda, T. Matsushima and Y. Yamada, Powder Technol., 2012, 217, 533.
42. N. Monzawa and K. Otsuki, Tectonophysics, 2003, 367, 127.
43. N. Keulen, R. Heilbronner, H. Stünitz, A.-M. Boullier and H. Ito, J. Struct. Geol., 2007, 29, 1282.
44. O. Ben-Nun, I. Einav and A. Tordesillas, Phys. Rev. Lett., 2010, 104, 108001.
45. D. J. Koeze, D. Vågberg, B. B. T. Tjoa and B. P. Tighe, Eur. Phys. Lett., 2016, 113, 54001.
46. D. S. Shimamoto and M. Yanagisawa, Phys. Rev. Res., 2023, 5, L012014.
47. O. Polana, M. Cabrera, M. Renouf, E. Azéma and N. Estrada, Phys. Rev. E, 2023, 107, L052901.
48. E. Bitzek, P. Koskinen, F. Gähler, M. Moseler and P. Gumbsch, Phys. Rev. Lett., 2006, 97, 170201.
49. B. P. Tighe, J. E. S. Socolar, D. G. Schaeffer, W. G. Mitchener and M. L. Huber, Phys. Rev. E, 2005, 72, 031306.
50. B. P. Tighe, A. R. T. van Eerd and T. J. H. Vlugt, Phys. Rev. Lett., 2008, 100, 238001.
51. A. R. T. van Eerd, W. G. Ellenbroek, M. van Hecke, J. H. Snoeijer and T. J. H. Vlugt, Phys. Rev. E, 2007, 75, 060302.
52. B. P. Tighe, J. H. Snoeijer, T. J. H. Vlugt and M. van Hecke, Soft Matter, 2010, 6, 2908.
53. E. Azéma and F. Radjai, Phys. Rev. E, 2012, 85, 031303.
54. C. S. O'Hern, S. A. Langer, A. J. Liu and S. R. Nagel, Phys. Rev. Lett., 2001, 86, 111.
55. M. Clusel, E. I. Corwin, A. O. N. Siemens and J. Brujić, Nature, 2009, 460, 611.
56. E. I. Corwin, M. Clusel, A. O. N. Siemens and J. Brujić, Soft Matter, 2010, 6, 2949.
57. C. P. Goodrich, A. J. Liu and S. R. Nagel, Phys. Rev. Lett., 2012, 109, 095704.
58. A. Donev, I. Cisse, D. Sachs, E. A. Variano, F. H. Stillinger, R. Connelly, S. Torquato and P. M. Chaikin, Science, 2004, 303, 990–993.
59. D. J. Koeze and B. P. Tighe, Phys. Rev. Lett., 2018, 121, 188002.
60. A. Lematre and C. Maloney, J. Stat. Phys., 2006, 123, 415–453.
61. E. DeGiuli, E. Lerner, C. Brito and M. Wyart, Proc. Natl. Acad. Sci. U. S. A., 2014, 111, 17054.

---

Footnotes

† Electronic supplementary information (ESI) available: Additional results and the theoretical expression of the bulk modulus. See DOI: https://doi.org/10.1039/d4sm01305k

‡ The power-law exponent for cumulative size distributions distributes around 2^42,43 and thus ν = 3 is typical of the size distributions.

---