Evolution of force networks in dense granular matter close to jamming

When dense granular systems are exposed to external forcing, they evolve on the time scale that is typically related to the externally imposed one (shear or compression rate, for example). This evolution could be characterized by observing temporal evolution of contact networks. However, it is not immediately clear whether the force networks, defined on contact networks by considering force interactions between the particles, evolve on a similar time scale. To analyze the evolution of these networks, we carry out discrete element simulations of a system of soft frictional disks exposed to compression that leads to jamming. By using the tools of computational topology, we show that close to jamming transition, the force networks evolve on the time scale which is much faster than the externally imposed one. The presentation will discuss the factors that determine this fast time scale.


Introduction
The interaction between particles that build dense granular matter (DGM) is characterized by the presence of a force network, mesoscopic structure that spontaneously form as a granular system is exposed to shear or compression. These networks are crucial for understanding the static and dynamic properties of DGM, and have been recently studied extensively, using the tools that vary from statistical methods [1][2][3] to network type of analysis [4][5][6] and application of topological methods [7][8][9][10][11][12].
Force networks exhibit complicated time dependent structures. Earlier approaches, used to analyze these structures, usually involved rather arbitrary choices of the force thresholds and often depended on ad hoc descriptions of the structures. To avoid these problems we define the force network [11] as a collection of sets {X θ } θ∈R such that for θ ≥ 0 the set X θ is the part of the contact network on which the force interactions between the particles exceed the value θ. In the series of papers [9][10][11][12], we make use of an algebraic topological technique that rigorously describes geometric structures of the sets {X θ } θ∈R over all force levels θ. This approach provides a quantitative description of the force networks. Moreover, it allows us to define a concept of distance between force networks, so we can quantitatively compare them.
In [12] we focused on comparing force networks in simulated granular systems built from soft, possibly inelastic and/or frictional particles exposed to steady slow compression. In that work we found that force networks' evolution was characterized by large distances between consecutively sampled states of the systems that were approaching jamming transition from below, and the evolution was much smoother for jammed systems. In particular close to jamming transition, we found rather dramatic evolution of the force networks in all considered systems (frictional or not, with elastic or inelastic particles).
The following unexpected findings obtained in [12] motivate the present peper. Distances between force networks at consecutive states of the system sampled at the rate used in [12] are comparable to the distances between force networks at the same packing fraction obtained from different realizations, see Figure 1. Sec. 2 provides a brief description of the plotted quantities while more in depth discussion of the figure is presented in Sec. 3. This result suggests that force networks evolve on a time scale that is much faster than the inverse sampling rate and therefore, the main topic of the present paper is the analysis of the system on extremely fast time scales. Furthermore, we discuss how the time scale characterizing evolution of force networks depend on two externally imposed time scales: the slow one imposed by external driving (compression in the case considered here) and the fast one imposed by the temporal discretization used in our simulations. We show that variation of these two scales provides further insight into the mechanisms governing force network evolution.

Methods
Here we give a brief overview of the techniques used; the reader is referred to [12] for more detailed description of the protocol and definitions of various quantities.
We consider a square domain in two spatial dimensions exposed to slow compression by the moving walls. The particles are polydisperse, soft, inelastic disks interacting by normal and tangential forces, with the latter ones modeled by Cundall-Strack type of interaction [13]. All quantities are expressed using the average particle diameter, d, as the length scale, the binary particle collision time τ c = π d/(2gk n ) as the time scale, and the average particle mass, m, as the mass scale. Here k n (in units of mg/d) is the spring constant used for modeling (linear) repulsion force between colliding particles, and g is gravity (used only for the purpose of scaling). For the initial configuration, particles (≈ 2000) are placed on a square lattice and given random initial velocities. Except if specified differently, the walls move at a uniform inward velocity v c = 2.5 · 10 −5 . We integrate Newton's equations of motion for both the translational and rotational degrees of freedom using a 4th order predictor-corrector method with time step δt = 1/50.
We use persistent homology [14] to quantitatively describe the structure of force networks. Briefly, each network is represented by two persistence diagrams that quantify the changes in geometry of the sets X θ for all force levels θ. The β 0 (β 1 ) persistence diagram encodes the appearance and disappearance of connected components (loops) by recording all values at which a connected component (loop) appears in a set X θ and all values of the force level at which the corresponding connected component (loop) disappears. For a rather brief discussion of persistence, we refer the reader to [12] and for a more indepth presentation to [11]. There are different types of well defined distances between the persistence diagrams. Each of these distances measures differences between the force levels at which distinct topological features appear and disappear. In particular, the distance between two β 0 persistence diagrams corresponding to different force networks measures differences between the force values at which distinct connected components/clusters are formed and merged together. In a loose sense, this corresponds to measuring differences between the force chains. The results reported in this paper do not depend on the particular choice of the distance and thus, for brevity, we will focus only on the so called Wasserstain distance, d W 1 , between the β 0 persistence diagrams (see [11] for a detailed definition).

