Getting hotter by heating less: How driven granular materials dissipate energy in excess

We investigate how the kinetic energy acquired by a dense granular system driven by an external vibration depends on the input energy. Our focus is on the dependence of the granular behavior on two main parameters: frequency and vibration amplitude. We find that there exists an optimal forcing frequency, at which the system reaches the maximal kinetic energy: if the input energy is increased beyond such a threshold, the system dissipates more and more energy and recovers a colder and more viscous state. Quite surprisingly, the nonmonotonic behavior is found for vibration amplitudes which are sufficiently small to keep the system always in contact with the driving oscillating plate. Studying dissipative properties of the system, we unveil a striking difference between this nonmonotonic behavior and a standard resonance mechanism. This feature is also observed at the microscopic scale of the single-grain dynamics and can be interpreted as an instance of negative specific heat. An analytically solvable model based on a generalized forced-damped oscillator well reproduces the observed phenomenology, illustrating the role of the competing effects of forcing and dissipation.


I. INTRODUCTION
The fascination of granular systems relies on their rich phenomenology, which eludes standard statistical physics [1]. Examples of surprising behaviors were reported in experiments and numerical simulations: Brazilnut effect [2], ratchet effect [3,4], spontaneous segregation of mixtures [5], nonequilibrium phase transitions [6], anomalous thermal convection [7,8], Kovacs-like memory effect [9], see also [10]. These features are ascribed to the dissipative nature of these systems, which generally cannot be treated via the introduction of effective parameters [11]. A fundamental open question is how the dissipation mechanisms relevant for different behaviors are related to the external energy injection. These mechanisms involve several scales, from particle-particle collisions to the interaction with boundaries. The relation between system kinetic energy K and input energy S involves the nonequilibrium response beyond the linear regime. In general, nonmonotonic behaviors have been observed in different models, from driven tracers in crowding environments (negative differential mobility) [12][13][14][15][16][17], to systems showing negative specific heat due to long-range interactions [18][19][20][21] or to the presence of baths at different temperatures [22,23]. Similar effects have been also observed in force-free cooling granular gases of aggregating particles [24], where the granular temperature can grow while the system energy decreases due to dissipative collisions.
These behaviors affect fluidization properties of the medium, with important consequences in industrial applications [25], where usually energy is fed via mechanical vibrations with frequency f and amplitude A. In some cases, a relation K ∼ S α , with S ∼ (Af ) 2 , has been derived [26][27][28][29]. Experimental studies focused on the spe-cific role of the forcing mechanisms, such as vibration amplitude, frequency or velocity [30][31][32][33][34][35][36][37][38][39][40][41][42][43][44][45][46][47]. In particular, an optimal frequency for energy transfer is found when the system is in a bouncing-bed state and A and f are varied keeping S fixed [30,31]. In this state the granular medium detaches from the driving plate and a resonant behavior is achieved via a synchronization between plate vibration and bed bouncing.
In this paper we consider vibration amplitudes small enough to keep the granular system always in contact with the oscillating plate. Our main result is that, even in this regime, if the input energy fed into the system is increased by increasing f (keeping A fixed), a nonmonotonic behavior is observed in both the kinetic energy of a driven vane immersed in the medium and in the kinetic energy of the granular medium itself. We show that there is an optimal frequency where the system reaches a maximum kinetic energy. Differently from previous results [30,31], our system is not in a bouncingbed or resonant state, and the granular internal energy is nonmonotonic with f , even if the energy input always increases. Our results can be interpreted as an example of negative specific heat, shedding light on the findings reported in [24]: Here, we investigate the complementary process, where the granular temperature decreases when the input energy increases.
The paper is structured as follows. In Sec. II we describe the experimental setup and the numerical model. In Sec. III we discuss the nonmonotonic behavior of the system energy as a function of the input energy, which represents the main result of the paper. Sec. IV is devoted to the study of the vane angular velocity, which shows a similar nonmonotonic behavior, related to the viscous properties of the granular medium. In Sec. V we present the analysis of the single-particle dynamics, which provides further details on the physical mechanisms underlying the observed phenomenology. The proposed theoretical model of a generalized driven-damped oscillator is described in Sec. VI. In Sec. VII we compare our findings with previous results, pointing out the main differences. Finally, in Sec. VIII we draw some conclusions on the effects and mechanisms observed in our system. Appendix A provides details on the numerical model.

