Anomalous flocking in nonpolar granular Brownian vibrators

Using Brownian vibrators, we investigated the structures and dynamics of quasi-2d granular materials, with packing fractions (ϕ) ranging from 0.111 to 0.832. Our observations revealed a remarkable large-scale flocking behavior in hard granular disk systems, encompassing four distinct phases: granular fluid, flocking fluid, poly-crystal, and crystal. Anomalous flocking emerges at ϕ = 0.317, coinciding with a peak in local density fluctuations, and ceased at ϕ = 0.713 as the system transitioned into a poly-crystal state. The poly-crystal and crystal phases resembled equilibrium hard disks, while the granular and flocking fluids differed significantly from equilibrium systems and previous experiments involving uniformly driven spheres. This disparity suggests that collective motion arises from a competition controlled by volume fraction, involving an active force and an effective attractive interaction resulting from inelastic particle collisions. Remarkably, these findings align with recent theoretical research on the flocking motion of spherical active particles without alignment mechanisms.

1 Supplementary Methods

The dynamics of active particles
From the theoretical work by Caprini et al. [1], the dynamics of active particles reads: Here the effective dissipative force |F d | ≡ mv 0 /τ I .When the system reaches a steady state, the rate of energy dissipation within the system equals the rate of energy input from the active force.Therefore, the magnitude of the active force |F a | = |F d | ≡ mv 0 /τ I = γ 0 v 0 , where γ 0 = m/τ I is the effective damping coefficient.
Our analysis focuses on the comparison of the dimensionless coefficients between the interparticle force F i and the active force γ 0 v 0 n i .First, by calculating the distribution of particle number density, we can determine the equivalent interaction force F i of inelastic collisions between particles.
The radial distribution function, also known as the pair correlation function, is defined as: where L is the system size, N is the total number of particles, |⃗ r jk | refers to the relative distance between two particles j and k, and ∆r is a small increment along the radial direction.Supplementary Fig. 1 The first peak of the pair correlation function in three systems.The figure of g(r = D) versus ϕ.Red circles refer to our experimental results of Brownian vibrators.Blue pluses are the experimental results with vibrated stainless steel spheres in Ref. [2].The black dotted line is the Carnahan-Starling equation.
The particle configurations of our system are different from those of equilibrium hard disks and the previous experiment using vibrated stainless steel spheres Ref. [2], as can be seen from the curve of g(r = D)(ϕ) ≡ g D (ϕ) shown in Supplementary Fig. 1 with D being the particle diameter.The black dotted line refers to the theoretical curve of Carnahan-Starling equation g CS D (ϕ) = [1 − 7/16ϕ]/(1 − ϕ) 2 , which is derived from the equation of state for non-attracting hard spheres [3].
Assuming that the inelastic collisions between particles can be effectively represented by an interparticle attractive potential, the radial distribution function is given by g(r) = exp(−V (r)/(kT )) • g eq (r), where V (r) represents the effective potential and kT ≡ 1 2 mv 2 T can be estimated from the thermal velocity of the particles, denoted as v T .The definition of the thermal velocity is provided below.The equilibrium radial distribution function, g eq (r), for purely repulsive hard spheres is obtained from the BGY theory as described in Ref. [4].The interparticle force is given by dr , and its maximum value is determined by: The thermal velocity, v T , can be obtained by calculating the ensemble-averaged root mean square (RMS) velocity across different packing fractions, ϕ, using the formula: where N ϕ represents the total number of particles at a given packing fraction ϕ, and v i denotes the RMS velocity of a specific particle i.Here the thermal velocity v T = 1.4 ± 0.1D • s −1 is weakly affected by ϕ [5].
Thus, we have completed the calculation of the equivalent attractive force of inelastic collisions at different packing fractions.Next, we will introduce the calculation details of the Péclet number P e at different packing fractions.
In our experiment, the persistence time, denoted as τ p , can be determined from the cutoff time corresponding to the initial segment of superdiffusion.This segment is identified by the condition where the local slope of the mean square displacement (MSD) curve, represented by k * = d(log(M SD)) d(log(t)) in double logarithmic coordinates, is greater than 1, as illustrated in Supplementary Fig. 2. Considering the fluctuations of the measured slope k * in the experiment, where k * = 1 ± 0.1 within the time scale of the random walk, we can determine τ p as τ p = max(t(k * > 1.1)).This maximum value of t is depicted in Supplementary Fig. 3.
To determine the root mean square of the active velocity, denoted as v 0 , singleparticle experiments are conducted.The estimation of v 0 is obtained as follows:  3 The persistence time τ p .The stair-like structure is due to the finite time resolution of 0.025 s.Since the experimental frame rate is 40 frames per second, the time resolution is 0.025s.
Here, the notation ⟨• • •⟩ represents the average over repeated experiments, and • • • indicates the time average over the particle trajectory.
The persistence length of particles, denoted as v 0 τ p , quantifies the distance over which a particle maintains directional persistence.It represents the range within which its motion remains coherent before being randomized by random forces.And the Péclet number P e ≡ v 0 /(D r D) is viewed as the ratio between persistence length and particle size which in this case is D=16 mm.Here the persistence time τ p = 1/D r [1].Therefore, theoretically, P e can also be calculated through the definition.However, due to the higher spatial resolution precision in the experiment (0.29mm per pixel) compared to the temporal resolution precision (average particle displacement of 0.6mm within 0.25 seconds), the P e presented in the main text is obtained from the corresponding mean square displacements over the persistent time.
Consider the motion of an individual Brownian particle on an aluminum alloy plate.When the vibrational driving is stopped, the time required for the particle with an average velocity v 0 to decelerate to 0 due to friction is τ f = v 0 /(µg), where µ is the friction coefficient and g is the acceleration due to gravity.µ has been measured to be approximately 0.420.This measurement was obtained by determining the angle of repose through an experimental setup involving a gradual and continuous tilting of a slope from its horizontal position.The angle of repose is identified as the point where the particle exhibits incipient instability and starts to slide.Then we obtain: In our vibration experiments, the motion of particles can be classified into two types: (a) collisions with the vibrating platform, and (b) a projectile-like leap motion where particles momentarily detach from the platform.The inertia time, denoted as τ I , of an individual vibrating particle is determined by two factors: the time required for the velocity to be completely dissipated through multiple collisions with the platform, and the cumulative duration of all the leap periods between these collisions.Mathematically, it can be expressed as: Here, τ c is the contact time for mutual collisions of spheres with radius R, and τ h = 38ms denotes the average duration of a single particle's leap, which is measured using a high-speed camera.The contact time τ c can be calculated using a formula derived from the reference [6][7][8]: In our experiment, a single collision between our Brownian vibrators and the platform can be effectively approximated as a collision between one of the 12 inclined legs and the platform.Considering that the legs of the particles are cylindrical with a radius of 0.5mm, we can treat the collision as a mutual collision between spheres with a radius R=0.5mm for calculation purposes.
Our Brownian vibrators are 3d-printed using CBY-01 resin.JS-UV-CBY-01 Light Yellow High-Temperature Hard Material is an ABS-Like Stereoscopic Modeling Resin with precise and durable properties, which is used in solid-state laser light curing molding.At a temperature of 25 °C, the densityρ = 1.14g • (mm) −3 , Poisson's ratio: ν = 0.39, Young's modulus: Y = 2GPa, the viscosity is 400 mPa • s.In our 3dprinting process, we utilized a curing depth of 0.135 mm and a construction layer thickness of 0.05 mm.The collision velocity, denoted as v c , can be estimated as , where τ h is the average duration of a single particle's leap and g represents the acceleration due to gravity.The maximum velocity of the vibrating platform along the z-axis is 46.8mm • s −1 .Since the collision velocity term v c in the equation above is raised to the power of 1/5, the correction for the platform motion's impact on the collision velocity can be neglected in the calculation.
The magnitude of the equivalent dissipative force is given by |F d | ≡ mv 0 /τ I .After reaching a steady state, the energy dissipated by the particle is equal to the energy input from the active force.Therefore, the magnitude of the active force where γ 0 is the equivalent drag factor, given by m/τ I , and m is the mass of the particle.
Here the magnitude of active force γ 0 v 0 is generated through the friction between the particle and the platform due to the unique structure of the inclined legs of the particles, which is different with the simulation [1].However, this distinction does not significantly impact the relevant analysis presented in the main text, as it primarily involves comparing the dimensionless factors associated with the second and third terms.
By comparing the maximum equivalent attractive force and the active force, the critical value can be evaluated.
The dimensionless coefficient, denoted as λ, exhibits an expected increase with an increase in ϕ.This trend is consistent with the fact that the effective attraction force also increases with ϕ due to the increasing collision frequency.For instance, at ϕ = 0.270, λ is approximately 0.72, while at ϕ = 0.317, λ reaches around 1.04.These values indicate that when the maximum effective attractive force surpasses the active forces, collective behavior starts to emerge according the theoretical predictions [1].

