Fluctuating hydrodynamics and mesoscopic effects of spatial correlations in dissipative systems with conserved momentum

Since the seminal paper of Einstein [1], it has been well known that the fluctuating behaviour of systems at the mesoscopic level reflects the hectic microscopic dynamics beneath. While the equilibrium behavior of mesoscopic fluctuations has been investigated and understood in detail [2, 3], a big effort is still being carried out to explore the fluctuating properties of non-equilibrium media [4]. These are known to lead, in great generality, to the emergence of spatial correlations and pattern self-organization [5, 6]. In this context, the crucial task of connecting microscopic and mesoscopic dynamics is considerably simplified when there exists a separation of scales, which makes it possible to introduce slow fields evolving under the so-called hydrodynamic equations [7].

An important class of systems exhibiting patterns includes two types of complex fluids: active matter [8, 9], such as bacteria or birds, and fluidized granular materials [10]. Interestingly, active and granular matter are often associated [11–14]. They are not only relevant for applied and biomedical sciences, but also offer fascinating challenges for kinetic theory [15]. Indeed, the lack of energy conservation in the microscopic dynamics makes them intrinsically out-of-equilibrium systems [16]. In granular and active fluids, the spectacular emergence of spatial patterns, particularly in vectorial fields such as momentum or orientation, is often understood in terms of hydrodynamic equations [17]. Furthermore, a relevant role is played by fluctuations, as an inevitable consequence of the relatively small number of their elementary constituents [18].

One of the most intensively investigated states in the realm of granular fluids is the homogeneous cooling state (HCS) [19, 20]. Therein, the granular temperature decays in time following Haff's law [21], whereas the system remains spatially homogeneous. Remarkably, the HCS is the reference state for the hydrodynamic description of granular fluids but it is unstable: for large enough inelasticity or system size, the scaled fluctuations of the transverse velocity increase (shear instability) and eventually density inhomogeneities arise (clustering instability) [22]. This instability for large system sizes makes it relevant to look into the finite size corrections to the physical quantities, like the cooling rate, since a 'thermodynamic limit' in which the system size is infinitely large cannot be taken without simultaneously making the inelasticity vanish. Notwithstanding, and to the best of our knowledge, these finite size corrections have only been investigated for a system of smooth inelastic hard spheres very close to the shear instability [23].

In deriving mesoscopic transport equations from microscopic rules, analytical results are needed and simple models are good candidates for this [24, 25]. In this paper, we study a 1D lattice model which implements two main ingredients of granular fluids: inelastic collisions and momentum conservation. Given the simplicity and appealing physical picture of the model, this novel approach may help to improve our current understanding of the complex behaviour of granular fluids. To start with, we recover some of the main features of granular fluids: more specifically, the large size and inelasticity shear instability. In addition, we are able to compute exactly the shape of velocity correlations, which allows us to, first, extend the known results of fluctuating hydrodynamics [26] and, second, obtain the finite size corrections to the cooling rate. The latter has a non-trivial dependence on the inelasticity and system size: for a given inelasticity, it changes sign at a certain value of the system size that is smaller than the one corresponding to the shear instability.

Finally, we would like to stress that momentum conservation is a physical constraint that certainly has a recognized role in the appearance of long-range spatial correlations [6, 27]. However, since it complicates the description and the derivation of exact results, it is rarely considered in its entirety. Here, starting from the microscopic dynamics, we are able to rigorously derive the mesoscopic equations that describe the average and fluctuating behaviour at the hydrodynamic scale, taking into account momentum conservation in full.

Fluctuating hydrodinamics in linear and nonlinear lattice diffusive models have been extensively studied in recent years, both in the conservative [28–31] and in the dissipative cases [32–34]. Inspired by [35], we consider a 1D lattice with N sites and given boundary conditions, either periodic or thermostatted, depending on the situation of interest. At a given time p, each site l possesses a velocity ${v}_{l,p}$ and the total energy of the system is ${E}_{p}={\displaystyle \sum }_{l=1}^{N}{v}_{l,p}^{2}$. In an elementary time step of the dynamics, with a probability discussed below, a pair of nearest neighbours $(l,l+1)$ collides inelastically and evolves following the rule ($0\lt \alpha \leqslant 1$)

$$\begin{eqnarray}&&{v}_{l,p+1}={v}_{l,p}-(1+\alpha ){{\rm{\Delta }}}_{l,p}/2,\qquad {v}_{l+1,p+1}={v}_{l+1,p}+(1+\alpha ){{\rm{\Delta }}}_{l,p}/2,\end{eqnarray} \tag{ 1 }$$

having defined the relative velocity

$$\begin{eqnarray}&&{{\rm{\Delta }}}_{l,p}={v}_{l,p}-{v}_{l+1,p}.\end{eqnarray} \tag{ 2 }$$