II. EXPERIMENTAL SETUP AND NUMERICAL MODEL
Our experimental setup is illustrated in Fig. 1 (see also Ref. [33]). The granular medium is made of N = 2600 steel spheres (diameter d = 4 mm, mass m = 0.27 g), contained in a conical-shaped floor cylinder (diameter 90 mm, minimum height 28.5 mm, maximum height 47.5 mm), which enhances the energy transfer from vertical to horizontal directions. A rigid lid (mass M top = 218 g) covers the system, to confine particles and to allow for a spectral analysis of the system oscillations. A vertical vibration is imposed to the coordinate of the container bottom z p (t) (1) A Plexiglas vane (height 15 mm, width 6 mm, length 35 mm) is suspended in the system. The vane can only rotate around the vertical axis, subjected to a constant torque T = 6 × 10 −3 Nm. The system is simulated by LAMMPS package [48]. The interaction among grains is described via the nonlinear Hertz-Mindlin model [50,51], see Appendix A for details. The relevant parameters are the stiffness of the nonlinear contact k n and the viscous damping coefficient γ n , accounting for dissipative interactions. The geometry of the system and numerical parameters are chosen to reproduce our experimental setup: The lid is made of 1773 granular particles glued together and the rotating vane is made of 4 × 10 particles glued together and overlapping by half a radius.