The hexatic order parameter
The hexatic order parameter, ψ j 6 , quantifies the similarity between the structure surrounding a particle j and a perfect hexagonal arrangement.It is defined as ψ j 6 = 1 nj nj k=1 e i6θ jk , where n j represents the number of nearest neighbors of particle j, and θ jk denotes the angle between the vector ⃗ r jk ≡ ⃗ r k − ⃗ r j and the x-axis.Notably, n j is not restricted to six; it can take different values, such as 2 for a linear structure, 3 for a trident structure, 4 for a center particle with four neighbors, and so on.
To illustrate, let's consider the case where n j = 2.When the center particle and its two neighbors form a straight line, ψ j 6 equals 1.Consequently, ψ j 6 characterizes the local particle-scale symmetry.Supplementary Fig. 5 The hexatic order parameter.(a-f) The absolute value of hexatic order parameter ψ j 6 at various ϕ ′ s.The color bar represents the absolute value of hexatic order parameter.
The 2d maps displaying the absolute value of the hexatic order parameter, ψ j 6 , at various ϕ values are presented in Supplementary Fig. 5.Each panel in the figure corresponds to the same particle configurations as shown in Supplementary Fig. 4.
The global mean hexatic order parameter, , is calculated by averaging the absolute value of the hexatic order parameter for all particles over different configurations.Here, N represents the total number of particles, and ⟨...⟩ denotes the average over different particle configurations.
In Supplementary Fig. 6, we compare the global mean ψ global 6 observed in our experiment with the results from Ref. [2], particularly showing a plateau for packing fractions ϕ ≤ 0.436.Moving on to Supplementary Fig. 7, we present the two-dimensional correlations of the local hexatic order, denoted as g 2D 6 (r, θ) =  Supplementary Fig. 8 The PDFs of particle displacement at ϕ = 0.270, 0.555, 0.690.The PDFs of u α /std(u α ) (α refers to the x or y component) at the volume fractions of ϕ = 0.270, 0.555, 0.690.Different colors and markers represent various time intervals.The black solid lines refer to Gaussian distributions.(a-c) The probability density functions (PDFs) of u x and u y , representing the x and y components of the particle displacement vector ⃗ u, are presented in Supplementary Fig. 8 and Supplementary Fig. 9, respectively.In these figures, the horizontal coordinates are normalized by the corresponding standard deviation of particle displacement.The solid black lines in the figures depict Gaussian distributions for comparison.
Notably, none of the PDF curves exhibit a Gaussian shape.As u x /std(u x ) or u y /std(u y ) varies above and below zero, each curve progressively deviates from the Gaussian distribution.However, in most cases, except for ϕ = 0.753, the Gaussian distributions can still offer an approximate description of the PDFs, particularly for the central part of the distributions.It is primarily in the tails of the PDFs, where the probabilities are less than 10 −3 , that the Gaussian approximation becomes less accurate.
In the case of the volume fraction ϕ = 0.753, the PDFs of u x /std(u x ) and u y /std(u y ) exhibit the most significant deviations from Gaussian behavior at long time intervals.This deviation is attributed to the coexistence of a large number of nearly perfect crystal particles and a relatively small number of grain boundary particles, which leads to distinct features in the PDF distributions.Supplementary Fig. 9 The PDFs of particle displacement at ϕ = 0.713, 0.753, 0.832.The PDFs of u α /std(u α ) (α refers to the x or y component) at the volume fractions of ϕ = 0.713, 0.753, 0.832.Different colors and markers represent various time intervals.The black solid lines refer to Gaussian distributions.(a-c) u x /std(u x ) (d-f)u y /std(u y ).
The departure of the PDFs from Gaussian behavior can be quantitatively assessed using the kurtosis, as depicted in Supplementary Fig. 10.The kurtosis is defined as follows: where u i represents a random variable, which can refer to either u x /std(u x ) or u y /std(u y ).It is important to note that the kurtosis of a Gaussian distribution equals 3.