Momentum is always conserved, ${v}_{l,p}+{v}_{l+1,p}={v}_{l,p+1}+{v}_{l+1,p+1}$, while energy, if $\alpha \ne 1$, is not: ${v}_{l,p+1}^{2}+{v}_{l+1,p+1}^{2}-{v}_{l,p}^{2}-{v}_{l+1,p}^{2}=({\alpha }^{2}-1){{\rm{\Delta }}}_{l,p}^{2}/2\lt 0$.

The definition of the model implies that there is no mass transport, particles are at fixed positions and they only exchange momentum and kinetic energy. We are also disregarding the so-called kinematic constraint in [35], namely a colliding pair is chosen independently of the sign of its relative velocity. This can be understood as the velocity of the particles representing not their motion along the lattice axis but rather along a transversal one: in fact, the hydrodynamics derived here replicates transport equations for granular gases in $d\gt 1$ restricted to the shear (transverse) velocity field; see figure 1.

Zoom In Zoom Out Reset image size

In the context of granular fluids, the model may be physically motivated as follows. We start from a $d\gt 1$ system that has been divided into 'slabs' that are perpendicular to the lattice direction. Specifically, each particle on the lattice represents one slab. In this sense, the parameter α that appears in the collision rule (1) should not be confused with the usual restitution coefficient defined in granular media, since here α stands for an effective inelasticity for the collisions between slabs. The connection with a 'real' granular fluid should be done at the level of the cooling rate that appears in the hydrodynamic equations; see section 3.

Now we write down the evolution equation for the velocities. At time p, the probability that the nearest-neighbours pair at sites $(l,l+1)$ collide is assumed to be ${P}_{l,p}\propto | {{\rm{\Delta }}}_{l,p}{| }^{\beta }$. Then,

$$\begin{eqnarray}&&{v}_{l,p+1}-{v}_{l,p}=-{j}_{l,p}+{j}_{l-1,p},\qquad {j}_{l,p}=(1+\alpha ){{\rm{\Delta }}}_{l,p}{\delta }_{{y}_{p},l}/2,\end{eqnarray} \tag{ 3 }$$

which is a discrete continuity equation for the (conserved) momentum, with ${j}_{l,p}$ being the momentum current, that is, the flux of momentum from site l to site $l+1$ at the time step p. Therein, ${\delta }_{{y}_{p},l}$ is Kronecker's δ and y_p is a random integer which selects the colliding pair with probability ${P}_{l,p}$. The evolution equation for the energy is obtained by squaring equation (3),

$$\begin{eqnarray}&&{v}_{l,p+1}^{2}-{v}_{l,p}^{2}=-{J}_{l,p}+{J}_{l-1,p}+{d}_{l,p}.\end{eqnarray} \tag{ 4 }$$

Again, we have defined an energy current ${J}_{l,p}$ from site l to site $l+1$ and the energy dissipation ${d}_{l,p}$ at site l as

$$\begin{eqnarray}&&{J}_{l,p}=({v}_{l,p}+{v}_{l+1,p}){j}_{l,p},\qquad {d}_{l,p}=\displaystyle \frac{{\alpha }^{2}-1}{4}({\delta }_{{y}_{p},l}{{\rm{\Delta }}}_{l,p}^{2}+{\delta }_{{y}_{p},l-1}{{\rm{\Delta }}}_{l-1,p}^{2})\lt 0.\end{eqnarray} \tag{ 5 }$$

The sink term ${d}_{l,p}$ only vanishes in the elastic case $\alpha =1$.

Stochastic simulations of the model are done as follows. At any Montecarlo step p, one site l is picked with probability ${P}_{l,p}\propto {| {{\rm{\Delta }}}_{l,p}| }^{\beta }$ and particles l and $l+1$ collide following the microscopic rules in equation (1). The simplest choice for ${P}_{l,p}$ corresponds to $\beta =0$; all the pairs are chosen with uniform probability, ${P}_{l,p}=1/L$, in which L is the number of pairs. This is often called in the literature the model of inelastic Maxwell molecules (MM) [35]. Note that L is basically equal to N but depends on the boundary conditions: for periodic boundary conditions, it is L = N, but if we consider the system coupled to two extra sites 0 and $N+1$, which introduce the appropriate boundary conditions, it is $L=N+1$. The periodic boundary conditions, sketched in figure 1, correspond to the free (undriven) evolution of the system and if l = N it is the pair $(N,1)$ that collides.

In the following, we discuss the hydrodynamic limit, fluctuations and correlations for a particular choice of ${P}_{l,p}$. Specifically, we consider the case of MM: such a choice is dictated by the will of simplifying the presentation and making clear the essential points. We postpone a more complete and general discussion to a more technical and detailed paper, in preparation. The theoretical results are compared to the numerical simulations described above. A large enough value of L, which is indicated in the figures, has been used to ensure the hydrodynamic limit, and we have averaged over 10⁵ realizations of the stochastic dynamics. Aside from MM, results for HS ($\beta =1$) are also shown in a few clearly marked cases.

