Master equation for the probability distribution functions of overlaps between particles in two dimensional granular packings
Abstract
How does the force chain network in a random granular material react to hydrostatic compression? We show that not only contacts, but also their opening and closing as well as interparticle gaps, i.e. virtual contacts, must be included for a comprehensive description of the system response involving the probability distribution functions (PDFs) of the extended force network. Considering overlaps/forces as stochastic variables, their conditional probability distributions (CPDs) are (numerically) determined and allow to formulate a Master equation for the PDFs. The insight one gets from the CPDs is striking: The mean change of contacts reflects non-affinity, while their fluctuations obey uncorrelated Gaussian statistics. In contrast, virtual contacts are mostly affine to isotropic (de)compression in average, but show multi-scale correlations with considerable probability of large “jumps” according to a stable distribution (cf. Lévy distribution), that allows for a generalized central limit theorem. Noting that all changes of the network during compression are similarly scaled by the distance from the jamming density, the evolution of the system is fully described by the Master equation.
pacs:
45.70.Cc,46.65.+g,61.43.-jThe statistical physics of soft matter, due to disorder, glassy dynamics, dissipation, etc. poses many difficulties for research and has wide practical applications. As one example, jammed granular assemblies are inherently out of equilibrium due to the dissipation of energy and the absence of temperature Jaeger et al. (1996). Our results on their enhanced statistical description will have wider consequences related to a generalization of the central limit theorem, which is relevant for completely different fields like self-similar models of, e.g. financial markets, photons in turbid media, or in spectroscopy, cf. pressure-broadened spectral lines ^{1}^{1}1 See http://en.wikipedia.org/wiki/Stable_distribution. Wikipedia: The Free Encyclopedia. Wikimedia Foundation, Inc. . Despite the non-equilibrium nature of granular materials, static properties of them around the rigidity transition density have been widely investigated as well as their relation to critical phenomena O’Hern et al. (2003); Majmudar et al. (2007); van Hecke (2010); Goodrich et al. (2012). It is challenging to explain these observations on static granular packings by statistical mechanics, e.g. based on the Edwards ensemble Edwards and Oakeshott (1989) or force ensemble Henkes and Chakraborty (2005); Tighe et al. (2010), where all configurations of particles are equally probable under a constant volume or force and torque balances of particles, respectively. Further studies Blumenfeld and Edwards (2003); Wang et al. (2012) and comparison Puckett and Daniels (2013) of these ideas lead us to the exploration into the statistical weight, i.e. the probability distribution functions (PDFs), of forces. Many theoretical investigations Metzger (2004); Henkes et al. (2007); Tighe et al. (2008) have been devoted to determine the functional forms of the PDFs obtained through experiments Erikson et al. (2002); Corwin et al. (2005); Desmond et al. (2013) and numerical simulations Silbert et al. (2002); Landry et al. (2003); Müller et al. (2010). There is, however, still much debate about the tails of the PDFs van Eerd et al. (2007) and how the closing/opening of contacts affects the PDFs at small forces Wyart (2012); Lerner et al. (2013); Charbonneau et al. (2012). Our novel probabilistic description links stochastic processes with the generalized central limit theorem and will allow for a micro-mechanic and statistical physics based macroscopic models of soft and granular matter.
In this Letter, we take a different approach from the previous works based on the mechanical equilibrium of static packings Edwards and Oakeshott (1989); Henkes and Chakraborty (2005); Tighe et al. (2010); Blumenfeld and Edwards (2003); Wang et al. (2012); Puckett and Daniels (2013); Metzger (2004); Henkes et al. (2007); Tighe et al. (2008). We regard overlaps between particles as stochastic variables and measure their conditional probability distributions (CPDs) by molecular dynamics (MD) simulations under isotropic (de)compression. The CPDs show striking features of not only contacts, but also closing/opening of contacts as well as interparticle gaps, i.e. virtual contacts. We introduce the Master equation for the PDFs of overlaps, where good agreement between its solutions and the MD simulations is established, as long as the increment of the area fraction is much smaller than the distance from the jamming point . Interestingly, the solutions are independent of history (the Markov property) and reversible if an effective increment , which governs the amplitude of non-affine response, is small enough. Our stochastic approach opens the door to a new physical description of static granular packings.
We use MD simulations of two-dimensional frictionless granular particles. The normal force between particles in contact is given by () with a spring constant , viscosity coefficient , overlap and relative speed in normal direction. A global damping force , proportional to the particle’s velocity , is also introduced to enhance the relaxation. We prepare static packings of the particles by a method similar to the one used in Ref. Ellenbroek et al. (2009). At first, we randomly distribute a binary mixture of particles with different radii, and (), and the same masses in a square periodic box, where each pair of particles is not in contact. Then, we adjust the size of the particles by rescaling every radius as until a target mean overlap is obtained, where , , and are time, an increment of time, and the mean overlap over all particles at time , respectively, and the initial value corresponds to . Here, we keep the mass constant and use a long length scale to grow the particles gently ^{2}^{2}2 The static packings prepared with longer length scales and give the same results concerning critical scaling of frictionless particles near the jamming point O’Hern et al. (2003); Majmudar et al. (2007); van Hecke (2010), while we cannot obtain the same results with a smaller value . . The particles lose their kinetic energies by means of inelastic contacts and global damping. We stop rescaling each radius when the accelerations of all particles drops below a small threshold and assume that the system is static Ellenbroek et al. (2009). We prepare samples for smaller systems () and samples for the largest system () by changing the initial random configuration of particles. However, we only report the results of and note that all results do not depend on the system size. In our simulations, the jamming point obeys the scaling van Hecke (2010), where we find , , and with the mean radius . Thus, the mean overlap is equivalent to the distance from the jamming point .
We apply isotropic compression to each packing by rescaling every radius by so that the area fraction changes from to . If an affine response of the packing is assumed Digby (1981), an overlap or interparticle gap , defined as the difference between the sum of the radii and the interparticle distance , changes to ^{3}^{3}3We neglect the higher order term .. However, the particles are randomly arranged and the force balance of each particle is broken by compression. Thus, the system relaxes to a new static state and the overlap or gap changes to a new value after relaxation, i.e. due to the non-affine response of overlaps. We introduce the Delaunay triangulation (DT) of static packings as sketched in Fig. 1, where the particles in contact are connected by red solid lines and the nearest neighbors without contact, i.e. a virtual contact, are connected by blue dashed lines. Since the DT is unique for each packing, virtual contacts are uniquely determined and the total number of contacts and virtual contacts are finite. Figures 1(a) and (b) show the DTs before compression and after compression and relaxation, respectively, where the overlaps between virtual contacts are defined as negative values and all kinds of transitions of overlaps are displayed. The overlaps and change to and , respectively; they do not change their signs and thus contacts are neither generated nor broken. We name these transitions contact-to-contact (CC) and virtual-to-virtual (VV), respectively. On the other hand, and change to and , respectively; a new contact is generated and an existing contact is broken, respectively. We name these cases virtual-to-contact (VC) and contact-to-virtual (CV), respectively ^{4}^{4}4In our simulations, the total number of contacts and virtual contacts are independent of the area fraction and always conserved under both compression and decompression. Although we do not take into account flips of Delaunay edges by large displacements of particles, we have not observed any flips if and the number of flipped edges are less than at most for . .
To characterize all transitions of overlaps, we plot them on scatter plots as shown in Figs. 2(i) and (ii), where and are scaled overlaps before compression and after relaxation, respectively. We also plot the affine response,
(1) |
with the effective increment ^{5}^{5}5 can be large, whereas is always small. and a coefficient , where is proportional to with an offset linearly dependent on . In these figures, the blue and red dots represent and , respectively. Though the overlaps after relaxation are not affine, the deviation of from is quite small if the increment is small and the system is far from the jamming point, i.e. for (Fig. 2(i)). However, if we increase or decrease , deviates visibly from the deterministic prediction and data points are more dispersed (Fig. 2(ii)). Figure 2(iii) is a sketch of the differences between affine and non-affine responses, where the mean values in (CC) , and (VV) , are described by linear fitting functions for
(2) | |||||
(3) |
respectively. The systematic deviation from is represented by the differences in slopes ( and ) and the offsets ( and ). We also introduce standard deviations of from and as and , respectively, which are almost independent of the initial values . Figure 2(iv) displays a double logarithmic plot of against , where all data collapse onto with . We also find , , , , and with , , , and , respectively, within the fitting range . Surprisingly, all parameters characterizing not only the mean values, but also the fluctuations are linearly scaled by Brunet et al. (2008). Because is comparable to the mean of (), virtual contacts behave affine in average except for their huge fluctuations around the mean value. In contrast, is always smaller than and intersects at independent of . This leads to small responses () of small overlaps (), which implies preferred tangential and hindered normal displacements as a sign of non-affine deformations Ellenbroek et al. (2009); Kruyt and Luding (2010). The scattered data in (VC) and (CV) are concentrated in the narrow regions (the inside of the dashed lines in Fig. 2(iii)), while linearly increases with in (VC) and there is no data in (CV).
Regarding the overlaps and as two stochastic variables, we introduce a conditional probability distribution (CPD) of the overlaps, , as
(4) |
satisfying a normalization condition van Kampen (2007), where and are the PDFs of overlaps before and after compression, respectively; gives a probability distribution of , which was before compression, and can be defined in (CC), (VV), (VC) and (CV). Figure 3(a) shows the CPDs in (CC), where all results with a wide range of collapse if we scale and by and , respectively. The solid line is a Gaussian distribution function
(5) |
with . Figure 3(b) displays the CPDs in (VV), where all results collapse as well, after the same scaling as for (CC). The solid line is here the stable distribution function Voit (2005)
(6) |
with and a dimensionless wave number , where the fitting parameters are given by and , i.e. the CPDs in (VV) are nearly Holtsmark distributions (). Figures 3(c) and (d) show the CPDs in (VC) and (CV) approximated by
(7) | |||||
(8) |
with and dimensionless length scales, and , respectively. The terms in parentheses on the right hand sides are required to satisfy the normalization conditions ^{6}^{6}6The normalization condition is given by . and well describe the dependence on (data are not shown). Note that the CPDs in (CC) and (VV) (Eqs. (5) and (6)) converge to a delta-function and those in (VC) and (CV) (Eqs. (7) and (8)) vanish in the limit of ^{7}^{7}7We rewrite and take , where we used . Because of , the CPDs in (VC) and (CV) vanish in the limit of ..
Now, we restrict to quite small values compared to and define an infinitesimal effective increment . Introducing a transition rate as , we rewrite Eq. (4) as the Master equation van Kampen (2007),
(9) |
Figures 4(a) and (b) display the numerical solutions of the Master equation under compression, where the CPDs are given by Eqs. (5)-(8). The increment of the area fraction is fixed to so that throughout the numerical integrations. Here, the initial condition is given by the PDF obtained through the MD simulations at with the area fraction before compression . The overlaps are scaled by the mean overlap at the initial state . Good agreement between the solutions (the red solid lines) and the PDFs from simulations (the open symbols) is established for small and thus we can confirm the Markov property of overlaps. Exchanging for , we can obtain the CPDs for decompression from to . The CPDs in (CC) and (VV) are given by replacing the linear fitting functions, and , with and , respectively, and the standard deviations, and , with and , respectively. The CPDs for decompression in (VC) and (CV) are given by those for compression in (CV) and (VC), respectively. Figures 4(c) and (d) show the numerical solutions of the Master equation under decompression with so that , where the initial condition is given by the PDF at with the area fraction before decompression . From these results, we can confirm the reversibility in the sense that the Master equation well describes the PDFs in both directions for small (even though the initial state before compression and the final state after decompression are not exactly the same). Note that a similar criterion for the reversibility is found in an experiment of sound propagation in sands Brunet et al. (2008).
In this study, we provide for the first time the Master equation for the PDFs of overlaps between particles. By considering the extended contact network, the Master equation is able to fully predict the evolution of the PDFs. The novel formulation has a huge impact from the application point of view as it provides a complete statistical description of the stress field in a granular assembly, basis for a macroscopic, continuum approach to large scale problems. Both the Markov property and reversibility are well retained only if the applied strain increment is infinitesimally small ^{8}^{8}8K. Saitoh, V. Magnanimo, and S. Luding, in preparation.. That is, by including contacts, virtual contacts and their mutual exchange into the description, the behavior is independent of history and the actual stress state is sufficient to predict the subsequent incremental response.
The CPDs show by themselves important properties. While the mean change of overlap responds in a non-affine way, astonishingly the fluctuations obey Gaussian statistics, indicating an independent stochastic evolution of overlaps/forces as reported, e.g. in Ref. Kruyt and Rothenburg (2002). In contrast, the CPDs in (VV), where no forces act on the virtual contacts, are nearly the Holtsmark distributions with much broader tails than the Gaussian, i.e. much larger jumps of negative overlaps can occur. This indicates a strongly correlated stochastic evolution over a wide range of length-scales. The probability for opening or closing contact transitions is exponentially decaying with distance from the zero, and causes the jump/gap in the complete PDFs, since just closed contacts are affected by repulsion (narrow exponential) and just opening contacts are free to open widely (wide exponential). Because both the Gaussian and Holtsmark distributions are members of the stable distribution family, fluctuations of contacts and virtual contacts in granular materials should obey the generalized central limit theorem Voit (2005), which has consequences for the statistical description of disordered systems in general. The strong deviation from an affine approximation Ellenbroek et al. (2009) for contacts and the enormous fluctuations of overlaps Henkes and Chakraborty (2005) for virtual contacts, as well as the opening and closing of contacts statistics, are all captured by linear scalings with the effective increment.
Clearly, there is the need of further studies on the physical origin of the surprising non-affine behavior and statistics of overlaps described above. The origin of the functional forms of the CPDs can give very interesting insights to the mechanics of granular materials, e.g. stochastic processes of overlaps. Now, analytic solutions or asymptotic solutions of the Master equation are important possible steps towards the understanding of the functional forms of the PDFs. The Master equation, however, requires the increment to be much smaller than . Thus, it can never reach , and the initial condition cannot be the PDF at . This means that the jamming transition is a singular limit of the Master equation. Our analysis can be easily extended to three dimensions and be examined and validated by experiments, e.g. by photoelastic tests Majmudar and Behringer (2005) or oedometer test of sands ^{9}^{9}9 Y.-H. Wang and Y. Gao, Granular Matter. DOI 10.1007/s10035-013-0457-1. .
We thank M. van Hecke, B. Tighe, H. Hayakawa, S. Yukawa, T. Hatano, H. Yoshino, and K. Kanazawa for fruitful discussions. This work was financially supported by the NWO-STW VICI grant 10828 and a part of numerical computation has been carried out at the Yukawa Institute Computer Facility, Kyoto, Japan.
References
- Jaeger et al. (1996) H. Jaeger, S. R. Nagel, and R. P. Behringer, Rev. Mod. Phys. 68, 1259 (1996).
- (2) See http://en.wikipedia.org/wiki/Stable_distribution. Wikipedia: The Free Encyclopedia. Wikimedia Foundation, Inc.
- O’Hern et al. (2003) C. S. O’Hern, L. E. Silbert, A. J. Liu, and S. R. Nagel, Phys. Rev. E 68, 011306 (2003).
- Majmudar et al. (2007) T. S. Majmudar, M. Sperl, S. Luding, and R. P. Behringer, Phys. Rev. Lett. 98, 058001 (2007).
- van Hecke (2010) M. van Hecke, J. Phys.: Condens. Matter 22, 033101 (2010).
- Goodrich et al. (2012) C. P. Goodrich, A. J. Liu, and S. R. Nagel, Phys. Rev. Lett. 109, 095704 (2012).
- Edwards and Oakeshott (1989) S. F. Edwards and R. B. S. Oakeshott, Physica A 157, 1080 (1989).
- Henkes and Chakraborty (2005) S. Henkes and B. Chakraborty, Phys. Rev. Lett. 95, 198002 (2005).
- Tighe et al. (2010) B. P. Tighe, J. H. Snoeijer, T. J. H. Vlugt, and M. van Hecke, Soft Matter 6, 2908 (2010).
- Blumenfeld and Edwards (2003) R. Blumenfeld and S. F. Edwards, Phys. Rev. Lett. 90, 114303 (2003).
- Wang et al. (2012) K. Wang, C. Song, P. Wang, and H. A. Makse, Phys. Rev. E 86, 011305 (2012).
- Puckett and Daniels (2013) J. G. Puckett and K. E. Daniels, Phys. Rev. Lett. 110, 058001 (2013).
- Metzger (2004) P. T. Metzger, Phys. Rev. E 70, 051303 (2004).
- Henkes et al. (2007) S. Henkes, C. S. O’Hern, and B. Chakraborty, Phys. Rev. Lett. 99, 038002 (2007).
- Tighe et al. (2008) B. P. Tighe, A. R. T. van Eerd, and T. J. H. Vlugt, Phys. Rev. Lett. 100, 238001 (2008).
- Erikson et al. (2002) J. M. Erikson, N. W. Mueggenburg, H. M. Jaeger, and S. R. Nagel, Phys. Rev. E 66, 040301(R) (2002).
- Corwin et al. (2005) E. I. Corwin, H. M. Jaeger, and S. R. Nagel, Nature 435, 1075 (2005).
- Desmond et al. (2013) K. W. Desmond, P. J. Young, D. Chen, and E. R. Weeks, Soft Matter 9, 3424 (2013).
- Silbert et al. (2002) L. E. Silbert, G. S. Grest, and J. W. Landry, Phys. Rev. E 66, 061303 (2002).
- Landry et al. (2003) J. W. Landry, G. S. Grest, L. E. Silbert, and S. J. Plimpton, Phys. Rev. E 67, 041303 (2003).
- Müller et al. (2010) M.-K. Müller, S. Luding, and T. Pöschel, Chem. Phys. 375, 600 (2010).
- van Eerd et al. (2007) A. R. T. van Eerd, W. G. Ellenbroek, M. van Hecke, J. H. Snoeijer, and T. J. H. Vlugt, Phys. Rev. E 75, 060302(R) (2007).
- Wyart (2012) M. Wyart, Phys. Rev. Lett. 109, 125502 (2012).
- Lerner et al. (2013) E. Lerner, G. Düring, and M. Wyart, Soft Matter 9, 8252 (2013).
- Charbonneau et al. (2012) P. Charbonneau, E. I. Corwin, G. Parisi, and F. Zamponi, Phys. Rev. Lett. 109, 205501 (2012).
- Ellenbroek et al. (2009) W. G. Ellenbroek, M. van Hecke, and W. van Saarloos, Phys. Rev. E 80, 061307 (2009).
- (27) The static packings prepared with longer length scales and give the same results concerning critical scaling of frictionless particles near the jamming point O’Hern et al. (2003); Majmudar et al. (2007); van Hecke (2010), while we cannot obtain the same results with a smaller value .
- Digby (1981) P. J. Digby, J. Appl. Mech. 48, 803 (1981).
- (29) We neglect the higher order term .
- (30) In our simulations, the total number of contacts and virtual contacts are independent of the area fraction and always conserved under both compression and decompression. Although we do not take into account flips of Delaunay edges by large displacements of particles, we have not observed any flips if and the number of flipped edges are less than at most for .
- (31) can be large, whereas is always small.
- Brunet et al. (2008) T. Brunet, X. Jia, and P. A. Johnson, Geophys. Res. Lett. 35, L19308 (2008).
- Kruyt and Luding (2010) O. D. N. P. Kruyt and S. Luding, Int. J. Sol. Struct. 47, 2234 (2010).
- van Kampen (2007) N. G. van Kampen, Stochastic Processes in Physics and Chemistry, 3rd edition (Elsevier B. V. Amsterdam, The Netherlands, 2007).
- Voit (2005) J. Voit, The Statistical Mechanics of Financial Markets, 3rd Edition (Springer-Verlag, Berlin, 2005).
- (36) The normalization condition is given by .
- (37) We rewrite and take , where we used . Because of , the CPDs in (VC) and (CV) vanish in the limit of .
- (38) K. Saitoh, V. Magnanimo, and S. Luding, in preparation.
- Kruyt and Rothenburg (2002) N. P. Kruyt and L. Rothenburg, Int. J. Solids Struct. 39, 571 (2002).
- Majmudar and Behringer (2005) T. S. Majmudar and R. P. Behringer, Nature 435, 1079 (2005).
- (41) Y.-H. Wang and Y. Gao, Granular Matter. DOI 10.1007/s10035-013-0457-1.