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 theory¹ or hydrodynamics², from which one derives that their frequency at long-wavelength is linear in the wavevector k, ω_L,T = v_L,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 experiments³. In liquids, the dynamics of longitudinal long-wavelength fluctuations have been experimentally ascertained⁴ 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 solids¹.

Leveraging on a simple viscoelastic model, Maxwell⁵ 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 proposed⁶ 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 (see⁷ for the complete history lesson) suggest that shear fluctuations in liquids obey the following telegrapher equation,

$${omega }_{T}^{2}+i{omega }_{T}{tau }_{M}^{-1}={v}_{T}^{2}{k}^{2}.$$

(1)

In Eq. (1), v_T is the transverse speed of sound related to the instantaneous shear modulus G_∞. By solving (1) (see⁸ for an extensive review), the dispersion of shear waves is obtained,

$${omega }_{T}=-frac{i}{2tau }+{v}_{T}sqrt{{k}^{2}-{k}_{g}^{2}}qquad ,{{mbox{with}}},qquad {k}_{g}^{-1}equiv 2{v}_{T}tau .$$

(2)

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 equations⁹. Above a critical wave-vector k_g, known as k-gap, the real part of the frequency becomes nonzero. Above k_g, 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, L_c ≡ 2π/k_g. 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 Frenkel⁶ based on a single particle picture), but also for large wave-vectors k ≫ k_g.

The k-gap is expected to appear at the melting temperature and to expand into the liquid phase as the temperature increases¹⁰. 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 liquids^{10,11,12,13,14,15,16,17} and other liquid systems^18,19,20, indicating the validity of the theoretical framework, even from a quantitative perspective²¹. Additionally, the similarity between vibrational modes at high frequency/wave-vector in liquids and solids^22,23,24, and the solid-like nature of confined liquids at low frequencies, as predicted by the concept of an elastic critical length L_c, have been experimentally verified^25,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 plasmas²⁷, which are systems of charged particles that interact strongly via Coulomb forces and exhibit liquid-like collective dynamics²⁸.

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 phase²⁹. 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 liquids³ 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 forces^{30,31,32,33,34,35}. However, due to the presence of permanent contacts and force chains³⁶, measuring the Hessian matrix directly in densely packed granular matter has proven to be a challenging task, even in the quasi-static limit^37,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 states^{41,42,43,44,45,46,47}. While transverse modes decouple⁴⁵, 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 challenging^48,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 systems^51,52. Additionally, the influence of 3D effects in quasi-2D horizontal systems has made it difficult to achieve homogeneous driving^53,54,55. While a few quasi-2D systems have achieved homogeneous driving^48,56,57, some lacked single-particle velocity Gaussian statistics^56,57. In contrast, others incorporated persistent unidirectional rotation at the single-particle level⁴⁸, 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 system^49,50. This system closely aligns with the active Brownian particles introduced in theoretical studies from the perspective of active matter⁵⁸. 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”⁵⁹.

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 plasma²⁷, 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 means⁴⁹. These features lead us to term the particle an “active Brownian particle”, closely mimicking the conditions studied in theoretical investigations of active matter systems⁵⁸. 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 driving⁴⁹. Unlike the mono-disperse systems studied earlier^49,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.

**Fig. 1: Active granular matter in the Lab.**

Unlike conventional polar particles studied previously^60,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 systems⁵⁰, where we observed global flocking, aligning with the theoretical investigation of active Brownian particles⁵⁸. 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 motion⁶³. 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 1d_s ≤ r ≤ 2d_s. 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 2d_s ≤ r ≤ 3d_s, and the third peak is in the range of 3d_s ≤ r ≤ 4d_s. 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 line⁶⁴.

**Fig. 2: Disappearance of medium-range order and the dynamical transition between a liquid-like to gas-like phase.**

To establish a connection between structure and dynamics, as achieved for supercritical fluids in ref. ⁶⁵, 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

$$Z(t)equiv langle {v}_{i}(0){v}_{i}(t)rangle =frac{1}{d}langle overrightarrow{v}(0)overrightarrow{v}(t)rangle ,$$

(3)

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 fluids⁶⁶. Despite the complete equivalence between the structural and dynamical criteria remains unproved, a direct connection between the dynamical crossover and thermodynamics has been demonstrated⁶⁷. 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 theory⁶, 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 connectivity⁶⁸, 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. ⁶⁶. We notice that there is now firm experimental evidence that the change of particle dynamics at the Frenkel line is seen in structural changes⁶⁹. 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,⁷). 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))⁷⁰. This equivalence has been already utilized in ref. ⁶⁶ 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.,⁷¹), Z(0) has been often associated to an effective granular temperature, T_e. Our results demonstrate therefore that T_e 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.

**Fig. 3: The role of the packing fraction in active granular matter.**

The radial pair correlation function is defined as,

$$g(r)=frac{1}{2pi rNrho }{sum }_{i=1}^{N}{sum }_{jne i}^{N}delta left(r-| {overrightarrow{r}}_{ij}| right),$$

