Low-frequency oscillations in narrow vibrated granular systems

We present simulations and a theoretical treatment of vertically vibrated granular media. The systems considered are confined in narrow quasi-two-dimensional and quasi-one-dimensional (column) geometries, where the vertical extension of the container is much larger than both horizontal lengths. The additional geometric constraint present in the column setup frustrates the convection state that is normally observed in wider geometries. This makes it possible to study collective oscillations of the grains with a characteristic frequency that is much lower than the frequency of energy injection. The frequency and amplitude of these oscillations are studied as a function of the energy input parameters and the size of the container. We observe that, in the quasi-two-dimensional setup, low-frequency oscillations are present even in the convective regime. This suggests that they may play a significant role in the transition from a density inverted state to convection. Two models are also presented; the first one, based on Cauchy's equations, is able to predict with high accuracy the frequency of the particles' collective motion. This first principles model requires a single input parameter, i.e, the centre of mass of the system. The model shows that a sufficient condition for the existence of the low-frequency mode is an inverted density profile with distinct low and high density regions, a condition that may apply to other systems too. The second, simpler model just assumes an harmonic oscillator like behaviour and, using thermodynamic arguments, is also able to reproduce the observed frequencies with high accuracy.


Introduction
Vibrated beds of granular materials present a wide range of different behaviours: phase separation [1,2], Faraday-like pattern formation instabilities [3,4], heap formation and convection [5,6], segregation [7,8], clustering [9] and periodic cluster expansions [10], as well as many others. These systems generally present a remarkable collection of distinct nonequilibrium inhomogeneous stable states for relatively small changes in the energy injection parameters. Hence, they are specially suited for the study of nonequilibrium phase transitions, as well as non-linear phenomena in general. Careful analysis of the microscopic mechanics behind the different transitions improves the comprehension of the complex dynamics present in driven granular systems. This gives further insight into when, and until what point, granular media behave like classical gases, fluids or solids, or whether they require an altogether different theoretical approach.
As can be seen by the aforementioned studies, the geometry of the system plays a fundamental role in determining the phenomena. Just by reducing the effective dimensionality of the system it becomes possible to observe behaviour not easily identifiable in fully three dimensional systems. The natural approach of study is then confining the grains to quasi-two-dimensional systems, where also particle-tracking methods become possible. Our study is inspired by one specific quasi-two-dimensional geometry that presents several distinct states in the energy injection parameter space: a vertical narrow box. That is, we focus on a vertically vibrated Hele-Shaw cell with the walls parallel to gravity, inside which the grains are located. The first reported classification of the different states present in this geometry was realised by Thomas et al. [11], in what would now be considered the low energy injection limit. Research then focused on the wave-like dynamics of the granular bed, and its variations with the frequency and the amplitude of oscillation [12,13]. It was with simulations that the energy input was considerably increased, and the existence of a density inverted state was first reported [14]. This state, named Leidenfrost after the analogous waterover-vapour phenomena [15], was then experimentally studied in depth by Eshuis et al. [16,17], as well as the buoyancy driven convection regime that is observed for higher energy inputs.
In the following simulational study the dimensionality of the vertical, narrow box is progressively reduced until both the width and depth are only five particle diameters wide, rendering the system effectively one-dimensional (see Figure 1). We refer to this last setup as the granular column. The first direct consequence of this further confinement is the frustration of the horizontally inhomogeneous states present in the wider geometries. Particularly, the suppression of convection makes it possible to directly observe the grains collectively oscillating at a much lower frequency than the energy injection frequency. Appropriately, we call these oscillations low-frequency oscillations (LFOs). Effective frequencies and amplitudes are defined and studied in the container length and energy injection parameter space. We then argue that LFOs are an essential feature of the dynamics of the narrow vibrated geometry, but it is only in Figure 1. Snapshots of three vertical narrow box systems with the same number of filling layers F = Nd 2 /l xly = 12, with N the total number of particles,d the (dimensional) diameter of the spherical particles, andl x andl y = 5d the (dimensional) width and depth of the container, respectively; and energy injection parameters, but different widthsl x . From left to right,l x = 100d,l x = 20d, andl x = 5d. The rightmost corresponds to the column geometry. Particles are coloured according to their kinetic energy.
the quasi-one-dimensional column setup that they can be easily isolated from the other collective grain movement of convection. Simulational measurements confirm this, as well as a continuum description of the system, which captures the correct frequency response for high energy inputs. The frequency behaviour is actually analogous to a forced harmonic oscillator, and is obtained mainly by considering a vibrated media with a high density region suspended over a low density one. This density inverted situation is indeed present, to different extents, in both the Leidenfrost and the convective regimes.

