Enhanced collective vibrations in granular materials

Abstract

Granular materials are defined as collections of macroscopic dissipative particles. Although these systems are ubiquitous in our lives, the nature and the causes of their non-trivial collective dynamics still remain elusive and have attracted significant interest in non-equilibrium physics. Here, we focus on the vibrational dynamics of granular materials. While the vibrational dynamics of random packings have been examined concerning the jamming transition, previous research has overlooked the role of contact dissipations. We conducted numerical and analytical investigations into the vibrational dynamics of random packings influenced by the normal dissipative force, which is the simplest model for contact dissipations. Our findings reveal that the kinetic energy per mode diverges in the low-frequency range, following the scaling law with the frequency ω_l, indicating that low-frequency modes experience strong excitation and that the equipartition of energy is violated. Additionally, the spatial structure factor of the velocity field displays the scaling law S_v(q) ∝ q⁻² with the wavenumber q, which signifies that the velocity field has an infinitely long range. We demonstrate that these phenomena arise from the effects of weaker damping on softer modes, where the particle displacements parallel to the contacts are minimal in the low-frequency modes, rendering normal dissipation ineffective at damping these modes.

---

1. Introduction

Granular materials, such as sand and flour, are defined as collections of macroscopic dissipative particles. While these systems are ubiquitous in our lives, their dynamics are highly complex, and many questions remain unanswered from a physics perspective.^1,2 Understanding their dynamics is crucial for engineering and other natural sciences as they provide guidelines for handling granular materials in industrial processes and also offer insights into various natural phenomena, such as snow avalanches and the formation of planetary rings.^3–5 A key feature of the dynamics of granular materials is the emergence of collective behavior.⁶ In dilute systems, granular particles in a free cooling state tend to form clusters, and the spatial heterogeneity results in non-trivial energy decay.^7,8 In dense systems, granular particles in sheared or vibrated states exhibit vortex structures in their velocity fields, which have been an important research focus.^9–17 Because granular materials are dissipative, these macroscopic dynamics emerge as non-equilibrium phenomena that respond to external manipulations of the systems.¹⁸

The physics of the jamming transition provides a useful framework for understanding the complex dynamics of dense granular materials. Random packings of athermal, frictionless particles serve as the central model in this area of research.¹⁹ This model achieves mechanical stability at a specific density, known as the jamming transition.^19–26 The jamming transition represents a form of critical phenomenon, where various physical properties exhibit power laws in relation to the proximity to the critical density.^19,23,27 Among these, the vibrational dynamics of random packings have been extensively studied. The vibrational modes of the packings are derived from the diagonalization of the dynamical matrix, which depends solely on the mass, inter-particle interactions, and the configurations of the constituent particles. In crystals, these vibrational modes manifest as plane waves, and the vibrational density of states adhere to the notable Debye law D(ω) ∝ ω^d−1, where d denotes the spatial dimension. In contrast, many vibrational modes of the random packings are spatially disordered, contributing to an excess density of states relative to the Debye law. D(ω) exhibits a plateau down to a characteristic frequency that reveals critical power-law behavior near the jamming transition. Furthermore, the spatially localized vibrational modes, referred to as quasi-localized modes, are known to emerge in the lowest frequency regime.^28–30 Importantly, the linear rheological properties of the packings are dictated by these vibrational modes,³¹ highlighting the significance of understanding their vibrational properties.

In parallel to these numerical and theoretical works on vibrational properties near jamming, the vibrational dynamics of dense granular materials have been studied experimentally. One method involves recording the covariance matrix of grain positions and diagonalizing it to obtain the vibrational modes. This analysis captures collective vibrations in the system.^32,33 Another method records the time correlation function of grain velocities and transforms it into the frequency domain. Through this approach, an increase in low-frequency vibrations near jamming was observed.^34,35 However, the connection between these experimental observations on granular materials and the above-mentioned understanding of vibrational modes of random packings remains unclear.^36,37 This arises because conventional vibrational mode analysis has neglected dissipation, which is present in experimental granular systems. Dissipation can alter the amplitude of vibrational mode excitations or even the vibrational modes themselves. To fully understand the vibrations of granular materials, it is essential to explicitly incorporate dissipation into vibrational mode analysis.^33,38,39

