Steady states of two-dimensional granular systems are unique, stable, and sometimes satisfy detailed balance

Understanding the structural evolution of granular systems is a long-standing problem. A recently proposed theory for such dynamics in two dimensions predicts that steady states of very dense systems satisfy detailed-balance. We analyse analytically and numerically the steady states of this theory in systems of arbitrary density and report the following. 1. We discover that all such dynamics almost certainly possess only one physical steady state, which may or may not satisfy detailed balance. 2. We show rigorously that, if a detailed balance solution is possible then it is unique. The above two results correct an erroneous conjecture in the literature. 3. We show rigorously that the detailed-balance solutions in very dense systems are globally stable, extending the local stability found for these solutions in the literature. 4. In view of recent experimental observations of robust detailed balance steady states in very dilute cyclically sheared systems, our results point to a self-organisation of process rates in dynamic granular systems.

1. Introduction Granular matter is ubiquitous in nature and plays a major role in our everyday life. Its near-indifference to thermal fluctuations has earned it a recognition as a new form of matter [1]. In spite of many decades of intensive theoretical, numerical, and experimental investigations into this form of matter, new aspects of its rich and complex behaviour are being discovered. The sensitivity of the large-scale behaviour and properties to the particle scale characteristics and structure has hindered the modelling of granular matter to date [2,3]. Consequently, of key significance is the modelling of granular dynamic evolution and the mechanically stable structures that such dynamics settle into. In particular, when the dynamics are quasistatic, the steady-state dynamics determine those stable structures and it is on this type of dynamics that we focus here.
Several methods have been proposed to describe and model the evolution of the underlying structure during quasistatic dynamics [4][5][6][7][8]. A general way to describe mechanically stable granular structures in two-dimensions (2D) is based on what is known as the cell order distribution (COD) [9][10][11]. A cell is the smallest void (loop) in the system, surrounded by particles in contact and its order is defined as the number of particles in contact surrounding it. During quasi-static dynamics, the COD evolves by intergranular contacts being made and broken, which split and merge cells, respectively. Such a process is shown schematically in Fig 1. The dynamic equa- 5,6 p 5,6 q FIG. 1: A sketch of cellular splitting and merging events: when the circled contact breaks, the two cells merge, 5 + 6 → 9 and when it is made the higher order cell splits, 9 → 5 + 6.
tions for the COD evolution are [12]: In these equations, Q k is the fraction of cells of order k, referred to in the following as k-cells. p i,j is the merging rate of i-and j-cells into an (i + j − 2)-cell, q i,j is the rate of the splitting of an (i + j − 2)-cell into an i-and a j-cell, and C is the highest possible cell order in the system. The factor 1/2 and the δ-functions ensure correct counting. The last term on the right hand side is needed because the total arXiv:2306.10526v2 [cond-mat.soft] 7 Jul 2023 number of cells changes with each merging or splitting event, which changes the fractions Q k . Rattlers, which are particles with one or no contact, were ignored in these equations, which is a good approximation for dense systems with low-order cells.
Including rattlers in the analyses that follow is possible, but the added complication does not add much insight and we disregard them. Wanjura et al. [12] found that, under some conditions, the steady-state cell order transitions of these far-from-equilibrium systems satisfy detailed balance, when the cell orders do not exceed six, a result corrected later to five [13]. The steady states of the evolution equations were also shown to be locally stable. Recent experiments on quasi-statically sheared 2D granular systems [14] have revealed a surprising observation -they always settled into steady states that satisfy detailed balance. Such robustness suggest that these steady states are not only stable, but also may be unavoidable. Moreover, these observations appear to contradict the paradigm that steady states of non-equilibrium dynamics cannot satisfy detailed balance [15]. Motivated by these experimental observations, we analyse here eqs. (1) in detail. We investigate the conditions for detailed balance and the properties of such steady states. We show that: (i) if a steady state satisfies detailed balance in systems where C = 6, 7 then it is the only possible steady state; (ii) there is strong numerical evidence that, in systems where C = 6, 7, or 8, only one physical steady state is possible, whether or not it satisfies detailed balance; (iii) the steady state in systems where C = 5 is not only locally but also globally stable. These findings provide a partial explanation for the observed convergence to detailed balanced steady states in [14].