III. OPTIMAL FORCING FREQUENCY AND ROLE OF DISSIPATION
In the main panel of Fig. 1 we report the vane kinetic energy contribution due to fluctuations, [49], where I = 353 g mm 2 is the momentum of inertia and Ω the angular velocity, measured in experiments as a function of f for different values of A. We also show a case with no lid to demonstrate the robustness of the observed behavior. In the inset we report results of numerical simulations showing that the model well reproduces the nonmonotonic behavior of the real system. In the following, we will be mainly interested in the total kinetic energy of the granular medium and in the dynamics of single grains, and therefore we will remove the vane in some simulations of the numerical model [54]. Numerical simulations allow us to investigate granular kinetic energy and single-grain motion, exploring a wider range of vibration frequencies, f ∈ (0, 1000] Hz, with respect to the experiments. As shown in Fig. 1, the vane kinetic energy is a nonmonotonic function of f (at fixed A): K v grows abruptly from zero to a finite value at a frequency threshold f 1 , related to the detachment condition [39], 2πf 1 = g/A. Then, after a maximum, K v starts to decrease, signaling that the granular medium leaves the state of maximal fluidization, hindering the motion of the vane by an increased effective viscosity. One can define a so-called friction recovery frequency, f 2 > f 1 , at which the system kinetic energy decays to about 1/2 of its maximum value. This behavior is related to dissipation mechanisms of the granular medium that depend on the viscoelastic properties of the material.
To better characterize the intrinsic behavior of the granular medium and its response to the external forcing, we focus on the total translational kinetic energy (the rotational kinetic energy being negligible, see Fig.  3a where v i is the grain velocity, in the absence of the suspended vane. Therefore, we introduce the adimensional quantity describing the energy input, S = (2πf ) 2 A 2 /(gd).
In Fig. 2 we report K(S), varying f at fixed A, for several values of γ n . The nonmonotonic behavior shows that the granular system reaches the maximal kinetic temperature at an optimal value of the input energy. When the input energy exceeds a certain threshold, related to f 2 , the system cools down because the dissipation effects increase. A key result is represented by the maximum position of the kinetic energy for different values of viscous damping coefficient γ n . We find that the peak shifts to the left when the dissipation in the system is increased. It is interesting to note that this shift doesn't occur if we vary k n signaling that the friction-recovery mechanism is not governed by the elastic contribution of the interaction (see Fig. 3(b)). Therefore the system can transfer more energy to the grain motion when it is vibrated at an optimal frequency. If the frequency increases, the overall external energy injected is larger but the dissipation mechanisms become dominant and the granular kinetic temperature decreases. This behavior can be interpreted as an instance of negative specific heat [22,24], occurring due to the subtle interplay between forcing and dissipation.
Conversely, if we increase S by increasing A at fixed f , we find a monotonic behavior. Namely, dissipative mechanisms are not activated and the kinetic energy keeps growing with the input energy as shown in Fig. 4. These features are a consequence of the permanent contact with the driving plate, which makes dissipation dominating at high frequencies. This is a novel phenomenon different from the resonant behavior in the bouncing-bed state clearly described in Ref. [30,31], where the optimal frequency is an increasing function of dissipation. Striking differences are also provided by the spectral analysis of the top plate oscillations which indicates that, in our system, the energy transfer is not maximized in the most coherent states (see Sec. VII C).

IV. VANE AVERAGE VELOCITY
Another crucial quantity in our experimental/numerical setup is the average angular velocity Ω of the vane in the steady state. While the variance is related to the kinetic energy of the system, the mean FIG. 3. Mean translational kinetic energy of the whole system K versus the driving frequency f at fixed A = 0.026 mm (simulations without the vane). (a) Comparison between the translational kinetic energy, the rotational one and the total (translational plus rotational). We show two curves reported in Fig. 2 as a function of f , instead of S. The rotational degrees of freedom follow the same nonmonotonic behavior but with smaller absolute values, so that their contribution does not affect the qualitative behavior of the total kinetic energy. (b) Study of the mean kinetic energy for different values of nonlinear stiffness: kn is varied at fixed γn = 2.9 × 10 7 (ms) −1 . The three curves have the same shape (the peak position don't change) but are vertically ordered by kn. This implies a larger kinetic energy and therefore less dissipation in the system with higher stiffness.
value of Ω is proportional to the inverse of the viscosity acting on the tracer during its motion. In the following we present the numerical study of Ω and K with particular attention on the role played by the parameters k n and γ n of the HM model in the characterization of its behavior. For an experimental study of this observable in the same setup see [33].
In Fig. 4b we show Ω as a function of S by increasing f at fixed A and vice-versa. We find a behavior very similar to the one of K reported in the previous section, the only qualitative difference is that Ω versus A reaches a constant value, while K versus A keeps growing also at larger amplitudes (compare with panel (a) of the same figure). To better understand the role of the interaction parameters we concentrate on Ω (f ) at con- stant amplitude. In Figs. 5(a) and (b) we vary γ n at fixed k n , while in Fig. 5(c) we vary k n at fixed γ n . We see that the parameter k n changes the typical shape of the curves near f 1 : we have a sudden peak when k n is high and a smoother behavior when it is small. In Fig.  5(b) we observe that at lower k n , reducing γ n clearly slows down the friction recovery. On the other hand, at higher k n ( Fig. 5(a)) we see the same effect but without a significant change in the height of the peak. This can be interpreted intuitively observing that a higher stiffness reduces the effect of the viscous damping during the collision. In Fig. 5(c) we see that also an increase of k n implies a growth of the recovery frequency, but this effect is less pronounced with respect to the one obtained decreasing γ n . We can conclude that f 2 is mostly ruled by the dissipation through γ n , in agreement with previous experimental results [33].
To sum up, a comparison between Figs. 5(a), 5(b) and 3(a) shows that the behavior of both Ω (a rheological response) and K (a thermodynamic observable) as a function of f is similarly affected by the change of γ n . On the contrary, from Figs. 5(c) and 3(b), it is clear that changing k n affects Ω more than K.

V. ANALYSIS OF THE SINGLE-PARTICLE DYNAMICS
The macroscopic features above described can be related to microscopic properties by investigating the single-grain dynamics. Typically, in dense systems, the diffusive motion of a single particle is hindered by the presence of many surrounding particles, inducing a cage effect [55][56][57]. Focusing on the Mean Squared Displacement (MSD) of a particle, one observes that, after a ballistic motion at short times, a plateau develops, signaling that the particle is trapped inside a cage. Then, at longer times, the particle manages to explore a larger region of the system and its MSD grows linearly in time. The relevant quantities are the cage size S c and the trapping time inside a cage τ c , both estimated from the first point of the MSD after the ballistic regime. This picture is fully supported by Fig. 6, where we show the MSD (averaged over about 20 particles randomly chosen in the system) in the horizontal XY plane, the plane relevant for the vane's motion. Rescaling time by the trapping time τ c and the MSD by the squared cage size, we find that curves collapse.
The nonmonotonic behaviors described above for K v and K, can be rationalized by the study of S c and τ c , as reported in the insets of Fig. 7. We find that τ c is a decreasing function of the frequency, τ c ∼ f −0.67 , whereas S c is a nonmonotonic function of S, varying f at fixed A, with a maximum at S 10 −2 . At low frequencies the cage size tends to be very small but the particles need a long time to be trapped. Conversely, at high frequencies cages are still small but the particles are trapped in a short time. Moreover at low frequencies, since the system is weakly perturbed by the input energy, a particle is able to explore larger regions for increasing S. However, for larger S (increased through f ), collisions become more frequent and the explored cage region decreases accordingly. If we estimate the average particle speed v XY on the XY plane as v XY ∼ S c /τ c , we find a nonmonotonic behavior, in agreement with the behavior of K (Fig. 7). This analysis shows that the single-particle dynamics reflects the same phenomenology occurring at the macroscopic scale.
Moreover, the nonmonotonic behavior of S c could be related to a change in the effective number of degrees of freedom in the system as the input energy is varied, in agreement with the explanation of negative specific heat suggested in Ref. [24] for cooling granular gases of aggregating particles. In that case, although the total energy of the system decreases due to dissipative collisions, the total number of clusters can decrease faster due to agglomeration, producing an increase of energy FIG. 9. Coarse-grained description of the single-particle dynamics. The springs k and the viscosity γ represent the cage that mediates the vibration of the shaker.
per particle.
We show in Fig. 8 the rescaled MSD on the xy-plane averaged over 20 particles and the relative cage size S c for different values of A. Here we see that the collapse is very good up to time t/τ c ∼ 1 (so the trapping time τ c does not vary with A) while S c (and consequently the speed v xy = S c /τ c ) grows monotonically with the driving amplitude. These last results show that the same phenomenology as a function of A at fixed frequency is present both at the macroscopic scale ( Ω and K) and at the single-particle one (v xy ).

VI. THE GENERALIZED DRIVEN-DAMPED OSCILLATOR: A SINGLE-PARTICLE PERSPECTIVE
Here we present a generalized model of a drivendamped oscillator that reproduces qualitatively the phenomenology of K studied in the simulations. In particular, our model shows: i) a nonmonotonic behavior of the energy as a function of the driving frequency f with a maximum in f * ; ii) an increasing behavior of the energy as function of the driving amplitude A; iii) the shift to the left of the frequency f * when the dissipation in the system is increased.
The model is obtained starting from an equation of motion for a generic particle in the system and then assuming that the energy of the whole system follows the same behavior of that of the single particle. This assumption is consistent with the previous analysis that shows a common phenomenology of single-particle quantities and macroscopic ones. In addition, the same assumption is done in PEPT experiments where the time-averaged behavior of a single tracked particle is considered representative of that of the global system [30].
A. Newton equation for a caged single particle The granular system is very dense in all regimes of shaking, so the short-time dynamics of the grains is expected to take place inside their cages formed by the surrounding particles. This short-time dynamics presents multiple relevant time-scales such as the one represented by the integration time dt (τ 1 ∼ 10 −5 s), that is a fraction of the duration of a single collision, and the one fixed by the inverse of the driving frequency τ d = 1/f ∼ 10 −3 − 10 −1 s. A third one is represented by the mean time between two collisions, reported in Fig. 7 (but also observed for the XY Z diffusion), which spans values in the range τ c ∼ 10 −4 − 10 −3 s. These time-scales are thus ordered as follows: τ 1 < τ c < τ d for all the cases presented in our study.
In order to estimate the mean kinetic energy over many collisions, we are interested in a time coarse-graining on a scale larger than τ c . At the same time, we want to write a differential equation for the coordinates of a particle in which the shaker dynamics appears as an external driving. In view of these arguments, we concentrate on time scales larger than τ c but smaller than τ d . To make our model as simple as possible we only consider the motion in the z direction. Therefore we consider a single particle of mass m confined in a one dimensional cage that consists in one spring of stiffness k with a resting length l 0 connected with the bottom of a vibrating box. The latter really represents the experimental/numerical container that oscillates following z p (t) = A cos(2πf t). The fact that the cage is actually made of fast vibrating particles also induces an effective viscosity with coefficient γ.
Our simplified coarse-grained description is sketched in Fig. 9, where we refer to the coordinate of the grain as z g . We remind here that the coefficients of the HM model used in the simulation (k n and γ n ) depend on the material properties and act on the fast time scale τ 1 . The way in which they are connected with the effective viscosity γ and stiffness k is not trivial. A reasonable value for k can be estimated considering that the modeled spring is actually made by a column of grains. Every grain with radius r and Young modulus Y can be thought of as a microscopic vertical spring with an elastic constantk given bỹ k = πY r/2 (for simplicity we are considering the grains as cylinders instead of spheres). Now the effective stiffness of the column results from the parallel of a mean numbern of such microscopic springs: k =k/n. For the parameters of our simulations and fixingn = 4.5 we find that k ∼ O(10 5 ) N/m. Regarding the effective viscosity, we will do an ansatz in the following. Looking at Fig. 9, it is straightforward to write for z g : where ξ(t) = z g − z p − l 0 , 2πf s = γ/(2m) and 2πf k = k/m ∼ O(10 4 ) s −1 . Now we come to the crucial hypothesis of our approach, i.e. the ansatz on f s . Some previous studies [30,33] suggest that the dissipation of energy due to interparticle collisions increases significantly when the driving frequency f grows. This is also visible in our data from Fig.  7, where we see that τ c becomes smaller for increasing f . Indeed, a reduction of the time between collisions means a growth of the number of dissipative events in the system (i.e. the collisions themselves). The simplest way to insert this dependence of the internal dissipation on the external driving is to take f s as an increasing function of f : f s = af α with a, α > 0. Bringing the variable z p (t) contained in ξ(t) to the right-hand side of Eq. (2) and adding gravity, we obtain the following equation: This is the equation for a driven-damped harmonic oscillator with characteristic frequency f k , viscous constant 4πaf α and external driving (2πf k ) 2 A cos(2πf t) that oscillates around the equilibrium position z eq g = l 0 − g/(2πf k ) 2 . The stationary solution of Eq. (3) is: Deriving z g (t), taking the square and then integrating over a period 1/f we can find the mean quadratic velocity of the particle: As showed in Fig. 10 for the specific value α = 2/3 ∼ 0.67 taken from the behavior of τ c (but it holds in general for α > 0), ż g 2 has a nonmonotonic behavior for f < f k and its maximum value shifts to the left as a is increased. In Fig. 10 we show the behavior of the mean total kinetic energy K = m eff ż g 2 /2 as a function of S, for several values of a. From a fitting procedure, we found 2πf k = 11905 s −1 , which is of the order k/m if the effective stiffness k is estimated by considering seriesparallel microscopic elastic constants relative to the grain material, as above illustrated. The prefactor m eff represents an effective mass and we find an optimal agreement with data for m eff = 912.6 g, which is of the order of the total mass of the system. Therefore, the analytical model predicts the general phenomenology of the 3D granular system, with a good quantitative agreement for high values of γ n .