(4)

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 = 2d_s. The first peak is at r = 1d_s, which indicates small-small particle pairing, the second indicates small-big particle pairing at r = 1.2d_s, and the third is big-big particle pairing at r = 1.4d_s, 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(ω) = α/(α² + ω²) (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.

**Fig. 4: From isotropic gas-like motion to collective liquid-like dynamics.**

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. ⁷². 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 g_L(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.

**Fig. 5: Emergence of wave-vector gapped shear waves in liquid granular matter.**

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 liquids^73,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. ⁷²).

**Fig. 6: Transverse dispersion relation under different driving frequency.**

Gapped shear waves

To analyze the vibrational dynamics of our system, we consider the experimental dynamical structure factor S_L,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 d_s 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) = v_Lk, where v_L is of the order 10² in units of d_s/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, S_T(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 k_g = 0.093 for ϕ = 0.811 and k_g = 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 v_T = 48.5 and v_T = 71.7 in the unit of (d_s/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 v_L ≈ 10v_th and v_T ≈ 5v_th near freezing⁷⁵, where v_th 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.

Table 1 Speed of sound

Full size table

The lifetime of local connectivity, or τ_LC, in classical liquids, is closely related to the Maxwell relaxation time, or τ, as per a research study⁶⁸. 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.

**Fig. 7: Longitudinal dispersion relation under different driving frequency.**

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 particles^49,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 system⁵⁰ by applying the theory of active Brownian particles⁵⁸, 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/s². The vibration amplitude A, defined as a/(2πf)², 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-dimensional⁴⁹.

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 v_x and v_y 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 systems⁵⁸.

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 theory⁷⁸.

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 packing³⁹ 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 glasses^37,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 boundary⁴⁹. 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,

$${C}_{ij}={langle n{(t)}_{i}n{(t)}_{j}rangle }_{t},$$

(5)

where n_i(t) is the displacement of ith degree of freedom at time t. The dynamical matrix can hence be calculated as,

$${D}_{ij}=frac{alpha ,{C}_{ij}^{-1}}{sqrt{{m}_{i}{m}_{j}}},$$

(6)

where m_i 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,

$${{boldsymbol{D}}}{{boldsymbol{u}}}={omega }^{2}{{boldsymbol{u}}}.$$

(7)

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. See⁷⁹ 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 liquids⁶⁸.

As an ulterior check of the validity of our criterion to fix α, we compute the mean particle kinetic energy, 〈E_kin〉 ≡ mZ(0) (see Fig. 3). Following the literature⁷¹, we can define an effective granular temperature T_e, and determine α using α = k_BT_e = 〈E_kin〉. 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

$${F}_{L,i}({{boldsymbol{k}}})={{boldsymbol{k}}}cdot {sum}_{j}{{{boldsymbol{u}}}}_{i}({{{boldsymbol{r}}}}_{j}){e}^{-i{{boldsymbol{k}}}cdot {{{boldsymbol{r}}}}_{j}},$$

(8)

and

$${F}_{Tz,i}({{boldsymbol{k}}})={{boldsymbol{z}}}cdot {{boldsymbol{k}}}times {sum}_{j}{{{boldsymbol{u}}}}_{i}({{{boldsymbol{r}}}}_{j}){e}^{-i{{boldsymbol{k}}}cdot {{{boldsymbol{r}}}}_{j}},$$

(9)

where the index i indicates the eigenvector field corresponding to eigenfrequency ω_i, the index j indicates the jth particle, r_j is the equilibrium position of jth particle and u_i,j is the eigenvector field corresponding to frequency ω_i and position r_j. 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,

$${C}_{T,L}({{boldsymbol{k}}},omega )={left| {sum}_{i}{F}_{(T,L),i}({{boldsymbol{k}}}){delta }_{omega ,{omega }_{i}}right| }^{2}.$$

(10)

The dynamical structure factor are given by,

$${S}_{L,T}({{boldsymbol{k}}},omega )propto frac{1}{{k}^{2}}{left| {sum}_{i}{F}_{(L,T),i}({{boldsymbol{k}}}){delta }_{omega ,{omega }_{i}}right| }^{2}.$$

(11)

Finally, the vibrational density of states (VDOS) can be obtained directly by counting the distribution of eigenfrequencies ω_i, which can be expressed as,

$$g(omega )={sum}_{i}delta (omega -{omega }_{i}).$$

(12)

Alternatively, one can integrate the current correlation function to obtain longitudinal and transverse components separately. The separated VDOS are given by,

$${g}_{L,T}(omega )={{mathcal{N}}}int{C}_{L,T}({{boldsymbol{k}}},omega )d{{boldsymbol{k}}}$$

(13)

where g(ω) = g_L(ω) + g_T(ω) 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.

$${S}_{T}(omega ,k)propto frac{{omega }^{2}}{{left({omega }^{2}-{Omega }_{T}^{2}right)}^{2}+{omega }^{2}{Gamma }_{T}^{2}}+frac{1}{pi }frac{gamma }{{omega }^{2}+{gamma }^{2}}+{S}_{{{rm{loc}}}}(omega ,k),$$

(14)

$${S}_{L}(omega ,k)propto frac{{omega }^{2}}{{left({omega }^{2}-{Omega }_{L}^{2}right)}^{2}+{omega }^{2}{Gamma }_{L}^{2}}.$$

(15)

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. S_loc(ω, k) is a term representing possible quasi-localized low-energy modes⁸⁰, 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.