# Rheology of dense vibrated granular flows: non-monotonic response controlled by granular temperature

Andrea Plati, Giuseppe Petrillo, Lucilla de Arcangelis, Andrea Gnoli, Andrea Puglisi, Alessandro Sarracino, E. Lippiello

► To cite this version:

Andrea Plati, Giuseppe Petrillo, Lucilla de Arcangelis, Andrea Gnoli, Andrea Puglisi, et al.. Rheology of dense vibrated granular flows: non-monotonic response controlled by granular temperature. 2026. [hal-05442832]

HAL Id: hal-05442832

https://hal.science/hal-05442832v1

Preprint submitted on 5 Jan 2026

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

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

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

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# Rheology of dense vibrated granular flows: non-monotonic response controlled by granular temperature

A. Plati,¹ G. Petrillo,² L. de Arcangelis,³ A. Gnoli,⁴,⁵,⁶ A. Puglisi,⁴,⁵,⁶ A. Sarracino,⁷ and E. Lippiello⁸

¹Université Paris-Saclay, CNRS, Laboratoire de Physique des Solides, 91405 Orsay, France

²The Earth Observatory of Singapore, Nanyang Technological University, Singapore

³Department of Mathematics and Physics, University of Campania "Luigi Vanvitelli", 81100 Caserta, Italy

⁴Institute for Complex Systems CNR, P.le Aldo Moro 2, 00185, Rome, Italy

⁵Department of Physics, University of Rome "La Sapienza", P.le Aldo Moro 5, 00185 Rome, Italy

⁶INFN, Sezione Roma2, Via della Ricerca Scientifica 1, 00133, Rome, Italy

⁷Department of Engineering, University of Campania "Luigi Vanvitelli", 81031 Aversa (CE), Italy

⁸Department of Mathematics and Physics, University of Campania "Luigi Vanvitelli", 81100 Caserta, Italy

(Dated: November 20, 2025)

We study the rheology of dense granular materials subjected to vertical vibration using numerical simulations of a stress-imposed vane rheometer. The effective viscosity increases with confining pressure, decreases with vibration amplitude, and exhibits a non-monotonic dependence on frequency: weakening is observed at intermediate frequencies but is lost at high frequencies. We show that the rheological response is governed by the balance between grain-scale agitation energy and the stabilizing effect of confinement. This framework reconciles previously observed trends in viscosity and friction weakening and emphasizes the central role of energy injection and dissipation in determining granular flow properties under vibration.

## I. INTRODUCTION

Granular materials are ubiquitous in natural environment and applications [1, 2]. Depending on the external forcing, they can exhibit both solid-like or fluid-like behaviour [3, 4]. The rheology of granular materials in shear cells or flowing over inclined planes has been extensively investigated [5–11]. At the same time, the role of vibration in determining granular flow properties has been highlighted in several theoretical, numerical and experimental studies [12–22]. External vibrations usually increase the flowability of a granular system and this can have dramatic consequences as in the case of earthquake-triggered landslides [18, 23] or can be exploited to facilitate the manipulation of granular systems in applications [24]. In laboratory experiments, vertical vibration is known to reduce to zero the yield stress and produce non-monotonic flow curves [13]. It may also induce a transition from thinning to collisional thickening [16].

The general setting for experiments or simulations to probe the rheology of granular materials under sinusoidal vibration considers the role of three main parameters: the vibration amplitude *A*, the vibration frequency *f* and the confining pressure *p*. Depending on the specific setup, different quantities can be used as proxies of the “resistance to flow” of the granular system. When the shear rate $\dot{\gamma}$ is imposed, one usually measures the macroscopic friction coefficient $\mu = \tau/p$ where $\tau$ is the resulting shear stress. Conversely, if the shear stress is imposed, the flow is characterized by the effective viscosity $\eta = \tau/\dot{\gamma}$. Non-vibrated granular systems have been extensively studied by focusing on the so-called flow curves $\mu(\dot{\gamma})$ and $\eta(\tau)$ [7, 10, 25], however - when vibrations are taken into account - the goal is to characterize the flowability of the material either through $\mu(\dot{\gamma}, A, f, p)$ or $\eta(\tau, A, f, p)$. To restrict the parameter space, it is sometimes useful to con-

sider $\dot{\gamma}$ and $\tau$ fixed, focusing on $\mu(A, f, p)$ and $\eta(A, f, p)$. In recent years, many experimental and numerical observations on these quantities have been reported [17, 19–21].

In Ref. [17], the effective viscosity of a dense, weakly fluidized granular system was investigated experimentally and numerically by varying the driving amplitude *A* and frequency *f* under a fixed confining load. It was found that $\eta$ decreases monotonically with *A*, whereas its dependence on *f* is non-monotonic: the viscosity decreases until a minimum at a frequency $f^*$ that shifts to lower values as the dissipative properties of the medium increase. A subsequent work [19] showed that these trends are linked to the behavior of the average kinetic energy of the system (i.e. the granular temperature, see Sec. IIA for its definition), which exhibits a maximum at $f^*$. In Fig. 1, we show the granular temperature as a function of *f* and *A* obtained extending the numerical dataset of Ref. [19]. It increases monotonically with *A* when the frequency is fixed and exhibits a non-monotonic trend with *f* when the amplitude is fixed, which correlates with the behavior of the effective viscosity obtained in experiments [17]. These results suggest that rheology is primarily controlled by the efficiency with which the medium absorbs energy from external vibration. As discussed in [19], while the non-monotonic dependence on *f* resembles a resonance phenomenon, the fact that $f^*$ depends on dissipation rules out a standard resonance, in which the peak position is fixed by the elastic properties of the system. The proposed interpretation instead attributes the non-monotonicity to the competition between the energy injected by the vibration mechanism, and the dissipation rate of the medium (which is related to the grain-grain collision rate). Indeed, both these competing effects are enhanced by the driving frequency, producing a non-monotonic trend.
<figure><img src="image_2.png" /><figcaption>FIG. 1. Granular temperature as a function of the oscillation amplitude *A* and frequency *f*. The granular temperature increases monotonically as a function of *A* at any fixed frequency. Meanwhile, for any fixed value of *A*, it exhibits a non-monotonic behavior, reaching a maximum at an intermediate frequency. This plot is obtained by extending the numerical results of Ref. [19] to a broader range of *A* and *f* values.</figcaption></figure>

