# Universal relaxation dynamics of sphere packings below jamming

Atsushi Ikeda, Takeshi Kawasaki, Ludovic Berthier, Kuniyasu Saitoh, Takahiro Hatano

► To cite this version:

Atsushi Ikeda, Takeshi Kawasaki, Ludovic Berthier, Kuniyasu Saitoh, Takahiro Hatano. Universal relaxation dynamics of sphere packings below jamming. Physical Review Letters, 2020, pp.058001. 10.1103/PhysRevLett.124.058001. hal-02591059

HAL Id: hal-02591059

https://hal.science/hal-02591059

Submitted on 2 Jun 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Universal relaxation dynamics of sphere packings below jamming

Atsushi Ikeda,¹,² Takeshi Kawasaki,³ Ludovic Berthier,⁴ Kuniyasu Saitoh,⁵ and Takahiro Hatano⁶

¹Graduate School of Arts and Sciences, The University of Tokyo, Tokyo 153-8902, Japan

²Research Center for Complex Systems Biology, Universal Biology Institute,
University of Tokyo, Komaba, Tokyo 153-8902, Japan

³Department of Physics, Nagoya University, Nagoya 464-8602, Japan

⁴Laboratoire Charles Coulomb (L2C), University of Montpellier, CNRS, Montpellier, France

⁵Research Alliance Center for Mathematical Sciences & WPI-Advanced
Institute for Materials Research, Tohoku University, Sendai 980-8577, Japan

⁶Department of Earth and Space Science, Osaka University, 560-0043 Osaka, Japan

(Dated: April 25, 2019)

We show that non-Brownian suspensions of repulsive spheres below jamming display a slow re-
laxational dynamics with a characteristic time scale that diverges at jamming. This slow time scale
is fully encoded in the structure of the unjammed packing and can be readily measured via the vi-
brational density of states. We show that the corresponding dynamic critical exponent is the same
for randomly generated and sheared packings. Our results show that a wide variety of physical
situations, from suspension rheology to algorithmic studies of the jamming transition are controlled
by a unique diverging timescale, with a universal critical exponent.

*Introduction.*— Disordered particle packings solidify when they are compressed [1–4]. In the absence of thermal motion [5–8], as in non-Brownian suspensions [9], emulsions of large droplets [10], foams [11], and granular materials [12], this phenomenon is a jamming transition [4]. Athermal frictionless spheres [13, 14] represent an idealized model to study the jamming criticality [15, 16]. At the critical packing fraction ϕ = ϕ<sub>J</sub>, the system is isostatic [14], the contact number is z = 2d where d is the spatial dimension, and the system is mechanically marginally stable [17]. In the jammed phase ϕ > ϕ<sub>J</sub>, the mechanical response [14, 18–22] and low-frequency vibrational modes [23–26] exhibit algebraic critical behaviour. Theoretical developments [27, 28], effective medium theory [29, 30], and an exact analytic solution in the d → ∞ mean-field limit [31–33] provide a good understanding of the jammed phase, including numerical values of critical exponents.

By contrast, current understanding of the unjammed
phase $\varphi < \varphi_J$ is more limited, in particular regard-
ing dynamics. An important physical observable in this
regime is the diverging shear viscosity of non-Brownian
suspensions, $\eta \propto (\varphi_J - \varphi)^{-\nu}$ [9, 34–36]. Even for ide-
alised frictionless spheres [34, 37], an exact mean-field
prediction for $\nu$ is not available, but several distinct
predictions exist: $\nu \approx 2.83$ in [39–41] and $\nu = 2$ in
[42]. Measured numerical values are spread in the range
$\nu \in [2, 4]$ [9, 34, 37, 38, 42–46]. We argue that physi-
cal situations where dynamics of unjammed packings is
relevant go well beyond steady shear rheology. In emul-
sions, nearly jammed packings are prepared by dispersing
the droplets in a solvent and performing centrifugation or
creaming under gravity [10, 47, 48], during which droplets
relax. Foams near jamming are also prepared by foam-
ing solutions and injecting the bubbles into the sample
chamber [49–51], during which bubbles can relax. These

