Exceptionally dense and resilient critically jammed polydisperse disk packings

Abstract

Understanding the way disordered particle packings transition between jammed (rigid) and unjammed (fluid) states is of both great practical importance and strong fundamental interest. The values of critical packing fraction (and other state variables) at the jamming transition are protocol dependent. Here, we demonstrate that this variability can be systematically traced to structural measures of packing, as well as to energy measures inside the jammed regime. A novel generalized simultaneous particle swap algorithm constructs overjammed states of desired energy, which upon decompression lead to predictable critical packing fractions. Thus, for a given set of particle sizes, states with extraordinarily high critical packing fractions can be found efficiently, which sustain substantial shear strain and preserve their special structure over the entire jammed domain. The close relation revealed here between the energy landscape of overjammed soft-particle packings and the behavior near the jamming transition points towards new ways of understanding and constructing disordered materials with exceptional properties.

---

1 Introduction

The jamming transition is a phenomenon observed in a broad range of soft matter, including foams, emulsions, colloidal aggregates, and granular materials, where the packing fraction ϕ of particles controls this rigidity transition^1–4 (Fig. 1a). At the jamming point, constituent particles form a percolating contact network, satisfying isostatic conditions,⁵ and the material behavior near the transition point has been described in the language of critical phenomena.^1,6–9 Numerous previous studies have focused on identifying this transition point and have settled on a critical packing fraction of “random close packing” ϕ_rcp ≈ 0.84 for frictionless disk packings in two dimensions (2D)^1,10 and ϕ_rcp ≈ 0.64 for frictionless sphere packings in three dimensions (3D).^11,12

Fig. 1 Continuous range of critical packing fraction in disk packings. (a) Schematic illustrations of jamming transition controlled by packing fraction, showing changes in contact numbers (color coded) with the isostatic condition satisfied at the jamming point. (b) Histograms of particle–particle distances r_ij show the emergence of overlaps between soft particles beyond the jamming point. Bars represent the accessible range of hard particle packings (ϕ ≤ ϕ_c) and soft particle packings (all ϕ). (c) Schematic of critically jammed structures for monodisperse (top), bidisperse (center), and polydisperse (bottom) frictionless disks. Monodisperse and polydisperse disks are colored by number of neighbors, bidisperse disks by size.

Recent research has shown, however, that the critical packing fraction at the jamming transition ϕ_c is not unique and strongly depends on the preparation protocol of the packing.^13–17 To understand mechanical properties of the material, it is crucial to understand and predict ϕ_c. Fig. 1c illustrates structural changes associated with different ϕ_c in 2D. The packing fraction of the hexagonal crystal ϕ_hex = 0.9069 … is proved to be the highest ϕ_c for monodisperse disks and also provides an upper bound on ϕ_c of bidisperse systems with sizes more similar than a threshold value,^18,19 which structurally corresponds to size segregation. In both mono- and bidisperse systems, lower ϕ_c indicates stronger positional disorder (increased defect density and smaller-scale segregation, respectively). Bidisperse systems with larger size contrast can avoid crystallization entirely and in many simulation protocols settle on a critical density ϕ_rcp ≈ 0.84, although lower-density random packings are possible with specialized protocols in monodisperse systems.²⁰

Systems with significant continuous size polydispersity are less often studied, but avoid long-range order very effectively. Again the vast majority of packing protocols yield ϕ_c ≈ ϕ_rcp. For a variety of given particle size distributions, we here show how to systematically construct packings with ϕ_c covering a wide range between ϕ_rcp close to ϕ_hex, without any long-range order (Fig. 1c). By connecting the concept of the energy landscape for jammed states to the critically jammed structures, our approach allows controlled access to distinct jamming configurations, including unusually dense states that exhibit extraordinary mechanical properties and structural robustness.

2 Results & discussion

2.1 Accessing metastable states by simultaneous swaps

In contrast to many common protocols of the Lubachevsky–Stillinger (LS) type,²¹ which approach the critical jamming point from unjammed (fluid) states by incrementally increasing packing fraction, we instead use polydisperse soft particles, generating overjammed states (ϕ > ϕ_c), which we relax by incrementally decreasing ϕ following^22,23 (Fig. 1b). While hard particle packings cannot access this overjammed regime as the bulk modulus diverges at the jamming point, soft particle packings can be stabilized within this regime under finite pressure. After decompression down to the jamming point, the resulting configurations can thus represent both hard and soft critically jammed frictionless disks.