More recently, Clark et al. [20] examined the friction coefficient $\mu(A, f, p)$ and observed a similar non-monotonic trend with frequency. However, rather than focusing on the frequency dependence, they analysed the behavior of the minimum friction coefficient $\mu_{min}(A, p)$ which was found to be controlled by the variable $A^2/p$. This scaling was motivated by the idea that friction weakening requires vibration amplitudes large enough to disrupt the underlying force network: while shaking at sufficiently high acceleration is enough to break individual contacts, an additional threshold in $A^2/p$ is necessary to destabilize mesoscale structures such as force chains.

In our study we go beyond this interpretation, based on the evidence that, regardless of whether the thresholds in shaking acceleration and $A^2/p$ are exceeded, friction weakening is always lost at sufficiently high $f$. Following the interpretation of [19], we highlight that the driving frequency governs the efficiency of energy transfer into the granular medium, which in turn controls the system's rheological properties. The centrality of energy transfer has also been evidenced in more recent studies, focusing on local rheological responses [21]. We therefore emphasize that any comprehensive description of friction weakening must include a detailed analysis of a wider frequency range.

In this article, we build on these findings to provide a unified framework for understanding the interplay between dissipation, energy injection, and rheological response in realistic vibrated granular systems. In Sec. II we describe our numerical setup which is a realistic *in silico* reproduction of the experimental system for vane rheometry studied in [17, 19]. In the same section we also provide an overview of the rheological response for different shaking parameters and confining pressures. This is done by discussing the behavior of the effective viscosity

$\eta(A, f, p)$. This analysis is complemented by the introduction of a minimal model (detailed in Appendix C), which accounts for the combined effect of shaking amplitude and confining pressure. In Sec. III, we show that the key variable controlling $\eta$ involves the ratio between the average kinetic energy of the grains and the confining pressure. There we also discuss a physical interpretation of these results based on the competition between two energy scales, the first related to grain motion and second related to the stability of local structures. In Sec. IV, we conclude our paper by discussing the link between our results and previous studies; we also provide a brief outlook for future developments of this work.

## II. EFFECTIVE VISCOSITY OF THE GRANULAR MEDIUM

Our numerical setup is the *in silico* reconstruction of the experimental apparatus used in [17, 19] and it has been shown to provide results in agreement with experiments in previous studies [19, 26]. We use the Discrete Element Method (DEM) [27] to perform numerical simulations of a monodispersed granular system with $N = 2600$ particles of radius $r = 2$ mm and mass $m = 0.27$ g confined in a cylinder with a conical-shaped base and covered by a circular lid of mass *M* and radius *R* = 45 mm, which matches the area of the cylinder's round face. The whole system feels gravity acceleration *g*. A rectangular-shaped vane, able to rotate around a fixed vertical axis, is immersed in the granular medium and subjected to a constant torque $T = 0.006$ Nm. The granular material is subjected to a sinusoidal vertical displacement $\Delta z(t) = A \sin(2\pi ft)$ which is imposed on the system's container. The torque value is chosen so that the vane does not move in the absence of vibration. We implement our numerical simulations through the LAMMPS package [28–30] which uses the Hertz-Mindlin contact model [30, 31] to describe grain-grain, grain-vane, grain-lid and grain-wall interactions. More technical details about our numerical methods are provided in Appendix A.

### A. Non-monotonic behavior as a function of the shaking frequency

Following a microrheology approach, we measure the angular velocity $\Omega(t)$ of the vane as a function of time and explore its behaviour varying the frequency $f$ and the amplitude $A$ of the imposed shaking. In view of this, we introduce the rescaled acceleration $\Gamma = (2\pi f)^2 A/g$.

We have checked that, for all the explored combinations of $A$ and $f$, our system always lies in a weakly fluidized regime which means that grains move inside the cage formed by the surrounding particles without experiencing substantial rearrangement except for that induced by the vane's motion. Once the vibration and the torque
<figure><img src="image_3.png" /><figcaption>FIG. 2. a) Effective viscosity as a function of the rescaled acceleration (varied through f) for three shaking amplitudes. Experimental data from [17]. b) Analogous results by numerical simulations. Each point is obtained by averaging over 3 statistically independent runs. c) Kinetic energy of the granular medium as a function of Γ (varied through f) for three shaking amplitudes (simulations performed without the rotating vane). d) Effective viscosity as a function of the granular temperature in numerical simulations. The decreasing trend is well approximated by the scaling $\eta \propto K^{\alpha}$ with $\alpha = -2$.</figcaption></figure>

