Experimental observation of gapped shear waves and liquid-like to gas-like dynamical crossover in active granular matter
Introduction
Collective modes are a direct macroscopic manifestation of coherent atomic motion and have a pivotal role in determining the thermodynamic, mechanical, and transport properties of physical systems. Phonons, collective lattice vibrations in solids, constitute an emblematic example as they determine most of the physics of solids at low energy, including their density of states, their heat capacity (Debye theory), and even possible superconducting instabilities (BCS theory). Phonon dynamics can be described using elasticity theory1 or hydrodynamics2, from which one derives that their frequency at long-wavelength is linear in the wavevector k, ωL,T = vL,Tk, with the transverse (T) and longitudinal (L) speeds of sound governed by the elastic moduli.
Because of the random atomic distribution and the absence of a fixed equilibrium reference frame, the fate of phonons and the vibrational properties of liquids represent a much harder challenge for both theory and experiments3. In liquids, the dynamics of longitudinal long-wavelength fluctuations have been experimentally ascertained4 to be qualitatively identical to that of solids, despite a smaller sound speed. On the contrary, the dynamics of transverse (or shear) long-wavelength fluctuations are radically different. Liquids have a vanishing static shear modulus and, at small wave-vector, they display a shear diffusion mode rather than propagating shear waves (transverse phonons) as in solids1.
Leveraging on a simple viscoelastic model, Maxwell5 proposed that shear stress in liquids has a characteristic exponential decay time τM ≡ η/G∞, where η is the shear viscosity and the instantaneous shear modulus G∞. This timescale is now known as the Maxwell relaxation time. Based on a more microscopic picture of liquid dynamics, Frenkel later proposed6 to identify such a timescale with the time of local particle re-arrangements, corresponding to hopping processes over potential barriers.
The emerging Maxwell-Frenkel picture of liquid dynamics (see7 for the complete history lesson) suggest that shear fluctuations in liquids obey the following telegrapher equation,
In Eq. (1), vT is the transverse speed of sound related to the instantaneous shear modulus G∞. By solving (1) (see8 for an extensive review), the dispersion of shear waves is obtained,
For long-wavelengths, one recovers the hydrodynamic shear diffusion mode with collective diffusion constant ({D}_{c}equiv {v}_{T}^{2}tau =eta /rho ) (with ρ the mass density of the system) predicted by Navier-Stokes equations9. Above a critical wave-vector kg, known as k-gap, the real part of the frequency becomes nonzero. Above kg, when the real part of ωT becomes larger than its imaginary part, propagating solid-like shear waves are then expected to emerge in liquids. This leads to a corresponding elastic-like response below a certain critical distance, Lc ≡ 2π/kg. In simpler terms, it means that liquids are predicted to exhibit solid-like transverse vibrational modes not only for high-frequency, ω ≫ 1/τ (as proposed initially by Frenkel6 based on a single particle picture), but also for large wave-vectors k ≫ kg.
The k-gap is expected to appear at the melting temperature and to expand into the liquid phase as the temperature increases10. It then reaches the maximum wave vector allowed at the edge with the gas phase. The existence of the k-gap, and its properties as described by equation (1), have been confirmed by several molecular simulations of classical liquids10,11,12,13,14,15,16,17 and other liquid systems18,19,20, indicating the validity of the theoretical framework, even from a quantitative perspective21. Additionally, the similarity between vibrational modes at high frequency/wave-vector in liquids and solids22,23,24, and the solid-like nature of confined liquids at low frequencies, as predicted by the concept of an elastic critical length Lc, have been experimentally verified25,26. However, due to the limitations of experimental scattering techniques at low k and ω, the k-gap has not been observed experimentally in classical liquids. The only recorded observation of the k-gap has been in complex dusty plasmas27, which are systems of charged particles that interact strongly via Coulomb forces and exhibit liquid-like collective dynamics28.
Granular materials differ from conventional thermal equilibrium systems, such as molecular gases, liquids, solids, colloidal liquids, or solids. They are made up of macroscopic particles that experience negligible thermal fluctuations compared to the typical energy scales of the system. Granular materials also have high dissipation due to inter-particle solid friction in the dense solid-like phase or inter-particle inelastic collisions in the fluid-like phase29. This means that a continuous external energy injection is needed to maintain the fluid-like phase of a granular system, making it a prototype of systems far from thermal equilibrium. However, this raises questions about whether the k-gap description, commonly used for molecular liquids or plasma, applies to a granular fluid. More in general, it remains unclear whether collective modes of a granular fluid are similar to those of classical liquids3 and, if so, what plays the role of thermodynamic variables such as temperature.
The densely packed configurations of granular materials often exhibit fluid-like behaviors when subjected to external forces30,31,32,33,34,35. However, due to the presence of permanent contacts and force chains36, measuring the Hessian matrix directly in densely packed granular matter has proven to be a challenging task, even in the quasi-static limit37,38,39,40. In contrast, loose granular matter, where collisions primarily govern particle interactions, has seen theoretical analyses of hydrodynamics and collective modes in granular fluids. This involves formulating transport equations for essential hydrodynamic quantities like mass, momentum, and heat, followed by a linear stability analysis of the homogeneous states41,42,43,44,45,46,47. While transverse modes decouple45, the longitudinal sector becomes intricate due to the non-conservation of energy, leading to significant modifications in the longitudinal channel. To the best of our knowledge, the discussion of the k-gap in the hydrodynamics of granular fluids has been absent, as it extends beyond the conventional long-wavelength hydrodynamic description. From an experimental standpoint, the primary advantage of granular fluids over traditional molecular liquids is the macroscopic size of the particles, which are measured in centimeters in this case. This larger size significantly facilitates the tracking of particle positions and dynamics using cameras. In contrast, such tracking is extremely challenging, if not impossible, with molecular liquids. In those cases, excitations can only be investigated through techniques like X-ray or inelastic neutron scattering, which are far more complicated, especially in the frequency and wave-vector ranges where the k-gap is expected to emerge.
If the k-gap description can apply to granular fluids, then a homogeneously driven granular fluid would be the most straightforward scenario to explore. However, creating a homogeneously driven experimental system has proven challenging48,49,50 due to the influence of gravity and anisotropic driving in three-dimensional (3D) systems or the implementation of boundary driving in quasi-two-dimensional (2D) vertical systems51,52. Additionally, the influence of 3D effects in quasi-2D horizontal systems has made it difficult to achieve homogeneous driving53,54,55. While a few quasi-2D systems have achieved homogeneous driving48,56,57, some lacked single-particle velocity Gaussian statistics56,57. In contrast, others incorporated persistent unidirectional rotation at the single-particle level48, introducing additional complexities in energy injection at the single-particle scale. Recently, Chen et al. designed an experimental system that achieves homogeneous driving, single-particle velocity and rotation statistics with Gaussian distributions of zero means in a quasi-2D system49,50. This system closely aligns with the active Brownian particles introduced in theoretical studies from the perspective of active matter58. We notice that our experimental setup can be considered as an active system since the external energy input, that maintains the system out of equilibrium, acts individually and independently on each “active particle”59.
Compared to a dusty plasma, a nonequilibrium system made up of micron-sized charged particles suspended in a plasma, an active granular system differs significantly in its interaction potential and driving mechanism. In dusty plasmas, the potential is governed by a Yukawa potential, while in active granular systems, interactions occur through inelastic collisions and solid friction. Additionally, the driving mechanisms differ: dusty plasma is driven by laser heating at the boundaries, while active granular systems experience homogeneous and random driving of individual particles. Our experimental findings add to the previous experimental observation of a k-gap dispersion in a dusty plasma27, demonstrating the universality of the k-gap phenomenon in the liquid phases of matter.
Active granular systems provide a novel platform for exploring the emergence of collective dynamics and showcasing a rich interplay of complex phases and phenomena. Our study focuses on bi-disperse active Brownian vibrators. Through measurements of the pair correlation functions, mean square displacements, velocity auto-correlation functions, vibrational density of states, and a detailed analysis of particle motion, we demonstrate that this active system exhibits both gas-like and liquid-like phases, depending on the packing fraction, despite pure hard-disk-like repulsive interactions. Within the granular liquid-like phase, we experimentally validate the existence of a k-gap in the dispersion of transverse excitations. This gap becomes more significant with a decrease in packing fraction and becomes ill-defined in the gas phase because of the disappearance of well-defined modes, aligning with theoretical expectations. Our results offer a direct experimental confirmation of the k-gap phenomenon, extending its relevance beyond classical thermal liquids to active granular systems, and reveal the existence of similarities between the physics of active granular matter and supercritical fluids.
Result and discussions
Granular Brownian vibrators
Figure 1 displays the top view of a layer of bi-disperse Brownian vibrators positioned on the surface of a shaker. This layer of particles is confined within a flower-shaped boundary to prevent the creep motion of particles near the boundary. Each Brownian vibrator has a flat, disk-shaped cap with twelve alternatively inclined legs below the cap. When a vertical sinusoidal vibration is applied to the supporting base, each single Brownian vibrator performs 2D Brownian motion. Previous studies by Chen et al. revealed that the translational and rotational velocities of a single Brownian vibrator follow Gaussian distributions with zero means49. These features lead us to term the particle an “active Brownian particle”, closely mimicking the conditions studied in theoretical investigations of active matter systems58. Moreover, in a collection of Brownian vibrators of the same size, the translational velocities of individual particles follow a Maxwell-Boltzmann distribution for low and intermediate speeds, but show high-energy tails that deviate significantly from the Maxwell-Boltzmann distribution for large speeds, which can be attributed to the inelastic collisions of particles and the homogeneous driving49. Unlike the mono-disperse systems studied earlier49,50, the present system is bi-disperse, which prevents crystallization at high packing fractions, as depicted in Fig. 1. More details about this setup are presented in the Methods. We emphasize that the frequency of the vertical sinusoidal vibration applied is of 100 Hz. In the rest of the manuscript, we will consider only time-scales which are parameterically longer than the driving frequency. In that regime, equilibrium thermodynamic and hydrodynamic concepts can still be applied.