Let us define, as usual, the following local averages over initial conditions and noise realizations: ${u}_{l,p}\equiv \langle {v}_{l,p}\rangle$, ${E}_{l,p}\equiv \langle {v}_{l,p}^{2}\rangle$ and ${T}_{l,p}\equiv {E}_{l,p}-{u}_{l,p}^{2}$. Their evolution is obtained under a series of assumptions. With the choice of MM, y_p is a uniform distributed random integer, namely $\langle {\delta }_{{y}_{p},l}\rangle =1/L$, with L being the number of nearest neighbour pairs. In addition, when considering the average dissipation at site l, there appear moments like $\langle {v}_{l,p}{v}_{l\pm 1,p}\rangle$. To the lowest order, we assume that neighbouring velocities are uncorrelated, that is, $\langle {v}_{l,p}{v}_{l\pm 1,p}\rangle ={u}_{l,p}{u}_{l\pm 1,p}$ (see appendix A).

Now we assume that ${u}_{l,p}$ and ${E}_{l,p}$ are smooth functions of space and time and introduce a continuum ('hydrodynamic') limit (HL) by defining macroscopic scales: ${\rm{\Delta }}x\sim {L}^{-1}$ and ${\rm{\Delta }}t\sim {L}^{-3}$. Each spatial derivative introduces thus a factor ${L}^{-1}$ in the continuum limit: therefore, the difference between the current terms in the balance equations is of the order of ${L}^{-3}$. On the other hand, the dissipation goes as $(1-{\alpha }^{2}){L}^{-1}$, which makes it useful to define the cooling rate as (see appendix A)

$$\begin{eqnarray}&&\nu =(1-{\alpha }^{2}){L}^{2}.\end{eqnarray} \tag{ 6 }$$

It is natural, on the scales defined by the HL, to define the mesoscopic fields $u(x,t)$, $E(x,t)$ and $T(x,t)$, as well as the average mesoscopic currents $j(x,t)={\mathrm{lim}}_{L\to \infty }$ ${L}^{2}\langle {j}_{l,p}\rangle =-{\partial }_{x}u(x,t)$, $J(x,t)={\mathrm{lim}}_{L\to \infty }$ ${L}^{2}\langle {J}_{l,p}\rangle =-{\partial }_{x}$ $[{u}^{2}(x,t)+T(x,t)]$ and the average mesoscopic dissipation of energy $d(x,t)={\mathrm{lim}}_{L\to \infty }\;{L}^{3}\langle {d}_{l,p}\rangle =-\nu T(x,t)$, which depends only on the granular temperature but not on the average local velocity $u(x,t)$, as expected on a physical basis. After computations detailed in appendix A, we get the HL of equations (3) and (4),which are ${\partial }_{t}u(x,t)=-{\partial }_{x}j(x,t)$ and ${\partial }_{t}E(x,t)=-{\partial }_{x}J(x,t)-\nu T(x,t)$ respectively. Then, we can write the average hydrodynamic evolution equations

$$\begin{eqnarray}&&{\partial }_{t}u={\partial }_{{xx}}u,\qquad {\partial }_{t}T=-\nu T+{\partial }_{{xx}}T+2{({\partial }_{x}u)}^{2}\end{eqnarray} \tag{ 7 }$$

over the length and time scales defined above. Note that, here, we have substituted $1+\alpha$ by 2, because ${\alpha }^{2}=1-\nu {L}^{-2}$, and we have already neglected ${L}^{-1}$ terms. These equations must be supplemented by appropriate boundary conditions for the situation of interest. Technical details are deferred to a later paper.

Let us consider the fluctuations of the microscopic currents and dissipation, that is, ${j}_{l,p}={\tilde{j}}_{l,p}+{\xi }_{l,p}$, ${J}_{l,p}={\tilde{J}}_{l,p}+{\eta }_{l,p}$, and ${d}_{l,p}={\tilde{d}}_{l,p}+{\theta }_{l,p}$. Tilde variables correspond to a partial average: they are averaged over the fast variables ${y}_{l,p}$ but not over the slow ones ${v}_{l,p}$. Thus, for example, ${\tilde{j}}_{l,p}=(1+\alpha ){{\rm{\Delta }}}_{l,p}/2L$. It is clear that this choice guarantees that all noises ${\xi }_{l,p}$, ${\eta }_{l,p}$ and ${\theta }_{l,p}$ have zero average. The noise correlations read $\langle \xi \xi ^{\prime} \rangle \sim {A}_{\xi }\delta (x-x^{\prime} )\delta (t-t^{\prime} )$, $\langle \eta \eta ^{\prime} \rangle \sim {A}_{\eta }\delta (x-x^{\prime} )\delta (t-t^{\prime} )$ and $\langle \theta \theta ^{\prime} \rangle \sim {A}_{\theta }\delta (x-x^{\prime} )\delta (t-t^{\prime} )$ with amplitudes ${A}_{\xi }=2{L}^{-1}T(x,t)$, ${A}_{\eta }=4{L}^{-1}T(x,t)[T(x,t)+2{u}^{2}(x,t)]$ and ${A}_{\theta }=3{L}^{-3}{\nu }^{2}\;{T}^{2}(x,t)$ (see appendix B). In the above relations, we have used the notation $\xi \equiv \xi (x,t)$ and $\xi ^{\prime} \equiv \xi (x^{\prime} ,t^{\prime} )$, and similar notations for $\eta ,{\eta }^{\prime },\theta ,\theta ^{\prime}$. Thus, the current noises are delta correlated in space and time, and their amplitudes scale as ${L}^{-1}$ with the system size L. On the other hand, the noise of the dissipation is subdominant with respect to the moment and energy currents, its amplitude being proportional to ${L}^{-3}$, and therefore it is usually negligible. Gaussianity for these noises can be easily demonstrated; see [32]. Interestingly, being in the presence of two fluctuating fields, correlations between different noises appear. Theoretical predictions for noise correlations, amplitudes and Gaussianity have been successfully tested in both MM and HS simulations; see appendices B and C.