In this work, we study the vibrational dynamics of a model of granular materials. We examine random packings of particles driven by external random forces as the simplest setup. In particular, we demonstrate the non-trivial roles played by dissipation in vibrational dynamics by comparing two distinct systems that share the same setups but the precise implementation of dissipative forces. The first system is the simplest model for granular materials, in which the effects of inelastic collisions are incorporated via normal dissipative forces. In the second system, dissipation is implemented by the Stokes drag force. Although this approach is well-suited to describing colloidal suspensions, we study it as a reference system. We solve the linearized equations of motion for these systems both analytically and numerically. Our findings indicate that the precise form of the dissipation significantly alters the vibrational dynamics. Notably, low-frequency vibrations are markedly amplified and long-wave-length collective motions arise in the case of the normal dissipative forces, whereas the law of equipartition of energy is satisfied under the Stokes drag.

This paper is organized as follows. We begin by introducing our models in Section 2. The simulation methods and physical observables are also presented in this section. Next, in Section 3, we analyze the excitations of vibrational modes using the analytical solutions of the linearized equations of motion. We demonstrate that the low-frequency modes are strongly excited in the case of normal dissipative force. In Section 4, we discuss the spatial properties of the vibrations. In particular, we show that collective motions arise due to the strong excitation of low-frequency modes.

2. Model and methods

2.1. Model

In this section, we describe our granular models. We consider a three-dimensional monodisperse system consisting of N particles with equal masses m. The system obeys the equation of motion:where r_i = (r_ix, r_iy, r_iz)^T is the position of particle i, and T denotes transpose operation. The vectors F^C_i, F^D_i, F^R_i denote the conservative force, dissipative force, and random force acting on particle i, respectively.

In our model, the particles interact through the harmonic potential:⁴⁰

where k is the spring constant, σ is the particle diameter, r_ij = |r_ij| = |r_i − r_j| is the distance between particles i and j, and Θ(x) is the Heaviside step function. Therefore, particles i and j interact only when they overlap r_ij < σ. Then, the conservative force acting on particle i is given by

For the dissipative force, we consider two different models. One is the normal dissipative force model:

∂_i denotes the set of particles in contact with particle i, and η_n is the coefficient of the normal dissipative force. v_i is the velocity of particle i and n_ij = r_ij/r_ij is the unit vector along the contact between particles i and j. Therefore, this dissipative force is proportional to the relative velocity parallel to the contact. This is the simplest model of inelastic collisions of granular materials. The other dissipative force model in our research is the Stokes drag model:

F^D_i = −η_sv_i, (5)

where η_s represents the coefficients of the Stokes drag. Although this is the simplest model of the dissipation of colloidal particles immersed in solvent, we study this model as well, for comparison.

For the external force, we consider the random white noise satisfying the following property:

〈F^R_i(t)F^RT_j(s)〉 = 2Bδ_ijδ(t − s)I₃, (6)

where B is a parameter that determines the strength of the white noise, δ_ij is the Kronecker's delta, δ(t) is the Dirac's delta function, and I₃ is the 3 × 3 identity matrix. 〈⋅〉 represents the average over noise. Note that this work focuses on the random force for its simplicity. In real granular materials, there is a wide variety of energy injections, such as vertical vibrations, shear deformations, and air fluidization. Some of them are characterized by the specific frequency and the resulting dynamics can depend on it. To avoid these complications, we focus on the random force as it has the simplest flat frequency characteristics.

2.2. Linearization of the equation of motion

We consider the model at the volume fraction of ϕ ≃ 0.74 which is higher than the jamming density. Our analyses focus on such a high density system to reveal general features of the vibrations of jammed granular materials without paying too much attention to the criticality of the jamming transition. We focus on a mechanically stable configuration, or in other words, an inherent structure of the model. The position of particle i in the inherent structure is denoted by R_i = (R_ix, R_iy, R_iz)^T. Hereafter, we introduce a compact notation of a 3N-dimensional vector like a = (a₁^T, a₂^T,…, a_N^T)^T. For example, the positions of particles in the inherent structure are denoted by R = (R₁^T, R₂^T,…, R_N^T)^T.

