Structural Evolution of Granular Systems: Theory

A general theory is developed for the evolution of the cell order distribution (COD) in planar granular systems. Dynamic equations are constructed and solved in closed form for several examples: systems under compression; dilation of very dense systems; and the general approach to steady state. We find that all the steady states are stable and for systems with CO$\,\leq 6$ they satisfy detailed balance. Illustrative numerical solutions of the evolution are shown. Our theoretical results are validated against an extensive simulation of a sheared system. The formalism can be readily extended to other structural characteristics, paving the way to a general theory of structural organisation of granular systems.

A general theory is developed for the evolution of the cell order distribution (COD) in planar granular systems. Dynamic equations are constructed and solved in closed form for several examples: systems under compression; dilation of very dense systems; and the general approach to steady state. We find that all the steady states are stable and for systems with CO ≤ 6 they satisfy detailed balance. Illustrative numerical solutions of the evolution are shown. Our theoretical results are validated against an extensive simulation of a sheared system. The formalism can be readily extended to other structural characteristics, paving the way to a general theory of structural organisation of granular systems. Introduction: Understanding and modelling selforganisation of dense granular matter (DGM) under external forces is essential to many natural phenomena and technological applications. Examples are: consolidation and failure of soils, packing of particulates in technological processes, initiation of avalanches, flow of slurries, and dense colloidal suspensions, to mention a few. This is also one of the most important problems in granular science [1-6] because both the dense flow dynamics and the large-scale properties of consolidated DGM depend strongly on the particle-scale structure. Such properties are: permeability to fluid flow [2], catalysis, heat exchange, and functionality of fuel cell electrodes [7]. In this paper we address this generic problem and develop a set of equations for the quasi-static structural evolution of two-dimensional (2D) systems of rigid grains.
Structural evolution of DGM is mediated by continual making and breaking of intergranular contacts. In turn, the temporal structural configuration determines the force transmission in the medium, which drives the evolution. The intergranular contact network can be regarded as a graph containing voids (or cells), each of which characterized by the number of grains enclosing it -its order. Developing a theory for the evolution of the cell order distribution (COD) is the specific aim here. This distribution has been argued [8] and shown [1, 2] to converge to a universal form, in which the intergranular friction is scaled away. Specifically, a set of evolution equations is constructed and solved in closed form, under some assumptions, both for very dense closed systems and for closed systems approaching a limit steady state. More general cases are solved numerically. We then compare the theoretical results with numerical simulations of sheared systems and show that there is good agreement. The evolution equations: When a 2D DGM undergoes slow quasi-static dynamics, it is possible to identify, at every moment, a contact network. In this graph the grain centres (however defined) are the vertices and edges connect grains in contact, with each edge corresponding FIG. 1. Making a contact between g and g , which occurs at rate q10,6, splits a 10-cell into two 6-cells (left to right). Conversely, breaking this contact, occurring at a rate p6,6, merges two 6-cells into a 10-cell.
to a specific contact. These make cells and we assign each a cell order defined as the number of grains surrounding it.
The structure evolves, through making and breaking of these contacts, which we call contact events (CEs). The former splits a cell into two smaller ones, while the latter merges two cells into a larger one. If a CE involves two rattler-free neighbour cells of orders i and j, the two processes satisfy i + j − − i + j − 2, exemplified in Fig. 1. Analogously, if the mother cell contains a rattler that participates in the event, the process sum rule is i + j − − i + j − 4. To model the dynamics of the COD, we define n k as the total number of k-cells in the system and their fraction of the total cell population, N c , as Q k = n k /N c . We also define the following rates, which we assume are system size-independent: p i,j = the merging rate of an i-and a j-cells into an i + j − 2-cell; q k,i = the splitting rate of a k = i + j − 2-cell into an i-and a k − i + 2-cells; r i,j = the merging rate of an i-and a j-cells, containing a rattler, into a k = i + j − 4-cell; and s k,i = the splitting rate of a k-cell, containing a rattler, into an i-and a (k − i + 4)-cell. The COD evolves via four basic CEs. arXiv:1904.06549v1 [cond-mat.soft] 13 Apr 2019 Two of creation: a k-cell is either a merger of j and k−j +2, or an offspring of a split large-order cell, and two of annihilation: either by splitting into two offsprings or merging to make a larger cell. Each of these processes has an equivalent when the combined cell contains a rattler, in which case the rates p j,k−j+2 and q k,j are replaced, respectively, by the rates r j,k−j+4 and s k,j .
If very large cells are rare then the occurrence of cells containing more than one rattler is very low and we ignore such events in the following. Then the evolution equations are: The 1/2 factor corrects for double counting and the terms containing δ-functions account for loss (gain) of two same order cells upon creation (annihilation) of a larger cell. The 'tilded' parameters in (1) are related directly to the size-independent rates: To simplify the following analysis, we ignore rattlers and set r i,j = s i,j = 0. Including the more realistic rattlerrelated terms is straightforward, but it would not add insight beyond the results presented below. From eqs.
(1), we note that k (k − 2)n k = E is a conserved quantity. This quantity has a physical interpretation: k n k ≡ N c and k kn k ≡ N z is twice the number of contacts, with z as the mean coordination number and N as the number of grains. Using Euler's topological expression in the plane to relate the numbers of vertices (=N grains), edges (=N z/2), and cells, N c : in which O( √ N ) are boundary terms. It follows that E = 2N +O( √ N ). Embedding the system on the surface of a sphere, E = 2(N − 1) exactly.
To make eqs.
(1) size-independent, we convert them for the fractions Q k . To this end, we define the deviation of the rattler-free steady state process i + j − − i + j − 2 from 'detailed balance', Similarly, is the deviation of the rattler-involving steady state process i + j − − i + j − 4 from a detailed balance-like behaviour. The sign of η i,j determines which direction of the process i+j i+j −2 is more frequent, with η i,j = 0 a detailed balance-like process. Generically, dilation or compression correspond to positive and negative η i,j , respectively.
Noting thatQ k = ṅ k − Q kṄc , the equations for the cell fractions becomė in which we useḋ Eqs. (5) are now conveniently system size-independent.
In practice, cell orders in realistic quasi-static systems cannot exceed an upper bound C and remain mechanically stable.
Focusing on closed systems with N ( 1) constant, we use Euler's relation again to derive a relation between the rates of change of N c and z: Exact solutions for 3-4 systems: The simplest non-trivial systems to study, realisable under highcompression processes, consist only of 3-and 4-cells. Eqs. (5) then reduce toQ The normalisation Q 3 + Q 4 = 1 makes the first two equations dependent. Eqs. (8) can be integrated straightforwardly to give Using (3), (8) and (9), we obtaiṅ This can be solved indirectly for Q 3 From this solution and eq. (8) we can obtain Q 4 and z. Examples of these solutions are shown in Fig. 2 for p 3,3 = 0.6 and q 4,3 = 0.4. We observe that the steady state is the same regardless of the different initial states and is determined only by the rates. As we show explicitly below, for C ≤ 6, this is not only a general feature but the steady state is also uniquely determined by the rates. This is illustrated for several 3-4-5 systems in Fig. 3 (b) and discussed in more detail below. The steady state: The rightmost sum in eqs. (5) describes the rate of change of the total number of cells and it vanishes at steady state. A direct calculation of eqs. (5) for very dense systems provides two general observations: for systems consisting of cell orders C ≤ 6, (i) there is only one possible steady state; (ii) this steady state satisfies η i,j = 0 for all i, j. Thus, the steady state of these systems satisfies detailed balance, which is very surprising in systems that are manifestly far from conventional equilibrium.
Moreover, while such a detailed balance-like solution is a valid steady state for any such dynamics, it is only one of the possible solutions for C > 6. This is because the number of relevant η processes in systems of cell orders up to C, is either (C − 3)(C − 1)/4, for C odd, or (C − 2) 2 /4, for C even. Having only C − 2 equations to determine the η i,j , these become under-determined for C > 6 leading to steady states that are strictly different from the above equilibria-like. For example, the steady states of a system with highest order C = 7 need only satisfy η 4,5 = η 3,6 = 0 and η 4,4 = −η 3,4 = −η 3,5 = η 3,3 . Therefore, if η 3,3 = 0, such steady states exist and they are not detailed balance-like. The approach to steady state: It is useful to understand the approach to the steady state, as well as its stability. Near this state, Q k (t) = Q s k + δQ k (t) is only slightly different from its steady state value, Q s k , . . , Q C ), and expanding the r.h.s of eq. (5) to linear order, we obtain for δ Q(t) in which the components of the constant matrix A are cumbersome combinations of p i,j , q i,j , and the steady state fractions Q s k . Denoting the ith eigenvector and eigenvalue of A, respectively, by v i and λ i , we have near the steady state Here, t 0 is the initial time and δ Q(t 0 ) the state at this time, which is presumed sufficiently close to steady state, such that Eq. (12) is a good approximation for all t ≥ t 0 . We first note that there can be no limit cycles, i.e. λ i is real for all i. This is because all 0 ≤ Q k , Q s k ≤ 1 are real for all k, which means that the eigenvectors composed of linear combination of them are real. Since all the components of the matrix A are also real, must be real for all i. Moreover, since the Q k s are normalised, at least one λ i must be non-zero before reaching steady state.
To show that the steady state is not unstable, i.e. λ i ≤ 0 for all i, we note that k δQ k = 0 and therefore k Q k remains normalised to this order. The existence of a λ i > 0 would lead to an indefinite growth of the terms involving e λit in (13). But since the Q k s are normalised at all times, then at some finite time, τ < ∞, this growth must stop abruptly. For this to happen, at least oneQ k s must diverge at t = τ . But since the r.h.s. of Eqs. (5) are all finite, this would violate the equations. Thus, there can arise no positive eigenvalue. Combined with the observation that there must be at least one nonzero eigenvalue, the conclusion is that the steady state is stable. Indeed, all the particular examples we studied numerically converged to a stable steady state, supporting this conclusion.
To determine the unique steady-state solution of the 3-4 system, we use the normalisation and detailed balancelike, η s 3,3 = 0, conditions to find with θ i,j ≡ q i+j−2,i /p i,j as rates ratio in expression (3).
represent the rate of approach to the steady state for p 3,3 = q 4,3 = 1.
The dynamics of 3-4-5 systems can be analysed similarly. Using η s 3,3 = η s 3,4 = 0 and the normalisation condition, the steady state solution satisfies and Using the positivity of the Q k , the uniqueness of the steady state can be established by studying the position of extrema of (16) [11] and it is determined by the rate fractions, θ i,j . The normalisation confines the dynamics to the plane Q 3 + Q 4 + Q 5 = 1 and the approach to steady state in this plane is illustrated in Fig. 3(b) for two different sets of rates, which give two distinct steady states. For C ≥ 6, the dynamics take place on the hypersurface k Q k = 1, on which there is at least one detailedbalance-like steady state. As discussed above, for C > 6, there may be additional steady states. Moreover, since all the steady states are stable then, if they exist, they must be distinctly separate on the hypersurface. Simulations: To test the theory, we carried out numerical simulations of quasi-static two-dimensional simple shear between parallel plates. We use a dissipative spheres model [12,13] and the implementation of the model was provided by the open source project LIGGGHTS [14]. Within this implementation, we used a Hertz force model for the grain-grain and grain-walls interactions, with: Young's Modulus 5 × 10 6 Pa, restitution coefficient 0.2, Poisson's ratio 0.5, and friction coefficient 0.5. The system consisted of N = 21, 690 spheres restricted to move in the plane z = 0 of four different diameters, 7mm (5402 grains), 9mm (5400), 11mm (6120) and 14mm (4770), respectively.
The system has a length of 1.65m and periodic boundary conditions in the x-direction. The initial state was generated by compressing the grains from an initially loose random distribution by a flat surface at constant pressure, P initial = 54.45Kg/m, in the y-direction until their total kinetic energy fell below 10 −12 J per particle. This state was regarded as mechanically stable. Then, maintaining the confining pressure, shear started in the x direction, with shear velocity v γ = 0.06m/s. No gravitational field was applied.
The time-step was 10 −6 s and the grain positions and velocities were saved every 200 time-steps for collecting detailed data on the contact network evolution. At each stop, cells and the grains surrounding them were identified, from which we tracked the evolution of the COD and obtained the rates of contact events, p ij and q ij . Cells of orders higher than 13 occurred very rarely and, therefore, were excluded from the analysis. Analysis of simulation data: The rates are determined from the numbers, N i,j , of the event i + j − → i + j − 2 and its reverse, respectively, during a time interval ∆t, and the total number of cells, N c : In general, the rates should depend on the force distribution and the shear rate, both of which change slightly during the simulation as the system dilates before reaching a steady state. Since the above analysis is for constant rates, this could complicate a direct comparison between the analytical solutions and the simulation data. Fortunately, the rates change little and their time dependence can be neglected. Low-statistics cells, e.g. for e ≥ 10, give large numerical errors in eqs. (18) and we set to zero the rates for which either Q i Q j or Q i+j−2 drops below δ ≈ 7.6 · 10 −4 , amounting to fewer than 15 cells, or there were fewer than 40 events in 0.1s. The threshold for non-zero rates was set at 60 events per 0.1s and the used rates spanned 60 − 1000 events per 0.1s. Interestingly, in spite of the long simulation time, the rates fluctuated significantly, 30% − 65%, perhaps because of the sensitively to small changes in Q k . Smaller fluctuations may be achieved in larger systems with more contact events. Fig. 4 shows the rate diagram, calculated from data collected during intervals of 0.1s and averaged over the entire evolution process. Such rate diagrams characterise uniquely a dynamic process. Using the calculated rates, shown in Fig. 4, we solved the evolution equations (5) numerically, with the initial COD determined by averaging over the first 0.1s in the simulation. While there was a reasonably good agreement between the solution and the simulation data when using the computed mean rates, we found that an even better agreement was achieved by correcting only q 6,4 by 10%. Note that this correction is much smaller than the aforementioned statistical fluctuations of the rates. The good agreement between our solution and the simulation data after this minor correction is shown in Fig. 5. Conclusion: To conclude, we have constructed master equations to describe the evolution of a key structural characteristic of granular media -the cell order distribution (COD) -from any initial state, given the cell merging and splitting rates. The equations yield a surprising result: the steady states of the non-equilibrium dynamics of granular systems, with cell order no higher than C = 6 satisfy detailed balance. This detailed balancelike steady state is unique and stable. Systems including cell orders of 7 and higher can converge to this steady state, as well as to other solutions that do not support detailed balance, all of which are stable fixed points of the equations.
We validated the theory by running a long simulation of a sheared system, determining the rates, and using those to solve the master equations. Indeed, the calculated solution agrees nicely with the simulated COD evolution, in spite of large fluctuations in the rates computed from the simulations. These fluctuations may be caused by 'clappers': Particle pairs that make and break their common contact repeatedly. We plan an experiment to study the relative contributions of clappers and non-clappers to the rate statistics.
In our analysis, we assumed constant rates throughout the dynamic process. This is unlikely to be the general case and extending the theory to time-dependent rates is the next step. Such time dependence is expected because the contact event rates is likely to be sensitive to the intergranular force distribution, which is positionand time-dependent. Indeed, it has been argued [2, 15] that the structural evolution of granular system is a selforganisation process coupled to the evolution of the intergranular forces. Therefore, these equations and their  extension form an important step towards a complete understanding of the structure-forces co-evolution. Another extension would be to open systems, when the number of particles is not conserved. This extension could be related to a grand-canonical statistical mechanical description of granular ensembles [16].  Structural Evolution of Granular Systems: Theory