Simulations
Simulations are performed using an event-driven molecular dynamics algorithm [18]. In this approach, particles move freely under the effect of gravity until an event take place, namely, a collision with another particle or a wall. The motion of the particles in between successive events does not have to be simulated: if their trajectory equation is known, time can be advanced in variable steps. This makes event-driven simulations considerably faster than usual soft-particle simulations, where time is advanced at constant steps, independent of particle interactions. However, the need of having an analytical expression for the particle motion is a strong condition that limits the possible interaction between particles. In the following, we consider the most common approach: perfect hard-spheres, which imply binary collisions and no overlap or longrange forces between particles. Collisions are then modelled by normal and tangential restitution coefficients, ǫ n = ǫ t = 0.95, and also static and dynamic friction coefficients, µ s = µ d = 0.1 [19]. In order to avoid inelastic collapse we use the TC model [20], with a constant value t c = 10 −6 (d/g) 1/2 , withd the (dimensional) diameter of the spheres, andg the (dimensional) gravity. (In the following, quantities without a tilde are dimensionless). That is, collisions between two entities are considered elastic if they occur more frequent than 10 −6 gravity timescale units. These values are known to reproduce complex behaviour observed in similar vibrated setups using stainless-steel spheres ofd = 1 mm tod = 5 mm in diameter [10,21].
The setup consists of an infinitely tall container of widthl x and depthl y inside which the grains can move. The boundaries of the container are considered solid, and have the same collision parameters as particle-particle collisions. The whole box (both the bottom and the side walls) is vertically vibrated in a bi-parabolic, quasi-sinusoidal way with a given frequencyω f and amplitudeÃ f . The use of a quadratic instead of a sine function gives a considerable speed advantage in simulations, as the prediction of collision times with the moving walls becomes substantially faster. Test simulations were performed using a sine function for exemplary cases, and no significant differences were observed [22]. Furthermore, we considerably varied the collision parameters and found essential phenomena to be robust, although friction was observed to play a significant role in the dynamics.
We now introduce dimensionless variables, which will be used for the rest of the text. The depth of the box is fixed, l y ≡l y /d = 5, and the horizontal width l x ≡l x /d is varied in the [5,100] region. N is always taken so that the number of filling layers F ≡ Nd 2 /l xly = N/l x l y = 12, which implies that N varies in the [300, 6000] range. Three different oscillation amplitudes are considered, A f ≡Ã f /d ∈ {0.4, 1.0, 4.0}. This allows us to compare with previous results, obtained for A f = 4.0, as also to extrapolate to lower amplitudes, where the vibrating bottom wall can be considered as a spatially fixed source of energy (i.e. a temperature boundary condition).
The dimensionless gravity timescale is given by t g ≡t/t g , witht g ≡ (d/g) 1/2 . Correspondingly, the dimensionless oscillation frequency ω f is scaled as ω f ≡ω ftg = ω f (d/g) 1/2 . Nevertheless, it is almost always more meaningful to measure time in periods of box oscillations,T = 2π/ω f , and thus we use t ≡t/T . In order to compare simulations with different energy injection parameters the dimensionless shaking strength is used, Finally, the mass scale m =m/m p is set to unity by takingm as the mass of one particle,m =m p .
Simulations are generally run for 10 5 T = 10 5 (2π/ω f ), unless otherwise stated. Particles are initially arranged in a strongly perturbed, low density crystalline state. We confirmed that this initial configuration has no influence on the steady dynamics by running a few simulations using the end state of the previous simulation as the initial configuration. Contrary to the experiments realised in [17], where the frequency of shaking is continuously increased, the energy injection parameters are kept fixed during any given simulation.

