{'fraction_non_alphanumeric': 0.03884621539446311, 'fraction_numerical': 0.01168186909905585, 'mean_word_length': 4.451908396946565, 'pattern_counts': {'":': 0, '<': 3, '<?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': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Fluctuating hydrodynamics for dilute granular gases: a Monte Carlo study

We investigate hydrodynamic noise in a dilute granular gas during the homogeneous cooling state, by means of a proper application of the Direct Simulation Monte Carlo (DSMC) algorithm. The DSMC includes a source of randomization which is not present in Molecular Dynamics (MD) for inelastic hard disks. Notwithstanding this difference, a fair quantitative agreement is found, including a violation of the fluctuation-dissipation relation for the noise amplitude of the same order observed in MD. This study suggests that deterministic collision dynamics is not an essential ingredient to reproduce, up to a good degree of approximation, hydrodynamic fluctuations in dilute granular gases.


I. INTRODUCTION
Granular gases represent an important benchmark for theories and methods of non-equilibrium statistical mechanics [1]. A granular gas is the idealized counterpart of fluidized granular media, i.e. collections of N ∼ 10 2 ÷ 10 3 macroscopic grains (typically spheres with diameters σ in the range 10 −4 ÷ 10 −3 m), contained in volumes V such that the packing fraction, or occupied volume, is smaller (often much smaller) than 10%: note that the small packing fraction ∼ N σ d /V is still compatible with a non-negligible collision rate, ∼ N σ d−1 /V , where d is the dimensionality of the system. Within these conditions, the main difference with respect to an ideal gas is given by inelastic collisions: kinetic energy is lost during an (hard core) interaction, ruling out many basic assumptions of equilibrium statistical mechanics and leading to several phenomena such as: breakdown of energy equipartition [2][3][4], spontaneous symmetry breaking (vortices and clustering) [5,6], "Maxwell demon" [7], "ratcheting" effects [8], and lot more. In both theoretical investigation and in real world applications, granular "particles" are often found in a small number, thousands is a common order of magnitude, as mentioned before. This makes fluctuations a fundamental ingredient of granular statistical descriptions [6]. While continuum theories, e.g. granular hydrodynamics and kinetic theory, describing local averages for density, velocity and energy and their evolution, have achieved an advanced degree of maturity [9], a good description of noise accompanying such averages is still an open problem. Near equilibrium, Einstein fluctuation theory, resumed by Landau and Lifshitz in their fluctuating hydrodynamics [10], based on equilibrium fluctuationdissipation relations (FDR) [11,12], simplifies the description of noise for transport processes. In granular gases, fluctuations do not satisfy, in general, FDR. The simultaneous presence of external random drivings and diluteness, usually guarantees the validity of FDR [13,14]. FDR is, instead, violated at packing fractions > 10% [15,16] and in a non-driven granular gas, even dilute, in the so-called Homogeneous Cooling State (HCS) [17,18].
A recent study on dilute inelastic hard particles in the HCS, with Molecular Dynamics (MD) simulations [18] and then an analytical treatment [19], have shown that fluctuations of the hydrodynamic transverse velocity field deviate from the classical Landau-Lifshitz theory for two reasons: 1) noise is not white, 2) it does not satisfy the FDR of the 2nd kind. The analytical study is based on projection operator formalism, starting from the "pseudo-Liouville" operator for the deterministic evolution which only includes free flights and instantaneous hard-core inelastic collisions at contact. The only statistical averaging is made with respect to an ensemble of initial conditions, while several hypothesis are introduced: among the others, spatial homogeneity and a special time-scaling (inherent to HCS), for one-particle and two-particle distributions, such that all time-dependence is contained in the "granular temperature" T g = m v 2 /d. The last part of this analysis makes use of the Molecular Chaos assumption, which is reasonable for dilute systems. Anyway, one wonders which ingredients of this sophisticated analysis are essential and which can be neglected. For instance, the use of Molecular Dynamics (MD) for simulations and of the rigorous manipulation of the pseudo-Liouville operator as starting point of the analytical treatment let us think that the deterministic nature of the hard spheres dynamics is crucial.
Here we present a different model, which in principle could produce different results, but in the end reveals to be in fair agreement with the previously cited study. It is a stochastic model, where collisions are treated randomly, in sharp contrast with deterministic dynamics employed in Molecular Dynamics simulations and described by the pseudo-Liouville operator. It consists of the classical "Direct Simulation Monte Carlo" (DSMC) algorithm [20] in its non-homogeneous form, i.e. placed on a spatial grid to measure modes of the transverse velocity field, with the use of "fictive particles" to guarantee Molecular Chaos (the algorithm is explained in the Appendix). The inelastic collision rule, relating post-collisional velocities to pre-collisional ones, and the time-rescaling are the only ingredients which are conserved with respect to the MD model. Our question is: are those ingredients sufficient to reproduce the hydrodynamic noise properties of the deterministic model? Similar studies have been performed for hydrodynamic fluctuations in the elastic case [21,22]. More recently, the DSMC approach has also been applied to fluctuations of the global energy fluctuations (a homogeneous, not spatially dependent quantity) in inelastic systems [23,24], obtaining a good agreement with the amplitude of fluctuations measured in MD simulations and explained in terms of projection operators applied to the deterministic dynamics [25]. A similar conclusion was available for that study: stochastic "Monte Carlo" treatment of collisions is sufficient to reproduce MD noise properties. A successive study [26] also showed how certain quantities (cumulants of fluctuations in the stationary state) may be obtained without resorting to projection techniques.
In the present study, in a sense, the observation is even more striking: the exchange of particles among different "copies", needed to reduce the incidence of re-collisions and finite size correlation effects, is a strong source of randomization. Surprisingly it does not affect the hydrodynamic fluctuations.
The organization of the paper is the following: in Section II we first describe the model, then in Section III we present the results, first verifying Molecular Chaos and then measuring the hydrodynamic noise in time-autocorrelation and amplitude. Finally we draw Conclusions. We have dedicated an Appendix to a detailed description of the DSMC algorithm used for this work.
where v ′ corresponds to post-collisional velocity, α is the coefficient of restitution and v ij = v i − v j is the relative velocity. The termσ is a random unit vector (see Appendix). The size of the system is L × L and the boundary conditions are periodic. The main hypothesis used here is that the system, starting in a uniform equilibrium state with granular temperature T g (0) = T 0 , evolves to the Homogeneous Cooling State (HCS), which is characterized by a single time-scale measured by the temperature T g (t): any other quantity depends on time only through T g (t). Apparently, as observed in many previous studies [27], in a homogeneous setting, Molecular Chaos is sufficient to guarantee this hypothesis.
The HCS is unstable against spatial fluctuations: this instability appears at scales larger than a critical length L c (α), therefore it can be avoided by taking the linear size of the system L < L c [28]. It is possible to analyze the effects of spatial fluctuations by deriving mesoscopic equations through a linearization around the HCS [5]. The resulting equations are generally coupled, but in the Fourier representation the transverse velocity field results decoupled from the other modes. We are interested, in particular, in the fluctuations of its largest mode, i.e. of the smallest wave-number |k m | = 2π/L. We choose the wave vector k m parallel to the axisx and the quantity above results to be where v y,j (t) is the component of the velocity of the particle j in the directionŷ and x j (t) is its coordinate along thê x axis. We want to verify the hypothesis that the fluctuations of U ⊥ (t) obey a linear Langevin equation, [5,18] where ν(t) is the kinematic viscosity which, in the HCS, is proportional to T g (see [9] for definitions), and v th (t) ≡ 2T g /m.
The last term in the Eq. (3) describes an internal noise corresponding to the rapid (microscopic) degrees of freedom of the system. This complex noise is assumed here to be white, Gaussian and with correlations given by with ξ * (t) the complex conjugate of ξ(t). Deviations from the white noise assumption are expected [19], but they are in general quite small. The most evident consequence of inelasticity is, instead, the variation of the dimensionless coefficient c in Eq. (3), which is equal to 1 at equilibrium (α = 1); its departure from c = 1 represents the fact that the fluctuations of U ⊥ (t) do not satisfy the FDR.
Since the system is cooling, the typical velocity of the particles becomes smaller and smaller, leading to increasing rounding errors and other numerical problems. These problems can be avoided with the use of a procedure which, in the HCS, maps the dynamics to a steady state by means of a time-rescaling [27,29]. Such procedure is explained in the next subsection.