Soft-particle packings in the overjammed regime can be modeled using interparticle elastic pair potentials in a simple yet generic manner (cf.Fig. 1b; see Methods). The entire system is then described (at zero temperature) by a multidimensional energy landscape whose degrees of freedom are the particle positions. An initial configuration relaxes to one of a huge number of mechanical equilibrium states, representing local minima of the energy landscape (metastable states, MS), which correspond to inherent structures in the parlance of the glass literature.^24,25 In recent work, the authors showed that deeply overjammed MS (ϕ = 1) with distinct energy levels exhibit distinct structural features as well as mechanical behavior, and that structural measures of the states can be utilized to efficiently construct low energy MS.²⁶ While this approach is effective in finding low energy MS, it is incapable of constructing MS at controlled higher energy levels, limiting our ability to efficiently explore the metastable energy landscape.

To access MS at different energy levels, we develop a novel algorithm, illustrated in Fig. 2a. The approach entails four steps: (1) generating an initial MS by annealing a random configuration, (2) converting the particle packing to a network system based on radical tessellation with a target network connectivity z_n, (3) annealing the network system with particle sizes as additional degrees of freedom, and (4) converting the network system back to a particle packing, reassigning original particle sizes in the order of the altered particle sizes, and subsequently annealing to a final MS (see Methods for details). The algorithm was evaluated for a variety of particle numbers (N = 144 to N = 1024), particle size distribution shapes (gamma, log-normal, and uniform area distributions), and size polydispersity (particle area coefficient of variation 0.2 ≤ c_A ≤ 0.5), and proved effective across all conditions. For clarity in presentation, we will concentrate on results for gamma distribution with c_A = 0.4 in the following sections.

Fig. 2 Novel algorithm for accessing metastable states for all MS energy levels. (a) Schematic of the generalized angle swap simulation protocol to generate the entire range of metastable state energies for overjammed packings. (b) Correlation between the average network coordination number z_n and the corresponding metastable state energy ε_r at ϕ = 1. ε_r,0 represents the energy level of monodisperse hexagonal packing at ϕ = 1 (see eqn (1)–(3)). (c) Correlation between z_n and mean contact number of the overjammed particle packing z.

This algorithm efficiently generates metastable states of overjammed packings for the entire energy range, with z_n as a control parameter (Fig. 2b). The energy is insensitive to z_n in the range 1 ≲ z_n ≲ 5, but sharply decreases from the highest values at z_n = 0 and again towards the lowest values at z_n = 6. The approach remains effective for a range of overjammed ϕ values (see ESI†). The mean number of particle contacts z in the final metastable state is similarly controlled by z_n (Fig. 2c), in agreement with the reported negative correlation between z and MS energy.²⁶Fig. 2b and c also demonstrates that the correlation is independent of the system size. The mean displacement of a particle between an initial MS and a final MS is far less than an average particle size in most cases, implying that local adjustment of particle positions and sizes is sufficient to lower energy significantly (ESI†). In particular, the contact topology of final metastable states faithfully retains the network topology for z_n ≳ 5.5.

2.2 Decompressing overjammed packings down to the critical jamming point

To investigate the relation between the metastable energy landscape and the jamming transition, MS of distinct energy levels at ϕ = 1 are quasi-statically decompressed by shrinking all particles proportionally until the system approaches the jamming transition point^22,23 (see Methods). The resulting critical jamming configurations are verified to fulfill the isostatic condition.

A main result of the present work is that the critical packing fraction obtained in this way shows a strong negative correlation with the energy of the overjammed metastable state it was generated from. Thus, a wide range of ϕ_c can be efficiently and systematically constructed (Fig. 3a and e). High energy MS at ϕ = 1 undergo the jamming transition at ϕ_c ≈ 0.84, agreeing with common values of ϕ_rcp.^1,27 As MS energy decreases, the resulting critical packing fraction continuously increases, and the lowest-energy MS transition at ϕ_c ≥ 0.89, significantly larger than typical values investigated before.^23,28,29 The correlation between MS energy and ϕ_c does not change for different system size (Fig. 3e), while the maximum ϕ_c continues to increase with N. This suggests that at least as wide a range of ϕ_c can be realized in the limit of large system size. This result is consistent with several previous studies finding that the jamming transition occurs at a continuous range of packing fractions and that it strongly depends on the preparation protocol.^13,16,30 The generalized swap algorithm presented here, however, provides a far more efficient way of constructing critically jammed configurations of exceptionally high packing fraction.