Phase Space
In order to validate our simulations, and explore further previous research, we first focus on the A f = 4.0 case, where the comparison with previous experiments and softparticle simulations undertaken by Eshuis et al. [17] is straightforward. Event-driven simulations are able to reproduce all previously observed states, as can be seen in the phase diagram in the {l x , S} space presented in Figure 2a. Furthermore, a quantitative comparison is possible by looking at the transition points in the l x = 100 case, where the experiments were realised. There is excellent quantitative agreement, within 5%, for the bouncing bed-undulations and the Leidenfrost-convection transition points, but a 30% error in the undulations-Leidenfrost one. The deviations could in part be explained by the nature of the transitions, as they are not sharp and are seen to present wide ranges of metastability. This makes it harder to define a precise transition point value, and motivates the use of transition regions, which we show in gray.
As can be seen from Figure 2a, for l x > 20, and as S is increased (keeping A f = 4.0), the system goes through a sequence of different non-equilibrium stable states: bouncing bed, bursts, undulations, Leidenfrost, convection and gaseous (S > 400, not shown). Some of these states disappear or appear as l x and A f are varied, but their relative order remains. In the following, we will focus on the Leidenfrost and convective states, where LFOs take place. For a description of the other states, we refer the reader to [11,17].
The density inverted Leidenfrost state owes its name to the analogous liquid-overvapour phenomena, where a thin layer of vapour over a hot surface significantly slows the evaporation of the droplet above it, by keeping it floating over the hot surface [23]. Figure 2b shows the packing fraction φ and the granular temperature T g as a function of z, for different amplitudes and frequencies of oscillation, all in the Leidenfrost state. The granular temperature is defined as twice the fluctuating kinetic energy per degree of freedom: where r i is the position of particle i, v i its velocity, and V the average velocity field. Indeed, a low temperature, high density region is suspended over a low density, high temperature one. Notice that the difference in density between the solid and gaseous regions is greater for higher ω f (red vs. black), but lower for higher A f (blue vs. red): these features will be relevant in our model discussion for the validity regions of a density profile approximation.
When S is further increased, the density of the solid region is seen to progressively decrease, leading to a buoyancy driven convective state (see Figure 1). Horizontal homogeneity is lost, leading to low density regions where particles go up and circulate around high density regions, where particles agglomerate and move mainly in the horizontal directions, towards the low density regions. The number of convection rolls (R) diminishes with increasing ω f , until the energy input is so high that particle motion is essentially uncorrelated and the system enters the gaseous state (S > 400, data not shown).
We now turn our attention to the lower amplitude regions. Figure 3 shows a phase diagram again in the {l x , S} parameter space, for different shaking amplitudes A f . As observed previously [17], and confirmed here for a wider range of parameters, the dimensionless shaking strength S is a better parameter than the dimensionless acceleration, Γ ≡Ã fωf 2 /g = A f ω 2 f for the characterisation of the Leidenfrost-convection transition. On the other hand, the transition points of bouncing bed-Leidenfrost (or undulations-Leidenfrost for A f = 4.0) vary significantly with S, but stay within 5% when compared in Γ. In general terms, the most significant influence of reducing A f is the disappearance of the bursts and undulations states; the large amplitude of the box oscillation plays a dominant role in the dynamics of these states.
We briefly remark that simulations were done until l x = 400, and no new states were observed, except for the coexistence of convection and Leidenfrost states for l x ≥ 200. The possibility of this coexistence provides new insight into the nature of the Leidenfrostconvection transition, and suggests further research.
If the length is reduced further, below the l x = 20 limit, the frequency needed to trigger convection progressively increases, until at l x ∼ 10 (a value slightly dependent on A f , see Figure 3) no convection was observed even for S = 10 4 . For A f = 4.0, undulations and bursts are also frustrated by the small size of the container. It is in this geometry that it becomes possible to observe the Leidenfrost state for higher S, where low-frequency oscillations (LFOs) can be directly observed and eventually, as S is increased, dominate the collective dynamics of the system.