When the random force is sufficiently weak, the model behaves as a solid and the particles move only around their positions in the inherent structure. In this situation, the prestress (or pressure) keeps the particles in contact with each other and the interparticle contacts never be broken by random forces even if the dissipative forces are absent. Hence, we linearize the equation of motion eqn (1) in terms of the displacements of the particles, u≡r − R, to obtain:

where Γ is the 3N × 3N damping matrix and is the 3N × 3N hessian matrix of the potential, which will be given below.

From the definition of the dissipative forces eqn (4) and (5), the damping matrix Γ is composed of the 3 × 3 matrix elements Γ_ij corresponding to the pairs of particles i and j. More specifically, the damping matrix for the normal dissipative force model eqn (4) is given by

where O₃ is the 3 × 3 zero matrix. The damping matrix for the Stokes drag model eqn (5) is given by

Similarly, the hessian matrix is composed of 3 × 3 matrix elements corresponding to the pairs of particles i and j. In general, the element depend on both the first derivative of the potential and the second derivative . Here, we employ a simple approximation , which is called the “unstressed” approximation. In this approximation, all the contacts between particles are replaced by the relaxed springs with the spring constant k. Then the element of the hessian matrix is given by:

In summary, the linearized equation of motion is controlled by the damping matrix Γ and the hessian matrix . These matrices are fully determined by the inherent structure. In the following sections, we solve this equations of motion for the normal dissipative force model (8) and the Stokes drag model (9) to discuss the vibrations of granular materials.

2.3. Observables

In this section, we introduce the observables to characterize the vibrations in our granular model. We first introduce the eigenequations of the hessian matrix:where e_l is the l-th eigenvector and λ_l the l-th eigenvalue (l = 1, 2,…, 3N). Physically, λ_l describes the effective spring constant for the mode e_l. Accordingly, we introduce the characteristic frequency of the l-th mode as For crystals and structural glasses where the dissipation does not play important role, these eigenvectors and eigenvalues carry all the relevant information of the vibrational dynamics of the system. In particular, thermal agitation excites these modes with the same energy, which is known as the law of equipartition of energy. However, in our granular model, the non-trivial dissipation is present, which may break the law of equipartition.

To quantify the excitations of the modes, we focus on the kinetic energy associated with each mode. To this end, we expand the displacement fields using the modes as

Here, C_l = e^T_lu is the displacement along the l-th mode, and also Ċ_l = dC_l/dt is the velocity along the l-th mode. Then the kinetic energy of the l-th mode is given byWe calculate for our granular models to examine the excitation of each mode.

Next, to quantify the spatial structure of the vibrational dynamics, we introduce a spatial correlation function of the velocities of particles:⁴¹

S_v(q) = 〈v_q·v_−q〉. (15)

Here, is the Fourier transform of the velocity field of the system, q is the wave vector, and v_i = du_i/dt is the velocity of the particle i. If S_v(q) → ∞ in the limit of large scale q → 0, it means that the spatial structure expands whole the system.

Finally, to quantify the alignment of the velocities of neighboring particles, we introduce the following parameter:

where N_i is the number of particles in the set ∂_i. This parameter quantifies the tendency of the alignment of velocities of neighboring particles. If Ψ_i = 1, particle i moves in the same direction as its neighboring particles. If Ψ_i = 0, particle i and its neighboring particles move in completely random directions.

2.4. Simulation setup