Here we focus our attention on the case of spatial periodic boundary conditions and an initial 'thermal condition': ${v}_{l,0}$ is a random Gaussian variable with zero average and unit variance, that is, ${T}_{l,0}\equiv T(x,0)=1$. Starting from this condition, the system is expected to typically fall into the so-called homogeneous cooling state (HCS), in which the velocity and temperature fields remain spatially uniform, and the temperature decays in time. Indeed, the spatially homogeneous solution of the average hydrodynamic equation (7) reads

$$\begin{eqnarray}&&u(x,t)=0,\qquad {T}_{\mathrm{HCS}}(x,t)=T(t=0){{\rm{e}}}^{-\nu t}.\end{eqnarray} \tag{ 8 }$$

The exponential decrease of the granular temperature is typical of MM, where the collision frequency is velocity-independent. It replaces the so-called Haff's law which was originally derived in the HS case, where ${T}_{\mathrm{HCS}}\sim {t}^{-2}$ because $\dot{T}\propto {T}^{3/2}$ [21].

The HCS is known to be unstable: it breaks down in too large or too inelastic systems [36]. In our model and in the hydrodynamic limit, this condition is expected to be replaced by a condition of large ν. The stability is studied by introducing rescaled fields $U(x,t)=u(x,t)/\sqrt{{T}_{\mathrm{HCS}}(t)}$ and $\tilde{T}=T(x,t)/{T}_{\mathrm{HCS}}(t)$ and by linearizing the hydrodinamic equations near the HCS, i.e. $T(x,t)={T}_{\mathrm{HCS}}(t)+\delta T(x,t)$ and $U(x,t)=\delta U(x,t)$. The analysis of linear equations becomes straightforward by space Fourier transforming, which gives

$$\begin{eqnarray}&&{\partial }_{t}\delta U(k,t)=\displaystyle \frac{\nu -2{k}^{2}}{2}\delta U(k,t),\qquad {\partial }_{t}\delta \tilde{T}(k,t)=-{k}^{2}\delta \tilde{T}(k,t).\end{eqnarray} \tag{ 9 }$$

Therefore, $\delta U$ is unstable for wave numbers that verify $\nu -2{k}^{2}\gt 0$. In the continuous variables we are using, the system size is 1, so that the minimum available wavenumber is ${k}_{\mathrm{min}}=2\pi$. Thus, there is no unstable mode for ν (lengths) below a certain threshold ${\nu }_{c}$ (L_c), with

$$\begin{eqnarray}&&{\nu }_{c}=8{\pi }^{2},\qquad {L}_{c}=2\pi \sqrt{2}{(1-{\alpha }^{2})}^{-1/2}.\end{eqnarray} \tag{ 10 }$$

On the other hand, for $\nu \gt {\nu }_{c}$ ($L\gt {L}_{c}$), the HCS is unstable and modes with wave numbers verifying $k\lt \sqrt{\nu /2}$ increase with time. This instability mechanism is identical to the one found in granular gases for shear modes [26]. Theoretical predictions and simulations results perfectly agree, as plotted in figure 2. It is important to stress that the amplification appears in the rescaled velocity $U(x,t)$ and not in the velocity $u(x,t)$. The same result is found and compares well with simulations in the HS case.

Zoom In Zoom Out Reset image size

Assuming space translation invariance, which is certainly valid in the HCS, we can write a hierarchy of equations for the spatial correlations defined as ${C}_{k,p}=\langle {v}_{j,p}{v}_{j+k,p}\rangle$ at a distance of k sites, at time p:

$$\begin{eqnarray}&&{C}_{0,p+1}={C}_{0,p}+({\alpha }^{2}-1){L}^{-1}({C}_{0,p}-{C}_{1,p}),\end{eqnarray} \tag{ 11 }$$