Low-Frequency Oscillations (LFOs)
Finally, we reach the column limit, where l x = l y = 5. In order to study LFOs the evolution of the vertical centre of mass of the particles is considered, z cm (t). Figure  4a shows z cm (t) for fixed A f = 1.0 and several different S. For comparison, nonstroboscopical and stroboscopical z cm (t) are shown for the S = 64 and S = 400 cases: the distinct high and low frequencies become immediately recognisable. The amplitude of the oscillations is seen to increase from the ∆z cm ≈ 1 to the ∆z cm ≈ 10 order, and present an appreciable regularity in time. While at S = 64 both oscillations are comparable in amplitude, and thus very hard to identify from direct observation, at S = 400 they have become clearly differentiable. Although LFOs are seen to be fairly chaotic (recall that there are only 300 particles in the column geometry, hence fluctuations play a leading role), we characterise them by a constant amplitude A 0 and a single frequency ω 0 , as an initial first order description.
First, let us focus on the frequency of the LFOs, ω 0 , which is clearly recognisable from the power spectra of z cm (t), presented in Figure 4b. The spectra are obtained by taking the discrete fast Fourier transform of z cm (t) over 20000T after an initial transient of 1000T , with a sampling rate of 0.05T . An average is then taken over 10 simulations with identical parameters but different initial conditions; although the shape and peaks are already recognisable from single simulations, the ensemble averaging reduces the noise considerably. The time window, the sampling rate and the transient time were varied and no significant differences were observed.
All spectra present two main features: the expected delta-like peak at ω f and its harmonics, and a broad peak one to two orders of magnitude lower, corresponding to the LFOs. The LFO frequency, ω 0 , is defined as the frequency of the maximum of this broad peak. After observing the different spectra it becomes evident that ω 0 depends on the energy injection parameters. Figure 5a shows ω 0 (S) for different l x and A f , remarkably scaling all LFO data. Notice that ω 0 decreases as S increases, i.e., the collective grain movement becomes slower as the shaking gets faster. The decay is faster than inverse linear, and can be fitted by a − 3 2 power with a 5% error (not shown). Let us also notice that the length of the container makes no discernible difference, as long as the system stays in the Leidenfrost state; the decreased data in the l x = 20 case are due to the Leidenfrost-convection transition. The collapse of the different amplitude curves is very good for A f = 0.4 and A f = 1.0, while for A f = 4.0 data slightly deviates. We interpret this decrease as the influence of the undulations state in the Leidenfrost regime; notice that for S ∼ 64 and A f = 4.0 the system is almost at the boundary between both states (see Figure 2a).
In order to quantify the relevance of the LFOs, we define the relative intensity of the ω 0 peak, I 0 , as the normalised distance from the low frequencies asymptotic limit to the maximum of the broad peak. Figure 5b shows I 0 (S) for different l x and A f . Although the dependency is not straightforward, it can be seen that LFOs become increasingly distinguishable from other movements until S ∼ 144, after which there is a decline, except for the highest amplitude case. Already at S = 25 oscillations should be discernible in the spectra as a peak twice as big as the low-frequency asymptotic limit. A f is seen to have a pronounced effect on I 0 ; higher amplitudes of oscillation lead to more pronounced LFO peaks.
Finally, we define the amplitude of the LFOs, A 0 , as the standard deviation of z cm (t): A 0 ≡ σ(z cm (t)). Data is considered only after t = 1000, to disregard transient states. Figure 5c shows A 0 (S) increasing in an almost linear way. The curves coincide, within their error, for A f = 0.4 and A f = 1.0, while for all other cases A 0 is consistently smaller. Nevertheless, S makes all curves comparable, further confirming its relevance for this system.