Supplementary Materials
The approach to steady state: In the section about the approach to steady state we claimed that k δQ k = 0. To establish this identity, we start from the full linear expansion of the evolution equations sufficiently close to steady state: Summing over the index k and using k Q k = 1, we obtain Renaming the index as k, the individual sums neatly cancel out: Furthermore, the cell fractions, Q k (t 0 ), must be normalised at all times and, in particular, close to steady state. It follows that k δQ k (t 0 ) = 0. Since kδQ k = 0, as the above calculation shows, we obtain that k δQ k (t) = 0 for all times t in the linearised regime as well.
Phase diagram: For C ≤ 6 the steady state can be uniquely expressed in terms of the fractions of the occurring rates for reaction and back-reaction, p i,j /q i+j−2,i . Therefore, multiplying all the rates by a constant factor does not affect the steady state, but only modifies the time by which the steady state is reached. This allows us to scale the rates and represent the steady state fractions, Q k and the mean coordination number, z, as contours in a phase diagram in a phase space spanned by the rate fractions. In particular, the phase diagram of a system, containing only cells of orders 3, 4 and 5, can be conveniently represented by a density plot, see Fig. 1. Fig. 1 shows phase diagrams of (a) the ultimate fraction of 3-cells and (b) the mean coordination number z, in the phase space spanned by the rates. Such phase diagrams demonstrate how the master equations can provide guidelines to design specific packing protocols. Clapping events: A particular phenomenon, observed extensively in all the simulations we ran, is the occurrence of 'clapping' events. These are contact events involving pairs of particles making and breaking the contact between them repeatedly. For simplicity, we defined non-clapping events in our analysis as events that occur only once during the system's dynamics. However, since simulations and experiments take place over finite duration, non-clapping events are defined practically as events that occur only once during the simulation/experiment. This provided us with a straightforward way to identify and distinguish between the two types of events in our simulations.
Evidently, where clapping occurs, the grains did not move away from one another, which may mean that the local configuration has not changed much. This implies that clapping events only add noise to the statistics of events and evolution of the COD. The probability that grains have moved away from one another and have come back together at a later time is expected to be very low and we have not observed any such event for the 500 clapping pairs, which we tracked explicitly during the entire simulation. A potential problem with the above definition is that events occurring close to the end of the simulation could be identified wrongly as non-clapping because the grains do not have the time to clap. This should not affect the statistics much because it means that half a clap has been counted as a non-clapping event and we just missed the other half of the clap. This would be significant only if the clap affects two cell types whose splitting or merging is rare. However, to avoid this potential error, reliable results should be obtained sufficiently early and far away from the end of the simulation. It is important to note that clapping does not contribute to the structural evolution on a coarse grained time scale, which the master equations come to describe. Clapping obeys 'detailed balance' in the sense that the number of 'reactions' equals the number of 'back-reactions' and η clap i,j = 0. Furthermore, the master equations,Q k = f ({η i,j }), are linear in the ηs, . This means that the contribution of the clapping events toQ in the master equations vanishes on a coarse-grained time scale that is sufficiently longer than the time between clapping events. We found in our simulations that most clapping pairs were active for some duration, went dormant, and became active again (see Fig. 2 below). Whether this is a result of the specific dynamics (e.g. shear, compression, etc.) is a question we intend to explore in future work.