Fig. 3 Construction of critical jammed states from overjammed MS. Representative jamming configurations for different critical packing fractions, particles in gray (a), color-coded by hexagonal orientation order (b), number of neighbors by radical tessellation topology (c), and angle order parameter (d). (e) The packing fraction ϕ_c of critically jammed states is strongly predicted by the energy level of corresponding overjammed packings. Symbols represent different sizes of individual simulation systems. (f) Hexagonal orientation order parameter 〈Ψ₆〉 and (g) coefficient of variation of radical tessellation topology c_n show negligible correlation with ϕ_c, while (h) the angle order parameter Θ shows excellent correlation.

It is important to emphasize that the critically jammed states represent an arrangement of particles with unchanged given sizes {A_i}, even though the algorithm transiently replaces them with altered disk areas . Had we constructed jammed states from the , even higher area fractions could have been obtained (see ESI†), due to the local size optimization. By contrast, we here present unusually dense packing solutions for the practical problem of jamming a given set of particles.

What sets high-ϕ_c polydisperse packings apart structurally? Many traditional measures of amorphous disorder fail to distinguish them from lower-ϕ_c structures. The average of the hexagonal orientation order parameter³¹ |Ψ_6,j| over all particles (see Methods), shown in Fig. 3b and f, increases very weakly with ϕ_c. Likewise, the coefficient of variation of the tessellation topology (see Methods) c_n, a successful indicator of disorder in other contexts,^32,33 is insensitive to ϕ_c (Fig. 3c and g). The radial distribution function g(r) also shows no discernible change across jamming configurations with distinct ϕ_c, making it an unsuitable measure for distinguishing them (ESI†).

By contrast, quantifying the local close-packing of disks directly provides a meaningful structural indicator. Fig. 3d and h confirm that the angle order parameter Θ^26,34,35 (see Methods), which was shown to be an excellent predictor of MS energy at ϕ = 1,²⁶ also correlates strongly with ϕ_c of critical jamming states. Hence, the angle order parameter is a valuable structural measure for the preparation of jamming configurations with a broad range of ϕ_c.

All findings are robust with respect to changing the size distribution shape: Fig. 4 compares results for gamma and log-normal distributions (two unimodal distributions), and by contrast for a uniform area distribution, all at c_A = 0.4. While the maximum ϕ_c for the uniform distribution is slightly lower, and the quantitative measures of disorder slightly different, all trends are the same and in particular the angle order parameter Θ predicts ϕ_c universally. Similar conclusions hold when c_A is varied in the range of interest (see ESI† for details), demonstrating again that this approach produces controlled-density critically jammed structures for many given size distributions. This is of interest as recent work shows that the phase behavior of glassy disk systems does depend on distribution shapes.³⁶

Fig. 4 Robustness of the algorithm across different size distribution shapes. (a) Relation between metastable state energy of overjammed packings and critical packing fraction for different area distribution shapes of the polydisperse systems. Only slight quantitative changes from Fig. 3 are found, confirming the conclusions for dependence on (b) hexagonal orientation order parameter, (c) coefficient of variation of the topology distribution, and (d) angle order parameter. N = 1024 for all data.

2.3 Mechanically resilient dense polydisperse packings