materials can also be probed rheologically by performing
step stress or strain experiments, monitoring the mechan-
ical response. In simulations, packings are often prepared
by relaxing a random configuration to mechanical equi-
librium [14]. In the unjammed phase $\varphi < \varphi_J$, relaxation
stops when particles are just touching, and $z < 2d$ [52].
Critical slowing down of relaxation algorithms near jam-
ming was reported [53].

Several pieces of information are known from the anal-
ysis of computer models. The diverging shear viscosity is
not associated to a slowing down of particle diffusion, as
would be the case in equilibrium due to Green-Kubo rela-
tions [54]. In fact, particle motion accelerates near jam-
ming [55, 56]. Instead, it was shown that the relaxation
dynamics after sudden shear cessation is characterised
by a time scale τ that diverges as ϕ_J is approached [45],
with the observation that η ∝ τ ∝ (Δz)⁻³.⁷ in d = 3,
where Δz = 2d − z [46]. Strikingly, ν seems to de-
pend on dimensionality [46], unlike other critical expo-
nent of jamming [14] and theoretical predictions [39–42].
In Ref. [37], a dynamical matrix was constructed for un-
jammed hard spheres packings under shear. Its spectrum
displays, in addition to more traditional signatures of soft
modes near jamming, an isolated low-frequency mode at
some ω_min ≈ Δz¹.⁵ for d = 3. This mode has no ana-
log above jamming, and it is directly connected with the
shear viscosity, η ≈ ω_min⁻² [37]. A theoretical argument
to explain this mode was also proposed, which predicts
ω_min ≈ Δz¹.⁷ [39–41]. Finally, Refs. [13, 57] studied the
relaxation dynamics of random packings of soft particles
below jamming after applying a step shear strain. The
ensuing relaxation dynamics slows down algebraically as
φ → φ<sub>J</sub>, with τ ∝ (φ<sub>J</sub> − φ)⁻³.³ for d = 3 [57].

Here, we suggest that all these forms of relaxation dy-
namics below jamming exhibit a critical slowing down
approaching jamming characterised by a unique diverg-
ing time scale $\tau$. Independently of the packing preparation, $\tau$ is fully encoded in the packing structure and can be readily accessed via the standard density of states, so that $\tau = 1/(2\omega_{\min}^2) \propto \eta$. This diverging time scale is characterised by a universal exponent which depends neither on the form of the repulsion potential nor on the preparation protocol, $\omega_{\min} \propto \Delta z^\alpha$, with $\alpha \approx 1.6$ in $d=3$.

**Model and methods.**— We consider a three dimensional 50:50 binary mixture of large and small soft spheres with diameter ratio 1.4. The potential energy of the system is given by $E = \sum_{i>j} \frac{\epsilon}{2}(1-r_{ij}/\sigma_{ij})^2\Theta(\sigma_{ij}-r_{ij})$, where $\Theta(x)$ is the Heaviside step function, $r_{ij}$ is the distance between particles $i$ and $j$, $\sigma_{ij} = (\sigma_i + \sigma_j)/2$ and $\sigma_i$ is the particle diameter of particle $i$. We study the athermal overdamped dynamics of the model,

$$\xi \frac{d\vec{r}_i}{dt} = -\frac{\partial E}{\partial \vec{r}_i}, \quad (1)$$

where $\vec{r}_i$ is the position of particle $i$, and $\xi$ is the viscous damping. We focus on the unjammed phase $\varphi < \varphi_J$, where the system relaxes at long times into a zero-energy configuration. We use the small particle diameter $\sigma$, $\epsilon$ and $\xi\sigma^2/\epsilon$ as units of length, energy and time.