LFO's in convective state
We now consider in detail the peculiar change of behaviour of ω 0 (S) and A 0 (S) for S ∼ 144 in the l x = 20 case. This is a sign of the Leidenfrost-convection transition, still present at this container length (see Figure 2a). During convection, z cm becomes a less relevant quantity, as there is no longer horizontal homogeneity. Nevertheless it is still possible to identify LFOs, even if the oscillations are entangled with the convective flow. The presence of LFOs in the convective regime should not be surprising if one notices that it also presents the essential feature of the Leidenfrost state: a high density, low temperature region suspended over a low density, highly agitated one, although there is an additional low density, highly convective zone above. Our model, derived in Section 3 below, suggests that when density inversion is present, LFOs exist. Figure 6 presents several different fields and snapshots that show that, indeed, density inversion is present in the convective regime, in addition to the horizontal inhomogeneity. All data is taken from the same simulation, and fields are time-averaged over 100T after an initial transient of 1000T , with data taken every 0.05T . The average velocity field, Figure 6a, clearly shows the presence of convective flow, with a small downwards band and a wider upwards region. Particles agglomerate at the bottom of the downwards flux side, as can be seen from the average density field (Figure 6b), and the two snapshots ( Figures  6d and 6e). This happens when downwards and upwards particles collide, leading to a high granular temperature region (Figure 6c). Note, then, that both sides correspond to low density, high temperature regions sustaining high density, lower temperature ones, although the density and temperature profiles vary considerably from left to right. The profile is more similar to the Leidenfrost case in the upwards flow region (left in the shown figures), as in the downwards flow region the high density area presents a comparable, although lower temperature to the low density region below.

Summary
Having possible experimental realisations in mind, the general picture is that LFOs are easier to observe for higher amplitude and frequencies of oscillation of the box, while keeping l x = l y small; it is at these configurations that LFOs have the highest amplitudes and better defined frequencies, as quantified by A 0 and I 0 , respectively. Let us now remember that at this limit we also observed the most clear phase separation in the Leidenfrost state, with distinct low and high density regions. In our model, presented next, the separation of the phases and the confinement of the system to a one-dimensional geometry implies the existence of LFOs, and the frequency is essentially determined by the ratio of the low and high densities.

Continuum model
After observing the collective movement of the particles in the column geometry, an oscillator-like description naturally comes to mind. The two coexisting frequencies observed in the spectra suggest a forced oscillator model, with clearly defined forcing and response frequencies. In the following we derive such frequency behaviour from a continuum description of the granular media. We begin by considering Cauchy's equations for mass and momentum conservation: where ρ corresponds to the material density, u = {u, v, w} is the velocity vector,σ the stress tensor and g the gravitational acceleration in the downwards direction, −ẑ. Furthermore, the material derivative is defined as D t ≡ ∂ t + u · ∇. We consider the same scaling as in simulations, with length scales in units of particle diametersd, time units given by gravityt g = (d/g) 1/2 , as alsoρ p , taken as the mass density of a single particle, ρ p =m p /Ṽ p , withṼ p = 1 6 πd 3 . As has been observed in simulations, the dynamics of the system in the column limit is effectively one-dimensional. This immediately suggest the consideration of ρ = ρ(z, t), u = w(z, t)ẑ andσ = σ zz (z, t). Substituting in (1) yields Furthermore, expanding (2), and using (3), one reaches a one-dimensional momentum conservation equation

Two phases approximation
In order to solve (4) it would be necessary to know both the density and the velocity profiles, ρ(z, t) and w(z, t). Our approach consists in eliminating the z-dependence from (4) by integrating in the vertical direction, and taking a first order approximation of the density profile ρ(z, t), and average values for the vertical velocity profile w(z, t). We begin by integrating (4) in the vertical direction with the bottom boundary, b(t), and top boundary, s(t), dependent on time, due to the movement of the bottom wall and the free surface at the top. The approximation of ρ(z) consists in dividing the system in two separate, constant density regions, inspired by the measured Leidenfrost state density profile. Let us remember that this approximation becomes increasingly better as S increases and A f decreases, as shown in Figure 2b. Consequently, a low density region is defined where ρ(z, t) = ρ g (t), for z < ξ; and a high density region where ρ(z, t) = ρ s (t), for z > ξ, with ξ = ξ(t) the position of the interface between the two regions. Figure 7 shows a schematic representation of this approximation, and the origin of its motivation. It then follows that the first integral in (5) can be expanded as Analogously, the third integral in (5) becomes with h g (t) ≡ ξ(t) − b(t) the height of the gaseous region, and h s (t) ≡ s(t) − ξ(t) the height of the solid region.
Notice that the second integral in (5), corresponding to the stress term, is a perfect integral, and thus only the stress boundary conditions are needed for its evaluation. We assume the stress through the system to be continuous in z, and thus it is not necessary to evaluate σ zz at the interface position ξ(t). Thus, from (5) we finally obtain:

Boundary conditions
It now becomes necessary to specify the boundary conditions. The shaking of the container implies that b(t) =Â fm sin(ω fm t), withÂ fm and ω fm the amplitude and frequency of energy injection in the model. At the top, s(t), we consider a free surface, and thus the kinematic boundary conditions are given by Furthermore, the stress at the bottom and top of the granular media are needed. The free surface at the top is straightforward: σ zz (z = s) = 0. At the bottom, on the other hand, we divide the stress contribution in two: mean (σ 0 b ) and fluctuating (σ b ) terms, where the mean term is straightforward: σ 0 b = Mg/η, with M the total mass of the system, M = Nm; and η the area of the base of the container, η = l x l y .
For the fluctuating part of the stress, σ b , we first consider the force applied to the granular medium by the moving bottom: with m b the mass being pushed by the bottom wall. In order to obtain m b , let us consider a moving platform of surface area η pushing an ideal, incompressible gas of density ρ g , in analogy to the moving box and the low density region observed in our system. Accordingly, the mass pushed by the box in time is given by dm p = ρ g ηv b dt.
Notice that this is valid for high S, where gravity effects on the dynamics of the particles can be ignored. Integrating, we directly get that m p = ρ g ηÂ fm sin(ω fm t). Substituting in (11): fm ω 2 fm (− sin 2 (ω fm t) + cos 2 (ω fm t)) = ρ g ηÂ 2 fm ω 2 fm cos(2ω fm t). (12) Notice that we have naturally obtainedÂ 2 fm ω 2 fm ≡ S m g as the amplitude of the force applied by the oscillating bottom, further suggesting that the shaking strength is the relevant parameter for the system in the high S limit. It then follows, from (11), that: Finally, substituting the stress boundary values in (8), we obtain:

Height averaging
The remaining two integrals in (14) involve the velocity profile, w = w(z, t), which varies in the vertical direction. In order to solve these integrals we height-average, that is, for a given quantity f (z), we consider its average valuē Notice that, from the first integral in (14), f would correspond to ∂ t w. Thus, before applying (15), we express the integral as a total time derivative. Considering that the boundaries are time dependent, it becomes necessary to use Leibniz integration rule, and thus the first integral in (14) can be expressed as Analogously, the second integral in (14) becomes, after using (10), Substituting (16) and (17) in (14), and using (15), we finally obtain: Based on the behaviour observed in simulations, we now assume that the high density region is incompressible. This implies that d t h s = 0, as also that the velocity Figure 7. From left to right, snapshots from simulations showing an LFO period, at phases 0, π/2, π, 3π/2 and 2π; the corresponding time averaged density profile, a representation of the two phases approximation made for the continuum equations, and finally a schematic representation of the model. The dashed line shows the position of the centre of mass, z cm , which in the model corresponds to the position of the interface between the two phases, ξ, which also corresponds to the position of the mass of a forced harmonic oscillator.
of the continuum media at the interface position is equivalent to the velocity of the interface, and hence to the velocity of the surface, that is, w(z = ξ) =w s = d t s = d t ξ.
Furthermore, we now use the fact that h g (t) = ξ(t) − b(t). Thus, substituting and dividing by ρ s , we obtain