Results
We start by recalling the results obtained for slow sampling rate of approximately 100 samples covering the packing fraction range [0.63, 0.9]. Since the evolution of force networks varies from realization to realization, and the distances between considered states of the system may be rather noisy, we considered the averages over a large number (20) of realizations [12]. Figure 1 shows that the distances between two consecutive states of a given realization are comparable to the distances between two states resulting from different realizations. This suggests that the memory of the system is lost on the time scale faster than the one defined by the slow sampling rate. To proceed, we consider much faster sampling rates. We also focus on presenting and analyzing the results from a single realization, since averaging the results over realizations obscures part of the information about force network evolution in a given system. Figure 2 shows d W 1 distance for four sampling rates r 1 , r 10 , r 100 and r 1000 ; here the index of a rate corresponds to the number of computational steps of duration δt between two consecutive sample points. The difference in the packing fraction, ρ, between two consecutive states sampled at the rate r i is ∆ρ = i · δρ where δρ ≈ 4 · 10 −8 is the change of ρ during one computational step. Note that r 1 corresponds to the fastest possible sampling rate, at which the state of the system is recorded at every time step. First, we note that, as expected, for faster sampling rates, there is a large number of points such that the distance is very small, suggesting that the force networks do not evolve much during a significant portion of the data shown. However, contrary to expectations, there is an increase of the peaks' sizes as the sampling rate is increased from r 1000 to r 100 and r 10 , followed by a gentle decrease as the rate is increased from r 10 to r 1 . This essentially means that one needs to sample at the rate comparable to r 10 (at which the peaks in d W 1 distance are the largest) to capture the corresponding large events present in networks' evolution. In turn, this implies that large changes between the force networks occur on the time scales shorter than τ c (capturing the system evolution with the sampling rate given by 1/τ c would correspond to r 50 ).
The number of data points shown in Fig. 2 is large; e.g., Fig. 2a shows about 10 6 data points. Since it is impractical to work with such a large number of data points, we proceed by carrying out a running average of the data shown in Fig. 2, and then proceed with considering a subset of the data extracted from the region very close to the jamming transition (that occurs at ρ ≈ 0.78 for the present system). Figure 3 shows the running averages for d W 1 distance. The averaging interval is given by 10 4 δρ ≈ 4 · 10 −4 . This averaging interval corresponds to the averaging time of about 200τ c which for a common granular system may range between a µsec and a msec. The actual number of points used for computing varies from 10 4 for the fastest sampling rate r 1 to 10 for the slowest sampling rate r 1000 .
Before jamming transition a typical distance between consecutive states grows as the sampling rate decreases, but at a rate that is typically much slower than expected. Note that the consecutive sampling rates differ by one order of magnitude. However, the differences between averaged distances is one order of magnitude only for the sampling rates r 1 and r 10 and becomes significantly smaller for slower sampling; the results are actually overlapping for r 100 and r 1000 . This finding suggests, consistently with the results shown in Fig. 2, that the time scale on which force networks evolve is comparable to τ c or even shorter. The fast time scale on which force networks evolve raises new questions. Recall that δt is the time interval needed for information about possible interactions to travel across a particle and hence δt is related to the speed of sound, v s , in the material. In particular, the approximate speed of sound, determined by the numerical time step, is given (in physical units) by d/(δtτ c ), and therefore, information about the features present in force networks that evolve on the time scale given by δt essentially propagates by the speed of sound. Thus, it is natural to ask whether the results may be influenced by δt. Figure 4 shows an influence of decreasing δt on d W 1 distance in a small interval based at ρ ≈ 0.78. Consequences of decreasing the δt are two fold. On one hand, as δt decreases, on average d W 1 decreases (the number of points at which d W 1 is very small increases as δt decreases). On the other hand, we notice that the peaks' heights do not necessarily decrease with δt, and more importantly, the peaks still relax on the scale corresponding to δt: for smaller δt's, the peaks relax faster (see the insets in Fig. 4). This finding is consistent with the general discussion above regarding the relevance of the speed of sound, and is showing clearly that in the considered system the force networks evolve on the time scale determined by the speed of sound in the material.
Another question to ask is related to the influence of the compression speed, v c , which defines a slow time scale in the problem. As v c is reduced, one may expect that the distance between the consecutive states of the system becomes smaller. As we show next, such trend is indeed observed in our results. Figure 5 shows d W 1 distance at the same ρ as in Fig. 4, but now for four (decreasing) v c 's. We see that, indeed, on average d W 1 decreases as v c is decreased, meaning that there are many δt intervals when the force networks essentially do not evolve. However, the height of the peaks does not necessarily decrease with v c ; we even find examples of larger peaks for slower compression. This finding can be understood by realizing that as the system is compressed slower, the particles have higher chance to rearrange, leading in some cases to large changes in the force network and increased distances between the consecutive states.

Conclusions
The force networks are found to evolve on the time scale that is much faster than the one introduced by externally imposed driving. As a consequence, an analysis of dynamics of force networks based on slow sampling rates may provide incomplete information about the networks' evolution, since slow sampling may miss large events that occur on faster time scales. In particular, we find that the time scale on which evolution occurs is specified by d/v s , where d is a typical particle diameter and v s the speed of sound in the material. It will be of interest to see to which degree such behavior can be seen in experimental systems. Furthermore, it should be noted that in the present work we focus on the systems very close to jamming transition; it remains to be seen whether such fast evolution as uncov-  . d W 1 distance between β 0 persistence diagrams corresponding to the consecutive force networks (one δt apart) for ρ ≈ 0.78 as the compression speed is reduced, with the part a) using the reference value, v c , and the parts b), c), and d) using v c /2, v c /4 and v c /8, respectively. δt is kept fixed. ered here may be seen further away from jamming, and/or for different type of external forcing.
The reported results rely on the use of computational homology that allows us to quantify the differences be-tween force networks as a granular system evolves. One particular strength of the computational homology approach is that it is system-independent and can be applied to any system, independently of the physical properties of the considered system, such as the type of interaction between the particles, or the number of dimensions of the physical space in which they live. Therefore, natural next step will be to consider particulate systems characterized by different physical properties (friction, polydispersity, cohesion, shape) or different geometry (2D versus 3D). Furthermore, the approach can be as well applied to experimental systems, such as those built from photoelastic particles, where detailed information about the force networks is available. Further analyses of such systems will be the subject of our future works.