are switched on, after a $t_0$-long initial transient, the vane reaches a stationary state with an average velocity $\langle\Omega\rangle = 1/T \int_{t_0}^{t_0+T} \Omega(t)dt$ which can be thought as a proxy of the shear rate. For all our simulations we have $t_0 = 13$ s and $T = 33$ s. We define the effective viscosity of the granular medium as $\eta = T/|\langle\Omega\rangle|$. The choice of this variable is motivated by our previous experiments [17] showing a non-monotonic behaviour of $\langle\Omega\rangle$ and in turn of $\eta$ as a function of the frequency. In Fig. 2a-b we provide a comparison between experimental and numerical data exhibiting the resulting non-monotonic dependence of $\eta$ as a function of $\Gamma$ (varied through $f$) for different $A$. As highlighted by the vertical lines in the plots, the decreasing branch of the $\eta(\Gamma)$ curves can be explained in terms of the fluidization threshold around $\Gamma = 1$ [17], where the shaking acceleration becomes equal to $g$. The physical mechanism underlying the viscosity strengthening at larger $\Gamma$ has been discussed in a previous work [19], where we showed that it is linked to the behavior of the average kinetic energy of the granular medium in the steady state $K = 1/T \int_{t_0}^{t_0+T} K_{\text{inst}}(t)dt$, where $K_{\text{inst}}(t) = \frac{m}{2N} \sum_i v_i^2(t)$ and $v_i(t)$ is the instantaneous speed of the grain $i$. The variable $K$ is usually called granular temperature (or directly proportional to it through the number of degrees of freedom which is constant here) and displays itself a non-monotonic dependence on the vibration frequency as shown in Fig. 2c. This phenomenon originates from the competing effects of injected vibration energy and (inelastic) grain-grain collision rate which both increase with $f$ [19]. A more detailed numerical analysis supporting this interpretation based on energy dissipation can be found in Appendix B. The behavior of $\eta$ as a function of the driving parameters $A$ and $f$ is then effectively

controlled by $K$, as highlighted in Fig. 2d. We point out here that to investigate the effect of the granular temperature on $\eta$ without considering trivial correlations due to the blade motion, we measured $K$ in simulations without the rotating blade.

## B. Effect of normal pressure

We now focus on the behaviour of $\eta$ as a function of $\Gamma$ (varied through $f$) for different couples of $\{A, p\}$ where $p = Mg/(\pi R^2)$ is the normal pressure due to the lid's weight. In Fig. 3a, we observe that all curves decrease up to the fluidization threshold and then reach a minimum $\eta_{\min}$ (highlighted with a larger mark) at $\Gamma = \Gamma_m$ before showing a gradual increase. As expected, increasing $A$ reduces the effective viscosity while raising $p$ has the opposite effect. In order to first focus on the most fluidized regime, we get rid for the moment of the $\Gamma$ dependence considering the behavior of $\eta_{\min}(A, p)$. In Fig. 3b, we show that $\eta_{\min}(A, p)$ is mainly controlled by the variable $A^2/p$. Remarkably, a similar data collapse has been previously observed also for the macroscopic friction coefficient in simulations of 2D sheared granular materials under vertical vibrations [20]. Our results generalize this observation to the effective viscosity of a 3D granular system probed through stress-imposed microrheology. The role of the variable $A^2/p$ was explained in Ref. [20] based on physical mechanisms involving the rupture of force chains in the bulk of the granular medium. In Appendix C we discuss a simple mechanical model which is able to reproduce the $A^2/p$ scaling and allows us to trace it back to the interplay between energy injection and non-linear elastic deformation. This shows that force chain
breaking is not necessary to explain the observed scaling.

### III. A UNIFIED PICTURE

Now we point out that the observation of a non-monotonic behavior of $\eta$ as a function of the shaking frequency (Sec. II A) is in striking contrast with the idea that the rheological response of the granular medium is mainly controlled by $A^2/p$. While the viscosity weakening branch can be understood by combining the acceleration threshold $\Gamma = 1$ and the $A^2/p$ scaling, the raising branch $\eta(\Gamma \ge \Gamma_m)$ is not. Indeed, simulations performed at $\Gamma \ge \Gamma_m$ (varied through $f$), with the same amplitude and pressure can manifest effective viscosities that vary over different orders of magnitude, see Fig. 3c. At the same time, we also verified that the kinetic energy of the system does not substantially depend on $p$ (not shown), thus the effect of $p$ on the effective viscosity (see Fig. 3a) can't be explained with a dependence of $K$ on $p$. Based on our numerical analysis, we define an energy scale related to the confining pressure $K_p = p\pi r^3$ (i.e. the work necessary to move of a distance $r$ against a constant force $p\pi r^2$) and propose the energy ratio $\bar{K} = K/K_p$ as the main variable controlling the effective viscosity of the system. In Fig. 3d, we plot $\eta(\Gamma \ge \Gamma_m)$ obtained for all the compatible frequencies, amplitudes and pressures varied in our analysis, showing a good collapse on the master curve $\eta \propto \bar{K}^\alpha$ with $\alpha \sim -2$.

The adimensionalized energy $\bar{K}$ can be interpreted as a competition between the typical agitation energy of the grains inside their cage $K$ and the typical energy (set by the confining pressure) needed to escape from such a cage $K_p$. As visible from the x-axis of Fig. 3d, our simulations are always in the regime $\bar{K} \ll 1$ where we expect spontaneous grain rearrangement to be extremely rare. The main idea behind a $\eta(\bar{K})$ dependence is that, even in this weakly fluidized regime, the larger $\bar{K}$, the easier it will be for the rotating blade to break grain cages in order to move.