We first obtained the inherent structure of the model in the same way as ref. 28. The rattler particles which have less than 4 contacting particles were removed iteratively. Starting from the inherent structures, we performed the molecular dynamics simulations of the linearized equation of motion eqn (7). We performed sufficiently long simulations to obtain the steady state of the system, and then performed product runs where we calculated the physical quantities defined in Section 2.3. To calculate the kinetic energies introduced by eqn (14), we diagonalized the hessian matrix to obtain the eigenvectors and the eigenfrequencies. We set the dissipative force parameter to which corresponds to typical restitution coefficient of granular particles ε ≃ 0.7.^42,43 The random force parameter was set to which justifies the linear approximation of the equation of motion. This corresponds to the situation in which the random forces are vanishingly weak but large enough for macroscopic motion to be clearly observed.

To integrate the equation of motion, we utilized the DPD-VV method,⁴⁴ which is an extension of the Velocity-Verlet method to particle systems with dissipation. The time step of the simulation was set to . The models with several different system sizes are considered: 1000 ≤ N ≤ 512 000. These system sizes enable us to observe the spatially heterogeneous vibrational dynamics. For each system size, we prepared 500 different steady-state configurations for N ≤ 125 000 and 250 configurations for N = 256 000, 512 000. The results are all averages over these samples. Note that we used the same inherent structures for both the Stokes drag model and the normal dissipative force model.

3. Kinetic energy on vibrational mode

In our models, the hessian matrix and the damping matrix Γ share the same matrix structure. This enables us to solve the equations of motion of our models. In the following, we calculate the kinetic energies using the analytical solutions of the equations of motion.

3.1. Stokes drag model

We first focus on the Stokes drag model eqn (9). Since Γe_l = η_se_l in this case, we can expand the equation of motion into each mode l (= 1, 2,…, 3N) as:where f_l is the random force for the l-th mode, defined by . This equation of motion can be solved by the variation of constants method (see Appendix A). Using this solution, we obtain the kinetic energy for each mode as: Therefore, all the modes are equally excited in the Stokes drag model. This is reasonable since this model can be seen as colloidal particles in a solvent, which are in the thermal equilibrium state.

To confirm this result, we performed MD simulations for the Stokes drag model (see Section 2.4 for details). The obtained kinetic energy are shown against the characteristic frequency ω_l as the blue symbols in Fig. 1. Clearly, the simulation result agrees with the theoretical result eqn (18) shown by the dashed line.

Fig. 1 The kinetic energy for each mode against the characteristic frequency. Symbols express the simulation results and lines represent the theoretical results eqn (18) and (21). Blue symbols are for the Stokes drag model and red symbols are for the normal dissipative force model. The law of equipartition of energy is violated for the later case.

3.2. Normal dissipative force model

Next, we consider the normal dissipative force model eqn (8). For this case, it is crucial to observe that the hessian matrix eqn (10) and the damping matrix eqn (8) have the same matrix structure and they are proportional . Therefore, the two matrices share the same eigenvectors and the following relation holds: This result means that the effective damping coefficient for the l-th mode is η_nλ_l/k. Since λ_l is the effective spring constant for the l-th mode, eqn (19) means that the damping becomes weaker as the mode is softer. This is natural due to the nature of the normal dissipative force. The hessian matrix eqn (10) describes the energy cost that comes from the relative displacements parallel to the contacts. If λ_l is small, such parallel displacements are small in the mode l. Since the normal dissipation damps the relative motions parallel to the contacts, the effective damping coefficient should be smaller when the parallel displacements are smaller. This effect, that is weaker damping for softer mode, is described by eqn (19).

Using eqn (19), the equation of motion can be decomposed into each mode as:

Compared to the Stokes drag force model eqn (17), the damping coefficient η_s is replaced by η_nλ_l/k. We can solve this equation by the same method, and we obtain the kinetic energy as The law of equipartition of energy is violated, and the low-frequency vibrations are strongly excited in the normal dissipative force model. The factor k/mω_l² comes from the effect of weaker damping for softer mode.

To confirm this result, we also performed MD simulations for the normal dissipative force model. The obtained kinetic energy are shown as the red symbols in Fig. 1. Clearly, the simulation results reproduce the theoretical result eqn (21) well.

4. Spatial structure of the velocity field

