Nonadditivity of Fluctuation-Induced Forces in Fluidized Granular Media
aa r X i v : . [ c ond - m a t . s o f t ] M a y Nonadditivity of Fluctuation-Induced Forces in Fluidized Granular Media
M. Reza Shaebani, ∗ Jalal Sarabadani, † and Dietrich E. Wolf Department of Theoretical Physics, University of Duisburg-Essen, 47048 Duisburg, Germany Department of Physics, University of Isfahan, Isfahan 81746, Iran (Dated: November 12, 2018)We investigate the effective long-range interactions between intruder particles immersed in a ran-domly driven granular fluid. The effective Casimir-like force between two intruders, induced bythe fluctuations of the hydrodynamic fields, can change its sign when varying the control parame-ters: the volume fraction, the distance between the intruders, and the restitution coefficient. Moreinterestingly, by inserting more intruders, we verify that the fluctuation-induced interaction is notpairwise additive. The simulation results are qualitatively consistent with the theoretical predictionsbased on mode coupling calculations. These results shed new light on the underlying mechanismsof collective behaviors in fluidized granular media.
PACS numbers: 45.70.Mg, 05.40.-a
Granular segregation has been extensively investigatedduring the last two decades aimed at revealing the under-lying complex dynamics [1, 2]. Besides the scientific in-terest, understanding the mechanisms of segregation is ofessential importance in geophysical [3] and industrial [4]processes. The behavior of granular mixtures, when me-chanically agitated, depends on a long list of grain, con-tainer, and external driving properties [2, 5]. The controlparameters can be tuned so that the demixing is initiated,reversed, or prevented [5–10]. While the phase behaviorof these systems is still a matter of debate, the nature ofparticle-particle interactions is known to play a crucialrole; two extreme limits can be distinguished: (i) the fullyfluidized regime where particles undergo only binary col-lisions, and (ii) the lasting contacts regime where durablefrictional contacts exist during a considerable part of theagitation cycle. While in the latter case the relevant pro-cesses are, e.g., reorganization, inertia, and convection[11], some studies reveal the existence of another mecha-nism in the fluidized regime: in the presence of intruderparticles, the hydrodynamic fields are modified especiallyin the inner regions between intruders, leading to effec-tive long-range interactions [8, 9, 12, 13]. Cattuto et al. [12] found that a pair of intruder particles experience aneffective force in a driven granular bed, originating fromthe modification of the pressure field fluctuations due tothe boundary conditions imposed by the intruders. SuchCasimir-like interactions are expected in thermal noisyenvironments confined by geometrical constraints [14].Most reports, so far, are about either binary mixtures[5, 6] or one or few intruder particles in a bed of smallerones [8, 10, 13]. An important question to address is howthe collective behavior is influenced by the number andarrangement of the intruders.In the present Letter, we study the effective interac-tions between immobile intruder particles immersed ina uniformly agitated granular fluid where all particlesundergo inelastic binary collisions (Fig. 1). We showthat the interaction between a pair of intruders exhibits a crossover from attraction to repulsion below a criticaldensity, as predicted in [12]. We here address the gen-eral conditions under which the transition happens, andpresent the phase diagram of the transition. Moreover,by comparing the behavior of two and multi intruder sys-tems, we find that the fluctuation-induced force is notderived from a pair-potential; inserting a new intruder af-fects the previously existing interactions in a non-trivialway, depending on the relative positions of the intruders.Such a feature together with the possible sign change ofthe forces make the multi-body interactions more com-plicated and may lead to a variety of collective behaviorssuch as segregation, clustering, or pattern formation. An-alytical calculations using the theory of randomly drivengranular fluids [15] confirm our findings.
Simulation method —
We consider a 2D granular fluidsimilar to the setup described in Refs.[12, 15, 16] bymeans of molecular dynamics simulations. We have areference system with two intruders A and B in which L /r =200, R /r =10, and D AB , /r =30 [see Fig. 1(a)]. (b)(a)(c) F A , B r / T ( × ) φ theorysimulation D AB =25D AB =35 α =0.7 α =0.9 FIG. 1: (color online) (a) Sketch of the simulation cell. (b) F A,B scaled by T /r vs φ . Comparison is made with the so-lution of Eq. (2) for the reference system (solid line), as wellas other values of the control parameters D AB and α (dashedlines). (c) Typical snapshots of dense (gray circles) and dilute(red crosses) states with φ =0 .
66 and 0 .
24, respectively. F A , B r / T ( × ) (a) D AB / r (b) T / T -1 -2 (c) L / r (× -2 ) F / h ( L ) L / r
FIG. 2: (a) F A,B vs D AB . The solid line is obtained via Eq. (2).(b) F A,B vs T . The dashed line indicates a linear relation.Values of F around 10 − T /r reflect the accuracy level of ourcalculation. (c) F A,B vs L . The dashed line corresponds to h ( L ). The inset shows that the deviation of the results from h ( L ) has no systematic dependence on L . φ =0 .
24 in all cases.
The reference volume fraction φ is 0 .
66, and the normalrestitution coefficient α is set to 0.8 for all collisions. Pe-riodic boundary conditions are applied in both directionsof the square-shaped cell to provide a spatially homoge-neous state. The system is coupled to an external heatbath that uniformly transfers energy into the system;the acceleration of each particle a i is perturbed instanta-neously by a random amount ξ i which can be consideredas a Gaussian white noise with zero mean and correla-tion h ξ ia ( t ) ξ jb ( t ′ ) i = ξ δ ij δ ab δ ( t − t ′ ), where a and b denoteCartesian components of the vectors, and ξ is the drivingstrength. The rate of the energy gain of a single particleaveraged over the uncorrelated noise source is ∂ t E = mξ [15]. Each particle also loses energy due to inelastic colli-sions at the mean-field rate of ∂ t E =( α − ωT / T is the granular temperature and ω is the collisionfrequency given by the Enskog theory [19]. Eventually,the system reaches a nonequilibrium stationary state bybalancing the energy input and the dissipation. Effective two-body interactions —
In the steady statewe measure the total force exerted by the granular fluidon each intruder along the x axis during the time inter-val τ ( τ ∼
250 collisions per particle). Due to the ob-served large fluctuations, the force is measured for morethan 10 consecutive time intervals τ . The probabil-ity distribution of the data is well fitted by a Gaussian[20] with the standard deviation σ =0 . T /r and thenonzero mean F = − . T /r , where T is the mean-fieldapproximation of the steady-state temperature deducedfrom the Enskog theory [15]. Using a similar analysisalong the y axis, we obtain zero force within the accu-racy of our measurements. F can be considered as themagnitude of the effective force F A,B that the intruder Bexerts on A, which is attractive in this case. We observethat, upon decreasing the volume fraction below a criti-cal value φ c ∼ .
57, the effective interaction F A,B becomesrepulsive [Fig. 1(b)], in agreement with the predictionof Ref. [12]. However, the transition is controlled notonly by φ , but also by D AB and α . One expects that,far from the transition region, increasing D AB decreasesthe magnitude of F A,B and it should eventually vanish at D AB = L/ ξ , it is shown that F is proportionalto the steady state temperature. Moreover, the resultsreveal the impact of dimensionality on the process: F increases slightly with L when we vary the system sizewhile other parameter values are kept fixed. The simu-lation results, shown in Fig. 2(c), can be well fitted by alogarithmic growth h ( L )= a ln( L )+ b (dashed line). Thisis contrary to what happens in three dimensional sys-tems, where the force is independent of L . Triple configurations and nonadditivity —
Next we ad-dress the interesting case of triple systems, where thethird intruder is located either on the x or y axis (Fig. 3).By choosing D AD /r =30 and φ =0 .
24, the effective force F A,BD exerted on particle A in the triple configuration(A,B,D) is compared to F A,B and F A,D obtained from thebinary systems (A,B) and (A,D), respectively. Note thatthe simulation is performed anew for each set of intrud-ers. Figure 3(a) shows that F A,BD is nearly the vectorsum of F A,B and F A,D . However, a comparison betweenthe sets (A,B,C), (A,B), and (A,C) in Figure 3(b) (with D AC /r =70) reveals that the force is definitely not pairwiseadditive in this case; F A,BC is even smaller than F A,B .In order to understand the mechanism behind thelong-range interactions and the transition, we draw at-tention to the fact that the fluctuating hydrodynamicfields, e.g. density [see Fig. 1(c)], are notably influencedby the geometric constraints, resulting in pressure im-balance around the intruders and effective interactionsbetween them. To establish a quantitative connectionbetween the effective force and the hydrodynamic fluc-tuations, we first employ mode coupling calculations[12, 15, 21] to evaluate the two-body interactions. Inthe nonequilibrium steady state, the hydrodynamic fields( p ( r ), T ( r ), n ( r )) fluctuate around their stationary val-ues ( p s , T s , n s ). The average pressure fluctuation p f ( r )in the presence of the boundary conditions imposed byintruders behaves analogously to the Casimir effect, i.e., p f ( r ) in the hatched region of Fig. 1(a) differs from thatof the cross-hatched region. Using the Verlet-Levesque ( a ) ( b )( c ) B D F =5.2 0.2
A,B
F =4.6 0.2
A,BC
F =4.9
A,BC ~ F = . . A , B D F = . . A , D F =1.3 0.2
A,C
F =5.2 0.2
A,B
A B C F / F A , B e D AC / r e F A,B + F
A,C F A,BC ...
FIG. 3: (color online) Comparison between the binary andtriple effective forces (scaled by 0 . T /r ) exerted on particleA in the presence of particles (a) B and D (b) B and C. (c) e , F A,BC and F A,B + F A,C versus D AC , the x position of particle C. attractionrepulsion α D AB /r0.30.40.50.6 φ AB /r0.10.30.50.7 φ α =0.7 α =0.8 α =0.9 9.5 -9.50.1 0.3 0.5 0.7 φ α D=25D=30D=50 1.0 -2.030 50 70 D AB /r0.70.80.9 α φ =0.50 φ =0.54 φ =0.56 FIG. 4: (color online) (a) Schematic phase diagram of the transition in the ( φ , D AB , α ) space. The surface corresponds to F A,B =0. (b-d) 2D profiles of the phase diagram where the solid lines mark the interface position. The dashed lines correspond tothe interface position for other values of the third control parameter. The color intensity reflects the magnitude of F (scaled by0 . T /r ), with blue (red) meaning repulsion (attraction). (e) F on particle A vs n , the length of the chain. The results of themode coupling, simulation, and pairwise summation of two-body forces are shown with solid line, full circles, and open circles,respectively. Inset: Mode coupling results (solid line) and their corresponding pairwise summations (dashed line) ( L =500 r ). equation of state for a hard disks system p ( n, T )= T f ( n )(with f ( n )= n (1 + φ / / (1 − φ ) , and n the number den-sity) [22], we expand the pressure up to second orderaround ( n s , T s ), and take the statistical average overthe random noise source: p f ( r )= f ′ ( n ) | ns h δn ( r ) δT ( r ) i + T s f ′′ ( n ) | ns h ( δn ( r )) i . By Fourier transforming δn ( r )and δT ( r ) one obtains p f ( r )= Z (cid:20) f ′ ( n ) | ns S nT ( k )+ 12 T s f ′′ ( n ) | ns S nn ( k ) (cid:21) d k , (1)where the integral is taken over the k vectors allowedat position r by the boundary conditions, and S ab ( k ) isthe pair structure factor defined as V − h δa ( k ) δb ( − k ) i .The detailed description of the structure factor calcu-lations will be reported elsewhere (see also Ref. [15]).Here we denote the integrand of Eq. (1) with g ( k , φ, α )and compare p f ( r ) for two surface points located on op-posite sides of intruder A with the same y coordinates.The related k vectors in the x direction are confined to D in ( y )= D AB − p ( R − y ) and D out ( y )= L − D AB − p ( R − y ) in the cross-hatched and hatched regions, respec-tively; Therefore the pressure difference between thesetwo points ∆ p ( y )= p (in) f − p (out) f has y dependence. Byintegrating over y , we arrive at the average pressure dif-ference between the gap and outside region:∆ p = Z R − R dy "Z π/r ∗ π/D in ( y ) dk x − Z π/r ∗ π/D out ( y ) dk x π/r ∗ π/L dk y g ( k x , k y , φ, α )2 R . (2)To ensure that the hydrodynamic description is valid,the integrals are restricted to the long wavelength range r ∗ =max(2 r, l ∗ ) (with l ∗ being the mean free path) andonly small inelasticities are considered. Moreover, D AB isalways chosen large enough (3 R ≤ D AB ) so that the short-range depletion forces [23] do not play a role. The ef-fective force F AB is calculated via Eq. (2) for differentvalues of φ or D AB and compared to the simulation re-sults in Figs. 1(b) and 2(a). The forces are of the sameorder of magnitude as those obtained from the simula-tions. The deviations can be attributed to the fact that the hydrodynamic fluctuations are correlated in the gapand outside regions. The correction due to this effect isproportional to ∂ g/∂φ which always has a sign oppositeto that of g , thus, Eq. (2) overestimates the magnitudeof the force. Using the mode coupling calculations, it isshown in Fig. 1(b) that the transition point is sensitive tothe choice of D AB and α . The set of control parametersfor sign switching of the Casimir force crucially dependson the physics of the system (see e.g. [24]). Figure 4 sum-marizes the calculations in a phase diagram in the ( φ , α , D AB ) space, which appears to be in remarkable accordwith the dynamical model. The phase diagram is notinfluenced by the choice of the steady-state temperature T s , while the magnitude of the force grows linearly with T s as expected for Casimir forces in thermal fluctuatingmedia [14]. Regarding the fact, that the leading termof g ( k , φ, α ) at small k is proportional to 1 /k [15], onealso finds from Eq. (1) that p f , and therefore the force,in the thermodynamic limit behaves as 1 /r in 3D whilediverges logarithmically as ln( L/r ) in 2D.
Multi-body effects —
In the triple system of Fig. 3(a),loosely speaking, because of the independence of k vec-tors in x and y directions one expects that the vector sumholds. This is in agreement with the simulation resultsin Fig. 3(a), neglecting the deviation due to hydrody-namic correlations. In the configuration of Fig. 3(b), theeffective force exerted on particle A results from the dif-ference between the range of available k modes on its leftand right sides. In the presence of particle C, the range ofallowed k modes decreases on the left side due to periodicboundary conditions, which causes a lowering of pressuredifference between both sides of A. The effective force F A,BC is thus smaller, compared to the binary interaction F A,B . The calculated force on particle A using Eq. (2) isshown with e F A,BC in Fig. 3(b); the agreement is satisfac-tory. We introduce a measure e = | F A,B + F A,C − F A,BC | / | F A,B | for quantifying the deviation from the case that the in-teractions are derived from a pair-potential. Figure 3(c)shows how e , obtained from the mode coupling calcula-tions, behaves when the position of particle C is varied FIG. 5: The trajectories of a test intruder subject to the effec-tive force field in a fixed square configuration of intruders, ob-tained from the mode coupling calculations for φ =0 .
24 (left),and 0 .
54 (right). Short range interactions are excluded (whitezones). Lighter colors mean stronger forces. The dashed linein the right figure indicates the border of an attractive sub-domain, where the trajectories of the test particle end up onthe fixed intruder. on the x axis. Figure 4(e) indicates that our analyticalapproach also provides reasonable estimates of the effec-tive force acting on particle A in a chain configuration,while the deviation from pairwise summation of two-bodyforces grows with the number of intruders, n .Finally, we investigate more complicated geometries bycalculating the effective force field of a fixed square struc-ture acting on a test intruder particle. To elucidate theimpact of sign switching on the force pattern we comparean intermediate density (around the transition zone) witha low density case. Figure 5 shows that the patterns areclearly different e.g. in terms of the number of equilibriumpoints. There exist also attractive subdomains of lengthscale ℓ around the fixed intruders in the right figure (e.g.the domain surrounded by the dashed line). Dependingon the choice of control parameters, ℓ ranges between 0(entirely repulsive patterns at low densities, as shown inthe left figure) and L (entirely attractive patterns at highdensities, not shown) leading to different types of collec-tive behavior such as segregation, clustering, and patternformation. Putting aside the pairwise summations, wecompare the mode coupling predictions with simulationresults at four selected points in Fig 5. The predicted TABLE I: Comparison between the theoretical (T) and sim-ulation (S) results for the selected points in Figs. 5(left) and5(right). Same force scales as Fig. 3. The maximum error barsfor the simulated F and θ are ± . ± ◦ , respectively. | ~F | ( A ) | ~F | ( B ) | ~F | ( C ) | ~F | ( D ) θ ( A ) θ ( B ) θ ( C ) θ ( D ) left-S 1 . . . . ◦ ◦ ◦ ◦ left-T 1 .
17 2 .
52 3 .
97 4 .
32 7 ◦ ◦ ◦ ◦ right-S 0 . . . . ◦ ◦ − ◦ − right-T 0 .
24 0 .
26 0 .
25 0 .
02 28 ◦ ◦ − ◦ ◦ forces are of the same order of magnitude as the simula-tion results (see table I), however, the errors of the forcesize | F | and its direction θ reach up to 30% and 15 ◦ , re-spectively. The differences are smaller at large distancesand also low densities. Indeed, the hydrodynamic corre-lations become more important in multi-body cases and,hence, one should include triplet or higher order struc-ture factors [25] in mode coupling calculations to properlytake the multi-body effects into account.In conclusion, we focus on the problem of long-rangefluctuation-induced forces between intruder particles im-mersed in an agitated fluid bed. The sign of the forcecan be reversed by tuning the control parameters. Fur-thermore, the multi-body interactions do not follow fromtwo-body force descriptions. A newly inserted intruder,depending on its position, may affect the pressure balancearound the other intruders. This suggests that the effec-tive force is not derived from a pair-potential in agree-ment with our simulation results. Our findings representa step forward in understanding the origin of collectivebehaviors in fluidized granular mixtures.We would like to thank I. Goldhirsch for helpful dis-cussions and J. T¨or¨ok and B. Farnudi for comments onthe manuscript. Computing time was provided by John-von-Neumann Institute of Computing (NIC) in J¨ulich. ∗ Electronic address: [email protected] † present address: Max Planck Institute for Polymer Re-search, D-55128 Mainz, Germany.[1] H. M. Jaeger, S. R. Nagel, and R. P. Behringer, Rev.Mod. Phys. , 1259 (1996).[2] A. Kudrolli, Rep. Prog. Phys. , 209 (2004).[3] O. Pouliquen, J. Delour, and S. B. Savage, Nature ,816 (1997); L. Hsu, W. E. Dietrich, and L. S. Sklar, J.Geophys. Res. , F02001 (2008).[4] J. C. Williams, Powder Technol. , 245 (1976); J. M.Ottino and D. V. Khakhar, Annu. Rev. Fluid Mech. ,55 (2000).[5] M. P. Ciamarra, M. D. De Vizia, A. Fierro, M. Tarzia, A.Coniglio, and M. Nicodemi, Phys. Rev. Lett. , 058001(2006).[6] D. C. Hong, P. V. Quinn, and S. Luding, Phys. Rev.Lett. , 3423 (2001); M. Tarzia, A. Fierro, M. Nicodemi,and A. Coniglio, Phys. Rev. Lett. , 198002 (2004); M.Tarzia, A. Fierro, M. Nicodemi, M. P. Ciamarra, and A.Coniglio, Phys. Rev. Lett. , 078001 (2005); K. M. Hilland Y. Fan, Phys. Rev. Lett. , 088001 (2008).[7] T. Shinbrot, Nature , 352 (2004).[8] D. A. Sanders, M. R. Swift, R. M. Bowley, and P. J.King, Phys. Rev. Lett. , 208002 (2004).[9] M. P. Ciamarra, A. Coniglio, and M. Nicodemi, Phys.Rev. Lett. , 038001 (2006); I. Zuriguel, J. F. Boudet,Y. Amarouchene, and H. Kellay, Phys. Rev. Lett. ,258002 (2005).[10] T. Schnautz, R. Brito, C. A. Kruelle, and I. Rehberg,Phys. Rev. Lett. , 028001 (2005).[11] E. Caglioti, A. Coniglio, H. J. Herrmann, V. Loreto, and M. Nicodemi, Europhys. Lett. , 591 (1998); T. Shin-brot and F. J. Muzzio, Phys. Rev. Lett. , 4365 (1998);T. Mullin, Phys. Rev. Lett. , 4741 (2000); G. Metcalfe,S. G. K. Tennakoon, L. Kondic, D. G. Schaeffer, and R.P. Behringer, Phys. Rev. E , 031302 (2002).[12] C. Cattuto, R. Brito, U. M. B. Marconi, F. Nori, and R.Soto, Phys. Rev. Lett. , 178001 (2006).[13] S. Aumaitre, C. A. Kruelle, and I. Rehberg, Phys. Rev.E , 041305 (2001).[14] M. Kardar and R. Golestanian, Rev. Mod. Phys. , 1233(1999).[15] T. P. C. van Noije, M. H. Ernst, E. Trizac, and I. Pago-nabarraga, Phys. Rev. E , 4326 (1999).[16] G. Peng and T. Ohta, Phys. Rev. E , 4737 (1998).[17] I. Goldhirsch and G. Zanetti, Phys. Rev. Lett. , 1619(1993).[18] T. P. C. van Noije, M. H. Ernst, and R. Brito, Phys. Rev.E , R4891 (1998). [19] S. Chapman and T. G. Cowling, The Mathematical The-ory of Non-uniform Gases (Cambridge University Press,Cambridge, 1970).[20] D. Bartolo, A. Ajdari, J. B. Fournier, and R. Golestanian,Phys. Rev. Lett. , 230601 (2002).[21] R. Brito and M. H. Ernst, Europhys. Lett. , 497 (1998).[22] L. Verlet and D. Levesque, Mol. Phys. , 969 (1982).[23] R. Roth, R. Evans, and S. Dietrich, Phys. Rev. E ,5360 (2000); C. N. Likos, Phys. Rep. , 267 (2001).[24] C. Hertlein, L. Helden, A. Gambassi, S. Dietrich, andC. Bechinger, Nature , 172 (2008); M. Levin, A. P.McCauley, A. W. Rodriguez, M. T. Homer Reid, and S.G. Johnson, Phys. Rev. Lett. , 090403 (2010).[25] P. Attard, J. Chem. Phys. , 3072 (1989); S. Jorge, E.Lomba, and J. L. F. Abascal, J. Chem. Phys.116