$$\begin{eqnarray}&&{C}_{1,p+1}={C}_{1,p}+(1-{\alpha }^{2}){L}^{-1}({C}_{0,p}-{C}_{1,p})+(1+\alpha ){L}^{-1}({C}_{2,p}-{C}_{1,p}),\end{eqnarray} \tag{ 12 }$$

$$\begin{eqnarray}&&{C}_{k,p+1}={C}_{k,p}+(1+\alpha ){L}^{-1}({C}_{k+1,p}+{C}_{k-1,p}-2{C}_{k,p}),\quad 2\leqslant k\leqslant \displaystyle \frac{L-1}{2},\end{eqnarray} \tag{ 13 }$$

$$\begin{eqnarray}&&{C}_{\frac{L+1}{2},p}={C}_{\frac{L-1}{2},p}.\end{eqnarray} \tag{ 14 }$$

A striking consequence of momentum conservation is the sum rule ${C}_{0,p}+2{\displaystyle \sum }_{k=1}^{(L-1)/2}{C}_{k,p}=0$. Then we expect that correlations are of the order $O({L}^{-1})$. For example, in the elastic limit $\alpha =1$, their equilibrium value is $\langle {v}_{j}{v}_{j+l}\rangle =-T{(L-1)}^{-1}$, $\forall l\ne 0$. We take equations (11) and (13) in the continuum limit, $x=(k-1)/L$ and $t=p/{L}^{3}$, and retain only terms up to $O({L}^{-1})$, obtaining

$$\begin{eqnarray}&&\displaystyle \frac{{\rm{d}}T(t)}{{\rm{d}}t}=-\nu [T(t)-{L}^{-1}\psi (t)]+O({L}^{-3})\end{eqnarray} \tag{ 15 }$$

$$\begin{eqnarray}&&{\partial }_{t}D(x,t)=2{\partial }_{{xx}}D(x,t)+O({L}^{-2}).\end{eqnarray} \tag{ 16 }$$

Here, we have introduced the notations $D(x,t)={LC}(x,t)$ and $\psi (t)={\mathrm{lim}}_{x\to 0}D(x,t)$. Expression (15) introduces a correction in the hydrodynamic average granular temperature, given by the nearest-neighbour-particle velocity correlation, whereas equation (16) is a diffusion equation for spatial correlations. Boundary conditions stem from equations (12) and (14), which give a reflecting boundary at $x=1/2$ and the sum rule for momentum conservation, $T(t)+2{\displaystyle \int }_{0}^{1/2}D(x,t)=0$. In the long time limit, we obtain the following scaled stationary solution

$$\begin{eqnarray}&&\overset{\sim }{D}(x)=-A\mathrm{cos}[\pi \sqrt{\displaystyle \frac{\nu }{{\nu }_{c}}}(1-2x)],\quad A=\displaystyle \frac{\pi \sqrt{\displaystyle \frac{\nu }{{\nu }_{c}}}}{\mathrm{sin}(\pi \sqrt{\displaystyle \frac{\nu }{{\nu }_{c}}})},\end{eqnarray} \tag{ 17 }$$

where $\overset{\sim }{D}(x)={\mathrm{lim}}_{t\to \infty }\;D(x,t)/{T}_{\mathrm{HCS}}(t)$. Note that the Fourier transform of $\overset{\sim }{D}(x)+\delta (x)$ ⁴ takes the form $S(k)=\frac{{k}^{2}}{{k}^{2}-\nu /2}$ with $k=2\pi n$ and n is a positive integer, which has been derived for the correlations of the velocity shear mode in $d\gt 1$ from Landau-like granular fluctuating hydrodynamics; see for instance [26]. This result reinforces the motivation of our 1D model as a simple picture for the shear mode of the velocities in $d\gt 1$.

In figure 3, we compare the theoretical prediction in equation (17) for the MM case with numerical results. Remarkably, such prediction compares well also with HS simulations, where the analytical computation appears to be more challenging. In conclusion, the mechanism that induces spatial correlations in the system seems to be independent of the particular interaction model.

Zoom In Zoom Out Reset image size

Equation (15) suggests that the Haff law has a finite size correction. We consider a perturbation around the HCS $T(t)={T}_{\mathrm{HCS}}(t)+{L}^{-1}\delta T(t)+...$ and $D(x,t)={D}_{\mathrm{HCS}}+{L}^{-1}\delta D(x,t)+...$. Making use of equation (15) and defining $\delta \tilde{T}(t)=\delta T(t)/{T}_{\mathrm{HCS}}(t)$, we obtain $\frac{{\rm{d}}}{{\rm{d}}t}\delta \tilde{T}=\nu {\tilde{\psi }}_{\mathrm{HCS}}$, with

$$\begin{eqnarray}&&{\tilde{\psi }}_{\mathrm{HCS}}=\displaystyle \frac{{\psi }_{\mathrm{HCS}}(t)}{{T}_{\mathrm{HCS}}(t)}=-\pi \sqrt{\displaystyle \frac{\nu }{{\nu }_{c}}}\mathrm{cot}(\pi \sqrt{\displaystyle \frac{\nu }{{\nu }_{c}}}).\end{eqnarray} \tag{ 18 }$$