We consider two types of initial configurations. First, we start from random initial configurations. We place $N$ spheres randomly in the box at packing fraction $\varphi$, and solve Eq. (1) until $E/N \le 10^{-18}$. We prepare 1000 independent initial configurations for each packing fraction. Second, we use steady-shear initial configurations. We perform overdamped simple shear simulations of the same model at finite shear rate $\dot{\gamma}$ to reach steady state. For these simulations, we use the algorithm developed in Ref. [44], where the pressure $p$ and the shear rate $\dot{\gamma}$ are fixed and $\varphi$ fluctuates during the simulation. We sample configurations from the simulations at $p = 10^{-5}$ and $\dot{\gamma} \ge 10^{-8}$, where the shear viscosity is Newtonian. After steady state is reached, we suddenly stop the shear and again solve Eq. (1) until $E/N \le 10^{-18}$. For each shear rate, which is specified by the average density $\langle \varphi \rangle$ in this paper, we obtain 50 independent initial configurations. We use standard periodic boundary conditions for the random case, and Lees-Edwards boundary conditions for the sheared case [58]. The number of particles is $N=3000$, unless otherwise noted.

**Relaxation dynamics.**— We first study the relaxation dynamics from random initial configurations. Figure 1(a) shows typical time evolutions of the potential energy $E(t)$ at various packing fractions. It displays a simple exponential relaxation at low density, but it exhibits two time regimes at higher density with a power-law decay $t^{-1}$ at short times, followed by an exponential decay $e^{-t/\tau}$ at long times. The final relaxation time $\tau$ increases rapidly as $\varphi_J$ is approached. We simultaneously measure the time evolution of the contact number $z(t)$ in Fig. 1(b) [60]. It decreases with time for $t < \tau$, which implies that contact breaking and restructuring occurs

<figure><img src="image_2.png" /><figcaption>FIG. 1: (a,b) Time evolution of the potential energy density $E/N$ and contact number $z$ starting from random initial configurations for various packing fraction $\varphi$. Symbols and lines are typical results for a single configuration with $N=3000$ and $10^6$, respectively. (c,d) The same as (a,b) with initial configuration drawn from the sheared steady state with $N=3000$. Dashed line indicates $E(t) \sim t^{-1}$ in (a) and $z=6$ in (b,d).</figcaption></figure>

on this time scale. On the other hand, $z(t)$ becomes nearly constant for $t > \tau$, which suggests that the final exponential decay of the energy takes place with a fixed contact network. Note that $z(t)$ approaches a finite value as $t \to \infty$, and particles are still touching in the final configurations, even though $\varphi < \varphi_J$. This is natural because in the overdamped dynamics, the particles' overlaps approach $0^+$ as $t \to \infty$, i.e. particles are just 'kissing' when they stop moving. The contact number $z(t \to \infty)$ continuously approaches the isostatic value $z=6$ as $\varphi \to \varphi_J$. The time fluctuations of the data seen at high density for $N=3000$ become very weak in much larger systems, $N=10^6$. In the following, we focus on the results for $N=3000$.

We also study the relaxation dynamics from steady-shear initial configurations, see Figs. 1(c,d). In contrast to the random case, the energy decays exponentially $E(t) \sim e^{-t/\tau}$ over the entire time domain, consistent with Ref. [45]. The contact number changes very little during the relaxation. This suggests the steady-shear initial configurations (i.e. instantaneous configurations in the steady-state shear flow) have almost the same contact networks as those of the final configurations, and they undergo exponential relaxation from the beginning.

**Vibrational density of states.**— We expect that the final exponential decay of the energy towards zero can be controlled by the potential energy landscape around the final configurations. Although the potential energy is nearly zero ($E/N \le 10^{-18}$) in the final configuration, particles still have well-defined contacts, and thus the
<figure><img src="image_3.png" /><figcaption>FIG. 2: The vibrational density of states $D(\omega)$ for final configurations obtained after relaxation from (a) random and (b) steady-sheared initial configuration. The average over initial configurations with the same density $\varphi$ for (a) and the same shear rate (specified by the average density $\langle\varphi\rangle$) for (b) broadens the isolated peak at $\omega_{\min}$.</figcaption></figure>

second derivative of the potential and the Hessian exist. Thus, we can calculate the vibrational density of states (vDOS) of the final configurations. After removing rattlers [61], we diagonalise the Hessian of each final configuration to obtain the eigenvalues $\lambda_k$ and eigenvectors $\vec{e}_k = [\vec{e}_{1,k}, \vec{e}_{2,k}, \dots, \vec{e}_{N,k}]$, where $k$ is the mode index. Since the system is hypostatic, we find $N\Delta z/2$ zero modes for each final configuration, as expected from the Maxwell criterion [17]. After removing these zero modes, the vDOSs is calculated as $D(\omega) = \frac{1}{N} \sum_k \delta(\omega - \sqrt{\lambda_k})$.