A. General steady states
The back-and-forth processes (i)+(j) ⇌ (i+j −2) are equivalent to chemical reactions in a multicomponent reactive system. Their net flux is η i,j = p i,j Q i Q j − q i,j Q i+j−2 and η i,j = 0 when they are balanced. At steady state, the sum of all the processes, i,jη i,j , vanishes by definition and eqs. (1) reduce to in which the bars indicate steady state values. Given C, there can be (C − 2) 2 /4 or (C − 2) 2 − 1 /4 processes, when C is even or odd, respectively. Focusing on the even case, for brevity, it is useful to rewrite eqs. (2) as Here, H is a (C − 2) × (C − 2) 2 /4 matrix and the vectorη's components are all the steady-state ηfluxes. Thus,η must exist within the null space of H. The normalisation constraint, kQ k = 1, reduces the number of independent first-order equations in (3) to C − 3. Thus, by Bézout's theorem [16] and, since η i,j are quadratic in the Q i s, the maximal number of solutions is 2 C−3 for any given set of rates, p i,j and q i,j .
Below, we show analytically that, at least up to C = 7, only one of these solutions is physical -the detailed balance steady state (when it exists), in which η i,j = 0 for all i, j. Indeed, an extensive numerical investigation over a wide range of parameters supports this conclusion -all other numerical solutions included either imaginary or negative Q k fractions.

B. The detailed-balance steady state is unique
The uniqueness of the detailed balance steady state was established for C < 6 [12,13]. We now extend this result to systems comprising arbitrary orders. Defining θ i,j = p i,j /q i,j , each balanced process satisfiesQ i+j−2 = θ i,jQiQj . It follows that The detailed-balance steady-state solution depends solely on the ratios θ i,j and onQ 3 . The value ofQ 3 can then be found from the normalisation condition: We note in passing that the values of θ i,j are not independent because cells of order k > 5 can be formed by more than one process [17]. An example of a system in which detailed balance is possible (indeed observed) is the experimental steady states observed in [11], which satisfy θ i,j = θ 0 for all i, j. This gives rise to an exponentially decaying COD. Imposing such a condition, we calculate explicitly the COD of the emerging detailed-balance steady state of two systems in which C = 7 and θ 0 = 0.5 and 2. These are shown in Fig. 2. Nevertheless, this dependence does not affect our following argument regarding the uniqueness of the detailed balance steady state. Since the rates p i,j and q i,j must be non-negative then θ i,j ≥ 0, for all i, j, and the left hand side of (5) is a monotonically increasing function ofQ 3 . Additionally, for C > 3, it vanishes atQ 3 = 0 and exceeds 1 atQ 3 = 1. Thus, only one solution exists in the range 0 <Q 3 < 1 and this is the only possible detailed-balance solution. This conclusion holds for both odd and even values of C. It follows that, if a system evolves into a detailed-balance steady state then that state is unique and independent of initial conditions. This goes some way to explain the robustness of the recently observed detailed-balance in [14].
The next question is whether or not there also exist steady states that do not satisfy detailed balance. There are (C − 2) 2 − 1 /4 or (C − 2) 2 /4 potential finite values of the fluxesη i,j to determine from eqs.
(2), for C odd or even, respectively. With only C − 3 independent eqs. available, in the case of C = 5 there are two eqs. for two unknowns. It follows that detailed balance is the only steady state for systems up to C = 5. Since such systems have only two processes, 3 + 3 ⇋ 4 and 3 + 4 ⇋ 5, each must be balanced separately, as no cycle is possible [13]. This detailed balance is then straightforward. However, sinceη i,j are underdetermined for C > 5, it was conjectured that, in addition to detailed-balance, such systems support an infinite number of other stable steady states [12,13]. To investigate this conjecture, we re-examine the steady state.
In the supplementary material, we use a similar analysis to show that the same conclusion holds for C = 7. This conclusion improves on the conjecture in [12]. We believe that this line of proof can be extended to C > 7, although not without substantial effort. However, a different approach is required for a general proof.
These results, combined with the extensive numerical investigations reported below, and the experimental observations in [14], lead us to conjecture that, when a detailed balance steady state exists, it is the only possible steady state for any value of C.