First order approximations
It now becomes relevant to consider the relative importance of each of these terms, in the region of phase space where simulations show that LFOs are present, that is, for S ≫ 1. First, we consider that ρ g /ρ s ∼ O(ǫ), a condition that holds better for S ≫ 1 and low A f , as shown in Figure 2b. On the other hand, ξ ∼ 10 ∼ O(1/ǫ), as can be seen from Figure 4. Furthermore, we measure from simulations that δ t ξ ∼ 0.2 ∼ O(ǫ) and δ tt ξ ∼ 0.1 ∼ O(ǫ), meaning that the dynamics of the LFOs are considerably lower than the typical velocity of grain diameters per gravity timescale, as can also be deduced by the previously obtained frequencies ω 0 . Let us also notice that h s ∼ h g ∼ 8 ∼ O(1/ǫ), again, from Figure 2b. Finally, from simulations we obtain thatw g ∼ 0.2 ∼ O(ǫ), and d twg ∼ 0.02 ∼ O(ǫ 2 ). Taking into account all these considerations, it becomes straightforward to see that the first term in (20) is O(ǫ 3 ), the second term is at most O(ǫ 2 ), the third is then O(ǫ), and the fourth term is O (1). Moreover, all terms on the right side are O(1). Thus, disregarding small terms in (20), after dividing by h s , we obtain where we have defined the mass of the solid region per unit base area η, m s ≡ h s ρ s , and the equivalent of the gaseous region, m g = h g ρ g . Equation (21) corresponds to a forced harmonic oscillator equation of the form: with natural frequency amplitude of forcing F 0 = 3 2 gρ g S m /m s , and constant C = 1 2 ρ g S m /m s + g(m g /m s − 1).

Model and simulations comparison
We have shown that, considering Cauchy's equations for continuum media, and making assumptions in concordance to the observed granular Leidenfrost state, the system becomes equivalent to a simple forced harmonic oscillator, expressed by (22). In this case, ξ is the displacement of the centre of mass around the equilibrium position at 0, ω 0m the natural frequency of the system, and F 0 and ω fm the amplitude and frequency of the forcing. The analogy of the forcing with the granular column is straightforward: ω 0m and A 0m would be equivalent to ω 0 and A 0 , respectively. Furthermore, we choose ξ to correspond to z cm , in order to directly compare with previous measurements. Notice that the natural frequency ω 0m does not explicitly depend on the forcing frequency ω fm , as can be seen in (23). The implicit dependence comes from the variation of ρ g and m s with ω fm , as observed in simulations, where, for fixed S, ρ g /ρ s increases with ω f , giving the correct inverse proportionality of ω 0m with ω fm . Therefore, in order to obtain a frequency from the model, only ρ g and m s need to be specified, which we measure from simulations.
Both quantities can be obtained from ρ(z), the density in the granular column as a function of height. In order to obtain an accurate average, we consider ρ * ≡ ρ(z − z cm (t), t), which makes all profiles directly comparable. This is analogous, in the model, to centering the profiles at the interface between the two distinct regions. It is then straightforward to compute ρ g as the average value of the density for z < z cm . On the other hand, m s we take as the total mass for z > z cm , taking care not to count particles that are in free flight above the solid region, as they do not have influence on the oscillator dynamics. This implies that although the center of our profiles is z cm , m g is not equal to m s . and A = 4.0 (c). All systems have l x = 5, F = 12. Simulation (black) corresponds to frequencies obtained from fast-Fourier transform of the simulation data, while continuum and thermodynamic/kinetic theory data points (blue and red, respectively) are obtained from models presented in sections 3 and 4, respectively, using data acquired from simulations.
The comparison between the frequencies obtained in simulations and from the model is presented in Figures 8a-c. For low amplitudes, A f ≤ 1.0, and high frequencies, S ≥ 144, the agreement between the frequencies is within the error bars. For lower S, or higher amplitudes, the assumption of two distinct phases, as also the approximation of ξ ∼ 1/ǫ, become less justified, resulting in the model consistently underpredicting the frequencies, with more than 50% disagreement at the point of the bouncing bed-Leidenfrost transition for A = 4.0. We believe that the prediction could be improved by considering more complex density profiles, as also by including terms of lower orders, although this exceeds the scope of our work. In general terms, the resulting onedimensional model turns out to be a remarkable well approximation for high ω f and low A f , showing that this many-particle, out-of-equilibrium system actually behaves as a regular forced harmonic oscillator when confined in a column, in the corresponding energy injection region.