In Fig. 2 we show the vDOS for both random and sheared configurations at various densities, averaged over independent final configurations. In both cases, the vDOS has three frequency regimes: a sharp edge at large frequency, a plateau at intermediate frequencies, and an isolated peak at low frequency. With increasing $\varphi$, the plateau and the isolated peak shift to lower frequency and become better separated. The high frequency edge and the intermediate plateau are well-known features of the vDOS seen above jamming, with an onset frequency for the plateau which scales linearly with the excess contact number [27]. Recently, the effective medium theory, which correctly predicts the vDOS of jammed configurations [29], was extended to the unjammed region $\varphi < \varphi_J$ and predicted $\omega^* \propto \Delta z$ [59]. Our data agree with this prediction. A novel feature of the reported vDOS is the isolated low-frequency mode. The corresponding peak appears broad in Fig. 2 because we performed an ensemble average, but we find a unique isolated mode (referred to as ILM) at frequency $\omega_{\min}$ in each individual configuration. This ILM was previously observed for sheared non-Brownian hard spheres [37]. The ILM we find should be similar to the one detected in sheared hard spheres, and our results show that it can be accessed from the standard analysis of vDOS of soft particles. To our

<figure><img src="image_4.png" /><figcaption>FIG. 3: Parametric plot of the relaxation time $\tau$ in the exponential regime of the relaxation against $1/(2\omega_{\min}^2)$ measured for each individual configuration. The dashed line indicates that $\tau = 1/(2\omega_{\min}^2)$ is exactly obeyed for both random (filled symbols) and sheared (open symbols) configurations.</figcaption></figure>

knowledge, the presence of the ILM has not been observed for random configurations before, and our results suggest that it will always be present in unjammed packings in a wide variety of preparation protocols.

We now connect the ILM to the relaxation time $\tau$ of the relaxation dynamics. To this end, we parametrically plot $\tau$ against $1/2\omega_{\min}^2$ for each individual configuration (the energy relaxation along a mode $k$ follows $e^{-2\lambda_k t}$) in Fig. 3, which confirms that $\tau = 1/2\omega_{\min}^2$ is precisely obeyed for configurations at different packing fractions and from different preparation protocol. We conclude that the final exponential relaxation of the energy in both random and sheared cases occurs along the ILM. We also find that the ILM dominates very strongly over all other modes for the sheared configurations [62], which presumably explains the absence of a power-law relaxation regime in that case.

**Critical behavior.**— We now focus on the critical behavior of $\omega_{\min}$ (and thus of $\tau$) near jamming. To avoid the unwanted effect of sample to sample fluctuations of $\varphi_J$ [14], we plot parametrically in Fig. 4 $\omega_{\min}$ against $\Delta z$, measured in each individual configuration. Clearly, the two data sets fully overlap without any rescaling over the entire range. Thus, the link between $\omega_{\min}$ and $\Delta z$ is universal, independent of the preparation protocol. This observation implies that the relaxation dynamics from a wide variety of initial configurations below jamming (quantified by $\tau$) is universal, and is directly linked to the shear viscosity $\eta$ of the suspension. We fit the data in Fig. 4 close to the jamming transition ($\Delta z < 0.3$) with a power law $\omega_{\min} \approx \Delta z^{\alpha}$, and find that $\alpha \approx 1.6$ works well. The measured exponent would decrease if we include data for larger $\Delta z$ in the fit. Of course, $\alpha$ could
<figure><img src="image_5.png" /><figcaption>FIG. 4: Parametric plot of the frequency of the slow mode, $\omega_{\min}$, against the distance to isostaticity $\Delta z = 6-z$ for each individual configuration. The dashed line indicates $\omega_{\min} \propto \Delta z^{1.6}$, valid at low frequency for both random (filled symbols) and sheared (open symbols) configurations.</figcaption></figure>