We find that the structures of high-ϕ_c critical states are remarkably resilient to material-scale deformation. To assess how materials respond to volumetric changes, jamming configurations at ϕ_c are re-compressed up to ϕ = 1 (Fig. 5a and b; see Methods for details), and particle center displacements ΔR are monitored. Since the compression protocol involves proportionally increasing particle sizes rather than reducing the periodic box size, these ΔR represent non-affine motion. The average over all particles (R₀ is the radius of a mean-area particle at ϕ = 1, and all rattlers are excluded) is a measure of structural rearrangement. We observe that structures unjamming at the highest ϕ_c have , while those unjamming at the lowest ϕ_c show substantial rearrangement (Fig. 5c). Thus, the particles in extraordinarily high-ϕ_c critically jammed states move almost perfectly affinely under uniform compression, maintaining their specific structure as the packing fraction changes. Very similar results are observed during the decompression process from ϕ = 1 to ϕ_c (ESI†), suggesting that structures with high ϕ_c reliably maintain their underlying structure during volumetric changes.

Fig. 5 Mechanical response of critical jamming states across ϕ_c. (a) Schematic of compressing critical jammed states, transitioning to overjammed states by uniformly growing particles. (b) Superposition of critical (blue dashed) and compressed (red solid) packing configurations for a typical low-ϕ_c critical state (left) and a typical high-ϕ_c critical state (right), indicating pronounced displacements ΔR of particle centers for the low-ϕ_c state and minimal displacements for the high ϕ_c state. (c) Averaged normalized particle displacement during compression from ϕ_c to ϕ = 1 (20 independent simulations of a fixed z_n for each data point). For configurations with the highest ϕ_c, the displacement is much reduced. (d) Schematic of cyclic pure shear deformation at ϕ_c. (e) Superposition of packing configurations before (blue dashed) and after (red solid) a γ = 0.04 shear cycle, for a typical low-ϕ_c critical state (left) and a typical high-ϕ_c critical state (right). Net displacement of individual particles vanishes in the latter case. (f) Mean particle displacement after a cycle as a function of applied strain γ for low-ϕ_c jammed states (yellow) and high-ϕ_c jammed states (purple). The latter states show very small displacement below a critical strain. N = 400 for all data.

As a complementary test of structural robustness, cyclic shear deformations at a range of strain amplitudes γ are applied to each critically jammed configuration and its evolution from an initial configuration is evaluated (Fig. 5d and e; see Methods for details). Fig. 5f illustrates that low-ϕ_c structures readily rearrange under cyclic shear with displacements roughly proportional to γ; they are only marginally stable. Structures with extraordinarily high ϕ_c, however, display qualitatively different behavior: they show negligible rearrangements up to substantially large γ values. Even after a finite displacement occurs, many such systems maintain this new configuration over a finite range of higher γ (ESI†). The high-ϕ_c jammed states have thus inherited the structural resilience of the low-energy overjammed MS from which they were constructed – those metastable states of low energy exhibit ultrastable characteristics, reminiscent of ultrastable glasses.^26,37–42

3 Conclusions

Here we have demonstrated an efficient way to construct packings of polydisperse hard disks with unusually high density. We find that these extraordinary structures are also exceptionally robust against shear and compression – indeed they inherit mechanical resilience from the overjammed MS they are constructed from because their relative particle positions are virtually unchanged from the overjammed states. Our results suggest that reducing the packing fraction of overjammed packings during the decompression process maintains the shape of the energy landscape around its deepest minima, mainly shifting the landscape down along the energy axis (Fig. 6), while higher minima do change their shapes and positions. As ϕ reaches the highest critical jamming fraction, the lowest minimum reaches zero energy and is lost. Further reduction of ϕ is necessary to unjam higher-lying minima, but their identity along the configuration axis is not maintained as faithfully as that of the lowest-energy states.

Fig. 6 Connecting the overjammed packing energy landscape to the jamming transition. Illustration of qualitative changes in the energy landscape upon decompression of an overjammed soft-disk system (ϕ decreases left to right). Particle packings are strictly overjammed in the high packing fraction limit and strictly unjammed in the low packing fraction limit, while a continuous range of packing fraction in between allows coexistence of overjammed, critically jammed, and unjammed states.

The jamming transition and its underlying structural features are closely connected to the understanding of the glass transition and glass dynamics, with low-lying minima on the energy landscape providing information about deep glassy states of low entropy. Very recent work on 2D glassy dynamics systems with polydispersity uses variants of simultaneous swap algorithms to efficiently equilibrate these systems to high packing fraction and high stability,⁴³ a direct indication that the present results have relevance in glassy systems and that the rare, disordered high-ϕ_c configurations may represent analogs of minimal entropy, ideal-glass states.^44,45 The robustness of our results against changes in the particle size distribution also hints at the universal character of these structures. Their further study, and the construction of analogous configurations in three dimensions, where studies of polydispersity effects are of great recent interest,^46,47 will allow systematic exploration of disordered states by packing density. Extraordinary structures being correlated with extraordinary material properties, this has potentially far-reaching applications in the design of specific mechanical, electrical, or thermal behavior in glasses and other amorphous materials.