B. The stationary representation
In the HCS, the granular temperature obeys the following equation: where ζ H is the cooling rate that results proportional to T g (t).
The stationary representation of the HCS [27] consists in introducing a new time scale τ defined by with ω 0 ant t 0 arbitrary constants, implying the definition of rescaled velocitiesṽ(τ ) = v(t)ω 0 t. It is easy to see [27] that observing the system on this new time scale is equivalent to apply a positive continuous drag to all particles ∂ τṽ (τ ) = ω 0ṽ (τ ). This naturally leads to define also the rescaled analogous of U ⊥ (t): We tune ω 0 in order to have T g (τ → ∞) = T 0 (so that the length of transients is reduced), a result which is obtained by taking ω 0 = ζ H (0)/2. In the steady state, we simplify the notation using v th = v th (0). With these choices, the Langevin equation (3) is mapped onto a new equation where λ 0 = L 2 /(N σ) ≡ 1/(nσ) is proportional to the mean free path andν = ν(t)/[λ 0 v th (t)] andk m = k m λ 0 are dimensionless rescaled viscosity and wave number, respectively. Ifνk 2 m − ω 0 < 0, which is equivalent to the condition of stability of shear modes in the HCS, the above equation leads to a dimensionless autocorrelation function given by with τ −1 0 ≡νk 2 m v th /λ 0 − ω 0 the characteristic time of decay. Based on Eq. (9), we can obtain the dimensionless kinematic viscosityν and the "FDR violation" coefficient c from a measure of C ⊥ (τ ):ν i.e.ν and c are obtained measuring the amplitude C ⊥ (0) and the decay time τ 0 .