Thermodynamic model
Remarkably, it is possible to obtain another accurate expression for ω 0 using a completely different approach, considering basic concepts from thermodynamics. Assuming a spring-like behaviour, the natural frequency of our medium is given by ω 2 0 = k/m s , with k the stiffness constant of the spring-like medium, and m s the mass sustained by the spring. We also know that k = ηB/h g , with B the bulk modulus, η the area of the spring, and h g it's height at rest. Assuming an adiabatic ideal gas, it is possible to relate the bulk modulus with the pressure, B = γP 0 , with γ the adiabatic index.
Notice that the ideal gas approximation is being used only for the gaseous region of the Leidenfrost regime, where densities are low and no significant correlation of the particles is observed. We then obtain: The pressure P 0 is taken as the force caused by the solid mass the spring sustains, m s , divided by the area of the container: P 0 = m s g/η. Thus, we finally reach: This significantly simple expression is remarkably accurate when compared with simulation measurements. Figures 8a-c also show ω 0 (S) for this model, taking h g to be the same as in the previous section, z cm . The adiabatic index is considered as γ = 1.67, the theoretical value for an ideal monoatomic gas. The agreement is again within error bars for high frequencies, and deviates considerably for lower frequencies, except for the A = 1.0 case, where low frequencies are also captured. Relating the two obtained LFOs frequencies, equations (23) and (24), one obtains γ = m g /m s ; the interpretation of this result remains a challenge.

Conclusions
A vertically vibrated bed of grains presents low-frequency oscillations (LFOs) due to the decoupling of the driving frequency and the dynamics of a high density region suspended by a lower density one. The relevance of these oscillations increases as the distinction between the two densities increases, that is, proportional to the frequency and inversely to the amplitude of oscillation of the system container. The LFO frequencies are inversely proportional to the driving frequency, and follow a common power law for a range of amplitudes. The amplitude of the oscillations, on the other hand, increases in an almost linear way with the frequency.
Event-driven simulations give an overall excellent qualitative and quantitative agreement with experiments and soft-particle simulations done in wider systems, although they show discrepancies in some critical transition values. We remark that the hard-sphere approximation can be meaningful even in systems with very highdensity regions, as present in the Leidenfrost state. The considerable speed advantage makes it extremely useful, and sometimes the only means to systematically study high dimensional parameter spaces.
Starting from Cauchy's equations for conservation of mass and momentum, integrating in the vertical direction and assuming two distinct low and high constant density regions, it is possible to reproduce the frequency behaviour observed in simulations. That is, a forced harmonic oscillator, with the natural frequency proportional to the ratio of the densities. This simple model is able to predict the natural LFO frequency for high excitation frequencies, where in fact the two phases are well separated. The non-linear terms, discarded in our analysis, should provide the necessary corrections for lower frequencies, as well as the consideration of a more realistic density profile. A second approach, using thermodynamic arguments, also gives a remarkably accurate expression for the frequencies, although in this case just a simple mass-spring system behaviour was assumed. The quantitative agreement of both models is nevertheless remarkable, taking into account the low number of particles involved, and the presence of very high and low density regions.
Further insight could be gained by appropriately coarse graining the granular medium in order to obtain stress fields, which would directly relate both models. A point of interest, not studied here, is how well do kinetic theory predictions hold in such a system, taking into account the reduced container size, the small number of particles and the presence of considerably different densities. Current work is being done on verifying the consistency of macroscopic fields obtained by theoretical arguments and coarse-grained simulational data.
We suggest that LFOs, here shown to be ubiquitous to vertically vibrated density inverted systems, could play a fundamental role in the Leidenfrost-convection transition. More specifically, LFOs could be the primary source of density fluctuations observed before convection is triggered, when one region of the system oscillates at a different phase than another. Understanding this will need further simulation and experimental research.