We point out that the ratio of energies defining $\bar{K}$ reminds of the competition between thermal agitation and activation energy barriers that characterizes the viscosity divergence in glass-forming supercooled liquids [32] with, however, the important difference that in these systems viscosity divergences are usually exponential or even steeper while in our case we observe a power law. This qualitative analogy with thermal systems is nevertheless conceptually useful to understand the relationship between rheological response and $\bar{K}(A, f, p)$ in our granular setup. A peculiarity of thermal systems is that temperature can be directly controlled by equilibrating the system with a thermal bath. Then, as a consequence of thermodynamic equilibrium, also the average kinetic energy of the system is fixed and proportional to the temperature. Instead, in a granular system, the average kinetic energy $K(A, f)$ depends non-trivially on the driving parameters. In other words, the relationship be-

tween injected and adsorbed energy is much more complex when thermodynamic equilibrium is violated, as in granular systems. Therefore, to understand how the driving parameters $A$ and $f$ control the rheology of a granular system it is important to first clarify how they act on its granular temperature $K$. Our numerical analysis shows that, since $\eta$ depends non-monotonously on $f$, there is no simple scaling in the form $\propto A^\beta f^\lambda$, with $\beta$ and $\lambda$ arbitrary exponents, that could predict the behavior of the effective viscosity. Panel c and d of Fig. 3 show in fact that $\eta$ depends on $A$ and $f$ through $K(A, f)$, which embodies the non-monotonic dependence on $f$.

### IV. CONCLUSIONS AND LINK WITH PREVIOUS STUDIES

In this paper we investigated the rheology of a dense vibrated granular system by means of realistic numerical DEM simulations of a stress imposed vane rheometry setup. We varied the vibration amplitude $A$ and frequency $f$ of the vertical sinusoidal shaking as well as the confining pressure $p$. We found that the variable controlling the effective viscosity is $\bar{K} = K/K_p$ where $K(A, f)$ is the typical agitation energy of a grain in the cage formed by the surrounding particles and $K_p \propto p$ is an energy scale set by the confining pressure that could be interpreted as the energy needed to escape from such a cage. Our proposed scaling takes into account the non-monotonic effect of $f$ on $\eta$ through the energy transfer mechanism discussed in Ref. [19] and Appendix B.

We now discuss the connection between our study and previous works on the rheology of vibrated granular materials in simplified geometries. We start referring to the numerical study reported in [20] which focuses on the macroscopic friction coefficient $\mu$ of a granular system subjected to a vertical sinusoidal shaking (with frequency $f$ and amplitude $A$) in a double plate geometry under a confining pressure $p$. Remarkably, similarly to the effective viscosity of our system, the frequency dependence of the friction coefficient $\mu(f)$ reported in this study exhibits a non-monotonic trend. When considering the full $\mu(A, f, p)$ dependence, the authors showed that the friction weakening behavior can be predicted by considering the region of control parameters identified by $\Gamma > 1$ and $A^2/p > a$ where $a$ is a threshold that depends on the restitution coefficient in the grain-grain interaction. However, the provided criterion only holds in the frequency domain where $\mu(f)$ decreases, thus not including in the analysis its non-monotonic behavior. We want to stress here that our study marks a step forward with respect to Ref. [20] for two key points. First, when the frequency dependence is not taken into account, we showed that the variable $A^2/p$ proposed for the friction-weakening behavior also controls the effective viscosity in our more realistic setup. Moreover, we have also proposed a minimal model for such a scaling in Appendix C. Second, when the full frequency dependence is taken into
<figure><img src="image_4.png" /><figcaption>FIG. 3. a) Effective viscosity as a function of Γ (varied through f) for two couples of amplitudes and pressures. b) Minimum effective viscosity as a function of A²/p for four pressures. c) Effective viscosity in the fluidized state (Γ > Γm) as a function of A²/p for four pressures and amplitudes. d) Effective viscosity in the fluidized state (Γ > Γm) as a function of K̄ ∝ K/p for four pressures and amplitudes. Diamond, square, upper- and lower-triangle markers correspond to A = 6, 14, 26, 58 µm respectively. Red, blue, green, yellow markers correspond to p = 31, 92, 308, 924 Pa respectively.</figcaption></figure>

account, we showed that it is not possible to predict the viscosity reduction by simply considering some thresholds on A²/p and Γ because it depends on the control parameters through K(A, f), which is the key variable embodying the non-monotonic dependence on f. We believe that the same idea could be used to generalize the analysis on the friction coefficient done in Ref. [20] by taking into account the non-monotonic µ(f) dependence. It is also important to point out that the relationship between granular temperature and rheology of vibrated granular materials was already introduced (although with a less systematic analysis) in our previous study [19] and later in Ref. [21] by Irmer et al. that showed how the granular temperature controls the local rheology in a double plate geometry. However, we point out that the latter study has been performed in a more idealized numerical setup without a direct experimental counterpart and without any specific mention of non-monotonic behaviors of rheological observables and granular temperature.

Finally, inspired by seminal studies on non-vibrated

granular flows [7], as an outlook for future development of this research line we highlight the possibility of studying the combined effect of vibration driving parameters and confining pressure also on the volume fraction of vibrated granular systems.

ACKNOWLEDGMENTS

A. Plati acknowledges funding from the Agence Nationale de la Recherche (ANR), France, grant ANR-21-CE06-0039. AG, AP and AS acknowledge funding from the Italian Ministero dell'Università e della Ricerca under the programme PRIN 2022 ("re-ranking of the final lists"), number 2022KWTEB7, cup B53C24006470006. EL acknowledges funding from the Italian Ministero dell'Università e della Ricerca under the programme PRIN 2022 PNRR P202247YKL, CUP B53D23033430001.

[1] B. Andreotti, Y. Forterre, and O. Pouliquen, *Granular media: between fluid and solid* (Cambridge University Press, 2013).