4 Methods

4.1 Soft particle packing

We model the interaction between particles i and j by a repulsive harmonic potential V(r_ij) (cf.¹) with normalized spring constant, namelyHere, r_ij is the distance between centers of particles i and j, σ_ij = σ_i + σ_j is an equilibrium distance equal to the sum of particle radii, and H(·) is the Heaviside step function. For circular particles, radii follow from areas A_i, whose distribution we fix. We stress that the results presented are independent of the particular functional form (1).

The total non-dimensional energy ε_r of a given metastable state and the reference energy ε_r,0 of a monodisperse hexagonal packing at packing fraction ϕ are then

We employ our custom Matlab code to generate the soft particle packing. An initial configuration is constructed at ϕ = 1 by choosing N particle sizes {A_i} randomly from a given distribution and placing the particles at random positions {R⃑_i} within a periodic box of fixed size. The initial configuration is then annealed to the nearest metastable state using the FIRE algorithm.⁴⁸ We use gamma, log-normal, and uniform area distributions and quantify their width by the coefficient of variation c_A. We restrict the width between c_A = 0.2 and c_A = 0.5 to avoid hexagonal ordering (crystallization) for low polydispersity and Apollonian fractal-like packings for high polydispersity. For results on distributions and polydispersities not discussed in the main text, see ESI.†

4.2 Simultaneous swap algorithm

The equilibrated MS overjammed packing is then translated to a network system of defined topology by finding all neighboring pairs by radical tessellation⁴⁹ along with their normalized distances, r_ij/σ_ij. Identifying particle positions with nodes of the network, network bonds are added sequentially in order of increasing normalized neighbor distance until the average coordination number of the network is larger than or equal to a target value z_n. To find equilibria of the network, we associate a harmonic potential energy with the bonds, with two major differences to the particle system: (1) the network topology is fixed while the contact topology of particle packings can change according to the change of particle positions; (2) node particle sizes, {A_i}, are considered as degrees of freedoms in addition to node positions. This latter treatment is reminiscent of recent work simultaneously optimizing particle positions and sizes for particle packings,⁴⁶ though that work does not enforce fixed network topology. The fixed network topology qualitatively changes the interaction from one-sided repulsive interactions between disks to two-sided repulsive-attractive interactions between nodes. This prevents the collapse of particle sizes to very small values, as these states are energetically unfavorable due to high attractive energy costs. We anneal the network system to the nearest metastable state, resulting in a new set of particle sizes . The initial particle sizes {A_i} are then reassigned to particle positions in the order of the sizes, effecting a simultaneous size swap for the entire particle system. This approach extends the conventional Monte-Carlo swap algorithm, widely used to expedite the search for equilibrium states in glass simulations.^50–52

4.3 Decompression and compression algorithm

For the decompression of overjammed packings, we use an initial decompression step size of Δϕ = 0.01 until the MS energy becomes zero, indicating that the configuration becomes unjammed. Iterative reduction of the step size then accurately finds the ϕ_c of the critical jamming state, defined as a configuration with ε_r between 10⁻¹⁶ and 2 × 10⁻¹⁶, for which we also verify that the isostatic condition is fulfilled and that any further decompression leads to loss of all contacts. In these critical states, the soft particles can be exchanged for hard disks without changing the structure.

For the compression of critically jammed states, we reverse the process and apply a compression step size (typically Δϕ = 0.01) until the packing fraction reaches the desired target value.

4.4 Hexagonal orientation order parameter

Using the numbers n_j of nearest neighbors of particle j in the radical tessellation of the packing, the hexagonal orientation order parameter³¹ is defined aswhere θ_jk is the angle of the line connecting the positions of disks j,k with respect to a reference direction.

4.5 Coefficient of variation of the tessellation topology