B. Beyond a simple resonance
Eq. (5) has a very simple form and its nonmonotonic behavior can be explained in terms of competitive mechanisms. From Eq. (5) we see that the mean quadratic velocity of the particle is proportional to the input energy of the shaker divided by an adimensional factor: can be considered the two competitive terms if f < f k . In this regime, u 2 is a decreasing function of f and enhances energy transfer, while u 1 (that contains the dissipation) increases with f and therefore has an opposite effect. In order to better understand the underlying mechanisms, we remind that for an ordinary driven-damped oscillator (namely with a viscous coefficient that does not depend on f ) the mean quadratic velocity is the same as Eq. (5), with the only difference that the dissipative term is proportional to the square of the rescaled driving frequency: u 1 = u 1 (f, 0). From this point of view, 1/u 2 can be interpreted as an energy gain factor: It grows before the resonant frequency (f = f k ) and then decreases. In the case of the usual damped oscillator, therefore, the nonmonotonic behavior is entirely explained by the nonmonotonic behavior of u 2 alone, and in fact dissipation does not change the peak position but only smoothes it. On the contrary, when α > 0 the nonmonotonic behavior has a different origin, coming from the competition between dissipation (u 1 ) and gain (u 2 ). This can be rigorously checked, deriving Eq. (5) to find the condition for the maximum: that turns out to be satisfied by f = f k only for α = 0. The competition is apparently present also in the ordinary driven-damped oscillator but in that case it is balanced by the (2πf ) 2 contained in V (f ) at the numerator in such a way that the nonmonotonic behavior of the energy can be explained only by the resonance. We finally note that another way to see competitive terms in Eq. (5) is to rewrite it as ż g 2 = (u 1 /V + u 2 /V ) −1 . In this form we have the inverse of a sum of two terms that, for f < f k , depend in opposite way and with different powers on f . This clearly gives raise to an extremal point also in the limit f f k , that is consistent with the values of the fitted parameters.
We stress that these mechanisms substantially differ from the standard resonance phenomenon, due to the presence of the damping term f s ∼ af α . Indeed, an increase of f induces the grains to adsorb the injected energy in a faster vibrating motion, experiencing a larger number of collisions (dissipative events) per unit of time. This means that at higher f , both the energy input and the energy output increase. The two phenomena compete and, because of their different functional dependencies upon f , an extremal point appears. We note that the maxima in the curves of Fig. 2 and 10 occur at frequencies much smaller than the fitted f k . Conversely, the ordinary resonance always occurs at the fixed characteristic frequency f k , independently of the viscous coefficient. For α > 0, the dissipation controls the position of the peak, reproducing the key feature of molecular dynamic simulations.