Top view of a layer of active Brownian vibrators. The packing fraction in this layer is ϕ = 0.822. In the upper left corner, one Brownian vibrator is shown. The total area of large and small particles has a fixed ratio of 1/1, and their diameters have a ratio of dl/ds = 1.4/1.
Unlike conventional polar particles studied previously60,61,62, the lack of particle-scale built-in asymmetry makes our experimental system unique and novel. The active force on each nonpolar particle results from collisions between the tilted legs and the vibrating bottom surface, causing the particle’s central axis to tilt slightly away from gravity, allowing only some legs–often represented as a single leg–to be propelled upon contact with the surface. These interactions produce minimal correlation in the contact angle, leading to random driving force directions while maintaining a nearly constant magnitude (see refs. 49,50 for details). In this system, there is no fluid flow from a surrounding solvent, and the main dissipation with the environment results from friction and inelastic collisions. Consequently, our system naturally belongs to the category of dry active matter.
In our present system of bi-disperse (nonpolar) Brownian vibrators, we have observed no global flocking within the several-hour experimental time window, which is likely due to the disorder introduced by bidispersity, in contrast to our previous monodisperse systems50, where we observed global flocking, aligning with the theoretical investigation of active Brownian particles58. Furthermore, we have not observed any phase separation, unlike the self-propelled binary colloids known as Quincke rollers, where significant demixing of small and large colloidal particles occurs after the system begins a global rotational collective motion63. On a microscopic level, the interactions between these Quincke rollers are influenced by electrostatic and hydrodynamical forces, which are quite different from the interactions observed in our nonpolar granular disks in a dry environment. While examining phase separation in nonpolar disks is an promising area for future research, it is currently beyond the scope of this study. Real-time videos of the particle motion are provided with from the Supplementary Movie 1 to 4 for different packing fraction to confirm these statements.
Structural and dynamical crossovers
We experimentally investigate our active granular system by measuring the pair correlation function g(r). In Panel (a) of Fig. 2, we show the results for different values of packing fraction ϕ, which is the fractional area occupied by the particles over the whole system. The first set of peaks corresponds to three peaks of g(r) within the range of 1ds ≤ r ≤ 2ds. It arises from the bi-dispersity of particle sizes, causing a single peak to split into three peaks. The second peak of g(r) is located within 2ds ≤ r ≤ 3ds, and the third peak is in the range of 3ds ≤ r ≤ 4ds. As we decrease the packing fraction, we observe a considerable decrease in the height of the first set of peaks, and the second and third peaks disappear. This observation implies that as we decrease the packing fraction, the medium-range order vanishes, and the system undergoes a structural crossover. This transition occurs at about ϕ = 0.618 and is known in the context of supercritical fluids as the Fisher-Widom line64.

