{'fraction_non_alphanumeric': 0.040161667730270155, 'fraction_numerical': 0.01076366730482876, 'mean_word_length': 4.421125461254612, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Development of a density inversion in driven granular gases

Granular materials fluidized by a rapidly vibrating bottom plate often develop a fascinating density inversion: a heavier layer of granulate supported by a lower-density region. We employ the Navier-Stokes granular hydrodynamics to follow a density inversion as it develops in time. Assuming a dilute low-Mach-number flow, we derive a reduced time-dependent model of the late stage of the dynamics. The model looks especially simple in the Lagrangian coordinates. The time-dependent solution describes the growth of a density peak at an intermediate height. A transient temperature minimum is predicted to develop in the region of the density peak. The temperature minimum disappears at later times, as the system approaches the steady state. At late times, the predictions of the low-Mach-number model are in good agreement with a numerical solution of the full hydrodynamic equations. At an early stage of the dynamics, pressure oscillations are predicted.


Introduction
An ensemble of hard spheres in motion, that collide inelastically and are characterized by a constant coefficient of normal restitution, represents the simplest model of granular gas [1]. It also provides an excellent example of complexity. Indeed, the properties of the granular gas as a whole, such as clustering [2,3,4,5], come out as non-trivial consequences of the simple algebraic laws that govern the velocities of the binary collisions of individual particles. Furthermore, being intrinsically far from thermal equilibrium, granular gas defies statistical mechanics, the cornerstone of the equilibrium science. As the result, formulating a universally applicable continuum theory of granular gas, even in the dilute limit, is not a simple task [6,7]. The use of the Boltzmann (or Enskog) equation, the starting point of a systematic derivation of the constitutive relations, is based on the Molecular Chaos hypothesis. This hypothesis is justified, in the dilute limit, for an ensemble of elastic hard spheres. Its use for inelastic hard spheres is, rigorously speaking, an uncontrolled assumption. Indeed, inelasticity of the particle collisions introduces inter-particle correlations which can invalidate the Molecular Chaos hypothesis. The correlations become more pronounced as the inelasticity of the collision increases. On the contrary, for nearly elastic collisions, 1 − r 2 ≪ 1 (where r is the coefficient of normal restitution) the correlations are small, and the Boltzmann equation can be used safely (again, in the dilute limit).
An additional crucial assumption, made in the process of the derivation of the hydrodynamics from the Boltzmann equation, is scale separation. The hydrodynamics is valid if the mean free path of the particles is much less than any characteristic length scale, and the mean time between two consecutive collisions is much less than any characteristic time scale, described hydrodynamically. This assumption (as well as the dilute gas assumption, if used) should be verified, in every specific system, after the hydrodynamic problem is solved and the characteristic length and time scales determined. One more criterion for the validity of granular hydrodynamics is the smallness of fluctuations, driven by the noise caused by the discrete nature of the particles. The accuracy of hydrodynamics improves with an increase of the total number of particles in the system N [8].
Once hydrodynamics is valid, it is highly rewarding using it. Hydrodynamics has a great predictive power and helps to identify important collective phenomena (shear flows and vortices, shocks, different modes of clustering flows, phase separation etc.) that are impossible to perceive in the language of individual particles. Once identified, these collective phenomena can then be investigated in experiment and simulations in more general (not necessarily hydrodynamic) formulations. A recent example: employing a variant of granular hydrodynamics, Livne et al. [9] predicted a novel phase separation instability in a system of inelastic hard spheres, driven by a rapidly vibrating plate at zero gravity. This system had been investigated earlier by many workers. However, the phase separation instability, and a plethora of interesting effects accompanying it, were overlooked. Undoubtedly, the reason for that was that hydrodynamics had not been fully appreciated and exploited [10].
In this paper we employ the Navier-Stokes granular hydrodynamics for a theoretical analysis of the process of formation of a density inversion in a dilute granular gas driven by a rapidly vibrating bottom plate in a gravity field. The phenomenon of a density inversion in the steady state of this system is well known in the granular community. On the contrary, the dynamics of the formation of the density inversion has never been investigated previously. As we will see, this dynamics is quite instructive. The hydrodynamic theory yields predictions of time-dependent quantities that can be checked in molecular dynamics simulations and experiment. Our additional motivation here is pedagogical. While dealing with this problem, we will present two useful hydrodynamic techniques. The first is reduction of the order of the hydrodynamic equations, based on a hierarchy of the length/time scales in the problem. The second technique is the use of the Lagrangian coordinates.
Here is a layout of the remainder of this paper. In Sec. 2 we will write down a set of hydrodynamic equations and boundary conditions that describe the dynamics of formation of a density inversion. By introducing scaled variables, we will delineate the two scaled governing parameters of the problem. In Sec. 3 we briefly review the properties of the steady state and the condition for a density inversion. In Sec. 4 we assume a low-Mach-number flow, derive a reduced hydrodynamic model and use it for an investigation of the late-time dynamics of the formation of the density inversion. Section 5 reports a numerical solution of the full hydrodynamic equations and compares the numerical results with the predictions of the low-Mach-number model. A discussion of the results is presented in Sec. 6.

The model problem and hydrodynamic equations
We will adopt the simplest possible model, in two dimensions, that exhibits a density inversion. Let N ≫ 1 identical nearly elastic hard disks of diameter d and mass m = 1 move without friction in a two-dimensional box of width L x and infinite height. The driving base is at y = 0, the side walls are at x = 0 and x = L x . The gravity acceleration g is in the negative y direction. The rapid vibrations of the base can be modeled, in a simplified way, by prescribing a constant granular temperature T = T 0 at y = 0 (we measure the granular temperature in the units of the velocity squared). This means that, upon collision with the base, the particle velocity is drawn from a Maxwell's distribution with temperature T 0 . The kinetic energy of the particles is being lost by inelastic hard-core collisions parameterized by a constant coefficient of normal restitution r. Particle collisions with the side walls are assumed elastic. Alternatively, one can specify periodic boundary conditions in the x-direction.
Introducing a hydrodynamic description, we deal with three coarse-grained fields: the particle number density n(r, t), mean flow velocity v(r, t) and granular temperature T (r, t). In addition, we make two important assumptions. First, we assume that that the collisions are nearly elastic: 1−r 2 ≪ 1. Second, we assume that the granular gas is dilute everywhere (including the density peak region, see below), that is n is small enough compared to the closepacking density n c . These assumptions enable us to use the classic version of the Navier-Stokes granular hydrodynamics [11], and to limit ourselves to the leading order terms in n/n c in the constitutive relations. The dilute gas condition will be checked a posteriori.
It has been shown recently that, when this system is large enough in the lateral direction, thermal convection may occur [12,13,14]. In this work we will assume that L x is sufficiently small, and the convection is suppressed by the viscosity and lateral heat conduction. Therefore, we assume flow solely in the vertical direction. Therefore, the governing equations are: Here P = −nT + η∂v/∂y is the vertical component of the stress tensor, q = −κ∂T /∂y is the heat flux, and I is the heat loss rate by the inelastic collisions. In the nearly elastic limit, the leading order contributions to the viscosity η and heat conductivity κ are given by their classical values for for elastic hard disks: The heat loss rate is the following: The boundary conditions for the problem include the fixed temperature and zero velocity v(y = 0, t) = 0 at the base. Now, as the total number of grains N is constant: the mass flux should vanish at infinity: Similarly, one should assume that the momentum flux at infinity is zero: as well as the energy flux: We refer the reader to the book of Landau and Lifshitz [15] for a detailed discussion of the divergent forms of the hydrodynamic equations, and for the expressions for the momentum and energy flux density used in Eqs. (8) and (9). What are the characteristic length/time scales of the problem? The fixed base temperature T 0 and gravity acceleration g define a characteristic macroscopic length scale λ = T 0 /g. The characteristic number density of the gas is therefore n 0 = N/(λL x ). Correspondingly, there are three characteristic time scales in the problem. The fast time scale of the problem, τ f = T 1/2 0 /g, is independent of the particle size. The heat conduction time scale is τ = (πλ/g) 1/2 (N d/L x ). The inelastic heat loss time scale τ l can be defined as τ −1 How do we know that the time scale τ f is indeed fast? Before we answer this question, let us first make a preliminary evaluation of the validity of hydrodynamics. We should demand that the characteristic mean free path of the particles ∼ (n 0 d) −1 be much less than the macroscopic length scale λ. This immediately follows L x /(N d) ≪ 1. Notice that the quantity N d/L x is of the order of the number of monolayers at rest (that is, when the system is not fluidized). Therefore, the validity of the hydrodynamics requires that the number of monolayers be much larger than 1. Now we compare the characteristic time scales τ f and τ introduced above. We see that the ratio τ f /τ is of the order of L x /(N d), which should be much less than unity if we want the hydrodynamics to be valid.
To make a full use of hydrodynamics, we introduce scaled variables. Let us measure time in the units of the heat conduction time τ , the coordinate in the units of λ, the temperature in the units of T 0 , the density in the units of n 0 and the mean flow velocity in the units of λ/τ . In the scaled variables, the governing equations (1) where ε = L x π 1/2 N d ≪ 1 and Λ = (1 − r 2 ) 1/2 2 ε are the two governing scaled parameters of the problem. Parameter Λ is of the order of the square root of the ratio between the heat conduction time scale and inelastic heat loss time scale.
The reduction of the four-dimensional space of governing parameters d, N, r and L x to the (ε, Λ) -plane is an immediate benefit from using hydrodynamics. Indeed, series of experiments or simulations with different d, N, r and L x , but with the same ε and Λ, should produce the same dynamics in the scaled variables. We will see in the following that, in the late stage of the dynamics (and in the steady state), the only relevant parameter is Λ.
Like in many other one-dimensional gasdynamic problems [17], it is convenient to make a transformation of variables from the Eulerian coordinate y to the Lagrangian mass coordinate m(y, t) = y 0 n(y ′ , t) dy ′ . The physical meaning of the Lagrangian mass coordinate is very simple: it is the mass content of the gas on the Eulerian interval (0, y) per unit length in the x-direction. Notice that the infinite interval 0 ≤ y < ∞ in the Eulerian coordinates is mapped, by this coordinate transformation, into a finite interval 0 ≤ m < 1, at all times. Going over to the Lagrangian coordinate m, we can rewrite Eqs. (10)-(12) as follows: ∂ ∂t The boundary conditions become T = 1 and v = 0 at m = 0, and n ∂T /∂m = n ∂v/∂m = 0 at m = 1.
Notice that the transformation from the Eulerian to Lagrangian coordinate involves the density n(y, t) that is unknown a priori. We should not worry about it. After having determined the density n(m, t) in the Lagrangian coordinates, one can go back, at any time t, to the Eulerian coordinate by using the relation

Steady state profiles and density inversion
Let us briefly review the steady state solution of the problem [13,14,18]. To obtain the steady state profiles n s (m) and T s (m) we simply put ∂/∂t = v = 0 in Eqs. (16)- (18). Then Eq. (16) is obeyed automatically, while Eqs. (17) and and respectively. As we see, the steady state is described by a single scaled parameter Λ. Integrating Eq. (20) with respect to m and using the boundary condition at m = 1 (that becomes simply nT = 0), we obtain n s T s = 1 − m. Now we express the density through the temperature: n s (m) = (1 − m)/T s (m). Using this relation in Eq. (21), we arrive at the following linear equation [13,14,18] for ω s (m) = T 1/2 s (m): where primes stand for the m-derivatives. This is a Bessel equation whose general solution is a linear combination of the modified Bessel functions of the first kind .
To transform back to the Eulerian coordinate y, we should evaluate the integral numerically. Figures 2 and 3 show the steady state temperature and density profiles in the Lagrangian and Eulerian coordinates for Λ = √ 10. Equation (24) enables one to find the critical (minimum) value of Λ for the density inversion. At small Λ the density decreases monotonically with m (and therefore with y). The density inversion is born, at m = 0, at Λ = Λ c . Therefore, to find Λ c , we demand that the derivative n ′ s (m) vanish at m = 0. Using Eq. (24), we can reduce this condition to the equation I 0 (Λ) = 2Λ I 1 (Λ), first obtained in Ref. [18]. Its solution yields Λ c ≃ 1.06569 . . . . At Λ > Λ c the density reaches its maximum in the Lagrangian point m * = 1 − Λ c /Λ > 0. Therefore, the position of the density peak in the Lagrangian coordinate grows with Λ monotonically. The respective (scaled) density maximum value is Using Eq. (25), we can compute the scaled Eulerian coordinate y * at which the maximum density maximum occurs: The dependence of y * on Λ is non-monotonic, see Fig. 1. Very close to the density inversion birth (at Λ > Λ c ) this dependence is linear: y * ≃ Λ/Λ c − 1, while at Λ ≫ Λ c one obtains y * ≃ (2Λ) −1 . Actually, the integral in Eq. (27) can be calculated analytically, in terms of a hypergeometric function.

Formation of a density inversion: Low-Mach-Number flow
Now let us return to the time-dependent problem. Can we exploit the presence of a small parameter ε 2 in Eqs. (17) and (18) to simplify the equations? One can expect that the respective terms cannot be neglected during the rapid initial stage of the dynamics. The duration of this stage is of the order of the acoustic time or, in our scaled units, of the order of ε. A comparison of different terms in the momentum equation (17) shows that the typical scaling of the velocity during this stage is ε −1 ≫ 1. In the physical units, the flow velocity here is comparable to, or larger than, the speed of sound. Correspondingly, at early times the ε 2 -terms are comparable to, or larger than, the other terms and therefore should be kept. The situation changes at later times, determined by the inelastic heat losses and heat conduction. Here the characteristic Mach number of the flow becomes small, and we can neglect, in the leading order, the ε 2 -terms in Eqs. (17) and (18). As the result, the momentum equation (17) is reduced to the hydrostatic condition (20), while in the energy equation (18) one can drop the viscous heating term. The continuity equation (16) remains the same as before. The boundary conditions at m = 0 do not change, but those at m = 1 are reduced to n ∂T /∂m = 0. The boundary condition n ∂v/∂m = 0 should be dropped altogether.
What is the physical meaning of the reduced model? A low-Mach-number flow proceeds on a time scale that is much longer than the time needed for establishing the (approximate) hydrostatic equilibrium, so the flow does not violate the hydrostatics. Similar reductions of nonlinear time-dependent gasdynamic problems, that employ a separation of time/length scales, have been successfully used in the context of the dynamics of an optically thin plasma that is cooled by radiation [19]. The presence of a heat loss term in the energy equation makes the latter problem similar to problems of granular dynamics.
The hydrostatic condition immediately yields p(m, t) ≡ n(m, t)T (m, t) = 1 − m, like in the steady state solution! Using it together with the continuity equation (16), we can rewrite the energy equation as an evolution equation for ω(m, t) = T 1/2 (m, t): This nonlinear parabolic equation should be solved with two boundary conditions: ω(0, t) = 1 and (1 − m) ∂ω/∂m → 0 as m → 1. Notice that the only relevant scaled parameter in the low-Mach-number problem is Λ. Once Eq. (28) is solved, we can find n(m, t) and v(m, t): To transform back to the Eulerian coordinates we use Eq. (19). Equation (28) was solved numerically with a MATLAB PDE solver. Figure 2 shows an example of numerical solution, at different times, for Λ = √ 10. The initial condition for the temperature, T (m, t = 0) = 1, corresponds to a barometric density profile, which looks especially simple in the Lagrangian coordinate: n(m, t = 0) = 1 − m. The solution shows cooling of the granular gas because of the inelastic heat losses (see Fig. 2 a). Simultaneously, a density peak develops at an intermediate height (Fig. 2 b). The material flux (Fig. 2 c) is everywhere negative: the material slowly falls from above until the steady-state with a zero mean velocity is reached. An important feature of the time-dependent solution is a transient temperature minimum that develops at intermediate times. It is caused by the fact that the inelastic heat losses are very small at large heights, where the gas density is small. Therefore, the gas cools down relatively fast at intermediate heights, which causes a temperature minimum there. Then a downward heat flux develops that removes the heat from above, finally creating a monotonically decreasing steady state temperature profile. One can see from Fig. 2 that the steady state temperature and density profiles coincide with the analytic solutions. The respective solutions in the Eulerian coordinates are shown in Fig. 3.

Formation of a density inversion: Early times
The low-Mach-number theory is invalid at early times, when the hydrostatics has not settled down yet, and the mean flow velocity can reach, or even exceed, the local speed of sound. To investigate these early times, and to check the validity of the low-Mach-number theory at later times, we solved numerically the full set of hydrodynamic equations (16)- (18) in the Lagrangian form.
The parameter ε was taken to be √ 5 · 10 −2 ≃ 0.02236. Equations (16) and (17) were solved with an explicit numerical scheme, based on the standard finite-difference method described in Ref. [20]. The viscous term in Eq. (17) was integrated explicitly, using a two-steps method, which is both stable and accurate. The energy equation (18) was treated implicitly, with a conservative scheme.
We started with simple initial conditions: the constant temperature T (m, t = 0) = 1, the barometric density profile n(m, t = 0) = 1 − m and a zero mean velocity. As expected, the late-time results show good agreement with the predictions of the low-Mach-number theory. Figure 2 shows that, at latest times, the full-hydrodynamics solutions coincide with those obtained with the low-Mach-number theory. At an intermediate time one can see a local increase of the temperature at high altitudes, which is predicted by the full hydrodynamics, but missed by the low-Mach-number theory. The reason for this transient temperature maximum is the relatively large viscous heating that occurs at high altitudes because of a steep velocity gradient there.
The early-time results of the full numerical solution, transformed back to the Eulerian coordinate, are shown in Fig. 4. One can see that the development of the density peak, and the formation of the transient temperature minimum, start already at early times. A more detailed analysis does show differences between the fast and slow theories. Indeed, a necessary condition for the validity of the low-Mach-number theory is the hydrostatic balance condition. This condition looks especially simple in the Lagrangian coordinate: p(m, t) = 1 − m. That is, when the hydrostatics sets in, the plots of the pressure versus m, at different times, should collapse into a single straight line. Figure 5 shows the pressure plots at different times. One can see that, at late times, the hydrostatic condition sets in as expected. At early times pressure oscillations are observed, and hydrostatics does not hold. Figure 6 shows the pressure history at the bottom plate m = 0, or y = 0. One can see that the pressure undergoes oscillations on the fast time scale of the system O(ε). The oscillations rapidly decay with time and, starting from t ≃ 0.4, the pressure stays almost constant which implies a low-Mach-number flow.

Discussion
We have employed granular hydrodynamics to describe a time-dependent granular flow that gives rise to a density inversion in a vibrofluidized granular system. We have predicted the formation of a transient temperature minimum at intermediate times. The temperature minimum ultimately disappears as the system approaches the steady state with a density inversion. We have also predicted pressure oscillations at early times, when the hydrostatic condition is not yet satisfied.
We hope that the predicted phenomena will be observed in molecular dynamics (MD) simulations and in experiments with vibrofluidized granular beds. The MD simulations should start from the barometric density profile, corresponding to the temperature of the bottom plate. A natural way to initiate the initial state is to first run the simulation with elastic particle collisions. After the barometric profile sets in, the collision rule should be switched into inelastic regime.
The hydrodynamic theory considered in this work is expected to be quantitatively accurate for nearly elastic collisions, 1 − r 2 ≪ 1, and in the dilute limit. The latter condition demands that the maximum number density of the particles be less than, say, 0.2 n c , where n c = 2/( √ 3d 2 ) is the hexagonal close packing density. The maximum density peak is achieved at the steady state. Therefore, using Eq. (26) and going back to the "physical" density, we can write the dilute regime criterion as An alternative form of the same criterion is We expect that the theory is still valid qualitatively when this criterion is not satisfied, but the maximum density is still less than about 0.5 n c . To achieve a better accuracy for such moderately dense flows, one can use the constitutive relations of Jenkins and Richman [11] that account for excluded volume effects.