[2] A. Puglisi, *Transport and fluctuations in granular fluids: From Boltzmann equation to hydrodynamics, diffusion and motor effects* (Springer, 2014).

[3] H. M. Jaeger, S. R. Nagel, and R. P. Behringer, *Reviews of modern physics* **68**, 1259 (1996).

[4] P. Eshuis, K. Van Der Weele, D. Van Der Meer, R. Bos, and D. Lohse, *Physics of Fluids* **19** (2007).

[5] R. A. Bagnold, *Proceedings of the Royal Society of London. Series A. Mathematical and Physical Sciences* **225**, 49 (1954).

[6] P. Jop, Y. Forterre, and O. Pouliquen, *Nature* **441**, 727 (2006).

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

[8] F. Da Cruz, S. Emam, M. Prochnow, J.-N. Roux, and F. Chevoir, *Phys. Rev. E* **72**, 021309 (2005).

[9] M. Bouzid, M. Trulsson, P. Claudin, E. Clément, and B. Andreotti, *Phys. Rev. Lett.* **111**, 238301 (2013).

[10] A. Fall, G. Ovarlez, D. Hautemayou, C. Mézière, J.-N.
Roux, and F. Chevoir, Journal of rheology **59**, 1065 (2015).

[11] S. Deboeuf and A. Fall, Journal of Rheology **67**, 909 (2023).

[12] H. Melosh, Nature **379**, 601 (1996).

[13] J. A. Dijksman, G. H. Wortel, L. T. Van Dellen, O. Dau-chot, and M. Van Hecke, Phys. Rev. Lett **107**, 108303 (2011).

[14] G. H. Wortel, J. A. Dijksman, and M. Van Hecke, Phys. Rev. E **89**, 012202 (2014).

[15] F. Giacco, L. Saggese, L. de Arcangelis, E. Lippiello, and M. Pica Ciamarra, Phys. Rev. Lett. **115**, 128001 (2015).

[16] A. Gnoli, A. Lasanta, A. Sarracino, and A. Puglisi, Scientific Reports **6**, 38604 (2016).

[17] A. Gnoli, L. De Arcangelis, F. Giacco, E. Lippiello, M. P. Ciamarra, A. Puglisi, and A. Sarracino, Phys. Rev. Lett **120**, 138001 (2018).

[18] J. Léopoldès, X. Jia, A. Tourin, and A. Mangeney, Phys. Rev. E **102**, 042901 (2020).

[19] A. Plati, L. de Arcangelis, A. Gnoli, E. Lippiello, A. Puglisi, and A. Sarracino, Phys. Rev. Research **3**, 013011 (2021).

[20] A. H. Clark, E. E. Brodsky, H. J. Nasrin, and S. E. Taylor, Phys. Rev. Lett. **130**, 118201 (2023).

[21] M. G. Irmer, E. E. Brodsky, and A. H. Clark, Phys. Rev. Lett. **134**, 048202 (2025).

[22] O. D'Angelo, M. Sperl, and W. T. Kranz, Phys. Rev. Lett. **134**, 148202 (2025).

[23] M. García-Rodríguez and J. Malpica, Natural Hazards and Earth System Sciences **10**, 1307 (2010).

[24] C. Wang, M. Moharekpour, Q. Liu, Z. Zhang, P. Liu, and M. Oeser, Construction and Building Materials **289**, 123164 (2021).

[25] F. Da Cruz, F. Chevoir, D. Bonn, and P. Coussot, Phys. Rev. E **66**, 051305 (2002).

[26] A. Plati, A. Baldassarri, A. Gnoli, G. Gradenigo, and A. Puglisi, Phys. Rev. Lett. **123**, 038002 (2019).

[27] P. A. Cundall and O. D. L. Strack, Géotechnique **29**, 47 (1979), https://doi.org/10.1680/geot.1979.29.1.47.

[28] S. Plimpton, Journal of Computational Physics **117**, 1 (1995).

[29] A. P. Thompson, H. M. Aktulga, R. Berger, D. S. Bolintineanu, W. M. Brown, P. S. Crozier, P. J. in 't Veld, A. Kohlmeyer, S. G. Moore, T. D. Nguyen, R. Shan, M. J. Stevens, J. Tranchida, C. Trott, and S. J. Plimpton, Computer Physics Communications **271**, 108171 (2022).

[30] https://docs.lammps.org/pair_granular.html.

[31] H. P. Zhang and H. A. Makse, Phys. Rev. E **72**, 011301 (2005).

[32] P. G. Debenedetti and F. H. Stillinger, Nature **410**, 259 (2001).

[33] I. I. Blekhman, *Vibrational mechanics: nonlinear dynamic effects, general approach, applications* (World Scientific, 2000).
## Appendix A: Details on numerical simulations

Our DEM simulations make use of the pair_style granular command of the LAMMPS package [28–30] to implement the Hertz–Mindlin [30, 31] model for grain–grain, grain–lid, and grain–vane contact dynamics. The lid is composed of 1773 granular particles glued together, and the rotating vane consists of 4 × 10 particles glued together and overlapping by half a radius. Grain–wall interactions are described by the same model, implemented with the wall/gran/region command. The Hertz–Mindlin model accounts for both the elastic and dissipative responses during mutual compression of the interacting objects and incorporates the dynamics of their relative translational and rotational motions (see [30] for the specific expressions of tangential and normal forces).