a The pair distribution functions g(r). b The normalized velocity auto-correlation functions (VACF). c The mean square displacement (MSD) of particles. d The second derivative of MSD with respect to time. Here, the values of MSD are properly normalized using 2Z(0), where Z(t) is defined in (3). In panels (a–d), the same set of packing fractions ϕ is chosen following the color scheme described in the legend in (d).
To establish a connection between structure and dynamics, as achieved for supercritical fluids in ref. 65, in Fig. 2b–d we present the experimental results for the velocity auto-correlation functions (VACF) C(t), the mean square displacements (MSD), and their corresponding second derivatives with respect to time (frac{{d}^{2}}{d{t}^{2}},{mbox{MSD}},(t)), as functions of time for various values of the packing fraction ϕ.
We begin by defining the unnormalized VACF Z(t) as
where the index i specifies the Cartesian component of the velocity (overrightarrow{v}) with i = x, y and d = 2 for our system. The statistical average 〈⋅〉 is first taken over different initial times ‘0’ for a given particle and then over all particles. The VACF C(t) is then defined as C(t) ≡ Z(t)/Z(0).
The VACF of a gas decreases continuously with time, while for liquids near melting point and solids, it shows a combination of an oscillatory and a decaying term. The presence of an oscillatory part in the C(t) can be identified by looking for the occurrence of a minimum or a change in its slope. For high packing fractions, C(t) displays a clear minimum below t ≈ 0.1s that gradually disappears as the packing fraction decreases, as shown in Fig. 2b. At a packing fraction ϕ = 0.618, the minimum in the VACF is no longer present, and C(t) becomes a continuously decreasing function, as expected in a gas. In the field of supercritical fluids, this dynamical transition determines the so-called Frenkel line that separates the rigid liquid phase, presenting oscillatory motion, and the non-rigid gas-like fluid phase. Evidence for the structural nature of the Frenkel line, hinting towards a possible equivalence with the Fisher-Widom line concept, has been reported in supercritical fluids66. Despite the complete equivalence between the structural and dynamical criteria remains unproved, a direct connection between the dynamical crossover and thermodynamics has been demonstrated67. Aware of these distinctions, in the rest of the manuscript, we will adopt the jargon rigid (liquid-like) and non-rigid (gas-like) states interchangeably. According to Frenkel’s theory6, this dynamical crossover corresponds microscopically to a situation where the jumping time between oscillatory motion around different local minima of the potential becomes comparable with the shortest vibration time. This dynamical crossover is also expected to coincide with the disappearance of collective shear waves at all frequencies in the liquid. Frenkel idea relies on a single particle picture, while the k-gap equation (1) describes collective dynamics. The time-scale τ appearing in Eq. (1) is therefore more correctly identified with the Maxwell relaxation time. In simple fluids the Maxwell time is very close to the lifetime of local connectivity68, that is another single particle concept that will be analyzed below. The critical packing fraction of ϕ = 0.618 is very close to the value at which medium-range order disappears in the pair correlation functions shown in panel (a). This suggests a significant link between structure and collective dynamics in active granular systems, similar in spirit to the results presented in ref. 66. We notice that there is now firm experimental evidence that the change of particle dynamics at the Frenkel line is seen in structural changes69. Despite a complete analysis of these structural changes being beyond the scope of the present work, in the Supplementary note 2, we provide a preliminary study of the experimental structure factor S(k) for four different packing fractions that confirms the gradual disappearance of structural order by decreasing the packing fraction. The relation between structural changes and dynamics in liquids is still poorly understood (see, for example,7). It would be necessary to explore this connection further in granular fluids.
Another substantial dynamical quantity, besides the VACF, is the MSD of particle motion. This is shown in Fig. 2c, where curves for the same set of packing fractions ϕ as in panel (b) are drawn. For ϕ = 0.618 and below (data not shown), MSD is quasi-ballistic for t < 0.1s, indicating underdamped particle dynamics. For t > 0.1s, the slope of MSD is close to one, showing the diffusive motion of a particle. However, for very large t, the slope of MSD deviates from one due to the finite system size as the length scale of MSD gradually approaches the system size. For ϕ > 0.618, the quasi-ballistic motions at small t are not clearly visible, and a subdiffusive regime at intermediate times emerges. The upper bound of this regime depends on the value of ϕ, beyond which the diffusive behaviors recover for the curves corresponding to ϕ = 0.773 and ϕ = 0.811. At ϕ = 0.822, within the entire observation time window of 0.025s ≤ t ≤ 1000s, the MSD shows subdiffusive behaviors, indicating the progressively more significant glassy dynamics with the increment of ϕ.
The MSD and the VACF are closely related to each other since (frac{{d}^{2}}{d{t}^{2}}{{rm{MSD}}}(t)=4Z(t))70. This equivalence has been already utilized in ref. 66 to investigate the “liquid-gas” transition. Upon comparing panels b and d of Fig. 2, we experimentally verified this equivalence, and the results show an excellent quantitative mutual agreement.
Granular systems are athermal in nature because of the macroscopic size of their constituents. At fixed activity, it is therefore necessary to understand which parameter plays an analogous role of temperature in thermal fluids, driving the system from a gas-like to a liquid-like phase. Both structural and dynamical observables suggest that the packing fraction ϕ of the Brownian vibrators plays such a role. In order to provide further qualitative evidence for this analogy, in Fig. 3, we plot the average velocity squared Z(0), that is proportional to the average kinetic energy, as a function of the packing fraction for small (red), large (blue) and all (gray) particles. We observe a clear anti-correlation between the two quantities, which is well fitted by a phenomenological function (Z(0)propto sqrt{{phi }_{c}-phi }), with ϕc ≈ 0.8318. In the literature (e.g.,71), Z(0) has been often associated to an effective granular temperature, Te. Our results demonstrate therefore that Te anti-correlates with ϕ, consistent with the dynamics experimentally observed both at the particle-level and collective scale. Temperature is clearly defined only in equilibrium thermodynamic systems. Defining temperature rigorously in non-equilibrium systems is often very challenging. Here, we introduce the term ‘temperature’ for our system only in a vague, intuitive sense, drawing an analogy with well-defined meanings for ordinary materials made of molecules in thermal equilibrium. It is necessary to notice that the scaling of the experimental data follows a mean-field behavior. We are not aware of any theoretical explanation of this phenomenon. This analysis suggests the existence of a critical packing fraction ϕc which might be connected to a jamming type transition in active granular systems. We leave the exploration of these two points for future research. Finally, we notice that the average velocity squared Z(0) for large particles is consistently larger than that of small particles, with this difference becoming more pronounced for small packing fraction. Given that the steady state velocity is a result of the balance between the energy injected and the energy loss caused by friction, this might be explained by the fact that larger particles experience stronger drag force induced by activity.

