Granular avalanche statistics in rotating drum with varied particle roughness

: We experimentally investigate the avalanche statistics of dry granular materials in a slowly rotating drum for ﬁve types of beads with varied surface roughness. For all beads, two distinct angles, i.e., repose angle θ r and maximal angle θ m , can be clearly deﬁned, and the avalanche size distributions P ( δθ ) are Gaussian-like. θ r , θ m , and the span in P ( δθ ) are all positively correlated with bead surface roughness. This observation thus contrasts with a power-law P ( δθ ) predicted by self-organized criticality, but is reminiscent of a ﬁrst-order phase transition. We speculate that both the inertia e ﬀ ect and the velocity-weakening mechanism during an avalanche process can enhance the ﬁrst-order features, which are however absent in plasticity of sheared amorphous solids. We also discuss the dependence between θ r and θ m for various particles, as well as the correlation between starting and stopping angles for an individual avalanche.


INTRODUCTION
Dry granular materials can form mechanically stable piles with a finite slope θ, as a simple manifestation of their particular properties distinct from ordinary phases of matter [1][2][3][4]. As θ of a pile slowly increases (by tilting or adding new particles to the top for instance), the surface granular particles start to flow, i.e., an avalanche occurs, and eventually the system will settle into a new θ. One would naturally relate this process to the well-acknowledged self-organized criticality (SOC) based on a sandpile model [5,6], in the sense that the granular pile self-tunes its θ fluctuating around a potential steady-state value, i.e., a critical point. Specifically for a granular pile in rotating drum, under quasistatic rotation, θ(t) as a function of time behaves as a jagged curve, formed by different linearly increasing intervals punctured by instant drops δθ (i.e., avalanches). In a SOC picture, the avalanche size distribution P(δθ) should be a power law. Although this was partly reproduced in experiments [7,8], more evidence has refuted SOC for real granular piles since θ(t) displays quasiperiodic feature [9][10][11]. Two crucial angles can be defined statistically, i.e., maximal angle θ m and repose angle θ r , before and after an avalanche respectively, which has been documented in early researches [12,13]. Then one obtains a finite peak ∆θ = ⟨δθ⟩ = θ m − θ r in P(δθ). Hence the avalanche c ⃝ The Author(s) 2023. Published by Science Press and EDP Sciences. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited. behavior in granular piles is reminiscent of a first-order phase transition instead of a second-order (SOC) [14]. Comparing the process of increasing θ from θ r to θ m with that of cooling a liquid in equilibrium, one can associate θ r with the freezing point which signals the appearance of nucleus, i.e., locally unstable beads that tend to flow, and θ m with the melting point when all nuclei grow to dominate the system, i.e., an avalanche occurs. Alternatively, θ m can also be associated with a spinodal point when the situation θ r < θ < θ m is considered to be metastable. The incompatibility between first-order and SOC behaviors has been rationalized in different approaches taking into account the flowing dynamics for the surface particles upon the pile during an avalanche [15,16].
The macroscopic measure θ(t) as well as the debate on first-order or SOC transition can offer insights into vastly different systems far beyond a granular pile. The Coulomb criterion defines the effective friction coefficient µ = tan θ, and thus the fact θ m > θ r can be associated with the hysteresis in solid-fluid transition for stressed granular materials, i.e., one needs a larger shear stress to trigger the flow than to maintain the steady flowing state [2,3]. Similar velocity-weakening phenomenon is also found in solid friction [17], landslides [18], etc. On the contrary, the stick-slip process of sheared amorphous solids typically possesses a power-law distribution in avalanche size (i.e., stress drop) [19]. A simple model has been proposed to unify the above-mentioned two facets of avalanche statistics [20]. Understanding this issue would benefit the studies on friction law, yielding/plasticity of amorphous materials, and different geophysical processes.
Here we experimentally probe the granular avalanches in a rotating drum, focusing on the macroscopic θ in the quasistatic regime. Despite the easy access of this experiment set-up [2,11,21,22], few studies have clarified the underlying connection between particle properties, macroscopic rheology, and avalanche statistics. Moreover, it has been claimed that microscopic friction should be responsible for the macroscopic hysteresis [23][24][25], but systematic experiments on θ controlled by friction are lacking to our knowledge. Utilizing the 3D-print technique, we are able to tune the surface roughness of granular beads up to the large friction limit [26,27], and hence obtain θ for different types of beads with broadly spanning friction. From extensive experiments, we determine two distinct θ r and θ m for all five beads. Combining with the data in previous literature, we further show a robust linear dependency between θ r and θ m , irrespective of specific particle properties like surface roughness or shape. We demonstrate that here the avalanche dynamics are more reminiscent of a first-order phase transition and no signature of SOC behaviors is found, since P(δθ) are Gaussian-like for all systems, i.e., θ evolves quasiperiodically as rotated. The detailed correlation between starting and stopping angles for each avalanche event is also studied, which exemplifies the inevitable inertia effect for a particle-level avalanche process.