VII. COMPARISON WITH PREVIOUS STUDIES
It is interesting to discuss some previous results [30,31] about energy transfer optimization in vibrated granular media to compare with those presented in this paper. In [30] the authors present experimental and numerical data for a vertically-shaken granular system in a 3D geometry. They find a behavior of the kinetic energy K versus the driving frequency f similar to the one found by us. Their study is done for increasing f and decreasing A keeping constant the shaker strength S = (2πf A) 2 /(gd). The main result is that there is a specific combination of f and A that, at a given fixed S, optimizes the energy transfer between the shaker and the granular bed. This nonmonotonic trend is rationalized by noting the following competitive effects: on the one hand, raising f makes the maximum rescaled acceleration Γ = (2πf ) 2 A/g larger enhancing the fluidization of the granular medium; on the other hand, decreasing A lowers the mean number of collisions with the shaker weakening the interaction with the external source of energy.
Regarding [31], the authors study the same setup but particularly concentrating on a 1D geometry, i.e. vertically-shaken columns of single grains. In this situation, when K is plotted against f at fixed S they find many peaks at the integer multiples of a specific frequency f * . They interpret this phenomenon as a resonant behavior: when the driving frequency is synchronized with the typical time of detachment (τ * = 1/f * ) of the granular column the energy transfer from the shaker to the system is optimized.
In the following, we comment on some crucial evidences about the differences between the results reported in our FIG. 11. Different regions of the parameter space investigated in our work with respect to Refs. [30,31]. study and those in [30,31].