The average velocity squared Z(0) for small (red), large (blue) and all (gray) particles as a function of the packing fraction ϕ. All the data are well fitted by a phenomenological function (Z(0)=asqrt{b-phi }), where b = 0.8318 for all curves. The error bar of averaged value here is obtained as the difference of minimum and averaged value.
The radial pair correlation function is defined as,
where ({overrightarrow{r}}_{ij}) is the vector between ith and jth particles, N is the total number of particles, and ρ is the particle density of system. As shown in Fig. 2 a, three peaks appear below r = 2ds. The first peak is at r = 1ds, which indicates small-small particle pairing, the second indicates small-big particle pairing at r = 1.2ds, and the third is big-big particle pairing at r = 1.4ds, which is also the diameter of the big particle.
From gas-like to liquid-like dynamics
To study collective motion, we examine displacement vectors and vibrational density of states (VDOS) g(ω), which can be obtained by diagonalizing the dynamical matrix computed from the displacement correlation matrix (see Methods).
In Fig. 4a, we observe that for the lowest packing fraction data, the VDOS decreases monotonically with frequency and can be accurately described by a Lorentzian line shape, g(ω) = α/(α2 + ω2) (dashed black line), at least for frequencies below 5 Hz. This line shape is indicative of purely Langevin diffusive dynamics. It is typical of a gas-like state, where particle collisions are almost uncorrelated and independent and can be described by kinetic theory. The Lorentzian fit becomes less accurate at higher frequencies, indicating that the low-packing fraction system is a dilute liquid rather than an ideal gas of free particles. Additionally, the longitudinal and transverse components of g(ω) are identical, confirming the emergent isotropy of the low-packing fraction phase. Due to the dilute packing, the constituent particles exhibit random, uncorrelated motion in both amplitude and direction. Individual particles’ motion is collisional, resulting in substantial displacements away from their initial positions.