Hence, the granular temperature follows

$$\begin{eqnarray}&&T(t)={T}_{\mathrm{HCS}}(t)[1+\displaystyle \frac{1}{L}{\tilde{\psi }}_{\mathrm{HCS}}\nu t+{\mathcal{O}}({L}^{-2})].\end{eqnarray} \tag{ 19 }$$

which is valid for not very long times (t not scaling with L). There is a critical dissipation value ${\nu }_{\psi }={\nu }_{c}/4=2{\pi }^{2}$ where ${\tilde{\psi }}_{\mathrm{HCS}}$ changes sign, and this determines a change of the time-derivative of $\delta \tilde{T}$. Thus, at finite (large) values of L, the temperature decays faster or slower than the Haff law if $\nu \lt {\nu }_{\psi }$ or $\nu \gt {\nu }_{\psi }$, respectively. In figure 4, we compare the predicted Haff law finite size effect with the simulation results, obtaining excellent agreement.

Zoom In Zoom Out Reset image size

In conclusion, we have discussed the rigorous hydrodynamic limit of a lattice model for granular fluids with momentum conservation and energy dissipation. Macroscopic equations reproduce the realistic evolution of the velocity shear mode, which is diffusive, as well as that of the temperature field, which includes heat diffusion, inelastic dissipation and viscous heating. A crucial phenomenon of inelastic fluids, that is, the shear instability of the homogeneous cooling state, is recovered.

The model allows us to derive the evolution of $\sim 1/L$ spatial correlations, which present non-trivial long-range extension due to momentum conservation and alter the temperature decay in an observable way. This opens new interesting paths of investigation, such as trying to relate the deviation from the Haff law found here with the renormalization of the cooling rate found in systems of smooth HS near the shear instability [23]. The appearance of long-range correlations under non-equilibrium conditions for conserved fields is a feature well expected on general grounds [5, 6], but rarely derived in full analytical detail.

Finally, we stress the importance of considering finite size effects in granular systems, since large sizes are rarely realized in experiments. In addition, the HCS is unstable for large L, and therefore finite size corrections cannot be disregarded by considering an arbitrarily large system. These facts, together with the scarcity of studies about finite size effects in granular matter, makes it worthwhile to further investigate this point in the near future.

Antonio Prados acknowledges the support of the Spanish Ministerio de Economía y Competitividad through grant FIS2014-53808-P. We would like to thank C A Plata, D Villamaina and A Baldassarri for useful discussions.

Average fields (over initial conditions and noise realizations) are defined as ${u}_{l,p}\equiv \langle {v}_{l,p}\rangle$, ${E}_{l,p}\equiv \langle {v}_{l,p}^{2}\rangle$ and ${T}_{l,p}\equiv {E}_{l,p}-{u}_{l,p}^{2}$. The microscopic equations for the evolution of averages at time p at site l are obtained by averaging equations (3) and (4) of the main text, obtaining

$$\begin{eqnarray}&&{u}_{l,p+1}-{u}_{l,p}=-\langle {j}_{l,p}\rangle +\langle {j}_{l-1,p}\rangle ,\end{eqnarray} \tag{ 20 }$$

$$\begin{eqnarray}&&{E}_{l,p+1}-{E}_{l,p}=-\langle {J}_{l,p}\rangle +\langle {J}_{l-1,p}\rangle +\langle {d}_{l,p}\rangle .\end{eqnarray} \tag{ 21 }$$

Averages of currents and dissipation can be computed assuming the local equilibrium approximation (LEA), which is explicitly stated as

$$\begin{eqnarray}&&{P}_{2}({v}_{l},{v}_{l+1};p)=\displaystyle \frac{1}{\sqrt{2\pi {T}_{l,p}}\sqrt{2\pi {T}_{l+1,p}}}{{\rm{e}}}^{-\frac{{({v}_{l}-{u}_{l,p})}^{2}}{2{T}_{l,p}}-\displaystyle \frac{{({v}_{l+1}-{u}_{l+1,p})}^{2}}{2{T}_{l+1,p}}}.\end{eqnarray} \tag{ 22 }$$

In the Maxwell molecules case ($\beta =0$), where one has $\langle {\delta }_{{y}_{p},l}\rangle =1/L$, computations using the LEA give⁵

$$\begin{eqnarray}&&\langle {j}_{l,p}\rangle =\displaystyle \frac{1+\alpha }{2L}({u}_{l,p}-{u}_{l+1,p}),\end{eqnarray} \tag{ 23a }$$

$$\begin{eqnarray}&&\langle {J}_{l,p}\rangle =\displaystyle \frac{1+\alpha }{2L}({T}_{l,p}-{T}_{l+1,p}+{u}_{l,p}^{2}-{u}_{l+1,p}^{2}),\end{eqnarray} \tag{ 23b }$$