D. Numerical investigation of non-detailed-balance steady states
The rates parameter space is infinitely large and most combinations of rates lead to non-detailedbalance solutions (with detailed balance only possible if the relation established in I C is satisfied, i.e. θ 3,3 θ 4,4 = θ 3,4 θ 3,5 ). To understand the nature of these solutions, we explored the parameter space numerically (see also supplementary material). Starting with C = 6, we tested the 4 8 = 65536 rate value combinations when each rate can assume any of the four values: 0.1, 0.5, 1.0, and 3.0. For each combination, we found all the solutions and noted the number of solutions. Unexpectedly, each combination gave rise to only one physical solution, with all the others containing either negative or complex values of someQ k .
To test the potential generality of this surprising observation, we solved numerically for the steady-state solutions in systems where C = 7 and 8. These have, respectively, 6 and 9 processes and 12 and 18 variable rates.
Owing to the required larger computational resources, we tested only 6 rate combinations for C = 7: , and in total 2 × 2 18 = 524288 different systems. We found that in none of these systems was there more than one physical steady state solution.
Based on these investigations, we conjecture that, for any choice of constant rates, there is only one physical solution regardless of the upper order, C. In Fig. 3, we show an example of the one non-detailedbalance solution when C = 6, for a set of rate parameters that also admits five other non-physical solutions (listed in the supplemental material). It can be observed in Fig. 3 that the difference between the breaking and making rates are the same for all rates. In particular, the vertical offset from the detailed balance line are the same and equal to A. A typical figure for C = 8, for a system that also does not satisfy DB, is shown in the supplementary material.

II. GLOBAL STABILITY FOR C = 5
A linear analysis of the steady states of eqs.
(2) has shown them to be asymptotically stable [12], a prediction that has been supported experimentally [14]. This, however, does not preclude possible limit cycles around the steady state away from the linear regime. We investigate next the global stability of the solution and show that, at least for the dense system comprising cells of orders 3−5, no such cycles exist.
Using the normalisation condition to eliminate Q 5 , the independent evolution equations can be written aṡ in which the rates q i,j and θ i,j are assumed for now to be constant (more on this condition below), R ≡ q 3,4 /q 3,3 , κ 3,3 ≡ θ 3,3 Q 2 3 −Q 4 , κ 3,4 ≡ θ 3,4 Q 3 Q 4 +Q 3 + Q 4 − 1, and time is scaled: t → t ′ ≡ q 3,3 t, such thaṫ The fractions are constrained by Q 3 , Q 4 ≥ 0 and Q 3 + Q 4 ≤ 1. We analyse eqs. (9) within this region of the Q 3 − Q 4 plane, using the theorems of Bendixson and Poincaré-Bendixson [18]. The former states that, in two-variables dynamical systems,ẋ = f (x, y) andẏ = g(x, y), if V (x, y) = ∂ x (f ) + ∂ y (g) is non-zero and has the same sign throughout a simplyconnected x − y domain, then no closed orbits can lie within that domain [18]. We define By the inequalities of the arithmetic and geometric means, Using (11) in (10), we establish that Since R > 0, V (Q 3 , Q 4 ) < 0 throughout the region to which the system is confined. Thus, according to the Bendixson theorem, this system has no limit cycles. Combining this result with the established uniqueness of the steady state [13] then, by the Poincaré-Bendixson theorem, the limit set contains only that steady state [18]. Thus, the detailed balance steady state is globally stable for any physical initial condition.