The previous section showed that the law of equipartition of energy is violated in the normal dissipative force model and the low-frequency vibrations are strongly excited due to the effect of weaker damping for softer modes. In this section, we show that this effect causes a collective behavior in the velocity field.

First, we use MD simulations to calculate Ψ_i that quantifies the similarity in velocities of the neighboring particles. Fig. 2 shows the snapshots of particles colored by the value of Ψ_i. In the Stokes drag model, Ψ_i of each particle is random and there is no noticeable spatial structure. On the other hand, in the normal dissipative force model, the particles with larger Ψ_i (colored in red) make clusters. This means that the particle tend to move in the same direction as the neighboring particles and the collective behavior emerges in the velocity field.

Fig. 2 Order parameter Ψ_i of the granular systems. Visualization of the spatial structure of the vibrations. Results for N = 16 000 are shown. (a) Stokes drag model and (b) normal dissipative force model. The color of the particles reflects the value of the order parameter, corresponding to the color bar on the right.

We quantify this collective behavior by calculating the spatial correlation function of velocities S_v(q). Fig. 3 shows the results. In the Stokes drag model, S_v(q) is constant for all q, which suggests a fully random structure of the velocity field. This is consistent with the thermal dynamics of this model as shown in the previous section. By contrast, in the normal dissipative force model, the correlation function increases and diverges at q → 0. In particular, the asymptotic behavior S_v(q) ∝ q⁻² appears in the low-q regime, which suggests the infinitely long-range correlation of the velocities of particles.

Fig. 3 Velocity correlation functions obtained from simulations. (a) Stokes drag model and (b) normal dissipative force model.

The asymptotic behavior S_v(q) ∝ q⁻² can be understood in terms of the effect of weaker damping for softer mode. Here, we focus on the fact that the eigenmodes of the “unstressed” hessian matrix can be approximated well by the plane waves in the low-frequency regime.²⁸ In this approximation, the mode l is characterized by the wave vector q and the polarization α = L, T₁, T₂, and we write

where s_q,α is the normalized polarization vector. Then, the Fourier component of the velocity field can be expressed as and the velocity correlation function is obtained as This expression should be valid for the low-q regime where the plane wave approximation works. In the Stokes drag model, eqn (18) leads to 〈Ċ_q,α²〉 = B/mη_s. Therefore, we obtain This result explains the simulation data for the Stokes drag model very well. In the normal dissipative force model, eqn (21) leads to 〈Ċ_q,α²〉 = Bk/m²η_nω_q,α². Within the plane wave approximation, the frequency follows the dispersion relation ω_q,α = c_αq, where c_L and c_T is the longitudinal and transverse sound velocity, respectively. Then, we obtain This result explains the simulation data for the normal dissipative force model where S_v(q) diverges as q⁻² in the low-q regime.

This calculation shows the effect of weaker damping for softer mode causes the infinitely long-range velocity correlation in the normal dissipative force model. This effect causes the strong excitation for the low-frequency modes as 〈Ċ_l²〉 ∝ ω_l⁻². In the low-frequency regime, the eigenmodes are approximated well by the spatially extended plane waves with the dispersion relation ω ∝ q, and then these modes are strongly excited as 〈Ċ_q,α²〉 ∝ q⁻². This leads to S_v(q) ∝ q⁻².

5. Conclusion and future works

In the present work, we examined the vibrational dynamics of random packings of particles agitated by random forces. We considered two types of models for dissipation: the Stokes drag model, which is appropriate for colloidal particles immersed in a solvent, and the normal dissipative force model, which describes dissipation due to contact between particles and is suitable for granular materials. By numerically and analytically solving the linearized equation of motion, we demonstrated that the kinetic energy per mode diverges in the low-frequency limit as in the normal dissipative force model. This sharply contrasts with the Stokes drag model, where the law of equipartition of energy applies effectively. This violation of the equipartition law in the normal dissipative force model results from the weaker damping effects for softer modes. Since the particle displacements parallel to the contacts are minimal in the low-frequency modes, normal dissipation is ineffective in damping these modes. Furthermore, we showed that the spatial structure factor of the velocity field exhibits the scaling law S_v(q) ∝ q⁻² in the normal dissipative force model, indicating that the velocity field is infinitely long-ranged. This also relates to the weaker damping effects for softer modes: low-frequency plane-wave modes are strongly excited, resulting in the emergence of long-range velocity correlations. These results establish that contact dissipation significantly alters the vibrational dynamics of random packings.