The experimental vibrational density of states (VDOS) for packing fraction ϕ = 0.618, ϕ = 0.773, ϕ = 0.811, and ϕ = 0.822 from panels (a) to (d). The VDOS is represented by gray, blue, and red lines, indicating the total VDOS and longitudinal and transverse components respectively. The VDOS is normalized by setting its area ∫g(ω)dω to 1, and the frequency is measured in Hz. The black dashed line in panel (a) shows the Lorentzian line shape α/(ω2 + α2) with α = 3.3. The top insets show the displacement configuration, with the intensity of the color indicating the amplitude of the single particle displacement, with darker shades representing more significant displacement. The displacement vectors have been enlarged ×1, ×2, ×3, and ×10 in panels (a) to (d), respectively. The bottom panels (e) and (f) show the trajectories of each particle, with the color changing from green to yellow indicating time evolution.
Upon increasing the packing fraction and entering the rigid liquid phase described above, the VDOS undergoes significant changes. A weak and broad peak emerges around 7 Hz in both the transverse and longitudinal components, indicating the emergence of strongly overdamped collective motion. The VDOS is no longer monotonic, and the transverse and longitudinal components begin to display a rich behavior, that can be possibly thought as the combination of a gas-like and a solid-like contribution as proposed in ref. 72. The correlated motion also appears at the level of the single particle displacement field. This field now presents geometric structures composed of vortex-like and string-like patterns and an increasing degree of heterogeneity with localized areas of large displacement (intense purple color) separated by more rigid regions can be observed (see Fig. 4b).
Moving further to ϕ = 0.811 (Fig. 4c) results in the disappearance of the Lorentzian gas-like contribution to the VDOS at low frequency. The total VDOS increases monotonically with frequency up to a peak located around ω = 15 Hz, which corresponds approximately to the average pseudo-Van Hove energy of the emergent collective longitudinal and transverse excitations, as shown in Fig. 5. The zero-frequency values of the total VDOS g(0) exhibit a substantial decrease with an increase in the packing fraction ϕ. This value correlates with the self-diffusion constant D, and a direct relation between D and g(0) can be derived for pure Langevin diffusion. This is consistent with the experimental data for the MSD presented in Fig. 2c. g(0) is largely dominated by the transverse component, and gL(0) is almost zero at ϕ = 0.811. Additionally, the system becomes strongly anisotropic as the longitudinal and transverse VDOS are significantly different, and higher-energy modes appear in the spectrum up to a frequency of approximately 35 Hz. The maximum frequency mentioned is determined by the highest acquisition rate of the cameras we use to track the position of particles in our granular matter system, which is 40 Hz. This maximum frequency has no intrinsic physical significance and is unrelated to the driving frequency.

The longitudinal (blue) and transverse (red) dynamical structure factors SL,T(ω, k) for packing fraction ϕ = 0.618, ϕ = 0.773, ϕ = 0.811 and ϕ = 0.822 from (a) to (d) on top and from (e) to (h) on bottom. ds is the diameter of small particles. In (c, d), the solid black lines show the fitting to the linear dispersion relation ({{rm{Re}}}(omega )={v}_{L}k). The fitted longitudinal speed of sound is vL = 101.2 and vL = 116.9 respectively in unit of (ds/s). The colored bullets indicate the peak position obtained by fitting SL(ω, k) at fixed values of the wave vector k. In (g, h), the numerical data are fitted with the k-gap formula ({{rm{Re}}}(omega )=sqrt{{v}_{T}^{2}({k}^{2}-{k}_{g}^{2})}) where kg is the size of the k-gap. We obtain a k-gap of 0.093 for ϕ = 0.811 (g) and 0.068 for ϕ = 0.822 (h) in units of ({d}_{s}^{-1}). The fitted transverse speed of sound is vT = 48.5 and vT = 71.7 respectively in unit of (ds/s). In all panels, the error bar indicates the linewidth of the corresponding peak in the dynamic structure factor. The dynamical structure factor is dimensionless and normalized by its maximum value, hence, the color bar goes from 0 to 1 in all panels.
As the packing fraction of the system increases to ϕ = 0.822, the value of g(0) becomes extremely small but still finite, indicating that the system is close to the solid phase but not quite there yet. This is later confirmed by the dispersion of shear waves. The longitudinal and transverse VDOS are linear in frequency up to about 8 Hz, as dictated by Debye’s law in 2D, and also commonly found in bulk liquids73,74. However, a sharper peak appears around 12 Hz, which is attributed to the flattening of the dispersion of the collective modes obtained from the dynamical structure factor (Fig. 5d, h). As a result, particle displacements become small, and granular particles move very little away from their initial positions. The dynamics and corresponding VDOS increasingly resemble those of a dense viscous liquid with enhanced solid-like elastic vibrations or a liquid near melting.
In summary, the analysis of the VDOS and the particle displacements reveal a continuous transition from a gas-like behavior typical of dilute liquids to a collective and viscoelastic motion characteristic of dense liquids by increasing the packing fraction. This perfectly aligns with the structural and dynamical transition between a gas-like liquid and a rigid liquid phase discussed in the previous section and displayed in Fig. 2. The study of collective modes performed in Fig. 5 confirms that ϕ plays the role of the inverse temperature of classical thermal liquids (see Fig. 3 and related discussion above). Indeed, the behavior of the VDOS shown in Fig. 4a–d is perfectly compatible with that found in liquids upon decreasing T (see, for example, Fig. 6 in ref. 72).

