Fluctuating hydrodynamics and correlation lengths in a driven granular fluid
Static and dynamical structure factors for shear and longitudinal modes of the velocity and density fields are computed for a granular system fluidized by a stochastic bath with friction. Analytical expressions are obtained through fluctuating hydrodynamics and are successfully compared with numerical simulations up to a volume fraction . Hydrodynamic noise is the sum of external noise due to the bath and internal one due to collisions. Only the latter is assumed to satisfy the fluctuation-dissipation relation with the average granular temperature.
Static velocity structure factors and display a general non-constant behavior with two plateaux at large and small , representing the granular temperature and the bath temperature respectively. From this behavior, two different velocity correlation lengths are measured, both increasing as the packing fraction is raised. This growth of spatial order is in agreement with the behaviour of dynamical structure factors, the decay of which becomes slower and slower at increasing density.
pacs:45.70.-n, 51.20.+d, 05.40.-a, 47.57.Gc
Granular media display a wide catalog of non-equilibrium phenomena . These materials are constituted by a number of elementary constituents, grains of typical diameter between and mm. The number is usually large enough to allow, or require, a statistical treatment. Unfortunately, interactions are non-conservative, resulting in the failure of equilibrium statistical mechanics. Kinetic theories, from Boltzmann equation to hydrodynamics , together with numerical simulations , are the best tools to describe those systems and to compare with real experiments, with the caveat of a proper adaptation to the peculiarity of granular interactions.
One of the debated points of granular kinetic theories is the way noise should be added to hydrodynamics in order to describe mesoscopic fluctuations [4, 5]. This is a general problem in non-equilibrium systems  (e.g. sheared fluids ), but here is even more pressing, given the rather small number of particles in a granular system: one has typically , even in experiments, so that fluctuations can hardly be neglected. Moreover, in granular systems the dynamics is non-conservative and therefore Fluctuation-Dissipation relations do not hold in general [8, 9, 10, 11], with exceptions in driven dilute cases [12, 13, 14]. In non-dilute systems, it is therefore difficult even to define a temperature, making tricky the modelization of fluctuations [15, 16].
A comprehensive study of the fluctuating hydrodynamics of a driven granular fluid is presented here. Static and dynamical structure factors are computed analytically in the framework of linearized hydrodynamics, and compared with extensive numerical simulations. A very good agreement is found between analytical and numerical results in a wide range of parameters, implying that, for this kind of model, fluctuating hydrodynamics is able to describe large scale fluctuations in a satisfactory manner.
The peculiarity of the model we have studied, when compared to others present in the literature [17, 18], is the prescription for the stochastic bath used to keep the system at stationarity. In particular our thermostat is able to equilibrate the system also when collisions are elastic . This happens because, in addition to a random driving, our thermostat acts on the particles also through a finite drag, in such a way that the temperature of the thermostat , different form the kinetic temperature of the fluid , is always well defined.
A remarkable feature of our model, related to the kind of thermostat we use, is the finite extent of velocity correlations in space. Indeed, the characteristic shape we find for transverse and longitudinal velocity structure factors allows us to define two non-equilibrium correlation lengths, and , which are known functions of the kinematic and longitudinal viscosities, respectively. This means that, instead of sampling trajectories of the system, out of equilibrium an average of static observables is enough to measure transport coefficients. The flattening of the velocity structure factors at equilibrium clearly results in our formulas from the vanishing of the velocity correlations amplitude, which is proportional to . In that case, access to transport coefficients is only possible through the study of the dynamical structure factors.
Finally, we study the behaviour of such coherence lengths at different packing fractions. We observe a significantly growth of the relative extent of correlations and with the packing fraction, where denotes the mean free path. This is in agreement with the slowing down of the dynamics also observed in dense granular fluids .
The paper is organized as follows. In Section 2 we discuss the model, the corresponding hydrodynamic equations and the specific forms of noise used. The dependence of the eigenvalues of the linearized hydrodynamic matrix is also studied in order to check the stability and bounds of the linear approximation. In Section 3 we present a complete study of the static and dynamical structure factors for the velocity and density modes, comparing analytical predictions and numerical results. In section 4 a discussion about the main findings is presented, with conclusions and perspectives for future work. In A and B are reported, respectively, formulas for the transport coefficients and some details on the noise terms.
2 Microscopic model and fluctuating hydrodynamics
We consider a system of inelastic hard spheres in dimensions with mass and diameter . Particles are contained in a volume , with the linear size of the system. We denote by the number density and by the occupied volume fraction, (in two dimensions ). The particles undergo binary instantaneous inelastic collisions when coming at contact, with the following rule
where () and () are the post and pre-collisional velocities of particle (particle ), respectively; is the restitution coefficient (in the elastic case ), and is the unit vector joining the centers of the colliding particles.
In order to maintain a stationary fluidized state, an external energy source is coupled to every particle in the form of a thermal bath . In particular, the motion of a particle with velocity is described by the following stochastic equation
Here is a drag coefficient (which defines the characteristic interaction time with the external bath, ), is a white noise with and (Greek indexes denote Cartesian coordinates), while represents the action of particle-particle inelastic collisions. The effect of the external energy source balances the energy lost in the collisions so that an out-of-equilibrium stationary state is attained . In particular, let us stress that such energy injection mechanism acts homogeneously across the whole system, differently from other mechanisms where the energy is directly supplied only to a part of the system, as for instance for fluids under shear, or for systems in conctact with vibrating walls.
The stationary state is characterized by two time scales and two energy scales: the time scales are and the mean free time between collisions ; the energy scales are the temperature of the thermostat and the granular temperature (equal sign holds only if or if ). In particular, in the dilute limit, the granular temperature satisfies the following equation in the non-equilibrium stationary state 
where and is the pair correlation function at contact. The model has a well-defined elastic limit , where the fluid equilibrates to the bath temperature, . The viscous drag term in Eq. (2) models the interaction between each particle and the thermostat. It is important to observe that is not related to the transport coefficients of the granular fluid and is fixed as a model parameter. As mentioned above, introduces a time scale in the system that rules the tendency of particle to relax toward equilibrium at temperature . The characteristic time of collisions, , in all our simulations will be kept much smaller than : for this reason is considered the microscopic time-scale of our system since it dictates the smallest scale of relaxation toward the non-equilibrium stationary state. In particular, in the coarse-grained hydrodynamic description to be introduced below, we will take care of comparing the characteristic decay time of different hydrodynamics modes with , to verify the presence of a sufficient separation of scales.
In the following we will present a thorough numerical analysis of model (2), using an event-driven molecular dynamics algorithm . In particular, we will consider periodic boundary conditions in dimensions. The fixed parameters of the simulations are , , and . The packing fraction is varied by changing the seize of the box, and we consider systems with . The simulation data on static structure factors are obtained for a system of particles, averaged over about 100 realizations, whereas the results on dynamical correlators are obtained with samples of particles, averaged over about 4000 realizations.
2.1 Linearized hydrodynamics
Because we are interested in the behavior of large-scale spatial correlations in our system, we introduce here the coarse-grained hydrodynamic fields and as follows:
In the above equations and are respectively the heat flux and the pressure tensor, see details in B, and . In the velocity equation the viscous drag term has been inserted, while in the temperature equation three terms have been added: the sink term , where is the collision frequency, takes into account the energy dissipated by inelastic collisions, while the terms represent the energy exchanged with the thermostat.
Eqs. (5) give a fair description of the mesoscopic degrees of freedom of a granular fluid as long as a proper separation of space and time scales is verified between those degrees of freedom and all the microscopic ones which are projected out. This condition is, of course, not always satisfied [26, 27, 28], but is not prevented in principle and is, indeed, realized in many experiments or simulations [29, 30, 31, 20, 32].
Eqs. (5) can be linearized around the stationary homogeneous state, where the hydrodynamic fields take the values , and . A system of linear differential equations for the fluctuations , with , can be considered, with the Fourier transform defined as
and with and respectively the shear and longitudinal modes, namely
being a unitary vector such that . The system in Eq. (5) in Fourier space becomes
with the dynamical matrix
where , is the thermal diffusion coefficient ( is the heat conductivity), while and are the kinematic and longitudinal viscosity respectively. Formulas for all parameters and transport coefficients are given in A. There we refer to the Enskog theory for dense elastic hard spheres (EHS) , which provides a good approximation, as observed in [18, 34]. The following sections are devoted to show how the viscosities and can be obtained as fit parameters of static and dynamical correlations. Such results will be compared with the dense EHS predictions, finding good agreement.
2.2 Spectrum of the hydrodynamic matrix and separation of time-scales
We analyze the eigenvalues of in order to study the linear stability of the model and to characterize the range of validity of time-scales separation required by hydrodynamics. From the expression in Eq. (9) we learn that the shear modes are decoupled from all the others, and the typical time-scale of their decay simply reads as . To obtain the typical time-scales for the fluctuations of the other hydrodynamic fields we solve, numerically, the equation , where is the identity matrix. The eigenvalues spectrum thus found is shown in Fig. 1, for parameters of the model and . All parameters are calculated according to the formulas reported in A.
The study of hydrodynamic eigenvalues shows that, even in the case of a quite inelastic and dense regime, a range of scales where mesoscopic relaxation times are larger than the microscopic times exists and a hydrodynamic description can be attempted. In particular, comparing the eigenvalues with the Enskog frequency (which will be verified, below, to be a good estimate of the real frequency in simulations), see Fig. 1, we find, for each eigenvalue, the interval of values of where is fulfilled. The subscript indicates eigenvalues related to heat (), sound () and shear () modes. The range is narrower for the case of the mode dominated by temperature fluctuations (denominated “heat mode” in the Figure, red thin curve), while it is larger for shear modes (green dot-dashed lines) and modes dominated by longitudinal velocities and density fluctuations (here referred to as “sound modes”, black and blue thick curves, continuous for real part and dashed for imaginary part).
A general observation is that eigenvalues never have a positive real part, i.e. no instabilities are found, thanks to the presence of the external bath. Moreover, for the whole range of studied, the two “sound modes” - as usual - are complex conjugate, i.e. they propagate with a -dependent velocity. A negative real part is always present, determining overall damping. At small values of one can identify a sound velocity by observing a linear relation ; instead the dispersion relation becomes strongly nonlinear at large . Surprisingly, the study of sound eigenvalues shows that they have bifurcations at very small and very large wave-numbers, becoming in both limits pure real numbers, i.e. losing their propagating behavior. Those bifurcations are due to the external damping, ruled by : indeed in the limit (keeping finite ), they disappear and the spectrum studied in [18, 20] is retrieved. Moreover, we find that the wavevector where the first bifurcation occurs, moves towards smaller values of when the dissipation is increased, namely when is decreased (at fixed ) or when is increased (at fixed ). Looking closely to the bifurcations, we see that one of the eigenvalues for “sound modes” approaches for (representing total number conservation) and the other tends to (as the shear one), see inset of Fig. 1. The bifurcation at large is perhaps non-physical, as it always falls out of the hydrodynamic range. For the eigenvalue of the heat mode we find that . In the numerical setup used below, we have , so that the purely exponential decay of sound modes with large waves is never observed.
2.3 Stochastic description with fluctuating hydrodynamics
In order to fully account for the spatial fluctuations of the hydrodynamic fields and for the decay in time of such fluctuations we must add some noise terms to the linearized hydrodynamic equations: the basic assumption under fluctuating hydrodynamics is the same as for average (deterministic) hydrodynamics, i.e. a good separation of scales between hydrodynamic fields and microscopic degrees of freedom. In the linearized hydrodynamic equations the small scale fluctuations have been projected out, but their feedback on large scale fluctuations can be recovered by a proper addition of noise terms to dynamical equations:
A derivation from first principles of the noise is beyond our scope. A kinetic theory with fluctuations has been recently proposed in [35, 36], for the homogeneous cooling regime, which is very different from our case. A similar treatment has been realized, only for the shear mode, in a driven case (random kicks without damping, i.e. ) , showing that noise can be safely assumed to be white, at difference with the cooling regime. In such a case, with an additional but reasonable assumption on the two-particle velocity autocorrelation functions (namely, that such functions have only components in the hydrodynamic subspace), it is found that the fluctuation-dissipation relation for the internal part of the noise is satisfied, as already assumed in [25, 18]. Following those previous studies, we will consider valid such an assumption. Notice that it does not imply that fluctuation-dissipation relations will be satisfied by the whole noise, which is composed by internal as well as external contributions. In summary we write
where the two sources of noise for the hydrodynamic fields fluctuations are put in evidence: the first is the external contribution coming directly from the thermal bath, namely the stochastic force of Eq. (2); the second is internal and enters through the constitutive equations for the heat flux and the pressure tensor. A detailed discussion on noises is presented in B. The external and internal noises are Gaussian with zero average. The variances of external noises can be obtained directly from Eqs. (2) and (4)
while the variances of the internal contributions are obtained by imposing the fluctuation-dissipation theorem (see B for details):
The internal and external noises are uncorrelated
The hydrodynamic analysis of model (2) consists, then, in solving the system of coupled linear Langevin equations (10). In particular we are interested in finding the explicit forms of the static and dynamical structure factors, respectively
3 Out-of-equilibrium correlations: Static and Dynamical structure factors
It is well known that spatially extended correlations develop in the non-equilibrium stationary state of a driven granular fluid [18, 20]. In particular, the velocity correlator , where the average is taken over noises, can be written as
with the cross terms . The two terms on the right of Eq. (18) can be studied separately.
3.1 Shear modes
According to the matrix (9) the shear modes are decoupled from the others in the linear approximation and their dynamics obeys a simple Langevin equation
where, from now on, we put in all formulas. The effect of internal and external noises is equal to a single complex noise with variance:
From this, the equal-time correlator can be easily calculated:
with . The meaning of the above equation is clear: the small length-scale physics, which depends on the inelastic collisions, is related to the granular temperature , while at large distances there are correlations with finite amplitude and extent . Notice that in the limit , keeping finite, from Eq. (21) one obtains the result of Ref. , where and long range correlations are observed. A similar power law behaviour is also observed in molecular fluids under shear . The different result we obtain in our model is due to the intrinsic cut-off introduced by the viscous drag . Indeed, the finite extent of correlations is even more clear when Eq. (21) is written in real space, yielding the spatial correlation function , which reads:
where is the 2nd kind modified Bessel function that, for large distances, decays exponentially
At equilibrium, namely when collisions are elastic and , equipartition between modes is perfectly fulfilled and the structure factor becomes flat, i.e. . Differently, in the granular case, where , equipartition breaks down and from Eq. (21) we have that for small and for large . We see here that out of equilibrium the quantity measures the range of static correlations of the vorticity field. The behaviour described above is in good agreement with experimental results obtained for driven granular fluids, as reported in  and, more recently, in .
From Eq. (19) we also find that fluctuations decay exponentially
with a characteristic time . Such a behavior is also observed for elastic fluids, with the only difference that in that case is constant. The length-scale can be therefore always connected to dynamical properties of the system. What is peculiar of the out-of-equilibrium regime is that the so defined also represents the extent of correlations of the vorticity field, thus establishing a remarkable link between static correlations and dynamical ones.
In Fig. 2 are reported the correlators for different values of and the same packing fraction measured in numerical simulations. Let us notice that the decay of is always exponential: we can therefore a-posteriori support the validity of linear hydrodynamics, from which the Langevin equations for shear modes is obtained. By interpolating with a parabola the characteristic time as a function of , see inset of Fig. (2), we obtain the shear viscosity . Let us stress the deep connection between statics and dynamics in the out-of-equilibrium regime: inserting the values of obtained from the dynamics into Eq. (21), we find curves that well superimpose the numerical data for , (see Fig. 3). The values of obtained from the dynamics of shear modes can be independently obtained as fit parameters of via Eq. (21). The values of obtained with the two different procedures are compatible, as can be seen from Tab. 1. They are also reasonably close to the dense EHS predictions, presented in the same Table.
3.2 Longitudinal modes
3.2.1 Static correlations
The same considerations discussed above for shear modes also hold for the other hydrodynamic modes, which are coupled each other. In order to study their behaviour we have to take into account all the elements of the dynamical matrix. In particular, the matrix of static structure factors with elements , is obtained solving the following linear system:
where the matrix of noises is such that:
where denotes a diagonal matrix with elements .
The expression of the longitudinal structure factor,
turns out to be the ratio between two even polynomial functions of the th order in :
In the above expression eight constants have been introduced, which depend in a complicate manner by all the parameters of the system. Let us focus here on the asymptotic behavior of at large and small values of . From Eq. (3) there follows the relation , which allows us to recast the series expansion around of the expression in Eq. (29) in the form
Up to the expression in Eq. (30) is equivalent to
namely a form analogous to the structure factor we found for the shear modes. In Fig. 4 we show the static structure factors for different packing fractions. Again, notice that in the limit with finite, the behaviour found in  is recovered. Albeit the above expression is in principle only valid for low values, it captures also the large limit, when fine oscillations are disregarded. Indeed, expanding Eq. (29) for large values we find:
Such a discussion shows that even for longitudinal modes the viscosity is related to a finite correlation length, measurable from static velocity correlations when the system is out of equilibrium with . The behavior of that length when the packing fraction is increased will be discussed in the last section.
3.2.2 Dynamical correlations
Dynamical correlations for longitudinal modes are less simple than those for the shear mode, since they are given by a superposition of different (real and imaginary) exponentials. The dynamical structure factors are obtained by solving the equation of motion (10), which, in the frequency domain, reads as
with the identity matrix,
The dynamical structure factors and take the explicit forms
|Sim||dense EHS||Sim||dense EHS|
|Sim||dense EHS||Sim||dense EHS|
|Sim||dense EHS||Sim||dense EHS|
All the coefficients appearing in these expressions can be evaluated within the dense EHS approximation, so that Eq. (41) can be used in order to obtain and from the numerical data. Clearly the writing of physical quantities, such as or itself, with dense EHS formulas leads to systematic errors. The amplitude of this error can be estimated for instance by comparing the theoretical prediction of dense EHS collision frequency with the collision frequency measured in simulations. In Table 2 are listed both dense EHS and numerical collision frequencies at different packing fractions. Again it is found that the dense EHS prediction is quite good.
In order to obtain the longitudinal viscosity and the thermal diffusion coefficient, we fit our numerical data for the longitudinal modes using Eq. (41), where all parameters but and are fixed to the dense EHS values and is the one measured in simulations. In Fig. 5 is shown for different values of and fixed , together with the best fit curves. The values of and so obtained, together with those computed within the dense EHS approximation, for different values of and , are reported in Table 3 and within errors are found independent of . This fact represents an a-posteriori check that we are in the regime of validity of linearized hydrodynamics. Indeed, in Eqs. (41) and (3.2.2) appears as -independent variable. Only at low packing fraction we observe a dependence on . This is perhaps due to diluteness, which implies too large mean free path or mean free time with respect to mesoscopic scales.
In Fig. 6 we also show the time decay of the dynamical structure factor. It can be appreciated, in the time domain, the superposition of different real and imaginary exponentials, which determines a mix of damping and propagation.
The dynamical structure factor at fixed and different packing fractions is reported in Fig. 7: it is remarkable the observation of a time-scale, individuated by the peak frequency of , which increases as the packing fraction is increased. Such behaviour is consistent with the observation, discussed in details below, of a growth of the correlation lengths defined above, together with the packing fraction.
4 Summary and conclusions: transport coefficients and non-equilibrium correlation lengths
In conclusion, we have studied both static and dynamical correlations for hydrodynamic fluctuations of the velocity and density fields: indeed we recall Eq. (41) which gives a direct relation between density and longitudinal velocity structure factors. A main comment concerns the good success of analytical predictions: comparison with simulations shows a fair agreement up to , with values for most of the parameters directly given in the dense EHS approximation. This signals a success of the three main ingredients: 1) dense EHS theory for transport coefficients, working even at quite high densities, 2) scale separation and granular hydrodynamics, and 3) prescription of Eqs. (12-14) for the hydrodynamic noise. The last ingredient is perhaps the most interesting, if one considers how hard is the characterization in simple terms of non-equilibrium systems. Moreover, the assumption of white internal noise is less obvious than in the case of simple random kicks without drag : in that case, the absence of drag () implies a relaxation time at large scales, , which diverges, making a solid base for assuming fast the relaxation of noise due to microscopic degrees of freedom. In our model, in principle, the drag could be large enough to make even large scales fast, making difficult to define hydrodynamic fields. Recent experiments show that, even if that could be a realistic situation, this does not change dramatically the qualitative behavior of structure factors .
One of the main results of our study is the presence of spatial order in the form of non-equilibrium velocity correlations. This should be related to the slowing-down of the dynamics with increasing packing fraction in granular systems [39, 20], and with the existence of a time-scale growing with the density, as observed in . As discussed in Sec. 3, in our model the competition between different relaxation mechanisms given by the kinematic and longitudinal viscosities and and the thermostat damping give place to a couple of length-scales characterizing non-equilibrium structure factors: and . The first trivial observation is that allowing such lengths diverge. It means that, according to the prediction of , the largest-scale correlations are always equal to the size of the system, which is the largest size available. In our model there is a cut-off on such correlations imposed by the viscous drag due to the interaction with the thermal bath. The existence of such a cut-off, which represents a fixed parameter at different packing fractions, allows us to show that by increasing the packing fraction the extent of correlations is effectively increased. While the absolute value of is only slightly enhanced when is increased, as also observed in experiments in , we stress that, in order to appreciate the physical meaning of such length at different densities, one has to compare it with the microscopic relevant spatial scale in the system, which is given by the mean free path of the particles and which also changes with the packing fraction. Then one finds that is remarkably increased at high densities, as can be seen in Fig. 8, and also in recent experiments . We may summarize the observed phenomenon saying that the higher is the packing fraction the higher must be the number of intermediate scattering events between two different particles in order to decorrelate their velocities. This scenario could be reflected in transitions of dynamical origin, at higher packing fractions [40, 41, 42, 39].
We thank U. M. B. Marconi and A. Vulpiani for helpful suggestions. The work of the authors is supported by the “Granular-Chaos” project, funded by the Italian MIUR under the FIRB-IDEAS grant number RBID08Z9JE.
Appendix A Dense EHS formulas for hydrodynamic coefficients
The definitions of the parameters entering the matrix are
In the dense EHS approximation the pressure can be written as
which in reads as
where has been made use of the relation between packing fraction and density , and of the definition with . Notice also that we have taken into account the correction due to the inelasticity, as given in . There are dense EHS formulas also for the dependence of the diffusion coefficients on the packing fraction and the granular temperature. In the following such formulas are written for a 2d system. In this case, we use the Verlet-Levesque approximation for the pair correlation function at contact: .
Appendix B Noises
In this appendix we present a detailed discussion of the noise terms appearing in the fluctuating hydrodynamic equations (10). Let us start from the external noises, which can be simply obtained from Eqs. (2) and (4). We find that the noise contributions to the equations for the velocity and temperature fields are, respectively
These are Gaussian noises with variances