III. NUMERICAL RESULTS
The aim of this section is to measureν and c by means of Equations (10) and (11) in DSMC simulations, comparing them to the results of MD. For this purpose a DSMC algorithm must be carefully devised to guarantee the validity of the HCS assumptions and to allow for the measure of W ⊥ (τ ). This is done in the next subsection.

A. Stationary representation for the inhomogeneous DSMC
The DSMC algorithm used here is explained in details in the Appendix. The essence of a DSMC algorithm, which has been originally designed as a tool to numerically solve Boltzmann equations, is a stochastic computation of collisions: pairs of particles and their orientation vectorσ are chosen randomly, with probabilities dictated by the Boltzmann collision integral. In homogeneous configurations, a pair of particles is chosen among all N particles of the system, disregarding spatial coordinates: the large number N guarantees Molecular Chaos. Here, to avoid shear mode instabilities, we keep the total size of the system under the critical size L c , expecting to observe a homogeneous regime. Nevertheless, the quantity under scrutiny, U ⊥ (t), requires that collisions are treated with a good spatial resolution, i.e. colliding particles must be close to each other to give the correct contribution to the variation of U ⊥ (t). For this reason, a inhomogeneous DSMC algorithm is necessary. The system is partitioned in m c non-overlapping cells of size l c : during the free-streaming step each particle can move from a cell to a neighboring one, but during the collision step the particle can collide with particles in the same cell only. The resolution for the measure of U ⊥ (t) is improved increasing m c ; at the same time, when m c increases, the average number of particles in each cell N c = N/m c decreases, threatening the Molecular Chaos assumption. In principle the perfect resolution would be achieved with l c ∼ σ, but this would result in a number of particles per cell N c ∼ nσ 2 which, for diluteness, is required to be much smaller than 1. We will see in the following (Fig. 2) that taking l c ∼ 15σ is sufficient to restore a good resolution for measuring U ⊥ (t). Such a choice, however, at the chosen packing fraction nσ 2 = 10 −2 gives N c ∼ 2.5, which is so low that the Molecular Chaos could be invalidated. A first check is to study the sensitivity of the rescaling procedure (7) (tested in [27] in the homogeneous DSMC simulations and in [30]) to the parameter m c (or equivalently N c ).
Having this aim, we have characterized the entire system analyzing its granular temperature T g as function of the rescaled time τ for different choices of the cell number m c at α and N fixed (see Fig. 1), comparing it with MD results. The time-rescaling is applied with a choice of ω 0 such that, in the steady state it is expected T g = T 0 . From the simulation data it results that a steady state is always reached, i.e. the granular temperature is constant in all cases. A homogeneous DSMC (i.e. m c = 1) gives an agreement with the expected value which is better, of a few percentages, with respect to MD. When m c is increased, however, the DSMC gives very bad results. This is due to the small value of N c , which results in strong finite size effects, for instance fake re-collisions which invalidate the Molecular Chaos hypothesis. The way of solving this problem is to simulate a larger number of particlesÑ = f N with f > 1, keeping the collision statistics per particle as that of a N -particles system (see the Appendix for a precise description). This is equivalent to simulate f copies of the original N -particles system and let particles of different copies swap at each time step. The effect of f > 1 is immediately seen in Fig. 1: the steady state granular temperature perfectly agrees with that measured in MD, as soon as f ∼ 10, even if m c = 225. Much larger values of f appear to slightly improve the value of T g . On the other side a too much large value of f can be very expensive in terms of cpu-time.
A crucial quantity, here, is the number of virtual particles per cellÑ c = f N c = f N/m c . If m c is increased (for instance in Fig. 1 it is changed from 225 to 400), f should be increased. Apparently it is sufficient to haveÑ c ∼ 25 to have a good agreement, for T g , between DSMC and MD. As discussed above, the necessity of a large number of cells comes from the sensitivity of U ⊥ (t) to the spatial localization of collisions. This is shown in Fig. 2, where the decay of C ⊥ (τ ) is displayed: the homogeneous DSMC (m c = 1) shows an exponential decay with a τ 0 completely different from the one observed in the MD simulations. This behavior is insensitive to changes of f from 1 to 10. In order to obtain the correct decay, it is necessary to use the inhomogeneous DSMC approach with m c ≫ 1. A partition of m c = 20 × 20 for N = 1000, in fact, reveals an exponential decay very close to the MD results: also in this case, a change of f from 10 to 100 is unimportant. When N is increased at fixed density (and mean free path), the size of cells must be kept constant, which implies m c ∝ N . If f is also kept constant,Ñ c is automatically preserved and we expect to have a good comparison with MD simulations in terms of Molecular Chaos, as well as decay of C ⊥ (τ ). This situation is fairly verified in Fig. 3, where the decay of C ⊥ (τ ) is shown to be always close to the MD results.
This approach is valid in general and, with proper values of f and m c , we can compare, for different values of α and N , the properties of the fluctuations of W ⊥ between DSMC and MD simulations. Guided by this preliminary analysis, we have set the value of f = 10 andÑ c = 25 for all the following discussion.