$$\begin{eqnarray}&&\langle {d}_{l,p}\rangle =\displaystyle \frac{{\alpha }^{2}-1}{2L}[{T}_{l,p}+\displaystyle \frac{{T}_{l+1,p}+{T}_{l-1,p}}{2}+{({u}_{l,p}-\displaystyle \frac{{u}_{l+1,p}+{u}_{l-1,p}}{2})}^{2}+{(\displaystyle \frac{{u}_{l+1,p}-{u}_{l-1,p}}{2})}^{2}].\end{eqnarray} \tag{ 23c }$$

The hydrodynamic limit consists of a change of spatial and time scales, from (l,p) to (x,t), related by size-dependent factors:

$$\begin{eqnarray}&&x=l/L,\qquad t=p/{L}^{3}.\end{eqnarray} \tag{ 24 }$$

This implies for a generic function ${f}_{l,p}$

$$\begin{eqnarray}&&{f}_{l+1,p}-{f}_{l,p}=\displaystyle \frac{1}{L}\displaystyle \frac{\partial }{\partial x}f(x,t)+{\mathcal{O}}(\displaystyle \frac{1}{{L}^{2}}),\end{eqnarray} \tag{ 25 }$$

$$\begin{eqnarray}&&{f}_{l,p+1}-{f}_{l,p}=\displaystyle \frac{1}{{L}^{3}}\displaystyle \frac{\partial }{\partial t}f(x,t)+{\mathcal{O}}(\displaystyle \frac{1}{{L}^{6}}),\end{eqnarray} \tag{ 26 }$$

which are introduced in equations (20) and (21) to get the final continuous equations in (x, t). Each discrete spatial derivative introduces a ${L}^{-1}$ factor in the HL. Then, the difference between the current terms in the balance equations is of the order of ${L}^{-3}$, because the average currents $\langle {j}_{l,p}\rangle$ and $\langle {J}_{l,p}\rangle$ are of the order of ${L}^{-2}$. Those terms, therefore, perfectly balance the $1/{L}^{3}$ dominant scaling on the left-hand side, i.e. the time-derivative. Since $\langle {d}_{l,p}\rangle$ is of the order of $(1-{\alpha }^{2})/L$, to match the scaling $1/{L}^{3}$ of the other terms, we define the cooling rate to be $\nu =(1-{\alpha }^{2}){L}^{2}$, which is assumed to be order 1 when the limit is taken. This choice automatically implies that when L increases one has that α approaches unity, a further reason to expect the validity of the LEA.

By retaining only the highest order terms in the equations, we get expression (7) of the main text. It is interesting to note that our expansion in terms of ${L}^{-1}$ is similar in spirit to the Chapman–Enskog expansion up to Navier–Stokes order, since we are keeping up to terms of the second order in the gradients (of the order of k², with k being the wave vector, in Fourier space). From a purely mathematical point of view, (7) becomes exact in the limit $L\to \infty$, but $\nu =(1-{\alpha }^{2}){L}^{2}$ of the order of unity, as stated in the previous paragraph. Interestingly, the dissipation field ${d}_{l,p}$ in equation (23c) admits an expansion in even powers of the gradients, as is also the case of granular fluids [37, 38]. However, in the above limit, the first terms in ${d}_{l,p}$ including the gradients are of the order of ${L}^{-2}$ as compared to the contribution $-\nu T$ at the Navier–Stokes order, that is, they would only be considered at the so-called Burnett order.

From a physical standpoint, equation (7) is approximately valid whenever the terms neglected upon writing it are negligible against the ones we have kept. Since the correlations $\langle {v}_{i}{v}_{i\pm 1}\rangle$ are expected to be of the order of ${L}^{-1}$ as compared to the granular temperature⁶ , we must impose that $L\gg 1$ and also $t\ll L$. On the other hand, the term proportional to the correlations in the evolution equation for the granular temperature is therefore of the order of $(1-{\alpha }^{2}){L}^{-1}$, which must be negligible against the second spatial derivative terms, of the order of ${L}^{-2}$. Then, $(1-{\alpha }^{2})L\ll 1$ must be further imposed when the correlations are neglected in equation (7). This condition, although less restrictive than $1-{\alpha }^{2}={\mathcal{O}}({L}^{-2})$, also reinforces the validity of the LEA. On the other hand, when the correlations are fully taken into account, as is the case of equations (15) and (16) of the main paper, the value of α is not restricted since the only assumption for writing them is that of homogeneity.

Noises with respect to averages appear in the currents ${j}_{l,p}={\tilde{j}}_{l,p}+{\xi }_{l,p}$, ${J}_{l,p}={\tilde{J}}_{l,p}+{\eta }_{l,p}$, and in the dissipation ${d}_{l,p}={\tilde{d}}_{l,p}+{\theta }_{l,p}$, with noises ${\xi }_{l,p}$, ${\eta }_{l,p}$ and ${\theta }_{l,p}$ defined to have zero average. The idea is that each term x is made of a $\tilde{x}$ part which is an average over the fast noise (that is, the collisions, which are counted by the fast stochastic variable ${y}_{l,p}$), but at fixed ${v}_{l,p}$ whose evolution is assumed to be slower than noise.