All this requires five parameters to be set, representing the stiffness and viscosity of both the normal ($k_n$, $\eta_n$) and tangential ($k_t$, $\eta_t$) forces, along with the tangential friction coefficient $\mu$. For all our simulations, we fix $k_n = 6.1 \times 10^7$ Pa, $\eta_n = 2.9 \times 10^7$ (m·s)$^{-1}$, $k_t = 7.0 \times 10^7$ Pa, $\eta_t = 2.3 \times 10^5$ (m·s)$^{-1}$, and $\mu = 0.5$ for grain–grain, grain–lid, grain–vane, and grain–wall interactions, which, for simplicity, are considered to be made of the same material. These numerical values are chosen within a range that, in previous studies [19, 26], proved to yield numerical results consistent with experimental observations. As shown in Ref. [19], these parameters influence the energy transfer mechanisms, but the observation of the non-monotonic response as a function of driving frequency remains robust against their variation.

Our simulations are initialized by randomly pouring the grains into the container and then letting the lid fall on them. During the pouring process, we impose a high-amplitude shaking to compact the system more effectively. After pouring, the specific A and f to be probed are set, and the torque on the vane is switched on.

## Appendix B: The role of energy dissipation

In order to support our interpretation based on energy dissipation for the non-monotonic behaviour of K, we studied the behaviour of the average work done by the dissipative forces in our system. Our analysis starts by considering that in a time interval $\Delta t$, the variation of the total internal energy in the system (i.e. kinetic plus potential) is equal to the work $W_{nc}$ done by the non-conservative forces on the system in that interval $U(t + \Delta t) - U(t) = W_{nc}$. Assuming a stationary state and taking the time average of this equality we obtain $\langle W_{nc} \rangle = 0$. In other words, different contributions to the work done by non-conservative forces must balance each other. In the setting used to study the granular temperature of the granular medium (i.e. without the rotating vane), we have the work done by the vibrating

<figure><img src="image_5.png" /><figcaption>FIG. 4. Ratio of the dissipative power to the granular temperature (main panel) and dissipative power (inset) as a function of $\Gamma$ varied through f. Simulations are in the same geometry as the ones discussed in the main text but removing both the vane and the top lid. The driving amplitude is fixed to $A = 26$ µm.</figcaption></figure>

container $W_w$, the work done by the confining lid $W_l$ and finally the work due to the dissipative components in the grain interactions $W_d$ (see Sec. A and [30]). Thus, assigning the plus sign to the work done on the system and the minus sign to the work done by the system, we have $W_{nc} = W_w + W_l - W_d$. The average work done by internal dissipative forces can then be expressed as $\langle W_d \rangle = \langle W_w \rangle + \langle W_l \rangle$. In the following, to simplify the analysis, we considered simulations without the lid (i.e. $\langle W_l \rangle = 0$ by definition); moreover, stationarity allows us to write down the work contributions during an interval $\Delta t$ in terms of power. All this leads to the central identity for the analysis that follows:

$$ \omega_d \Delta t = \omega_w \Delta t, \qquad (B1) $$

where $\omega_d$ and $\omega_w$ are, respectively, the power associated to the dissipative interactions and the one associated to the vibrating walls. Practically, from our numerical data, we measured the instantaneous force $F_i^w$ between the walls and each grain i at contact for each simulation step $dt$, along with the displacement $dS_i$ travelled by the grain in the same $dt$. We point out that, to further simplify this numerical measure, we switched off the tangential interactions between the walls and the grains in the simulations performed for this analysis. Then, using Eq. (B1) we obtain the per-particle dissipation power

$$ \bar{\omega}_d = (Ndt)^{-1} \sum_i F_i^w \cdot dS_i. \qquad (B2) $$

This procedure is much easier than evaluating the dissipation rate calculating the work done by the dissipative contact forces for each couple of interacting grains. Nevertheless, we also performed some auxiliary simulations in a simplified geometry with only a couple of grains to check the consistency of our approach (not shown).

In Fig. 4, we plot the ratio $\bar{\omega}_d/K$ as a function of $\Gamma$ varied through f and observe an increasing monotonic trend. This quantity represents the rate at which energy is dissipated per unit of granular temperature at a given driving
frequency. Its monotonic increase tells us that maintaining a stationary state with a granular temperature K at higher driving frequencies implies a higher dissipation rate than at lower frequencies. In other words, our analysis provides a direct confirmation of our main argument for explaining the non-monotonic trend of the granular temperature as a function of f discussed in the main text. Increasing the driving frequency increases the amount of energy available for absorption by the granular medium (i.e. the walls have a higher kinetic energy). However, it also increases the dissipation rate, which makes the energy absorption process less efficient. In the inset of Fig. 4 we also plot $\bar{\omega}_d$ as a function of f showing that it is decreasing at large f. This tells us that the monotonic increase in $\bar{\omega}_d/K$ is due to the fact that the decay of dissipative power at large f is slower than that of the granular temperature.

## Appendix C: A minimal two-block model for the $\bar{A}^2/\bar{p}$ scaling

Here we present a generalized model of a driven-damped oscillator that qualitatively reproduces the main features of the frictional behavior observed in this study. The model is sketched in Fig. 7.

It extends the model introduced in [19] to account for the kinetic energy measured in DEM simulations, which shows a non-monotonic dependence on the driving frequency with a maximum at $f^*$, an overall increase with driving amplitude, and a shift of the peak frequency toward lower values as dissipation grows.

The model is derived from the equation of motion of a single particle, assuming that the collective dynamics of the granular system can be approximated by the behavior of an individual grain. Since the medium remains dense across all shaking regimes, particle motion at short timescales is confined within cages formed by surrounding grains. The key parameter is the mean collision time $\tau_c$, which decreases as the driving frequency increases, leading to more frequent dissipative events.