The transverse dynamical structure factors under different driving frequencies f = 100, 110, 120, 130 Hz from (a) to (d). The packing fraction is fixed to ϕ = 0.822, and the acceleration is fixed to a = 2.5g.
Gapped shear waves
To analyze the vibrational dynamics of our system, we consider the experimental dynamical structure factor SL,T(k, ω) (details in Methods). Figure 5 shows a color map of the obtained dispersion relations of the longitudinal and transverse parts for different packing fractions ϕ ranging from the gas-like phase ϕ = 0.618 to the dense liquid-like phase ϕ = 0.822. At the low packing fraction ϕ = 0.618, no distinct collective excitation is observed in both the longitudinal and transverse sectors (panels (a) and (e)). Instead, an incoherent diffusive signal, which extends from zero frequency to approximately 10 Hz roughly independent of the wave vector, dominates the response. Moreover, the longitudinal and transverse components are almost identical, consistent with the emergent isotropy at low packing fractions and the results for the VDOS in Fig. 4a. This is also compatible with the structural analysis and the dynamical crossover observed in the VACF shown in Fig. 2.
At packing fraction of ϕ = 0.773, panels (b) and (f), the system is in a dilute liquid phase. This phase displays only weak indications of collective motion, with the exception of a strong diffusive signal around ω = k = 0. There are also faint indications of collective modes in the longitudinal and transverse sectors. However, these excitations are strongly overdamped, which makes it difficult to identify their energy using wave-vector cuts precisely. Despite this, there appears to be a k-gap emerging around (k=0.08{d}_{s}^{-1}), where ds is the diameter of the small particles. This can be anticipated despite the challenges in identifying the energy of the excitations due to the broadening of the color map in panels (b) and (f).
As we move to higher packing fraction data, ϕ = 0.811, 0.822 shown in panels (c–g) and (d–h), the fingerprints of liquid-like collective modes become more pronounced. By fitting the dynamical structure factor, we have extracted the frequency Ω(k) and the line-width Γ(k) of the lowest collective modes in both the longitudinal and transverse sectors. Details regarding the fitting process can be found in the Methods section, and additional figures can be found in the Supplementary note 3. In the longitudinal sector, as shown in panels (c) and (d), we observe a mode that disperses linearly at large wavelengths, with ΩL(k) = vLk, where vL is of the order 102 in units of ds/s. The speed increases slowly as the packing fraction increases, confirming that ϕ plays the role of inverse temperature in classical thermal liquids. At (kapprox 0.1,{d}_{s}^{-1}), the dispersion bends down towards a constant pseudo-Van-hove plateau. Additionally, the longitudinal sound mode has a linewidth that increases with the wave vector k and becomes overdamped for small wavelengths.
The transverse component of the dynamical structure factor in the dense liquid phase at a high packing fraction is shown in panels (g) and (h). Apart from a strong diffusive signal near the origin, ST(k, ω) presents a distinctive k-gap dispersion for the transverse shear waves, which is confirmed by using cuts at constant k (red symbols). The value of the k-gap is approximately kg = 0.093 for ϕ = 0.811 and kg = 0.068 for ϕ = 0.822 in units of ({d}_{s}^{-1}). As the packing fraction ϕ increases, the value of the k-gap becomes smaller, which is equivalent to decreasing the temperature in classical liquid systems. We find that the frequency of the gapped shear waves is well fitted by the theoretical formula ({Omega }_{T}(k)=sqrt{{v}_{T}^{2}({k}^{2}-{k}_{g}^{2})}), with a transverse speed of sound of the order of vT = 48.5 and vT = 71.7 in the unit of (ds/s) for ϕ = 0.811 and ϕ = 0.822, respectively. The linewidth of the shear waves becomes larger with the wave vector k, similar to longitudinal waves.
In simple fluids with steep interactions, it has been found that vL ≈ 10vth and vT ≈ 5vth near freezing75, where vth is the thermal velocity of the system. Using our experimental data, we have estimated the ratio between the sound velocities and the thermal velocity computed from Z(0) (see Table 1). We found that both ratios are one order of magnitude larger than in simple fluids. We speculate that this difference is a direct consequence of the activity of the system and of the strong deviations from local thermal equilibrium.
The lifetime of local connectivity, or τLC, in classical liquids, is closely related to the Maxwell relaxation time, or τ, as per a research study68. For a value of ϕ = 0.822, our analysis yields a value of τLC as 0.3 seconds, which is in proximity to the Maxwell relaxation time of 0.11 seconds, calculated from the size of k-gap (as shown in panel (d) of Fig. 5). This similarity in values confirms the validity of our data analysis method for obtaining the dynamical structure factor. Also, it suggests that classical liquids in the supercritical state and granular fluids share similarities at the smallest particle level.
To understand how various parameters in the experimental setup is also necessary, in addition to the packing fraction ϕ, influence the k-gap behavior. We conducted additional experiments to investigate this by increasing the driving frequency from 100 Hz to 130 Hz while maintaining the highest packing fraction of ϕ = 0.822. We want to emphasize that the acceleration of the vibrator was held constant. This means that by changing the frequency f, the amplitude of the oscillations was also altered. In Fig. 6, we present the transverse spectrum for four different driving frequencies: f = 100, 110, 120, 130 Hz, with a constant packing fraction of ϕ = 0.822 and acceleration a = 2.5g.
We observe that the position of the k-gap is largely independent of the driving frequency, at least in the range we considered. In contrast, the dispersion of shear waves is sensitive to the driving frequency and exhibits a complex dependence on the frequency f. A similar pattern is seen in the dispersion of longitudinal phonons, as illustrated in Fig. 7. Further experimental investigation is necessary to thoroughly analyze how collective motion varies with driving frequency and amplitude.