EXPERIMENTAL DETAILS
As shown in Figure 1A surface can be recorded when rotated. Figure 1B shows an example of the binarized image for the surface profile, from which the angle θ is measured by fitting a straight line to the surface margin. We confirm that this linear profile holds for all Ω covered here. As illustrated in Figures 1C-1G, we consider spherical beads with five different surface properties, namely steel, ABS, glass, 3D-printed plastic (3DP), and 3D-printed plastic with bumpy surface (BUMP) [26]. They all possess nearly 50:50 bi-dispersity with d = 5/6 mm. Specifically, the BUMP beads in Figure 1G can mimic particles possessing very large friction. Smaller glass beads with d = 2-4 mm are also included to investigate the finite size effect.

QUASISTATIC REGIME
As the drum angular velocity Ω increases, the granular dynamics will crossover from discrete avalanches to continuous flow [28][29][30]. Figures 2A-2C show the evolution of θ as a function of accumulated rotational angle A = Ωt for three typical Ω. We clearly observe the distinction between ( Figure 2A) the jagged curve and ( Figure 2C) the fluctuation in continuous flow, as well as ( Figure 2B) the intermediate state with moderate Ω.
Only for the small-Ω quasistatic regime, discrete avalanches of size δθ can be clearly defined by the instant drop in θ, i.e., δθ = θ 2 − θ 1 , as illustrated in Figure 2A. Statistically, the angles before and after an avalanche respectively define θ m and θ r , i.e., ∆θ = ⟨δθ⟩ = ⟨θ 2 ⟩ − ⟨θ 1 Notice in Figure 2 we sample θ with a constant ∆A = 0.1 • for all Ω, and thus both the average angle Θ = ⟨θ⟩ (in the quasistatic case Θ = (θ m + θ r )/2 by definition) and its standard deviation σ(θ) can fairly quantify the dynamic change as Ω grows. Figures 3A and 3B respectively show Θ and σ(θ) as functions of Ω for glass beads in four different sizes, in which the plateaus in Θ and σ(θ) end at different Ω. This is due to the intermediate state ( Figure 2B) when the jagged curve (i.e., avalanche regime) is occasionally interrupted by noises, also found in refs. [11,30]. For this state, Θ varies mildly while σ(θ) decays sharply. Hence the plateau in σ(θ) is the proper indicator for the quasistatic regime, which indicates the persistence of sequential avalanches. In Figure 3B, the locations of Ω when σ decays are varied for different systems. This can be rationalized by the dimensionless Froude number Fr, originally defined by the ratio between two characteristic velocities: v 1 = ΩD and v 2 = √ gD, respectively associated with the external rotation and gravity, i.e., Fr = v 1 /v 2 = Ω √ D/g. Here we note that v 1 should be replaced by the free surface velocity scale v * 1 . Since the depth for the surface flow in unit of d would not vary considerably when different beads are used, the downward flux scales as v * 1 d, which simply balances the upward flux caused by rotation ΩD 2 . This mass conservation gives v * 1 = ΩD 2 /d and thus the modified Froude number is Fr * = Fr · D/d. Figure  3C presents the normalized fluctuation σ/σ 0 versus Fr * , in which a clear transition at Fr * = 0.3 is observed, compatible with ref. [11]. As a physical interpretation, the boundary of quasistatic regime means that the time scale for particle restabilization under gravity is comparable to that for lifting a new particle to the top, i.e., v * 1 and v 2 are comparable.