To keep the model simple, we restrict the analysis to vertical motion and describe a particle confined in a one-dimensional cage, represented by a spring attached to the base of a vibrating box oscillating with amplitude A and frequency f. The cage, composed of neighboring grains, also induces an effective viscosity. Numerical simulations show that $\tau_c$ becomes smaller for increasing the driving frequency f. This effect has been modeled by assuming a viscosity $\gamma$ which is an increasing function of f: $\gamma = \gamma_0 f^\alpha$ with $\alpha \simeq 0.67$ estimated from numerical simulations.

To study the effect of friction, we consider a similar model in which a second particle is placed on top of the system to mimic the role of the upper plate over the granular material. A repulsive spring of stiffness $k_2$ is then included in the model to describe the interaction between the grains and the top plate.

Accordingly, because of the interaction with the top

plate, the particle is subject to a spring-like force of the form

$$F_{\text{spring}} = -k_2 |z - z_t|^\beta \operatorname{sign}(z - z_t), \quad (\text{C1})$$

where z is the position of the particle and $z_t$ the position of the top plate. Here $\beta = 1$ or $\beta = 3/2$, depending on whether a Hookean or a Hertzian spring model is considered, respectively. A similar force acts on the particle because of the interaction via a spring of stiffness $k_1$ with the bottom plate, whose position $z_p$ presents regular oscillations $z_p(t) = A \cos(2\pi ft)$.

In addition, the particle experiences a viscous force

$$F_{\text{viscous}} = -\gamma (\dot{z} - \dot{z}_p), \quad (\text{C2})$$

where $\gamma$ is the viscous damping constant.

The top plate is subject to two forces: an upward force from the interaction with the particle,

$$f_b = k_2 |z - z_t|^\beta \operatorname{sign}(z - z_t), \quad (\text{C3})$$

and a downward force $F_P$ proportional to the external pressure. For simplicity, we consider a unitary mass $m = 1$ for both the particle and the top plate.

The key quantity is the difference between the top plate position $z_t$ and the particle position z. When $z_t - z > l_0$, with $l_0$ the rest length of the spring, there is no interaction between the top plate and the grains. This corresponds to a zero-friction scenario. Conversely, friction arises when the top plate is “attached” to the particle.

The mean value of friction coefficient can thus be estimated as the percentage of contact time:

$$\langle \mu \rangle = \frac{\tau_{\text{contact}}}{\tau_{\text{tot}}}, \quad (\text{C4})$$

where $\tau_{\text{tot}}$ is the total observation time and $\tau_{\text{contact}} \le \tau_{\text{tot}}$ is the fraction of time with $z_t - z \le l_0$. When the two particles are always in contact, $\langle \mu \rangle$ reaches its maximum value, while $\langle \mu \rangle = 0$ corresponds to the absence of contact.

To evaluate $\langle \mu \rangle$, it is necessary to solve the coupled equations for the dynamics of $z(t)$ and $z_t(t)$. This is a non-trivial task even in the Hookean case, and therefore we present below results of the numerical solution of the equations.

## 1. Results

In Fig. 5 we show the average friction $\langle \mu \rangle$ as a function of the dimensionless parameter $\Gamma = \frac{A f^2}{F_P}$ varies through f. Each curve corresponds to fixed values of A and $F_P$, while the frequency f is varied. Results are reported for the Hertzian case, but qualitatively similar behavior is found for the Hookean case.

We observe a clear non-monotonic trend: friction decreases with increasing $\Gamma$, reaches a minimum $\mu_{\min}$ at $\Gamma > 1$, and then increases again. Both the position of
<figure><img src="image_6.png" /><figcaption>FIG. 5. Average friction $\langle \mu \rangle$ as a function of $\Gamma = Af^2/F_P$ in the Hertzian case. Different curves correspond to different values of A and $F_P$. The equations of motion are solved with $k_1 = k_2 = 1$, and $\gamma = 5$.</figcaption></figure>

the minimum and the value of $\mu_{min}$ depend explicitly on A and $F_P$.

In Fig. 6 we analyze the dependence of the minimum friction $\mu_{min}$ on the driving amplitude A, for different confining pressures $F_P$, comparing Hertzian and Hookean models. As shown in Fig. 6a, at small amplitudes $\mu_{min} \simeq 1$, while for A larger than a characteristic value (which increases with $F_P$), $\mu_{min}$ decreases monotonically.

A striking difference emerges between the two contact models. For the Hookean case (Fig. 6b), rescaling A by $F_P$ leads to a perfect collapse of the data for different $F_P$. Conversely, in the Hertzian case (Fig. 6c), the collapse is recovered only when A is rescaled by $F_P^{1/2}$.

## 2. Analytical interpretation of the different collapse of $\mu_{min}$ between Hookean and Hertzian cases

Figure 6 highlights the different scaling behaviors of the minimum friction $\mu_{min}$ as a function of A and $F_P$ for Hookean and Hertzian interactions. In the Hookean case, the scaling is exact: rescaling A and $F_P$ by the same factor yields identical results, i.e., the data are strictly invariant under proportional changes of A and $F_P$. In contrast, for the Hertzian case the collapse obtained with $A/F_P^2$ is only approximate and limited to a certain range of parameters, with deviations emerging outside this regime. The origin of this difference can be clarified by analyzing a simplified version of the model, where the stiffness of the bottom spring is taken in the limit $k_1 \to \infty$. In this regime, the vertical position $z_1$ of the lower block is constrained to follow the bottom plate, i.e. $z_1 = z_p + \ell_0$, where $z_p$ is the vertical position of the bottom plate and $\ell_0$ is the rest length of the lower spring. Thus, the dynamics reduce to those of the top