The average curl of the particle displacement field
We focus on the particle displacement field to quantify the large-scale collective motion.First, the particle density field can be coarse-grained with the kernel function Ψ(⃗ r) = e −(2⃗ r/D) 2 : Supplementary Fig. 10 The kurtosis of the PDFs of particle displacement.
The kurtosis of the PDFs of u x /std(u x ) (a) and u y /std(u y ) (b) for a series of time intervals and different packing fractions ϕ ′ s.
Next, the particle displacement field can be coarse-grained within a time interval ∆t: Here, the displacement of the particle i is ⃗ u i (∆t) = (⃗ r i (t 0 + ∆t) − ⃗ r i (t 0 )), and the collective behavior is evident in Fig. 1 (e) in the main text for ∆t = 100s.In a quasi-two-dimensional system, the average curl of the displacement field Ω is defined as: Here, the angle brackets represent the spatial average followed by a time average over different t 0 .

Important timescales in the experiment
Our experiment involves four crucial timescales: (1) The time interval of stochastic driving, estimated by the vibration frequency of the electromagnetic shaker, is approximately 1/f = 0.01 seconds.
(2) The persistence time, denoted as τ p , is obtained from the particle's mean square displacement (MSD) curve and is around 0.25 seconds for a single particle system.
(3) The particle's mean free time, denoted as t c , is illustrated in Supplementary Fig. 12.Here, t c (ϕ) is defined as the time when the MSD of particles at a given packing fraction ϕ reaches the mean separation distance of neighboring particles, which is equal to ( π/(2 √ 3ϕ) − 1)D.In the granular fluid phase (0.111 ≤ ϕ ≤ 0.270) and the flocking fluid phase (0.317 ≤ ϕ ≤ 0.690), which are the most interesting regimes in terms of physics, the mean free time between particle-particle collisions is much larger than 0.25 seconds.This indicates that stochastic driving randomizes the particle's velocities well before two adjacent particles interact.
(4) The period of collective motion around the system's center is approximately 4000 seconds.In the flocking fluid phases, the MSD peaks around 2000 seconds correspond to half a rotation period of the outermost particles' collective motion.These particles contribute the most to the MSD, as their displacement follows an approximately semicircular path from one side of the system to the other, covering a distance of about 30 particle diameters.Once the displacement reaches this magnitude, the MSD starts to decrease.