A. A different granular phase
Previous studies consider a range of parameters which is completely different from the one considered by us. More precisely, in our experiments/simulations we consider f ∈ [20, 1300] Hz and A = 0.014, 0.026, 0.039, 0.053 mm, so that S ∈ [0.0003, 1.15]; on the contrary Refs. [30,31] investigate the ranges f ∈ [5,80] Hz and A ∈ [0.6, 8.0] mm, corresponding to S ∈ [0.9, 11.5], see Figure 11. They consider amplitudes A much larger than those studied in our system. For those values of A the system detaches from the driving plate reaching the usually called bouncing bed state, where it is possible to define the time of free flight τ * responsible for the aforementioned resonant behavior. We also stress that in [30] the granular medium is in a dilute phase where every grain explores all the accessible space whereas in [31] the grains are arranged in a dense state and each particle in the column remains in the same position with respect to the others.
In our article we show that a vibrofluidized state exists at much smaller values of A (see Fig. 11) which are sufficiently small to keep the granular medium always in contact with the bottom of the container. This is supported by the collision rate as a function of f , plotted in Fig. 12. In our system, indeed, the collision rate is an almost linearly increasing function of f , in sharp contrast with the explanation of the phenomenon reported in Refs. [30,31]. More precisely we have checked (visually in the experiments and quantitatively in the simulations) that, for all the driving parameters we explored, the grains are always arranged in a dense packing where they vibrate around an almost fixed position experiencing rare and slow rearrangements. As a consequence, in the considered range of parameters, the notion of time of free flight is meaningless and the mechanism responsible for the nonmonotonic behavior of K has nothing to deal with a resonant mechanism. Indeed, at variance with the bouncing-bed state, the permanent contact allows the oscillating driving plate to continuously transfer kinetic energy to the granular medium as the driving frequency is increased. In our case, the nonmonotonicity of K originates from a competition between driving energy and dissipation, which both increase with frequency, but with different powers of f . In this sense, taking the point of view of the authors of [30,31], we observe a maximum of the internal granular energy even if the input energy is increasing.