B. Analysis of the fluctuations
In the previous section we made the hypothesis that the fluctuations of the slowest mode of the transverse velocity field U ⊥ (t) follow the Langevin equation (3), or (8) in the rescaled representation: this implies that its correlation function follows the exponential decay in Eq. (9), with characteristic time τ 0 . The density is the same as in Fig. 1.
In order to apply the correct thermostatv = ω 0ṽ , and to obtain the kinematic viscosityν from Eq.(9), it is necessary to know the value of ω 0 = ζ H (0)/2. This is computed, in the homogeneous case, giving the theoretical expressions [27] ζ where the coefficient a 2 (α) is We have verified that the prediction in Eq. (12) is good also for the inhomogeneous DSMC and for MD. Eq. (9) is expected to be valid in the large N limit: for this reason we have performed simulations for different values of N (1000, 1500, 2000, 2500 and 3000), fixing the particle density n of the system, and measuring the time decay τ 0 . The data obtained using (10) and (12) in DSMC and MD, are shown in Fig. 4. The dimensionless viscositỹ ν appears to be a linear function of N −1 , therefore a best linear fit allows to extrapolate the asymptotic value for N → ∞. The agreement between MD and DSMC data is good in all cases, even at not too high values of N .
Our second question is if the inhomogeneous DSMC can reproduce the MD results for the noise term in the Langevin equation (3) or (8), in particular the quantityν ′ ≡ cν that is the viscosity term appearing in the noise. Analogously to Fig. 4, we have analyzed its dependence versus N −1 . Also in this case we have found a linear law, shown in Fig. 5: a small discrepancy is observed, for this quantity, between DSMC and MD data, when N is increased. The disagreement results more evident if the collisions are more inelastic. Our interpretation of this small difference is that the DSMC dynamics at small time-scales is slightly different from MD: indeed, this disagreement is less evident at large time-scales, for instance in the measure ofν. If we consider also that the DSMC and MD results about the correlation of the energy fluctuations are very similar [23], we tend to conclude that this small discrepancy could be due to the use of small cells (large m c or equivalently smallÑ c ), which is necessary in the measures of C ⊥ (0). A more general discussion of the cause of this discrepancy is given at the end of this section.
Notwithstanding this slight discrepancy, the value of c =ν ′ /ν results significantly larger than 1 also in DSMC. The inhomogeneous DSMC is then able to reproduce the breakdown of the FDR as already found with the MD approach [18]. This is evident from Fig. 6, where theν andν ′ are plotted as function of α for different values or the particle number and for N → ∞. The extrapolated results for N → ∞, forν are in good agreement with the MD data and with the theoretical analysis of [19], while those forν ′ still underestimate them, as shown in the fourth panel of Fig.  6.
The ratio c is shown in Fig. 7: it clearly displays strong deviations from 1 (violations of FDR), which reach 12% at α = 0.8 also in DSMC. The dot-dashed curve represents the result for c from the theoretical analysis of [19] with the additional assumption of white noise. We note a close comparison with the data from DSMC, suggesting that the effect of this protocol is an effective time-decorrelation of noise. Such a conclusion should be taken with care, in view of other simulations, performed with different choices of f and m c : these results are not shown here, since have a quite high statistical uncertainty due to the large computational time required (increasing f and m c results in an increase of fictive particlesÑ ). They suggest that, when f and m c are increased keepingÑ c = f N c constant, the value of c stays substantially unchanged. On the other side ifÑ c is increased (e.g. by increasing f at constant m c ) the value of c slightly increases and seems to slowly tend toward the corresponding MD value, which is also compatible with the theoretical analysis with a colored noise. The presence of time-correlation in the noise, even in DSMC, is in agreement with the results for C ⊥ (τ ) shown in Fig. 2 and 3, which are not straight exponentials, but some bending can be observed at small τ .
Apart from these small discrepancies, the main fact is that the inhomogeneous DSMC is able to reproduce the violation of FDR, which was not obvious. In fact, it is reasonable and well verified that average values (e.g. transport coefficients) obtained in DSMC agree with dilute MD simulations; anyway it is less trivial to observe good agreement for fluctuations, which -in non-equilibrium situations -can be much more sensitive to the detailed mechanisms of the dynamics.