The spatial connected correlation function of the displacements
The spatial correlation function of the instantaneous velocities, as depicted in Supplementary Fig. 13(a), demonstrates that the velocity field lacks polar order in both the granular fluid phase and the flocking fluid phase.However, for dt=10s and 100s in Supplementary Fig. 13(c) and (d), the spatial correlation function of the displacements shows long-range spatial correlations at ϕ = 0.555, 0.713 and 0.832, confirming that there is no collective motion at short time scales, and the flocking fluid phase can only be observed at time scales greater than 10s.
Here the spatial connected correlation function of the displacements is defined as: , δu i and δu j are the displacement of the particles at positions i and j in the time interval of dt.
Supplementary Fig.4The Voronoi diagrams.(a-f) Voronoi diagrams at various ϕ ′ s.The color bar represents the order of the Voronoi cells.Supplementary Fig.4illustrates the Voronoi diagrams at different values of ϕ.In the granular fluid phase (ϕ = 0.270), the Voronoi cell order, represented by the number of edges, ranges from 4 to 9 and exhibits a random spatial distribution.In the flocking fluid phase (ϕ = 0.555), the order ranges from 4 to 8, and there is no noticeable localized ortho-hexagonal aggregation.In comparison, in the poly-crystalline phase (ϕ = 0.713), the majority of Voronoi cells are of order six, and the grain boundaries consist of pairs of grains with orders 5 and 7.This observation holds true for ϕ = 0.753 as well.With the exception of defects near boundaries, the crystalline region at ϕ = 0.832 comprises Voronoi cells of order six that span the entire system.By utilizing the area S v of Voronoi polygons, we can calculate the local packing fraction ϕ local = πD 2 4Sv and analyze its distribution and fluctuation.
et al. experimental resultsSupplementary Fig.6The global mean of ψ 6 at various ϕ ′ s.Red circles refer to Brownian vibrators.Blue pluses are the results of Ref.[2].
Supplementary Fig.13The spatial connected correlation function of the displacements.(a) The time interval dt=0.025s, while it's the same as the connected correlation function of the instantaneous velocities; (b-d) dt=1s,10s,100s. )