Discrete Element Method (DEM) simulation:
The main text describes simple shear simulations. To study another dynamic process, we analysed a long-time DEM simulation of a bi-axial compression. The simulated system consisted of 22, 381 discs, whose diameters were distributed log-normally [1, 2]. The discs interaction was harmonic, with tangential and normal spring constants k t and k n , respectively, obeying k t /k n = 1/4. The mean overlap at the contacts is estimated as d/D = σ c /k n = 10 −5 , where σ c is the initial confining pressure. This overlap is much smaller than the one used for the simple shear simulation. The restitution coefficient was set to 0.98 -much larger than in the simple shear simulations. The initial configuration was prepared with a very low inter-particles friction. At the start of the simulation, the particles were assigned a friction coefficient of µ = 0.5, the same as in the simple shear simulation. The system was periodic in both dimensions and was subjected to a gradual vertical compression keeping the lateral confinement pressure constant. This resulted in a quasi-static dilation with a diagonal shear. The y-axis shows the inverse time between subsequent claps. It can be observed that some clapping pairs, e.g. the blue and orange pairs, are active for a while, lay dormant and become active again at a later stage.
The particle stiffness and the high restitution coefficient enhance clapping, making it a good system to study this phenomenon. We chose to trace 448 individual clapping pairs throughout the simulation, as tracing all pairs would have been too time-consuming. We found that a pair can indulge in repeated clapping, disengage for a while, and become active again later. This can be seen in Fig. 2, in which we show the inverse time between claps as function of the time in the simulation for 10 arbitrarily picked pairs. Clapping dominated the contact events in this simulation: only 0.7 % of all events were nonclapping, suggesting that the structure is hardly changed within large regions. In this respect, a dilation process is quite different from simple shearing in that it does not mix the system quite as well.
The wide distribution of the time intervals between claps of a specific pair is evident from the very long tail of the histogram in Fig. 3. The clapping seem also to be spatially uncorrelated, occuring uniformly across the system, as shown in Fig. 4.
Unfortunately, the dilation simulation does not reach a steady state in which the rates remain constant. The most likely reason is that the boundary forces change continually with time, which changes the statistics of the intergranular contact force distribution, in turn affecting the rates. Indeed, the rates appear to decay with time exponentially, ∼ e −t/τ , with τ ranging from 0.2 to 2.2, depending on the specific rates. Since the theory presented is limited to constant rates, we have not pursued a detailed modelling of the COD evolution of this simulation. The continual decay of rates, in spite of the long time simulation, further suggests that a steady state may not be reached before the dilation disrupts mechanical equilibrium. This highlights the fact that not every process reaches necessarily a steady state. * ccw45@cam.ac.uk † rbb11@cam.ac.uk