B. A different role of dissipation
In the bouncing bed state the optimal frequency f * is obtained when the vibration period τ * is synchronized with the flight time of the granular bed leading to Eq.(5) of Ref. [31] where g is the gravitational acceleration, the restitution coefficient (i.e. the fraction of energy lost during a collision), N the number of grains and β a fitting parameter that depends on the details of the experimental setup. Eq. (7) clearly shows that the resonant frequency f * = 1/τ * increases when the restitution coefficient is decreased: the authors rationalize this observation with the fact that more dissipated energy in collisions implies shorter flights, i.e. smaller τ * or higher resonant frequencies f * . In our case, conversely, our peak frequency decreases when the friction coefficient γ n of the normal force between the grains (and between grains and borders) is increased. This is the evidence of a very different mechanism at work: when we increase γ n , i.e. dissipation per single collision, the range of frequencies where dissipation wins the competition is larger, and this corresponds to a shift of the peak to the left. In Fig. 13, we show the spectra of the z-coordinate (measured with a laser sensor, see [33] for details) of the top plate for different values of the driving frequency f and A = 0.053 mm. This figure sheds light on the relation between chaotic motion and energy transfer in our system. We see that the most coherent motion (pronounced peaks on the driving frequency and its harmonics) is obtained for f = 53 Hz (panel (a)), before the energy maximum. At the driving frequency where our maximum in the energy occurs (panel (b), f = 100 Hz), we see that the peaks on the harmonics are less pronounced and more broadband. For increasing driving f (decreasing transferred energies), we see a more and more chaotic spectrum with the appearance of distinct peaks at noninteger multiple of f (panel (c)) and the disappearance of the peaks in the harmonics (panel (d)). Our scenario, therefore, is that of an increasing chaoticity -with f -of the granular dynamics, irrespective of the energy transfer, i.e. both before and after the energy maximum. This is quite different from the one observed in [30,31] where the energy transfer to the system is optimized at the frequency such that the granular system is more regular and shows the most defined peaks in its spectrum.  TABLE I. Numerical values for the material properties and the coefficients of the viscoelastic interaction. We use the index a for the grain-grain, grain-vane and grain-lid interactions while index b is used for the grain-wall ones.