To obtain the correlations of noise, we exploit a series of conditions. Explicit calculations are discussed here for the case of the momentum current noise ${\xi }_{l,p}$. It is clear that the definition ${\tilde{j}}_{l,p}=(1+\alpha ){{\rm{\Delta }}}_{l,p}/2L$ corresponds to the above prescription for the noise. First, it is straightforward that $\langle {\xi }_{l,p}{\xi }_{l^{\prime} ,p^{\prime} }\rangle =0$ for $p\ne p^{\prime}$, because y_p and ${y}_{{p}^{\prime }}$ are independent random numbers. Second, we take into account that $\langle {\delta }_{{y}_{p},l}{\delta }_{{y}_{p},l^{\prime} }\rangle ={\delta }_{l,l^{\prime} }\langle {\delta }_{{y}_{p},l}\rangle ={\delta }_{l,l^{\prime} }/L,$ and the fact that all the other contributions are of the order of ${L}^{-2}$. Thus, for $p=p^{\prime}$ we have $\langle {\xi }_{l,p}{\xi }_{l^{\prime} ,p}\rangle ={(1+\alpha )}^{2}\langle {{\rm{\Delta }}}_{l,p}^{2}\rangle {\delta }_{l,l^{\prime} }/4L+O({L}^{-2})$. At this point, the quasi-elasticity of the microscopic dynamics makes it possible to (i) substitute $(1+\alpha )/2$ by 1 and (ii) calculate $\langle {{\rm{\Delta }}}_{l,p}^{2}\rangle$ by using the LEA, to obtain

$$\begin{eqnarray}&&\langle {\xi }_{l,p}{\xi }_{l^{\prime} ,p^{\prime} }\rangle \sim \displaystyle \frac{1}{L}\;2{T}_{l,p}\;{\delta }_{l,l^{\prime} }\;{\delta }_{p,p^{\prime} }.\end{eqnarray} \tag{ 27 }$$

In the large size system, ${j}_{l,p}$ scales as ${L}^{-2}$ (see [32]). Therefore, the mesoscopic noise of the momentum current is defined as $\xi (x,t)={\mathrm{lim}}_{L\to \infty }\;{L}^{2}{\xi }_{l,p}$, and $j(x,t)=\tilde{j}(x,t)+\xi (x,t)$, in which, again, $\tilde{j}(x,t)={\mathrm{lim}}_{L\to \infty }\;{L}^{2}{\tilde{j}}_{l,p}$. Going to the continuous limit, remembering equation (24), and taking into account that ${\delta }_{l,l^{\prime} }/{\rm{\Delta }}x\sim \delta (x-x^{\prime} )$ and ${\delta }_{p,p^{\prime} }/{\rm{\Delta }}t\sim \delta (t-t^{\prime} )$ we get the noise amplitude of the velocity current in the main text. Identical considerations lead to the amplitude for the energy current noise. For the fluctuations of dissipation, the dissipation term is split again as ${d}_{l,p}={\tilde{d}}_{l,p}+{\theta }_{l,p},$ with $\langle {d}_{l,p}\rangle =\langle {\tilde{d}}_{l,p}\rangle$. We know from the dissipation current definition that $\langle {\theta }_{l,p}{\theta }_{l^{\prime} ,p^{\prime} }\rangle =0$ for $p=p^{\prime}$. Making use of the LEA and in the large size system ${d}_{l,p}$ scales as ${L}^{-3}$ and it is expected that the noise should have the same scaling. Going to the continuous limit and taking in account equations (24) and (27), the result in the main text is recovered.

The cross correlations between different noises are straightforwardly obtained, along similar lines, with the result

$$\begin{eqnarray}\langle \xi (x,t)\eta (x,t)\rangle & = & \langle \eta (x,t)\xi (x,t)\rangle =\displaystyle \frac{4\;T(x,t)u(x,t)}{L},\\ \langle \xi (x,t)\theta (x,t)\rangle & = & \langle \theta (x,t)\xi (x,t)\rangle =0,\\ \langle \eta (x,t)\theta (x,t)\rangle & = & \langle \theta (x,t)\eta (x,t)\rangle =0.\end{eqnarray} \tag{ 28 }$$

A comparison for the amplitudes of noise for the velocity and energy currents is shown in figure C1 . A case with the MM interaction ($\beta =0$) is considered. The simulations are performed with periodic boundary conditions, therefore without energy injection, and starting with a non-homogeneous initial condition. The initial mesoscopic velocity profile and homogeneous granular temperature are $u(x,0)={u}_{0}\mathrm{sin}(2\pi x)$ and $T(x,t)\equiv {T}_{0}$, respectively, with ${u}_{0}={T}_{0}=1$.

Zoom In Zoom Out Reset image size