Our results suggest that the vibrations in granular materials are more collective than those in the thermal glass system. As mentioned in the introduction section, the vibrational dynamics of granular materials were studied experimentally, and the collective vibrations along with an increase in low-frequency modes were observed.^32–35 It would be interesting to quantitatively analyze these experimental results based on those obtained in this work. However, for this purpose, our model remains too simplistic. The present study considered only the normal dissipative force, while in real experiments, there must be situations where both the normal dissipative force and the Stokes drag are involved. Besides, the tangential dissipative force should also be present in granular materials. To account for this in our model, we need to treat the case where the damping matrix Γ and the hessian matrix do not share the same matrix structure. In this case, the vibrational modes are no longer the eigenvectors of the hessian matrix . Furthermore, the dynamics of real-world granular materials are affected by the friction force arising from the contacts between particles. Understanding these effects on vibrational dynamics will be the focus of future studies on granular materials.

Data availability

The data that support the findings of this study are available from the corresponding author upon reasonable request. Source data are provided with this paper.

Conflicts of interest

There are no conflicts of interest to declare.

Appendix

A. Derivation of the kinetic energy

In both the Stokes drag model and the normal dissipative force model, the linearized equations of motion for each mode have the following formwhere m, η, λ are constants and f is the random noise satisfying 〈f²〉 = 2B. This Langevin equation can be solved using the standard method of variation of constants.

We first transform this equation to the matrix form:

where a(t) = (C(t), Ċ(t))^T, F(t) = (0, f(t)/m)^T, andThe covariance matrix of F is given by Since this is the first order differential equation, the formal solution can be written as Because the kinetic energy is given by 〈Ċ²〉, we need to evaluate the covariance matrix of a. This can be calculated as
where † denotes the conjugate transpose. By diagonalizing Θ, the matrix e^tΘ and e^tΘ† can be calculated in a lengthy but straightforward way. The obtained result is In the case of the Stokes drag model, the parameter is given by η = η_s and then In the case of the normal dissipative force model, the parameter is given by η = η_nλ_l/k and then

Acknowledgements

This work was supported by JSPS KAKENHI Grant Numbers JP20H01868, JP20H00128, JP22K03543, JP23H04495, and JP24H00192. In this research work, we used the mdx: a platform for building data-empowered society.

References