VIII. CONCLUSION
Counter-intuitive phenomena, forbidden by standard thermodynamic arguments, can occur in nonequilibrium systems. Negative differential mobility and negative specific heat are typical examples [16,17,22,24]. This kind of nonmonotonic behavior in the current-force relation is due to the combination of competing mechanisms. Here we have unveiled a similar effect in a complex many-body system with dissipative interactions: The interplay between external forcing and internal dissipation in granular media can result in a nonmonotonic behavior of the system kinetic energy as a function of the input energy, representing an instance of negative specific heat. Moreover, our analysis explains the important role played by the vibration frequency, triggering specific dissipation mechanisms. These arise at different scales, from rheological behavior to single-particle motion. The observed phenomenology may have a deep impact application in several fields related to granular matter physics. Our molecular dynamics simulations are performed through LAMMPS package [48] and the system geom-etry reproduces the experimental setup of Ref. [33]. The values of Ω and K reported in the main text are the result of a time average over 3 × 10 7 time steps (78 equivalent seconds) in the steady state. Regarding the interaction between the grains, we used the Hertz-Mindlin (HM) model [50][51][52] to solve the dynamics during the contact. This viscoelastic model takes into account both the elastic and the dissipative response to the mutual compression. In addition, it includes in the dynamics not only the relative translational motion but also the rotational one. The model is described by the following equations: ξ ij (t ) ds(t ).
These equations are written for two particles with radius R i , R j , center positions r i , r j translational velocities v i , v j and rotational velocities ω i , ω j . The relative velocity is defined as g ij = (˙ r i − ω i × R i n) − (˙ r j + ω j × R j n) where n = ( r i − r j ) / | r i − r j |; we call g N ij and g T ij the two projections, respectively normal and tangential, to the surface of contact. The instantaneous normal compression is represented by ξ ij (t) = R i + R j − | r i − r j | and its derivative isξ ij (t) = | g N ij |. During the contact, the particles are subjected to a normal force F N ij and a tangential one F T ij ; both these components have an elastic and a dissipative contribution respectively characterized by the coefficients k n , k t , γ n , γ t . In the normal force F N ij we can see an elastic term that derives from the hertzian theory of contact mechanics [53] characterized by a nonlinear dependence on the displacement. The HM model is used to describe the interactions between all elements of the simulation (the box, the vane and the lid) considering the flat surfaces of the box as spheres with infinite mass and radius. Further discussions on the physical meaning and estimation of the model parameters and on the choice of the simulation time step dt can be found in the Supplemental Material of [56]. Here we report the numerical values used for the present study (Tab. I).