REPOSE ANGLE AND MAXIMAL ANGLE
In the following, we focus on the measurements for the quasistatic regime. We emphasize that the distinction between θ r and θ m is not unique for the rotating drum set-up, but reflects the universal features of granular materials. By slowly adding grains to the peak of a heap in which no obstacle exists against the downward flow, clear intermittent avalanches with θ r and θ m can still be found [29,31] Θ is positively correlated with the particle size in Figure 3A, reflecting the finite size effect in accordance with refs. [11,21]. The ratio d/D naturally quantifies the finite size, given that d/W is sufficiently small. For the quasistatic regime in Figure 3A, we plot Θ as a function of d/D in Figure 4A, compared with other results for the similar glass beads [11,21]. Also included in the inset is the associated angle difference ∆θ = θ m −θ r . Different experiments share a similar trend, i.e., both Θ and ∆θ can be extrapolated to finite values as d/D → 0. We also notice that ∆θ is relatively robust against different experimental realizations. Physically, the surface properties of glass beads intrinsically determine both the effective friction and the angle hysteresis, which remain finite to the large system size limit. It is then interesting to do a quantitative survey using different types of beads.
At the fixed Ω = 0.4 • /s, we investigate the quasistatic avalanches for five types of beads in identical size 5/6 mm (see Figure 1). Their associated θ r and θ m are listed in Table 1. These angles would possibly shift to lower values if smaller beads are applied due to the finite size effect displayed in Figure 4A. But, considering the identical d/D, the order in θ r or θ m already indicates the levels of surface roughness for different types of beads. Figure 4B shows θ m versus θ r for these five beads, along with the data of small glass beads in Figure 4A and granular particles varied in shape or surface properties collected from previous literature [11,13,21,31]. A decent linear relation between θ r and θ m emerges irrespective of particle properties and experimental details, which has been noted in ref. [21]. Relatively, particles with notable asphericity or surface roughness possess large θ r and θ m . Furthermore, this relation is weakly ruined by the finite size effect, i.e., our data for different beads of d = 5/6 mm, d/D ≈ 0.03 (solid cubes) slightly deviate from the master curve, while those for small glass beads d/D < 0.02 (empty cubes) do not. After checking d/D for all the data in Figure 4B, we choose [θ r , θ m ] with d/D < 0.02 to conduct the linear fit (straight line), i.e., θ m = 1.256θ r − 3.699. This interesting observation indicates that one can inspect both the shape and surface properties of granular particles in a unified approach.  [13] Liu et al. [21] Borzsonyi et al. [31] Balmforth et al. [11] Figure 4C replots θ m as a function of ∆θ. θ m (or θ r ) are positively correlated with ∆θ for our five beads (solid cubes), supporting the view that both the effective friction coefficient and the angle hysteresis are simply monotone increasing functions of microscopic friction coefficient (or surface roughness) for different beads. This can also be inferred from the linear fit in Figure 4B, i.e., by simple algebra one obtains θ m = 4.906∆θ + 14.449 (grey line in Figure 4C). According to refs. [23][24][25], microscopic friction is responsible for a finite angle hysteresis ∆θ, and thus the intercept in Figure 4C refers to the finite effective friction when particles are microscopically frictionless. Such reference agrees with that a macroscopic friction exists even for a granular solid composed of frictionless particles due to its particulate nature. However, the value θ r = θ m ≈ 14 • is larger than the prediction ∼ 6 • obtained from the simulation of quasistatically sheared frictionless beads [23]. This accurate comparison might suffer from including those aspherical particles for the fit in Figure 4B. Also, the experimental uncertainty in ∆θ is not negligible.