1. H. M. Jaeger and S. R. Nagel, Science, 1992, 255, 1523.
2. P. G. de Gennes, Rev. Mod. Phys., 1999, 71, S374.
3. N. Borderies, P. Goldreich and S. Tremaine, Icarus, 1985, 63, 406.
4. H. M. Jaeger, S. R. Nagel and R. P. Behringer, Rev. Mod. Phys., 1996, 68, 1259.
5. A. Daerr and S. Douady, Nature, 1999, 399, 241.
6. I. S. Aranson and L. S. Tsimring, Rev. Mod. Phys., 2006, 78, 641.
7. I. Goldhirsch and G. Zanetti, Phys. Rev. Lett., 1993, 70, 1619.
8. S. N. Pathak, Z. Jabeen, D. Das and R. Rajesh, Phys. Rev. Lett., 2014, 112, 038001.
9. J. A. C. Gallas, H. J. Herrmann and S. Sokołowski, Phys. Rev. Lett., 1992, 69, 1371.
10. S. L. Conway, T. Shinbrot and B. J. Glasser, Nature, 2004, 431, 433.
11. P. Jop, Y. Forterre and O. Pouliquen, Nature, 2006, 441, 727.
12. S. Abedi, A. L. Rechenmacher and A. D. Orlando, Granular Matter, 2012, 14, 695.
13. A. Fall, G. Ovarlez, D. Hautemayou, C. Mézière, J.-N. Roux and F. Chevoir, J. Rheol., 2015, 59, 1065.
14. K. P. Krishnaraj and P. R. Nott, Nat. Commun., 2016, 7, 10630.
15. C. R. K. Windows-Yule, New J. Phys., 2016, 18, 033005.
16. P. Kharel and P. Rognon, Phys. Rev. Lett., 2017, 119, 178001.
17. J. Kozicki and J. Tejchman, Granular Matter, 2018, 20, 48.
18. K. Ahmad and I. J. Smalley, Powder Technol., 1973, 8, 69.
19. C. S. O'Hern, L. E. Silbert, A. J. Liu and S. R. Nagel, Phys. Rev. E, 2003, 68, 011306.
20. H. M. Princen and A. D. Kiss, J. Colloid Interface Sci., 1986, 112, 427.
21. F. Bolton and D. Weaire, Phys. Rev. Lett., 1990, 65, 3449.
22. T. S. Majmudar, M. Sperl, S. Luding and R. P. Behringer, Phys. Rev. Lett., 2007, 98, 058001.
23. M. van Hecke, J. Phys.: Condens. Matter, 2009, 22, 033101.
24. M. Clusel, E. I. Corwin, A. O. N. Siemens and J. Brujić, Nature, 2009, 460, 611.
25. E. Brown, N. Rodenberg, J. Amend, A. Mozeika, E. Steltz, M. R. Zakin, H. Lipson and H. M. Jaeger, Proc. Natl. Acad. Sci. U. S. A., 2010, 107, 18809.
26. R. P. Behringer and B. Chakraborty, Rep. Prog. Phys., 2018, 82, 012601.
27. C. S. O'Hern, S. A. Langer, A. J. Liu and S. R. Nagel, Phys. Rev. Lett., 2002, 88, 075507.
28. H. Mizuno, H. Shiba and A. Ikeda, Proc. Natl. Acad. Sci. U. S. A., 2017, 114, E9767.
29. G. Kapteijns, E. Bouchbinder and E. Lerner, Phys. Rev. Lett., 2018, 121, 055501.
30. M. Shimada, H. Mizuno and A. Ikeda, Phys. Rev. E, 2018, 97, 022609.
31. A. Lemaître and C. Maloney, J. Stat. Phys., 2006, 123, 415.
32. C. Brito, O. Dauchot, G. Biroli and J.-P. Bouchaud, Soft Matter, 2010, 6, 3013.
33. S. Henkes, C. Brito and O. Dauchot, Soft Matter, 2012, 8, 6092.
34. E. T. Owens and K. E. Daniels, Soft Matter, 2013, 9, 1214.
35. S. A. Blue, S. C. Wright and E. T. Owens, Phys. Rev. E, 2024, 110, 014902.
36. E. Somfai, M. van Hecke, W. G. Ellenbroek, K. Shundyak and W. van Saarloos, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2007, 75, 020301.
37. K. Saitoh, R. K. Shrivastava and S. Luding, Phys. Rev. E, 2019, 99, 012906.
38. J. Chattoraj, O. Gendelman, M. Pica Ciamarra and I. Procaccia, Phys. Rev. Lett., 2019, 123, 098003.
39. D. Ishima, K. Saitoh, M. Otsuki and H. Hayakawa, Phys. Rev. E, 2023, 107, 034904.
40. D. J. Durian, Phys. Rev. Lett., 1995, 75, 4780.
41. S. Henkes, K. Kostanjevec, J. M. Collinson, R. Sknepnek and E. Bertin, Nat. Commun., 2020, 11, 1405.
42. J. M. Lifshitz and H. Kolsky, J. Mech. Phys. Solids, 1964, 12, 35.
43. C. S. Sandeep, K. Senetakis, D. Cheung, C. E. Choi, Y. Wang, M. R. Coop and C. W. W. Ng, Can. Geotech. J., 2021, 58, 35.
44. G. Besold, I. Vattulainen, M. Karttunen and J. M. Polson, Phys. Rev. E, 2000, 62, R7611.

---