IV. CONCLUSIONS
In this work we have obtained two main results. First, we have individuated which are the basic ingredients -in a Direct Simulation Monte Carlo algorithm -to correctly measure the time-decay of the autocorrelation function of fluctuations of a hydrodynamic field, keeping valid the assumption of Molecular Chaos: the volume must be divided in small cells for the purpose of collision computation, even if the system is assumed to be spatially homogeneous; the good resolution to appreciate the "real" time-decay of C ⊥ (τ ) is achieved only when N c ∼ 2 which is too small to avoid fake recollisions: for this reason, a number (order ∼ 10) of "virtual" copies of the system is necessary to restore Molecular Chaos. This procedure is similar to that adopted in [31]. The study of the transverse velocity field could also be done by dividing the system in slides, as in [21,24]. Such a procedure can be used to analyze U ⊥ (τ ), while we have no knowledge of its adaptation to other space-dependent fields. Note that -when DSMC is used to study the  Fig. 6). The circles and squares correspond to DSMC and MD data, respectively. The solid is the theoretical prediction [19]. The value of reference c = 1 (dashed line) is obtained if FDR are satisfied. The dot-dashed curve is the result of the theoretical analysis of [19] with the additional assumption of white noise.
homogeneous cooling state, i.e. keeping the system size below the critical size for instability, L < L c -one usually does not divide the system in cells, and still obtains excellent agreement for one-point observables. Previous study of cooling granular gases with inhomogeneous DSMC [30,31] where in fact proposed to analyze the departure from the HCS.
Our second result concerns a quantitative agreement with Molecular Dynamics results, which is very good for the exponential tail of the autocorrelation of fluctuations. The amplitude of fluctuations, on the other side, shows small discrepancies with respect to MD results, which become more evident as α is reduced and N increases. In the extrapolated N → ∞ limit, these discrepancies can reach (for α = 0.8) an underestimation of order 10%: nevertheless it is still possible to appreciate a ∼ 12% of deviation from the validity of FDR, which is one of the most peculiar property of granular hydrodynamic fluctuations in the HCS. In DSMC simulations, deviations from a perfect exponential decay of C ⊥ (τ ) are also visible, at small times, but these are only qualitatively similar to those observed and predicted in [19].
Hydrodynamic noise is always important in granular systems, which -in terms of number of elementary constituents -are much smaller than molecular fluids (10 4 instead of 10 20 particles). Future research is needed to understand the properties of granular hydrodynamic noise in steady state models: even if we have shown that these properties do not depend on the detailed collisional mechanism, it seems that they depend on the choice of energy driving protocol [32]. This is just an alternative representation of the algorithm, which underlines the randomizing mechanism -intrinsic of DSMC and not present in MD.