The longitudinal dynamical structure factor under different driving frequencies f = 100, 110, 120, 130 Hz from (a) to (d). The packing fraction is fixed to ϕ = 0.822 and the acceleration to a = 2.5g.
Although we cannot conduct experiments with frictionless granular particles for direct comparison, we believe that interparticle friction is less significant than the interactions resulting from inelastic collisions between particles. Several pieces of evidence from our previous study of monodisperse particles49,50 support this hypothesis.
First, we observed no correlations between individual particles’ translational and rotational degrees of freedom. This suggests that particle rotation due to interparticle friction is largely independent of their translational motions. Second, we successfully explained the flocking behavior observed in our previous monodisperse system50 by applying the theory of active Brownian particles58, which completely disregards interparticle friction, as the theory assumes.
Third, we employed vertical driving instead of the more conventional shear driving, reducing interparticle friction’s impact. Lastly, even in our bidisperse systems at the highest packing fraction, the system remains unjammed, indicating that the dominant interactions are due to interparticle collisions rather than contact forces associated with permanent contacts, as seen in jammed solids. Therefore, we anticipate that the results will be comparable to those of frictionless particles.
Conclusion
Our study demonstrates that, despite significant differences in the size of composed particles and the absence of a classical thermodynamic description, granular matter on the Brownian vibrator exhibits various similarities with classical liquids in the supercritical region both in the collective and particle-level motions. These alike aspects include pair correlation functions, velocity autocorrelation, mean square displacement on particle levels, vibrational density of states, and dispersion relation of collective excitations. With an increase in the packing fraction ϕ, distinct phases from gas-like to condensed states (liquid or solid) have been observed on both particle and collective motion levels, suggesting the role of ϕ as an effective inverse temperature variable, that is corroborated by the anti-correlation between the average kinetic energy and the packing fraction. However, it is needed to notice that the analogy of temperature is only vague and intuitive and is limited to our specific system. Furthermore, our work experimentally revealed the emergence of the k-gap, as predicted by viscoelasticity theory, in granular fluids. This provides a direct link between classical liquid and active granular matter and experimentally confirms the phenomenon of k-gap and its theoretical premises.
Methods
Experimental setup
Our experimental system comprises a horizontal layer of granular particles driven vertically by a sinusoidal oscillation with a fixed frequency f and amplitude A induced by an electromagnetic shaker.
Unless indicated otherwise, the frequency f is set at 100 Hz, and the maximum acceleration a is 2.5 times the gravitational acceleration g = 9.8m/s2. The vibration amplitude A, defined as a/(2πf)2, is 0.062 mm, significantly smaller than the particle’s vertical dimension (6 mm). Therefore, we can neglect the vertical displacement of a particle, treating the system as quasi-two-dimensional49.
The upper left corner in Fig. 1 depicts a granular particle as a disk-shaped body with 12 alternately inclined supporting legs. The small particle’s disk has a diameter d of 16 mm and a thickness of 3 mm. The legs, with a height of 3 mm, are inclined inward by 18. 4°, and alternately deviated from the mid-axis plane by ± 38. 5°.
For a single particle, the distributions of rotational velocity and translational velocity components vx and vy follow Gaussian distributions with negligible mean values compared to their standard deviations, typically less than ten percent. This indicates that the motion of a single particle is both random and isotropic, leading us to term the particle an Active Brownian Particle (ABP), closely mimicking conditions studied in theoretical investigations of active matter systems58.
To prevent crystallization, we utilize bi-disperse particles with a size ratio measured in terms of disk diameters of 1: 1.4 and a number ratio of 2: 1 for small and large particles. These parameters maintain an approximately equal area ratio between small and large particles across a range of packing fractions. The packing fraction ϕ is defined as the ratio between the area occupied by all particles and the confining area of the particle layer.
These ratios are derived from jamming studies of binary disks, particularly the research conducted by O’hern et al.76,77. Their simulations employed a size ratio of 1:1.4 because a large disk can have up to seven nearest neighbor small disks, while a small disk can have up to five nearest neighbor large disks. This configuration aligns with the theoretical concept of the 7-5 defect roles in the KTHNY theory78.
In the original O’hern algorithm, a 1:1 number ratio was utilized. However, it was later discovered that a 2:1 number ratio between small and large disks, which corresponds to an area ratio of 1:1, yields even better results. A recent experimental study on the vibrational density of states in jammed 2D granular packing39 analyzed how the density of states changed when the number ratio varied from a hexagonal crystal (1:0) to the most disordered binary mixture (2:1). The density of states for a system with a 2:1 ratio of small to large disks exhibited several features remarkably similar to those found in molecular glasses37,39.
The particles are placed on top of an aluminum plate (60 cm × 60 cm) and confined within a flower-shaped boundary designed to suppress creep particle motions along the boundary49. We initiate the experiment with all particles randomly and uniformly placed on the base plate. After applying vibration for two hours, we achieve an initial state of particle configuration. Subsequently, the particle layer undergoes continuous vertical vibration, while a Basler CCD camera (acA2040-180kc) records particle motion at 40 frames/s for at least an hour.
Displacement correlation matrix and dynamical structure factor
The displacement correlation matrix C is defined as,
where ni(t) is the displacement of ith degree of freedom at time t. The dynamical matrix can hence be calculated as,
where mi is the mass of the ith degree of freedom, and α is a dimensionful parameter which will be later specified. Diagonalizing the matrix D, one obtains the eigenvalues κi and the eigenfrequencies ({omega }_{i}=sqrt{{kappa }_{i}}). The eigenvector fields u are then defined by solving the eigenvalue problem,
Up to this point, the precise numerical values of the obtained eigenfrequencies ωi are only determined up to an unknown energy scale α. In athermal systems, such a scale cannot be determined using the temperature T, since the latter is not well defined. See79 for an extensive discussion on this issue. In order to fix the value of α, we resorted to a more phenomenological approach. More precisely, we have rescaled all the eigenfrequencies by setting the value of the largest one to coincide with the highest observable frequency of our instrument, namely 40Hz. In doing so, the corresponding eigenfrequencies have now the correct physical dimension. We validate a posteriori our hypothesis by noticing that the time scale obtained from k-gap τ is compatible with the average local connectivity time, as observed in classical thermal liquids68.
As an ulterior check of the validity of our criterion to fix α, we compute the mean particle kinetic energy, 〈Ekin〉 ≡ mZ(0) (see Fig. 3). Following the literature71, we can define an effective granular temperature Te, and determine α using α = kBTe = 〈Ekin〉. By adopting this method, we find a normalization for α that is compatible, apart from an ({{mathcal{O}}}(1)) prefactor, with the previous estimates. This confirms the validity of our previous arguments. Finally, we emphasize that the value of α, apart from fixing the physical dimension of our frequencies, affects only the values of the frequencies but not the qualitative physical trends discussed in the main text, as for example the shape of the VDOS or of the dispersion of collective excitations. At the same time, the value of the k-gap does not depend on the determination of α.
To separate the longitudinal and transverse components properties, we Fourier transform the eigenvector fields
and
where the index i indicates the eigenvector field corresponding to eigenfrequency ωi, the index j indicates the jth particle, rj is the equilibrium position of jth particle and ui,j is the eigenvector field corresponding to frequency ωi and position rj. T, L stand respectively for transverse and longitudinal. Finally, z is the spatial coordinate perpendicular to the 2D (x, y) plane. The current correlation function can be obtained as,
The dynamical structure factor are given by,
Finally, the vibrational density of states (VDOS) can be obtained directly by counting the distribution of eigenfrequencies ωi, which can be expressed as,
Alternatively, one can integrate the current correlation function to obtain longitudinal and transverse components separately. The separated VDOS are given by,
where g(ω) = gL(ω) + gT(ω) and ({{mathcal{N}}}) is just a normalization factor to ensure that ∫g(ω)dω = 2N where N is the number of particles. These two methods give same results for the total VDOS.
Data analysis
The following fitting functions have been used to analyze the experimentally obtained dynamical structure factor (see additional figures in the Supplementary note 3). The transverse and longitudinal components of the dynamical structure factor are respectively fitted using.
The first terms are simply a damped harmonic oscillator with energy ΩL,T(k) and linewidth ΓL,T(k). The second term is a quasi-elastic contribution modeled with a Lorentzian function. The dispersion relations ΩT,L(k) are shown in Fig. 5 as colored bullets and the corresponding linewidths ΓT,L by the related errors bars. Sloc(ω, k) is a term representing possible quasi-localized low-energy modes80, and it is not relevant for the present discussion. More details on the fitting procedure and the raw data for the dynamical structure factors can be found in the Supplementary note 3.
Responses