<figure><img src="image_7.png" /><figcaption>FIG. 6. Analysis of $\mu_{min}$ as a function of deformation amplitude and pressure. (a) Trend of $\mu_{min}$ as a function of deformation amplitude (A) for two extreme pressure values: $F_P = 0.1$ (circles) and $F_P = 0.001$ (left triangles). Filled symbols represent the Hertzian contact model, while empty symbols represent the Hookean contact model. (b) Normalization of the Hookean data, showing the collapse of curves for various pressures as a function $A/F_P$. (c) Normalization of the Hertzian data, showing the collapse of curves as a function of $A^2/F_P$.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>FIG. 7. Schematic representation of the two-block model.</figcaption></figure>

block, whose vertical position $z_2$ satisfies the equation

$$m\ddot{z}_2 + \gamma\dot{z}_2 + k_2(|z_2 - z_1|)^\beta \operatorname{sign}(z_2 - z_1) = N - F_P, \quad (C5)$$

where $N = k_2\ell_0^\beta$ represents the reaction force ensuring that, in the absence of confining pressure ($F_P = 0$), the equilibrium position of the top plate is $z_2 = z_1 + \ell_0$.

Assuming that the bottom plate oscillates as $z_p = A\cos(\omega t)$, and introducing the relative distance between the top block and the bottom plate as $z = z_2 - z_p$, Eq. (C5) can be rewritten as

$$m\ddot{z} + \gamma\dot{z} + k_2(|z - \ell_0|)^\beta \operatorname{sign}(z - \ell_0) = N - F_P + mA\omega^2 \cos(\omega t) + \gamma A\omega \sin(\omega t). \quad (C6)$$
Introducing dimensionless variables

$$y = \frac{z}{\ell_0}, \quad A' = \frac{A}{\ell_0}, \quad F'_P = \frac{F_P}{N}, \quad \omega' = \omega t_0, \quad \gamma' = \frac{\gamma \ell_0}{N t_0},$$

with the characteristic time scale $t_0 = \sqrt{\frac{m\ell_0}{N}}$, Eq. (C6) can be expressed in non-dimensional form as

$$\ddot{y} + \gamma' \dot{y} + (1+y)^\beta \text{sign}(1+y) = 1 - F'_P + A'f(\omega'), \quad (C7)$$

where derivatives are taken with respect to the scaled time $t' = t/t_0$, and

$$f(\omega') = \omega'^2 \cos(\omega't') + \gamma' \omega' \sin(\omega't').$$

a. Hooke case $\beta = 1$

Setting, $\beta = 1$ in Eq. (C7), after rescaling the variable as $v = y/F'_P$, we obtain

$$\ddot{v} + \gamma' \dot{v} + v = 1 + \frac{A}{F'_P} f(\omega'), \quad (C8)$$

which clearly shows that the evolution of $v$ depends only on the ratio $A/F'_P$. This implies that, upon rescaling the driving amplitude $A$ by the external pressure $F_P$, the explicit dependence on $F_P$ vanishes. Such a scaling argument explains the exact data collapse observed in Fig. 6(b).

b. Hertz case $\beta = 3/2$

In this case, we focus on small oscillations around the steady-state solution, i.e. $y \ll 1$. Expanding Eq. (C7)

up to second order in $y$, we obtain

$$\ddot{y} + \gamma' \dot{y} + \beta y + \lambda y^2 = -F'_P + A'f(\omega'), \quad (C9)$$

with $\lambda = 3/8$.

We next decompose $y(t) = q(t) + \psi(t)$, where $q(t) = \langle y(t) \rangle$ denotes the average of $y(t)$ over one period $T = 2\pi/\omega$, and $\psi(t)$ represents the fast fluctuating component. Substituting into Eq. (C9) and averaging over one period yields

$$\ddot{q} + \gamma' \dot{q} + \beta q(1 + q/4) + \lambda \langle \psi^2 \rangle = -F'_P. \quad (C10)$$

The evolution equation for $\psi(t)$ can be obtained by subtracting Eq. (C10) from Eq. (C9). Following the so-called inertial approximation [33], we approximate it as

$$\ddot{\psi} = A'f(\omega'), \quad (C11)$$

which gives $\langle \psi^2 \rangle = \frac{\omega'^2}{2} (A'^2 \omega'^2 + \gamma'^2)$.

Substituting this result into Eq. (C10), we obtain

$$\ddot{q} + \gamma' \dot{q} + \beta q(1 + q/4) = F'_{P} - \lambda \frac{\omega'^{2}}{2} (A'^{2} \omega'^{2} + \gamma'^{2}). \quad (C12)$$

Since this equation holds for $y \ll 1$, we introduce the rescaled variable $v = q/F'_P$. At leading order, we obtain the oscillatory equation

$$\ddot{v} + \gamma' \dot{v} + \beta v = 1 - \lambda \frac{A'^2 \omega'^4}{2F'_P} - \lambda \frac{\gamma'^2 \omega'^2}{F'_P}. \quad (C13)$$

When $A\omega^2$ (the acceleration associated with the oscillating bottom plate) is much larger than $\gamma\omega$ (the dissipative contribution of the medium), the oscillations of the upper block are primarily governed by the ratio $A^2/F_P$. This condition is expected to hold near the frequency corresponding to the minimum in the friction coefficient. Therefore, Eq. (C13) provides a theoretical justification for the scaling observed in Fig. 6c. In summary, our analysis demonstrates that the $A^2/F_P$ scaling naturally arises from non-linear effects in the effective interaction between the granular layer and the upper plate.