III. CONCLUSIONS AND FUTURE WORK
To conclude, we studied, both theoretically and experimentally, the nature of the detailed-balance steady states into which the non-equilibrium dynamics of granular matter has been found to settle. We have proven that, for any maximum cell order, C, there can only be one detailed-balance steady state. We have also shown rigorously that, if a detailedbalance steady state solution exists up to C = 7 then it is the only solution.
Intriguingly, by solving the steady state equations numerically for 614400 systems up to C = 8, we found that there is always only one physical solution, in which Q k is real and positive for all k. We conjecture that the evolution equations (1) yield only one physical solution for the steady state, which may or may not satisfy detailed balance. This conjecture is supported by clear experimental observations of detailed-balance steady states for systems with C > 10 [14].
Next, we used the theorems of Bendixson and Poincaré-Bendixson to show that the detailed balance solution of the dynamics of systems with C = 5 is globally stable and no periodic orbits exist. This may explain the robustness of such solutions observed experimentally [14].
It should be noted that, in our discussion, the steady-state rates p i,j and q i,j are constant, but of arbitrary functional form, subject to the condition of detailed balance. However, they need not be constant and evolve during the approach to the steady state.
Another intriguing implication of our results is the following. We have found that most systems with C > 5, settle into steady states that do not necessarily satisfy detailed balance. For example, solving and finding the only solution when C = 6, p 3,5 = p 4,4 = q 3,3 = q 3,4 = q 3,5 = 1, and q 4,4 = 0.5, we find that not all η i,j = 0, namely, no detailedbalance. Yet, the experiments of Sun et al. [14] reveal that, in a range of quasi-statically cyclically sheared 2D granular systems, the steady states always satisfy detailed-balance. This suggests that the cell breaking and merging rates in those experiments were neither arbitrary nor time-independent. Rather, they must have evolved as the granular systems self-organised into values that satisfy detailed balance. This seems the most plausible explanation for the emergence of detailed balance in those experiments and could reconcile the discrepancy between those observations and what appears to be a violation of the paradigmatic Klein principle [15].
Rate equations similar to (1) have been used for modelling many evolution processes. In particular, in models of aggregation-fragmentation (AF). Nevertheless, there are some similarities and differences between the granular dynamics we study here and those models. The similarities are that our cell order fractions, Q k , are analogous to aggregate sizes and all models have transition rates. Additionally, many models of AF assume detailed balance, implicitly or explicitly [19][20][21][22][23]. One minor difference is that the normalisation of our cell order fraction takes into consideration the changing total number of cells, which which gives rise to the last term in our eq. (1). There are, however, more significant differences. One is that, unlike in most of AF models, for example, we do not assume the mathematical forms of the steady-state rates, p i,j and q i,j , subject to the condition that they satisfy detailed balance. This makes our analysis more general and more applicable than the studies that make such assumptions. Another difference is that most studies of AF dynamics simplify the analysis by allowing the size of aggregates to tend to infinity, which is often not physical. In contrast, our analysis applies to arbitrarily finite highest order, C. The third difference is quite fundamental. As mentioned, many models of AF processes assume detailed-balanced steady states. While this phenomenon is well established for systems in equilibrium, the current belief in the community is that it cannot be satisfied in out-ofequilibrium steady states. Indeed, our motivation to study steady states of sheared granular system is the surprisingly strong experimental evidence in [14] of persistent detailed-balanced in the steady states of their granular dynamics. To the best of our knowledge no such evidence exists for the non-equilibrium systems, to which some AF models presume to apply. This also means that our results apply only to AF systems either in strict equilibrium or where detailed balance has been established experimentally.