In the radical tessellation of the packing, as in any 2D tessellation without borders, Euler's theorem ensures that the mean number of neighbors is equal to six. Therefore, the coefficient of variation of the tessellation topology is

4.6 Angle order parameter

The angle order parameter averages the absolute differences |θ_jkl − θ_jkl,0| for all triplets of particles that are mutual neighbors by radical tessellation, where θ_jkl,0 is the angle of locally sterically optimal packing, with all three particles in contact, and θ_jkl is the actual angle within jammed particle packings, i.e.,see ESI† for more details.

4.7 Cyclic shear algorithm

For the cyclic shear deformation, pure shear strain is first applied quasi-statically with an increment of 0.001 until a target strain γ ≤ 0.1 is reached. Then the system is quasi-statically sheared back to a state of zero strain. If a configuration unjams during shear, it is excluded from the analysis. For each target strain γ, the mean particle displacement is determined.

Author contributions

SK and SH designed research; SK performed all simulations and analysis; SK and SH wrote the paper.

Data availability

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

Conflicts of interest

There are no conflicts of interests to declare.

Acknowledgements

The authors are grateful for helpful discussions with Ludovic Berthier, Eric Corwin, Varda Hagh, Robert B. Kusner, and Vinothan N. Manoharan.

Notes and references

1. C. S. O’Hern, L. E. Silbert, A. J. Liu and S. R. Nagel, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2003, 68, 011306.
2. M. van Hecke, J. Phys.: Condens. Matter, 2009, 22, 033101.
3. D. Bonn, M. M. Denn, L. Berthier, T. Divoux and S. Manneville, Rev. Mod. Phys., 2017, 89, 15.
4. R. P. Behringer and B. Chakraborty, Rep. Prog. Phys., 2018, 82, 012601.
5. W. G. Ellenbroek, V. F. Hagh, A. Kumar, M. F. Thorpe and M. van Hecke, Phys. Rev. Lett., 2015, 114, 135501.
6. W. G. Ellenbroek, E. Somfai, M. van Hecke and W. van Saarloos, Phys. Rev. Lett., 2006, 97, 258001.
7. A. J. Liu and S. R. Nagel, Annu. Rev. Condens. Matter Phys., 2010, 1, 347–369.
8. P. Olsson and S. Teitel, Phys. Rev. Lett., 2007, 99, 178001.
9. J. D. Sartor, S. A. Ridout and E. I. Corwin, Phys. Rev. Lett., 2021, 126, 048001.
10. A. Zaccone, Phys. Rev. Lett., 2022, 128, 028002.
11. J. G. Berryman, Phys. Rev. A: At., Mol., Opt. Phys., 1983, 27, 1053–1061.
12. G. D. Scott and D. M. Kilgour, J. Phys. D: Appl. Phys., 1969, 2, 863.
13. P. Chaudhuri, L. Berthier and S. Sastry, Phys. Rev. Lett., 2010, 104, 165701.
14. A. B. Hopkins, F. H. Stillinger and S. Torquato, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2013, 88, 022205.
15. C. F. Schreck, C. S. OHern and L. E. Silbert, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2011, 84, 011305.
16. M. Ozawa, L. Berthier and D. Coslovich, SciPost Phys., 2017, 3, 027.
17. H. Brouwers, Soft Matter, 2023, 19, 8465–8471.
18. G. F. Tóth, Acta Mathematica Academiae Scientiarum Hungarica, 1972, vol. 23, pp. 263–270.
19. A. Heppes, Discrete & Computational Geometry, 2003, vol. 30, pp. 241–262.
20. S. Atkinson, F. H. Stillinger and S. Torquato, Proc. Natl. Acad. Sci. U. S. A., 2014, 111, 18436–18441.
21. B. D. Lubachevsky and F. H. Stillinger, J. Stat. Phys., 1990, 60, 561–583.
22. N. Xu, J. Blawzdziewicz and C. S. O'Hern, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2005, 71, 061306.
23. K. W. Desmond and E. R. Weeks, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2009, 80, 051305.
24. F. H. Stillinger and T. A. Weber, Science, 1984, 225, 983–989.
25. P. G. Debenedetti and F. H. Stillinger, Nature, 2001, 410, 259–267.
26. S. Kim and S. Hilgenfeldt, Phys. Rev. Lett., 2022, 129, 168001.
27. F. Bolton and D. Weaire, Phys. Rev. Lett., 1990, 65, 3449–3451.
28. C. Voivret, F. Radjai, J.-Y. Delenne and M. S. El Youssoufi, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2007, 76, 021301.
29. X. Du and E. R. Weeks, Phys. Rev. E, 2024, 109, 034605.
30. T. Bertrand, R. P. Behringer, B. Chakraborty, C. S. O'Hern and M. D. Shattuck, Phys. Rev. E, 2016, 93, 012901.
31. P. M. Chaikin, T. C. Lubensky and T. A. Witten, Principles of condensed matter physics, Cambridge University Press, Cambridge, 1995, vol. 10.
32. A.-K. Classen, K. I. Anderson, E. Marois and S. Eaton, Dev. Cell, 2005, 9, 805–817.
33. S. Kim, Y. Wang and S. Hilgenfeldt, Phys. Rev. Lett., 2018, 120, 248001.
34. H. Tong and H. Tanaka, Phys. Rev. X, 2018, 8, 011041.
35. H. Tanaka, H. Tong, R. Shi and J. Russo, Nat. Rev. Phys., 2019, 1, 333–348.
36. P. Sampedro Ruiz and R. Ni, J. Chem. Phys., 2020, 153, 174501.
37. S. F. Swallen, K. L. Kearns, M. K. Mapes, Y. S. Kim, R. J. McMahon, M. D. Ediger, T. Wu, L. Yu and S. Satija, Science, 2007, 315, 353–356.
38. Y. Guo, A. Morozov, D. Schneider, J. W. Chung, C. Zhang, M. Waldmann, N. Yao, G. Fytas, C. B. Arnold and R. D. Priestley, Nat. Mater., 2012, 11, 337–343.
39. S. Singh, M. D. Ediger and J. J. de Pablo, Nat. Mater., 2013, 12, 139–144.
40. M. Ozawa, L. Berthier, G. Biroli, A. Rosso and G. Tarjus, Proc. Natl. Acad. Sci. U. S. A., 2018, 115, 6656–6661.
41. W.-T. Yeh, M. Ozawa, K. Miyazaki, T. Kawasaki and L. Berthier, Phys. Rev. Lett., 2020, 124, 225502.
42. T. Yanagishima, J. Russo, R. P. A. Dullens and H. Tanaka, Phys. Rev. Lett., 2021, 127, 215501.
43. F. Ghimenti, L. Berthier and F. van Wijland, arXiv, 2024, preprint, arXiv:2402.06585 DOI:10.48550/arXiv.2402.06585.
44. C. P. Royall, F. Turci, S. Tatsumi, J. Russo and J. Robinson, J. Phys.: Condens. Matter, 2018, 30, 363001.
45. V. Bolton-Lum, R. C. Dennis, P. Morse and E. Corwin, arXiv, 2024, preprint, arXiv:2404.07492 DOI:10.48550/arXiv.2404.07492.
46. V. F. Hagh, S. R. Nagel, A. J. Liu, M. L. Manning and E. I. Corwin, Proc. Natl. Acad. Sci. U. S. A., 2022, 119, e2117622119.
47. C. Anzivino, M. Casiulis, T. Zhang, A. S. Moussa, S. Martiniani and A. Zaccone, J. Chem. Phys., 2023, 158, 044901.
48. E. Bitzek, P. Koskinen, F. Gähler, M. Moseler and P. Gumbsch, Phys. Rev. Lett., 2006, 97, 170201.
49. A. Okabe, B. Boots, K. Sugihara and D. Kendall, Spatial Tessellations: Concepts and Applications of Voronoi Diagrams, Wiley & Sons, 1992.
50. T. S. Grigera and G. Parisi, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2001, 63, 045102.
51. A. Ninarello, L. Berthier and D. Coslovich, Phys. Rev. X, 2017, 7, 021039.
52. C. Brito, E. Lerner and M. Wyart, Phys. Rev. X, 2018, 8, 031050.

---

Footnote

† Electronic supplementary information (ESI) available. See DOI: https://doi.org/10.1039/d4sm00426d

---