AVALANCHE STATISTICS
Beyond the sample-averaged measurement, i.e., θ r and θ m , as discussed above, here we carefully inspect the statistics of individual avalanche events. Figure 5A presents the avalanche size distribution P(δθ) for our five beads, where the five numbers in legend refer to the total numbers of included avalanches for different beads. These five curves share a similar shape, that peaks at a finite δθ, i.e., approximately the average ∆θ = ⟨δθ⟩ considering the probable non-Gaussianity. In Figure 5B, the distributions of the normalized variables ε = (δθ − ∆θ)/σ(δθ) for five systems nearly collapse on the normal distribution (curve), and only mild deviations especially at small-ε regime are observed. One probable avenue to inspect the small-δθ regime in Figure 5A is to subtract the Gaussian part from the whole distribution. The remaining part could be associated with SOC behaviors. Considering that the resolution in δθ is estimated by d/D [14], i.e., about 1.7 • , and that δθ spans smaller than a decade, avalanche dynamics here displays no meaningful features of SOC. In other words, recalling Figure 2A, the Gaussian-like P(δθ) indicates a quasiperiodic evolution in θ with the quasistatic rotation. Combining with the [θ r , θ m ] data reported above, we obtain a coherent picture that a bead with larger roughness is associated with larger θ r , θ m , ∆θ, as well as a wider span in P(δθ).
As pointed out in ref. [14], the distinguishable θ r and θ m are reminiscent of the freezing and spinodal points for the supercooled (metastable) branch in a phase diagram. If this view holds accurately, the instant starting and stopping angles θ 2 and θ 1 should both be narrowly distributed, and their intercorrelation during a single avalanche should be negligible. Figures 6A and 6B respectively show P(θ 1 ) and P(θ 2 ), in which one recognizes that both angles possess Gaussian-like distributions with similar span for each bead. And, P(θ 1 ) and P(θ 2 ) merely show tiny overlaps for any bead, suggesting that our system size is thermodynamically reliable. Like ε in Figure 5B, in Figure 6C, the normalized variables ε 1 for θ 1 (or ε 2 for θ 2 ) are well characterized by the normal distribution. Notice that mild deviations are found for θ 1 (solid points), especially for BUMP.
We further check the correlation between θ 1 and θ 2 for each avalanche event, by calculating the conditional average ⟨θ 1 |θ 2 ⟩ as a function of θ 2 . In order to include all beads, we extract the averages θ r and θ m from the results, as shown in Figure 6D (y-label is simplified as θ 1 − θ r ). We find that θ 1 are negatively correlated with θ 2 , i.e., a larger instant θ 2 tends to be followed by a smaller θ 1 . The consequent δθ = θ 2 − θ 1 is positively correlated with θ 2 , similarly shown in Figure 6E. This consistency can also be proven by the relations θ 1 − θ r ≈ −0.4(θ 2 − θ m ) and δθ − ∆θ ≈ 1.4(θ 2 − θ m ), marked by the straight lines in Figures 6D and  6E. It is reasonable to speculate that the weak correlation between θ 1 and θ 2 is caused by the specific flowing dynamics, since a larger θ 2 accompanies with stronger inertia effect such that the system rests at a smaller θ 1 after an avalanche. Accordingly, the span (or standard deviation) of P(δθ) is also larger than that of P(θ 1 ) or P(θ 2 ) because of the above correlation, shown in Figure 6F. Also included are the results for glass beads with d = 2 mm (solid points), and one observes that all three standard deviations approach zero as ∆θ → 0. Thus, as bead roughness declines and system size is enlarged, distributions for θ 1 , θ 2 , and δθ gradually evolve to delta function with zero average, supplementing

DISCUSSION
In this work, we experimentally investigate the avalanche statistics in a quasistatically rotating drum for five types of beads with varied surface roughness (or friction coefficient as the ideal equivalence). The evolution of granular pile slope θ allows one to define the starting and stopping angles θ 2 , θ 1 , before and after an avalanche event, respectively. In the average meaning, both the repose and maximal angles (θ r = ⟨θ 1 ⟩ and θ m = ⟨θ 2 ⟩), as well as their difference ∆θ which quantifies the hysteresis, are all positively correlated with the bead roughness. One could further infer the extreme situation for frictionless granular particles, possessing no hysteresis but a finite slope (i.e., macroscopic effective friction). The avalanche size (i.e., δθ = θ 2 − θ 1 ) distributions P(δθ) are Gaussian-like for all beads, indicating that θ evolves quasiperiodically and thus the avalanche dynamics in rotating drum is a first-order transition rather than SOC irrespective of specific bead roughness. We also show the weak correlation between θ 1 and θ 2 (also δθ), i.e., a larger θ 2 tends to be followed by a smaller θ 1 , suggesting the inertia effect. The original sandpile model for SOC [5] has not included the inertia effect involved in a real granular avalanche. Also, the velocity-weakening in either microscopic or macroscopic (effective) friction [15] should be taken into account for an initiated flow. Therefore, for our rotating drum, once θ m is exceeded, an avalanche occurs and both the above factors (they are not necessarily independent) lead to a restabilized angle smaller than θ m , i.e., θ r which even strongly depends on θ m . This contrasts with a single critical point as expected based on SOC. However, by slowly adding particles to a pile, SOC behaviors can still be reproduced since now avalanche is achieved not catastrophically but in a quasistatic manner [8], i.e., a clear θ r is forbidden. Alternatively, introducing mechanical perturbation would destabilize a pile before θ m is reached, also resulting in a SOC domain [32]. All these considerations indicate that the competition between SOC and first-order transition is considerably affected by the driving protocols for granular avalanches on the particle level. On the contrary, the stress avalanche statistics in quasistatically sheared amorphous solids display notable agreement with the SOC prediction [19]. This indicates that avalanches on particle and contact levels should be treated differently. A unified framework to explain the above rich phenomena is still a hot topic [20]. Any progress in this direction will benefit the understanding of granular constitutive law, plasticity of amorphous solids, and related geophysical processes like landslides and earthquakes.

Author contributions
Y.W. conceived the concept. A.P. carried out the experiments. A.P., Y.Y., and Y.W. analyzed the data and wrote the paper.