{'fraction_non_alphanumeric': 0.03764423626185308, 'fraction_numerical': 0.010967668228036102, 'mean_word_length': 4.623835528429168, 'pattern_counts': {'":': 1, '<': 0, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Close-packed granular clusters: hydrostatics and persistent Gaussian fluctuations

Dense granular clusters often behave like macro-particles. We address this interesting phenomenon in a model system of inelastically colliding hard disks inside a circular box, driven by a thermal wall at zero gravity. Molecular dynamics simulations show a close-packed cluster of an almost circular shape, weakly fluctuating in space and isolated from the driving wall by a low-density gas. The density profile of the system agrees very well with the azimuthally symmetric solution of granular hydrostatic equations employing constitutive relations by Grossman et al., whereas the widely used Enskog-type constitutive relations show poor accuracy. We find that fluctuations of the center of mass of the system are Gaussian. This suggests an effective Langevin description in terms of a macro-particle, confined by a harmonic potential and driven by a delta-correlated noise. Surprisingly, the fluctuations persist when increasing the number of particles in the system.


Introduction
Despite much progress in the last two decades, modeling of granular flow remains a challenge for physicists and engineers [6; 15]. The simplest case for a first-principle modeling seems to be a granular gas, or rapid granular flow: a flow dominated by binary particle collisions [5; 11]. Here one can use the model of inelastically colliding hard spheres to develop kinetic and hydrodynamic descriptions. The simplest version of this model involves binary collisions with a constant coefficient of normal restitution ε: where primed quantities stand for post-collisional velocities, and e i j ≡ (r i − r j ) / r i − r j . Assuming molecular chaos allows to use the Boltzmann or Enskog kinetic equation, properly generalized to account for the inelasticity of the particle collisions. Systematic derivations of the hydrodynamic equations from the Boltzmann or Enskog equation [3; 21; 26] use an expansion in the Knudsen number and, in some versions of theory, in inelasticity [26]. A finite inelasticity immediately brings complications [11]. Correlations between particles, developing already at a moderate inelasticity, may invalidate the molecular chaos hypothesis [2; 25; 27]. The normal stress difference and deviations of the particle velocity distribution from the Maxwell distribution also become important for inelastic collisions. As a result, the Navier-Stokes granular hydrodynamics should not be expected to be quantitatively accurate beyond the limit of small inelasticity.
A further set of difficulties for hydrodynamics arise at large densities. In mono-disperse systems strong correlations appear, already for elastic hard spheres, at the disorder-order transition. As multiple thermodynamic phases may coexist there [7], a general continuum theory of hard sphere fluids, which has yet to be developed, must include, in addition to the hydrodynamic fields, an order-parameter field. Furthermore, even on a specified branch of the thermodynamic phase diagram, we do not have first-principle constitutive relations (CRs): the equation of state and transport coefficients. Last but not least, kinetic theory of a finite-density gas of elastic hard spheres in two dimensions has notorious difficulties related to the long-time tail in the velocity pair autocorrelation function [28].
It has been argued recently that discrete particle noise may play a dramatic role in rapid granular flow [1; 12; 23]. A new challenge for theory is a quantitative account of this noise. A promising approach at small or moderate densities is "Fluctuating Granular Hydrodynamics": a Langevintype theory that takes into account the discrete particle noise by adding delta-correlated noise terms in the momentum and energy equations, in the spirit of the Fluctuating Hydrodynamics by Landau and Lifshitz [18]. There is a recent progress in this area in the case of small densities [4]. For high densities such a theory is presently beyond reach.
This work addresses granular hydrodynamics and fluctuations in a simple two-dimensional granular system under conditions when existing hydrodynamic descriptions [3; 21; 26] break down because of large density, not large inelasticity. In view of the difficulties mentioned above, attempts of a first-principle description of hydrodynamics and fluctuations should give way here to more practical, empiric or semi-empiric, approaches. One such approach to a hydrodynamic (or, rather, hydrostatic, as no mean flow is present) description was suggested in 1997 by Grossman et al. [13]. The present work puts it into a test in an extreme case when macro-particles (granular clusters with the maximum density close to the hexagonal close packing) form. The model system we are dealing with was first introduced by Esipov and Pöschel [9]. It is an assembly of N ≫ 1 identical disks of mass m, diameter d and coefficient of normal restitution ε, placed inside a circular box of radius R at zero gravity. The circular wall of the box is kept at constant temperature T 0 . We measure, using molecular dynamics (MD) simulations, the radial density profiles of the system, including the close-packed part. Furthermore, we solve numerically a set of granular hydrostatic equations which employ the CRs by Grossman et al. and show that, in a wide range of parameters, there is good agreement between the two. We also show that, for the same setting, the Enskog-type CRs [14] perform poorly. Finally, we investigate in some detail the fluctuations of the macro-particle's position, by measuring the radial probability distribution function (PDF) of the center of mass of the system. These fluctuations turn out to be Gaussian, which suggests an effective Langevin description of the system in terms of a macro-particle, confined by a harmonic potential and driven by a white noise. Surprisingly, the fluctuations persist as the number of particles in the system is increased.

MD simulations
We employed a standard event-driven algorithm and thermal wall implementation [24]. We put m = d = T 0 = 1 and fixed R = 100 and ε = 0.888, while the total number of particles N served as the control parameter in the simulations. We were mostly interested in a hydrodynamic (low Knudsen number) regime, when the mean free path is small compared to the system size. This requires N ≫ R/d = 100. For N in the range of a few hundred, we observed a dilute granular gas with an increased density in the center of the box. The clustering in the center becomes more pronounced as N grows. The clustering can be easily explained in the hydrodynamics language: because of the inelastic collisions the granular temperature goes down as one moves away from the wall toward the center of the box. Combined with the constancy of the pressure throughout the system, this causes an increased particle density at the center. As N increases, the particle density in the center approaches the hexagonal close packing value n c = 2/( √ 3d 2 ). Figure 1 shows snapshots of the system for three different, but sufficiently large, values of N. A perfect hexagonal packing is apparent. Movies of these simulations show that the cluster position fluctuates around the center of the box, while the cluster shape fluctuates around a circular shape.
Our diagnostics focused on two radial distributions: the number density of the particles n(r) and the PDF of a ra-dial position of the center of mass of the system P(r m ), see below. As we are interested in steady-state distributions, we disregarded initial transients.

Hydrostatic theory
As the fluctuations are relatively weak, it is natural to start with a purely hydrodynamic description. In the absence of time-dependence and for a zero mean flow this is essentially a hydrostatic theory. It operates only with (time-independent) granular density n(r), temperature T (r) and pressure p(r). The energy input at the thermal wall is balanced by dissipation due to inter-particle collisions, and one can employ the momentum and energy balance equations: Here κ is the thermal conductivity and I is the rate of energy loss by collisions. Note that the heat flux, entering the thermal balance in Eq. (2), does not include an inelastic term, proportional to the density gradient [3; 21; 26]. In the nearly elastic limit 1−ε ≪ 1, that we are interested in, this term can be neglected. The boundary condition at the thermal wall is where r and φ are polar coordinates with the origin at the center of the box.
To proceed from here, we need CRs: an equation of state p = p (n, T ) and relations for κ and I in terms of n and T . As we attempt to describe close-packed clusters, the standard techniques, based on the Boltzmann or Enskog equations, are inapplicable. Grossman et al. [13] suggested a set of semi-empiric CRs in two dimensions, that are valid for all densities, all the way to hexagonal close packing. Their approach ignores possible phase coexistence beyond the disorder-order transition and assumes that the whole system is on the thermodynamic branch extending to the hexagonal close packing. Grossman et al. employed free volume arguments in the vicinity of the close packing, and suggested an interpolation between the hexagonal-packing limit and the well-known low-density relations. The resulting CRs [13] read and Here l is the mean free path, which is given by an interpolation formula [13]: where a = 1 − (3/8)  (2) in rectangular geometry. We adopted the same values of α and γ for the circular geometry. Employing Eqs. (4) and (5) we can reduce Eqs.
(2) to a single equation for the rescaled inverse density z(r, φ ) ≡ n c /n(r, φ ). In the rescaled coordinates r/R → r the circle's radius is 1, and the governing equation becomes where and is the hydrodynamic inelasticity parameter introduced in Ref. [20]. As the total number of particles N is fixed, z −1 (r, φ ) satisfies a normalization condition: where is the average area fraction of the particles. Azymuthally-symmetric solutions of Eq. (7), z(r, φ ) = Z(r), are described by the following ordinary differential equation: while the normalization condition (10) becomes Equations (12) and (13), together with the obvious boundary condition form a complete set. The hydrostatic density profile is completely determined by two scaled parameters Λ and f , and can be found numerically as we will show shortly.
Are the azimuthally symmetric states stable with respect to small perturbations? We performed marginal stability analysis to find out whether there are steady-state solutions with broken azimuthal symmetry, z(r, φ ) that bifurcate continuously from an azimuthally symmetric solution Z(r). Under additional assumption that the possible instability of the azimuthally symmetric state is purely growing (that is, not oscillatory), the marginal stability analysis yields the instability borders. The marginal stability analysis goes along the same lines as that developed for the rectangular geometry [16; 17; 20; 19; 23]. Let us search a steady-state close to an azymuthally-symmetric state: As z(r, φ + 2π) = z(r, φ ), k must be an integer which can be chosen to be non-negative. Substituting Eq. (15) into Eq. (7) and linearizing around the azimuthally symmetric state Z(r), we obtain a linear eigenvalue problem, where k = k( f ,Λ ) plays the role of the eigenvalue: where ζ k (r) = Φ k (r) F(Z). Equation (16) is complemented by the boundary conditions ζ k (0) = 0 and ζ k (1) = 0 (17) and can be solved numerically. Let us ignore for a moment the quantization of the eigenvalue k and, while looking for k, assume that it is a (positive) real number. In that case, a numerical solution yields, for a fixed Λ , a curve k = k( f ), see Fig. 2 for two examples. At a fixed k, the azimuthally symmetric state is unstable within an interval of area fractions. The foots of one such curve corresponds to the (hypothetical) case when k tends to zero (that is, the azimuthal wavelength tends to infinity). The instability interval becomes narrower when k is increased, and it shrinks to a point at a maximum k max , signalling that a density modulation with a sufficiently short azimuthal wavelength should be stable for all f . For example, when Λ = 5 × 10 4 , the marginal stability curve k = k( f ) has its maximum at k max ≈ 0.38 (see Fig. 2) that is less than unity. Going back to the physical case, where k is a positive integer, we see that all the values for k, determined from the eigenvalue problem, are unphysical. That is, there are no solutions to the eigenvalue problem that would satisfy the boundary conditions and the quantization condition k = 1, 2, 3 . . . . We observed a similar behavior for values of Λ up to 10 8 . Though k max increases with Λ , the increase is extremely slow: slower than logarithmical (see the inset of Fig. 2). These numerical results strongly indicate that the azimuthally symmetric states are stable with respect to small perturbations. This is in marked contrast with the presence of bifurcating states with broken symmetry in similar settings of a granular gas driven by a thermal wall, but in rectangular [16; 17; 20; 19; 23] and annular [8] geometries. The absence of the bifurcating states with broken symmetry in the circular geometry gives a natural explanation to the persistence of circle-shaped cluster shapes as observed in our MD simulations. Now let us compare the azimuthally symmetric density profiles, found from our hydrostatic calculations, with the results of MD simulations. As noted above, in the MD simulations the coefficient of normal restitution ε = 0.888 was fixed, while the number of particles N was varied. Accordingly, the scaled parameter Λ = 9980 was fixed but f was varied. Figure 3 shows a comparison of the hydrostatic radial density profiles with the radial density profiles obtained in MD simulations, for four different values of N. The MD profiles were averaged over 1000 snapshots. It can be seen from Fig. 3 that for N = 1716 the theory overestimates the density in the center of the box. This is expected as, for relatively small N, the maximum density is considerably less than the close-packing density, and the accuracy of the CRs by Grossman et al. is not as good. For larger N the agreement rapidly improves. We also computed the density profiles using another set of semi-empiric CRs: those obtained in the spirit of Enskog theory [14]. One can clearly see from Fig. 3 that the Enskog-type CRs predict unphysically high densities in the cluster, and is also less accurate at intermediate densities. Figure 4 compares, at different N, the maximum radial densities predicted by hydrostatic theory with the CRs by Grossman et al. and those obtained in the MD simulations. The agreement is very good for the densities approaching the close packing density.

Fluctuations
Now we turn to stationary fluctuations of the radial coordinate r m (t) of the center of mass of the system. The radial PDF P(r m ,t) is normalized by the condition where we have returned to the dimensional coordinate. Typical MD results are presented in Fig. 5, which shows log P(r m ) versus (r m /R) 2 , for four values of N. 1 The observed straight lines clearly indicate a Gaussian. This finding suggests a Langevin description of the macro-particle. Consider a macroparticle, performing an over-damped motion in a confining harmonic potential U(r m ) = kr 2 m /2 and driven by a (deltacorrelated) discrete-particle noise η(t). The Langevin equation for this problem reads where h is the damping rate, η(t)η(t ′ ) = 2hΓ δ (t −t ′ ), and Γ is an effective magnitude of the discrete particle noise. In the limit Γ → 0, the (deterministic) steady-state solution of Eq. (19) is r m =ṙ m = 0: the macro-particle at rest, located at the center of the box. At Γ > 0, the steady state PDF is (see, e.g. Ref. [10]) where σ 2 = Γ /k, and the normalization constant is computed for σ ≪ R. The variance σ 2 is the ratio of Γ (a characteristic of the discrete noise) and k (a macroscopic quantity). Unfortunately, the present state of theory does not enable us to calculate either of these two quantities. An important insight, however, can be achieved from the N-dependence of σ , obtained by MD. In analogy with equilibrium systems, one might expect the relative magnitude of fluctuations to decrease with increasing N. Surprisingly, this is not what we observed, see the bottom panel of Fig. 4. One can see that σ (N) approaches a plateau, that is fluctuations persist at large N. 2

Discussion
We employed a simple model system [9] to investigate the density profiles and fluctuation properties of dense clusters emerging in driven granular gases. The density profiles, obtained in our MD simulations, are well described by hydrostatic equations which employ the constitutive relations suggested by Grossman et al. [13]. The good performance of these constitutive relations is in agreement with previous results in rectangular geometry, with and without gravity [13; 22]. Marginal stability analysis yields a natural explanation to the circular cluster shape, observed in the MD simulations. Finally, the observed Gaussian fluctuations of the center of mass suggest an effective Langevin description in terms of a macro-particle in a confining potential of hydrodynamic nature, driven by discrete particle noise. The fluctuations persist as the number of particles in the system is increased, and this surprising finding awaits a proper theoretical interpretation.