increase a little further if data for even smaller $\Delta z$ could be included. A more precise numerical determination of $\alpha$ is a worthwhile goal for future work, and the vDOS approach used here seems the easiest route for this task.

**Nature of isolated low-frequency mode.**— Finally, we show that the ILM is qualitatively different from any other mode in the spectrum. Let us consider the evolution of the potential energy obtained by deforming the system along each eigenvector, $\vec{r}_i \equiv \vec{r}_i^0 + x\vec{e}_{i,k}$, where $\vec{r}_i^0$ is the position of particle $i$ in a relaxed configuration, and $x$ is the amplitude of the deformation. We can then track $E = E(x)$ for each mode. We find that $E(x)$ has a parabola shape for all the modes except for the ILM. Fig. 5 shows that $E(x)$ for the ILM is finite for $x < 0$, but zero for $x > 0$. Namely, the potential energy along the ILM is a one-sided harmonic potential. In addition, the contact number jumps discontinuously from $z(x > 0) = 0$ to a finite value for $x < 0$, which is again different for other modes. This one-sided potential energy landscape connects the mechanical vacuum of unjammed configurations with zero potential energy to the slightly compressed sphere packings obtained when the system is undergoing its relaxational dynamics. The curvature of the potential is given by $\omega_{\min}^2$, which thus controls the long-time relaxation dynamics. The direct connection to the shear viscosity stems from the numerical observation that under simple shear, particle displacements are essentially occurring along the ILM. It is remarkable that a single mode captures essentially the entire relaxational physics, rather than a more complex spectrum associated to more complex time dependences as is typically the case in disordered materials [3].

**Conclusion.**— We studied the relaxation dynamics of soft particles in the unjammed phase. We considered

<figure><img src="image_6.png" /><figcaption>FIG. 5: The one-sided potential energy landscape along the mode with the lowest frequency $\omega_{\min}$ measured for a single configuration at $\varphi = 0.64$. It follows the harmonic behavior $E = \frac{1}{2}\omega_{\min}^2 x^2$ (indicated by the dashed line) for $x < 0$. The discontinuous change of the contact number at $x = 0$ is also shown.</figcaption></figure>

two types of initial configurations, random and steady sheared. Despite exponential differences at short time, the relaxation is always exponential at large times, $E \sim \exp(-t/\tau)$. All final configurations display an isolated low-frequency mode at frequency $\omega_{\min}$, and the relation $\tau = 1/(2\omega_{\min}^2)$ is accurately obeyed for each individual configuration, independently of their preparation protocol. The final relaxation process thus occurs along the isolated low-frequency mode, which represents the softest mode that can simultaneously and collectively remove all overlaps from the final packing. Our main quantitative result is that the late stage of the relaxation dynamics is universal, $\omega_{\min} \propto \Delta z^{\alpha}$, with $\alpha \approx 1.6$ in $d=3$.

Our results, on the one hand, establish the universality of the relaxation of athermal particles near the jamming transition. Because the relaxation from two extreme types of initial configurations exhibit the same slowing down, we expect the same in many other physical situations, such as creaming emulsions and injecting foams. It would be interesting to study other protocols. We are currently studying dynamics after a step shear strain [57] more quantitatively. It was reported that during the relaxation dynamics after a step strain [57], the stress exhibits a power law decay $\propto t^{-1/2}$ at short times, followed by an exponential decay. This is similar to our findings with random initial conditions, since stress (and pressure) scale as $E^{1/2}$ for harmonic spheres. On the other hand, our results open a new way to study the viscosity divergence at the jamming transition without applying shear. The exact mean-field solution was recently derived for several properties of jammed configuration, but this approach is still unable to predict the exponents $\alpha$ and $\nu$ discussed in our work. Our results suggest that
these critical exponents can be theoretically analysed by studying the simpler problem of the quench dynamics in the unjammed phase without applying a shear flow.

We thank Hideyuki Mizuno, Masanari Shimada, Matthieu Wyart, and Douglas Durian for useful discussions and comments. This work was supported by KAKENHI grants No. 16H04034, 16H04025, 16H06018, 16H06478, 17H04853, 18H01188, 18H05225, and 18K13464, and a Grant from the Simons Foundation (# 454933, L. B.). Some calculations were performed using the Research Center for Computational Science, Okazaki, Japan.

[1] R. G. Larson, *The Structure and Rheology of Complex Fluids* (Oxford University Press, New York, 1999).

[2] P. Coussot, *Rheometry of Pastes, Suspensions, and Granular Materials* (Wiley, New York, 2005).

[3] L. Berthier, and G. Biroli, Rev. Mod. Phys. **83**, 587 (2011).

[4] A. J. Liu, and S. R. Nagel, Nature **396**, 21 (1998).

[5] A. Ikeda, L. Berthier, and P. Sollich, Phys. Rev. Lett. **109**, 018301 (2012).

[6] P. Olsson, and S. Teitel, Phys. Rev. E **88**, 010301 (2013).

[7] A. Ikeda, L. Berthier, and P. Sollich, Soft Matter **9** 7669 (2013).

[8] M. Dinkgreve, M. A. J. Michels, T. G. Mason, and D. Bonn, Phys. Rev. Lett. 121, 228001 (2018).

[9] F. Boyer, E. Guazzelli, and O. Pouliquen, Phys. Rev. Lett. **107**, 188301 (2011).

[10] T. G. Mason, J. Bibette, D. A. Weitz, Phys. Rev. Lett. **75**, 2051 (1995).

[11] F. Bolton, and D. Weaire, Phys. Rev. Lett. **65**, 3449 (1990).

[12] H. M. Jaeger, S. R. Nagel, and R. P. Behringer, Rev. Mod. Phys. **68**, 1259 (1996).

[13] D. J. Durian, Phys. Rev. Lett. **75**, 4780 (1995); Phys. Rev. E **55**, 1739 (1997).

[14] C. S. O'Hern, L. E. Silbert, A. J. Liu, and S. R. Nagel, Phys. Rev. E **68**, 011306 (2003).

[15] A. J. Liu, S. R. Nagel, W. van Saarloos, and M. Wyart, in *Dynamical Heterogeneities in Glasses, Colloids and Granular Materials*, edited by L. Berthier, G. Biroli, J.-P. Bouchaud, L. Cipelletti, and W. van Saarloos (Oxford University Press, Oxford, 2011).

[16] M. van Hecke, J. Phys. Condens. Matter **22**, 033101 (2010).

[17] S. Alexander, Phys. Rep. **296**, 65 (1998).

[18] W. G. Ellenbroek, E. Somfai, M. van Hecke, W. van Saarloos, Phys. Rev. Lett. **97**, 258001 (2006).

[19] B. P. Tighe, Phys. Rev. Lett. **107**, 158303 (2011).

[20] E. Lerner, E. DeGiuli, G. During, and M. Wyart, Soft Matter, **10**, 5085 (2014).

[21] K. Karimi, and C. E. Maloney, Phys. Rev. E **92**, 022208 (2015).

[22] H. Mizuno, L. E. Silbert, and M. Sperl, Phys. Rev. Lett. **116**, 068302 (2016).

[23] L. E. Silbert, A. J. Liu, and S. R. Nagel, Phys. Rev. Lett. **95**, 098301 (2005).

[24] A. Ikeda, L. Berthier, and G. Biroli, J. Chem. Phys. **138**, 12A507 (2013).

[25] P. Charbonneau, E. I. Corwin, G. Parisi, A. Poncet, and F. Zamponi, Phys. Rev. Lett. **117**, 045503 (2016).

[26] H. Mizuno, H. Shiba, and A. Ikeda, Proc. Natl. Acad. Sci. U. S. A. **114**, E9767 (2017).

[27] M. Wyart, L. E. Silbert, S. R. Nagel, and T. A. Witten, Phys. Rev. E **72**, 051306 (2005).

[28] L. Yan, E. DeGiuli, and M. Wyart, EPL (Europhysics Letters), **114**, 26003 (2016).

[29] M. Wyart, EPL (Europhysics Letters) **89**, 64001 (2010).

[30] E. DeGiuli, A. Laversanne-Finot, G. Düring, E. Lerner, and M. Wyart, Soft Matter, **10**, 5628 (2014).

[31] G. Parisi, F. Zamponi, Rev. Mod. Phys. **82**, 789 (2010).

[32] P. Charbonneau, J. Kurchan, G. Parisi, P. Urbani, and F. Zamponi, J. Stat. Mec. **2014**, P10009 (2014).

[33] S. Franz, G. Parisi, P. Urbani, and F. Zamponi, Proc. Natl. Acad. Sci. U. S. A. **112**, 14539 (2015).

[34] P. Olsson and S. Teitel, Phys. Rev. Lett. **99**, 178001 (2007).

[35] Y. Forterre, and O. Pouliquen, Annu. Rev. Fluid Mech. **40**, 1 (2008).

[36] T. Hatano, J. Phys. Soc. Jpn **77**, 123002 (2008).

[37] E. Lerner, G. Düring, and M. Wyart, Proc. Natl. Acad. Sci. U. S. A., **109**, 4798 (2012).

[38] 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. **105**, 175701 (2010).

[39] E. Lerner, G. Düring, and M. Wyart, EPL (Europhysics Letters) **99**, 58003 (2012).

[40] G. Düring, E. Lerner, and M. Wyart, Phys. Rev. E **89**, 022305 (2014).

[41] E. DeGiuli, G. Düring, E. Lerner, and M. Wyart Phys. Rev. E **91**, 062206 (2015).

[42] K. Suzuki, H. Hayakawa, Phys. Rev. Lett. **115**, 098001 (2015).

[43] B. Andreotti, J.-L. Barrat, and C. Heussinger, Phys. Rev. Lett. **109**, 105901 (2012).

[44] T. Kawasaki, D. Coslovich, A. Ikeda, and L. Berthier, Phys. Rev. E **91**, 012203 (2015).

[45] P. Olsson, Phys. Rev. E **91** 062209 (2015).

[46] P. Olsson, Phys. Rev. Lett. **122**, 108003 (2019).

[47] J. Brujic, S. F. Edwards, D. V. Grinev, I. Hopkinson, D. Brujic, and H. A. Makse, Faraday Discuss. **123**, 207 (2003).

[48] I. Jorjadze, L.-L. Pontani, and J. Brujic, Phys. Rev. Lett. **110**, 048302 (2013).

[49] B. Hertzhzt, S. Kakadjian, and M. Moan, Colloid Surf. A **263**, 153 (2005).

[50] G. Katgert, M. E. Mobius, M. van Hecke, Phys. Rev. Lett. **101**, 058301 (2008).

[51] R. Kurita, Y. Furuta, N. Yanagisawa, N. Oikawa, Phys. Rev. E, **95**, 062613 (2017).

[52] P. K. Morse and E. I. Corwin, Phys. Rev. Lett. **112**, 115701 (2014).

[53] S. Atkinson, G. Zhang, A. B. Hopkins, and S. Torquato, Phys. Rev. E **94**, 012902 (2016).

[54] J. P. Hansen and I. R. McDonald, *Theory of Simple Liquids*, (Elsevier, Amsterdam, 1986).

[55] P. Olsson, Phys. Rev. E **81**, 040301 (2010).

[56] C. Heussinger, L. Berthier, and J.-L. Barrat, EPL (Europhysics Letters) **90**, 20005 (2010).

[57] T. Hatano, Phys. Rev. E **79**, 050301 (2009).

[58] M. Allen and D. Tildesley, *Computer Simulation of Liquids* (Oxford University Press, Oxford, 1987).
[59] G. Düring, E. Lerner, and M. Wyart, Soft Matter, **9**, 146 (2013).

[60] Time evolution of the contact numbers in Figs. 1(b,d) are calculated without removing rattlers. The contact numbers shown in Figs. 4, 5, where we performed quantitative analysis, are obtained after removing rattlers.

[61] We remove the particles whose contact numbers are less than 3. We confirmed that, with this choice, the rattler

removal does not affect the vDOS. Note that the vDOS changes if we remove the particles whose contact numbers are less than 4.

[62] We calculate the overlap between the overall displacements of particles during the relaxation and the eigenvectors of the vibrational modes. We find that the ILM has exceptionally large overlap among all the other modes.