#  Dynamic constraints predict the relaxation of granular materials

Qinghao Mao  Department of Physics, University of Chicago, Chicago, 60637, IL,
USA  Yujie Wang  State Key Laboratory of Geohazard Prevention and
Geoenvironment Protection, Chengdu University of Technology, Chengdu 610059,
China  Department of Physics, College of Mathematics and Physics, Chengdu
University of Technology, Chengdu 610059, China  Walter Kob  Department of
Physics, University of Montpellier and CNRS, F-34095 Montpellier, France
${\dagger}$<EMAIL_ADDRESS>
${\ddagger}$<EMAIL_ADDRESS>
∗<EMAIL_ADDRESS>
Abstract Granular materials such as sand, powders, and food grains are
ubiquitous in civil engineering, geoscience, agriculture, and medicine. While
the influence of friction between the grains on the static structure of these
systems is well understood, its impact on the dynamics is an open problem.
Here we use particle-based simulations of a granular pack under cyclic shear
and discover that the relaxation time of the system is a non-monotonic
function of friction. By introducing the concept of dynamic constraints, we
reveal that this re-entrant dynamics is due to the competition between
increasing frictional coupling and a concurrent change in the structure of the
granular pack. Our theoretical approach, which unifies the dynamics of
friction-less systems with frictional ones, is applicable to other systems
that have a complex free energy landscape and a dynamics which involves time-
dependent constraints, thus setting the stage for a description of the dynamic
behavior of a large class of complex systems.

Due to the importance of granular materials for our daily life, they have been
in the focus of interest of a multitude of studies [1, 2, 3, 4, 5]. Critical
insight into their static properties has come from various approaches,
including Maxwell’s method which connects mechanical constraints with the
rigidity of a structure [6, 7, 8, 9, 10], force chain networks [11, 12, 13,
14], as well as the Edwards ensemble [15, 16, 17, 18], which allows a
statistical mechanics description of the structure. Despite these successes,
our understanding of the microscopic dynamics of granular packs, which
features limit cycles, caging, convection, and a strong dependence on the
history and driving protocol [19, 20, 21, 22, 23, 24, 25, 26, 27], remains
unsatisfactory [28]. This lack of insight is in contrast to the case of
thermal glass-formers, which have a structure similar to the one of granular
systems but whose relaxation dynamics is described reasonably well by multiple
approaches [29, 30, 31, 32]. One reason for this striking difference is that
particle friction renders granular systems dynamically constrained and
dissipative [33, 19, 34, 35], hence requiring non-equilibrium approaches [2].

Theoretical studies of model glass-formers have revealed that dynamical
slowing down can be caused by the presence of constraints [36, 37, 38]. These
findings and the fact that the rigidity of granular systems can be explained
by constraint counting [39, 40, 41] hint that this approach can be generalized
to describe the dynamics of driven granular materials. Here we probe this
dynamics and demonstrate that dynamic constraints indeed allow to rationalize
it. We argue that this approach is quite general and hence can be applied to a
multitude of disordered many-body systems [42, 43, 44, 45].

System: We investigate a 2:1 mixture of 2D frictional Hertzian particles
having a diameter $d_{s}=\text{1~{}cm}$ (our unit of length) and
$d_{b}=\text{1.4~{}cm}$, and a density $\rho=1\text{g}\cdot\text{cm}^{-2}$.
The tangential component of the force, ${\bf F}_{t}$, is restricted by the
Coulomb cone $|{\bf F}_{t}|\leq\mu|{\bf F}_{n}|$, where $\mu$ is the friction
coefficient and ${\bf F}_{n}$ is the normal force at the contact. Details of
the model and the simulations are given in the Methods. The total number of
particles is $N=10^{4}$, and they undergo cyclic shear at constant pressure
with a strain amplitude of 0.05 in a box of size 120 cm $\times$ 120 cm
consisting of four rigid walls made of particles, see Inset of Fig. 1b. In the
following, time will be given in number of shear cycles.

Re-entrant dynamics: Figure 1 shows the time dependence of the translational
and rotational mean squared displacements, TMSD and RMSD, respectively. The
TMSD demonstrates that at small $\mu$ the dynamics has a fast, diffusion-like
behavior while an increase of $\mu$ to 0.014 leads to a motion in which the
particles are temporarily caged before they finally start to diffuse. For
somewhat larger frictions the TMSD is small and basically constant, i.e., the
particles are completely caged. This $\mu-$dependence is thus reminiscent of
the slowing down of thermal glass-forming systems if the temperature is
lowered [29], except that here it is the frictional coupling between the
particles that is the control parameter. Surprisingly one finds that if
friction is increased beyond $\mu\approx 0.2$, the dynamics speeds up again
and at the highest $\mu$ the TMSD becomes hyperdiffusive. This re-entrant
phenomenon, which is in qualitative agreement with the results of Ref. [23]
for a 3D system, hints that the relaxation dynamics is controlled by two
competing processes that depend on $\mu$ and the nature of which will be
elucidated below.

The RMSD shows that at very small $\mu$ the rotational motion of the particles
is slow and ballistic, i.e., increases quadratically with time, a result that
is reasonable since the low friction and the slow shearing makes it hard for a
particle to change its angular velocity. With growing $\mu$ the RMSD at short
times increases quickly and its slope in this double logarithmic
representation decreases to 1, indicating a rotational diffusive dynamics. For
$\mu$ larger than $\approx 2\cdot 10^{-3}$ this dynamics slows down again
while its time-dependence turns back to ballistic. Finally, when $\mu$ is
larger than $\approx 0.2$, the rotational movement speeds up again and becomes
diffusive. Hence we conclude that also the rotational degrees of freedom show
a re-entrant behavior as a function of friction.

Since the time dependence of the MSDs does not show for all $\mu$ a diffusive
behavior, the diffusion constant cannot be used to quantify the dynamics. To
this aim we therefore monitor MSD($t=1$) as well as the long-time slope of the
MSD in a double-logarithmic plot. Figure 2(a) shows that for $\mu\lesssim
2\cdot 10^{-3}$ the TMSD has a slope of 1.0 and then drops quickly to a value
that is basically zero, i.e., the dynamics changes from diffusive to caged.
This caging regime extends from around $\mu_{\rm min}=2\cdot 10^{-2}$ to
$\mu_{\rm max}=0.2$ and will be referred to as “glass-zone” since the TDOF are
completely frozen. For $\mu>\mu_{\rm max}$ the slope shoots back up to $1.0$
and remains constant. The exponent for the RMSD decreases from 2.0 at $\mu=0$
(ballistic motion) to a value slightly above $1$ at $\mu\approx 10^{-3}$. If
$\mu$ is increased further, the slope rises quickly and becomes 2.0 in the
glass-zone, demonstrating that some of the particles undergo free rotation
[22, 46]. For $\mu\geq\mu_{\rm max}$, the slope decreases again and becomes
1.0, signaling diffusive rotation.

The short time dynamics is quantified by the value of the MSD at $t=1$, Fig.
2(b). At small friction the TMSD($t=1$) is large and stays constant for
$\mu\leq 10^{-3}$, i.e., in this range of friction the short time
translational dynamics is independent of $\mu$ and the system is similar to 2D
glass-formers [47]. Upon approaching $\mu_{\rm min}$, TMSD$(t=1)$ starts to
drop and reaches a minimum at around $\mu=0.15$. This decrease, which is
directly related to a shrinking cage-size, is compatible with a power-law with
an exponent of -2.0. Beyond $\mu_{\rm max}$, TMSD$(t=1)$ increases by two
orders of magnitude, signaling the acceleration of the dynamics, before
reaching a plateau. The RMSD$(t=1)$ shows qualitatively the same trend,
although the $\mu$-dependence is much weaker. Note that the first signs that
friction influences significantly the dynamics are seen at around
$\mu=10^{-3}$. This surprisingly small value implies that the dynamics of real
granular systems cannot be described theoretically if friction is not taken
into account.

Fig. 2(c) presents the $\mu$-dependence of different characteristic times
$\tau$ (defined in the caption). All times increase exponentially quickly upon
approaching $\mu_{\rm min}$ (left inset of panel (c)), beyond which $\tau$
cannot be determined reliably anymore. In contrast to this we find for
$\mu\geq\mu_{\rm max}$ that $\tau(\mu)$ is given by a power-law with a
critical point at $\mu_{\rm max}$, right inset of panel (c), hinting the
relation to the isostaticity point of the pack [10]. The presence of these two
different growth laws indicates that the mechanisms leading to the increase of
$\tau$ are distinct.

Dynamic constraints: For many-body systems the slowing down of the dynamics is
often related to a change in structural quantities like the packing fraction
or coordination number [29, 30]. Fig. 3(a) shows that the cycle averaged
coordination number, $\overline{Z}$, is slightly above 4.0 if $\mu$ is small,
while for large friction it slightly exceeds 3.0, thus compatible with the
values expected for a 2-dimensional static pack of particles with infinite
elastic moduli [8]. That $\overline{Z}$ (as well as the packing fraction, see
Extended Data Fig. 2) decreases simply monotonically with $\mu$ indicates that
static quantities like $\overline{Z}$ are not able to rationalize the re-
entrance dynamics and one needs also to consider the dynamics on the contact
level. Hence we probe the probabilities $\eta$ that during a cycle a contact
i) breaks, ii) persists and slides, or iii) is completed frictional locked,
see cartoons in Fig. 3(b). That panel demonstrates that $\eta_{\text{\rm
break}}$ is low and basically independent of $\mu$, except for a small rise at
the largest $\mu$. $\eta_{\text{slide}}$ decreases monotonically to zero with
increasing friction, while $\eta_{\text{lock}}$ increases monotonically and
eventually saturates, which demonstrates that with increasing $\mu$ the
contact dynamic crosses over from sliding to rolling.

Maxwell’s idea which relates the number of constraints to the static
mechanical rigidity of a structure can be generalized to the dynamic case by
defining for each contact in the pack a cycle-averaged number of constraints
as follows. During a cycle, one monitors the contact between particles $i$ and
$j$ and defines $C^{ij}(t)$, the associated number of constraints for this
contact by distinguishing the three cases introduced above: i) The contact
breaks, which corresponds to a normal force $|{\bf F}_{n}^{ij}|=0$ and hence
the contact is not constrained, i.e., $C^{ij}(t)=0$; ii) The contact has
$|{\bf F}_{n}^{ij}|>0$ with a tangential force $|{\bf F}_{t}^{ij}|>\mu|{\bf
F}_{n}^{ij}|$, thus is sliding, giving a constraint of 1; iii) If $|{\bf
F}_{t}^{ij}|<\mu|{\bf F}_{n}^{ij}|$ one has slipless rotation which
corresponds to two constraints. The cycle-averaged dynamic constraint of the
contact, $\overline{C^{ij}}$, is thus given by

$\overline{C^{ij}}=\frac{1}{T}\int_{0}^{T}C^{ij}(t)dt,\qquad{\rm where}\quad
C^{ij}(t)=\left\\{\begin{aligned} &0\quad{\rm if}\quad|{\bf
F}_{n}^{ij}(t)|=0\\\ &1\quad{\rm if}\quad|{\bf F}_{n}^{ij}(t)|>0\text{ and
}|{\bf F}_{t}^{ij}(t)|=\mu|{\bf F}_{n}^{ij}(t)|\\\ &2\quad{\rm if}\quad|{\bf
F}_{n}^{ij}(t)|>0\text{ and }|{\bf F}_{t}^{ij}(t)|<\mu|{\bf
F}_{n}^{ij}(t)|.\\\ \end{aligned}\right.$ (1)

($T$ is the cycle time and details on the calculation are in the Methods.) The
number of constraints acting on particle $i$ is given by
$\overline{C^{i}}=\sum_{j}\overline{C^{ij}}$ and this quantity will indicate
whether or not the motion of the particle is blocked. The system average of
$\overline{C^{i}}$, $\overline{C^{pp}}=N^{-1}\sum_{i=1}^{N}\overline{C^{i}}$,
gives the average number of constraints per particle and shows, Fig. 3(a),
that for small friction the number of constraints is just above 2.0, rises
gently towards 4.0 if $\mu$ approaches $\mu_{\rm min}$, before decreasing
again at the largest $\mu$. The $\mu$-dependence of $\overline{C^{pp}}$ can be
rationalized by making the mean-field approximation
$\overline{C^{pp}}\approx\overline{Z}\cdot\overline{C}/2$, where
$\overline{C}$ is the system- and time-averaged number of constraints per
contact, shown in Fig. 3(a). ($\overline{C^{pp}}$ is compared with mean-field
in the SI.) One observes that at small friction $\overline{C^{pp}}$ grows
because $\overline{C}$ increases, while $\overline{Z}$ stays constant. The
product reaches a maximum at around $\mu_{\rm min}$ and stays constant up to
around $\mu_{\rm max}$ since the slight decrease of $\overline{Z}$ is
compensated by the increase of $\overline{C}$. For $\mu>\mu_{\rm max}$,
$\overline{Z}$ drops quickly, i.e., the structure becomes more open, while
$\overline{C}$ remains constant, resulting in a quick drop of
$\overline{C^{pp}}$, ensuing an acceleration of the dynamics. This result is
thus strong evidence that the non-monotonic behavior of the relaxation
dynamics can be understood directly from the $\mu$-dependence of the number of
dynamic constraints on the particles.

The PDF of $\overline{C^{i}}$, Fig. 3(c), allows to get insight on its
average, $\overline{C^{pp}}$. For $\mu=0$ there are two narrow peaks at 4 and
5, corresponding to particles that have, respectively, four and five
frictionless contacts, i.e., $\overline{C^{ij}}=1$. For somewhat larger $\mu$,
they are absorbed into the broader peak which shifts to the right, i.e., most
of the contacts are at least temporarily locked. Upon approaching the glass-
region, $\mu=0.02$, the PDF starts to show again several narrow peaks,
$\overline{C^{ij}}\approx 7.8$ and 9.6, and an increase of $\mu$ to 0.3 makes
these two peaks sharp and move to 8 and 10, respectively, and three others
appear, at 4, 6, and 12, corresponding to particles that have 2, 3,…,6
neighbors with locked contacts, i.e., $\overline{C^{ij}}=2$. Note that the
peak at 4 is the last to emerge, at around $\mu_{\rm max}$ (Inset). This
demonstrates that for $\mu>\mu_{\rm max}$ a typical local structure of the
system is a linear chain in which each particle has exactly 2 permanently
locked contacts, see cartoon in Fig. 3(c), indicating that spatial
correlations of constrained contacts play an important role for the relaxation
dynamics and below we will see that this is indeed the case.

Spatial structure of constraints: We now probe the nature of the cooperative
motion, often observed in other glassy materials [29, 48], in our system and
relate it to the local structure and the constraints. For this we compare the
spatial map of $\overline{C^{ij}}$ with the one of $D_{\rm min}^{2}$, the non-
affine displacement of a particle with respect to its nearest neighbor shell
[49].

The upper left parts of Fig. 4 show the maps of $\overline{C^{ij}}$ for
increasing $\mu$, and Fig. 5(a) presents the cluster analysis. We find that
this structure transforms from a percolating $\overline{C^{ij}}\leq 1$ network
at $\mu=10^{-5}$, Fig. 4(a), to a percolating $1<\overline{C^{ij}}\leq 1.99$
one when $\mu$ approaches $\mu_{\rm min}$, Fig. 4(b), while contacts
$\overline{C^{ij}}>1.99$ are scarce. This demonstrates that the dynamic
slowing down is not related to completely locked contacts but instead to the
percolating network with intermediate $\overline{C^{ij}}$. Even deep in the
glass zone, Fig. 4(c), the fully locked contacts do not yet percolate, but
they do so at the highest $\mu$, Fig. 4(d). This percolating network consists
of chain-like structures pictured in Fig. 3(c), and therefore can deform more
easily in a collective manner.

The network of $\overline{C^{ij}}$ determines how the system can relax,
described by $D_{\rm min}^{2}$ in the lower right part of the panels. At small
$\mu$, Fig. 4(a), $D_{\rm min}^{2}$ is spatially very heterogeneous, since the
percolating network $\overline{C^{ij}}\leq 1$ constraints the dynamics only
weakly, thus giving rise to the heterogeneous dynamics found in other glass-
forming systems. For $\mu=10^{-2}$, Fig. 4(b), the dynamical heterogeneities
increase strongly, quantified via $\chi_{4}$ in Extended Data Fig. 3, and the
average magnitude of $D_{\rm min}^{2}$ decreases, since the dynamics is now
confined by the percolating $1<\overline{C^{ij}}\leq 1.99$ network. Inside the
glass zone, Fig. 4(c), the constraints are so strong that most of the
particles cannot change anymore their neighbors and hence the zones with
significant $D_{\rm min}^{2}$ become small islands. If $\mu$ exceeds $\mu_{\rm
max}$, Fig. 4(d), zones with large $D_{\rm min}^{2}$ reappear, but in contrast
to the case of small $\mu$, they form a multitude of small clusters, which are
embedded in the highly constrained contact network.

Dynamic constraints in configuration space: The influence of the constraints
on the dynamics can be understood via the propagation of the system in
configuration space, illustrated by the cartoons on the upper right corner of
the panels in Fig. 4. At very low $\mu$, the system can freely explore the
landscape between blue stable areas since there are only very few constraints.
With increasing $\mu<\mu_{\rm min}$, the exploration is hindered by barriers
from increasing frictional constraints, red lines, therefore the motion
becomes more collective and slower.

Inside the glass region, the system is likely to occupy places surrounded by
many more barriers, therefore the relaxation ceases. To accommodate the
strain, the system forms narrow micro-shear bands, mainly oriented
horizontally and vertically, in which the contacts slide because they have a
small $\overline{C^{ij}}$, Fig. 5(b). The particles inside these bands have a
large $D_{\rm min}^{2}$ at maximal strain, but most motion is reversible.
Thus, although trapped near a local minimum of the free energy landscape, the
system has self-organized its structure permitting it to have a high
elasticity.

When $\mu>\mu_{max}$, the number of stable structures increases, since
friction lowers the minimal coordination number, allowing new relaxation
channels to circumvent the barriers. This is accomplished by the deformation
of the open highly-constrained network, which allows the particles enclosed to
undergo a non-affine motion to relax. Some of these cells in the network
become unstable during the cycle, i.e., the local bubble bursts, leading
ultimately to the relaxation of the entire network.

The evolution of the structure in configuration space demonstrates the dual
role of friction: Increasing $\mu$ results in a higher $\overline{C}$ which
means more dynamic barriers, thus slower relaxation. At around $\mu_{\rm max}$
friction allows the creation of many stable states with low $\overline{Z}$,
facilitating relaxation via new pathways of deformation and reconstruction of
the network.

Outlook: Our results demonstrate that the number of dynamic constraints is a
good indicator of the relaxation dynamics, and provide more insight than the
often studied force chain network, see Extended Data Figs. 4, 5. This implies
that in frictional systems the dynamics is governed by the distance of a
contact from its Coulomb boundary and not by the magnitude of its force, a
conclusion that holds for all materials, irrespective of their friction
coefficients (hydrogels, grains, gears particles,..), thus opening the door to
tune the relaxation dynamics by controlling the friction coefficient via
methods like surface chemical treatment [50]. Our approach allows us to treat
the dynamics of granular material and frictionless hard-sphere glass formers
within the same framework, which will permit to unify the physics of various
disordered systems. The approach of dynamic constraints to study dynamics is
very general and hence it can be adapted to other soft matter systems like
suspensions, deformable cells, and biological tissues [41, 43, 44]. As many
systems (in and out of equilibrium) have complex free energy landscapes, and
often have a dynamics subject to time-varying constraints, applications in
domains like the traveling salesman problem, financial markets, and ecology
network [45, 42], are expected to be feasible.

Methods:

Model and Simulations: The contact force between particles is Hertzian, and
the tangential component is subject to a Coulomb cutoff [51, 52]. Thus the
total force between two particles $i$ and $j$ is given by

${\bf F}(r)=\sqrt{R_{i}R_{j}/(R_{i}+R_{j})}\sqrt{\delta}\left[(k_{n}\delta-
m_{\rm eff}\gamma_{n}v_{n}){\bf n}-(k_{t}\Delta s_{t}+m_{\rm
eff}\gamma_{t}v_{t}){\bf t}\right].$

Here ${\bf n}$ and ${\bf t}$ are, respectively, the normal and tangent vectors
at the contact, and the overlap between the particles is given by
$\delta=R_{i}+R_{j}-r$, where $R_{i}$ is the radius of the particle. $v_{n}$
and $v_{t}$ are the relative speed of the particles, projected on ${\bf n}$
and ${\bf t}$, respectively, and $m_{\rm eff}$ is the effective mass. $\Delta
s_{t}$ is the displacement between two particles in the tangential direction.
The absolute value of $k_{t}\Delta s_{t}$ is subject to cutoff at $\mu
k_{n}\delta$. We parameters we use are $k_{n}=2\cdot
10^{8}\text{g}\cdot\text{cm}^{-1}\cdot\text{s}^{-2}$,
$\gamma_{n}=11,900\text{cm}^{-1}\cdot\text{s}^{-1}$, $k_{t}=\frac{2}{7}k_{n}$,
and $\gamma_{t}=\frac{1}{2}\gamma_{n}$. Simulations were carried out using the
LAMMPS software [53] with a time step size of $4\cdot 10^{-5}$s. The
simulation shear box, shown in Inset of Fig. 1(b), consists of four amorphous
walls to avoid heterogeneous crystallization of the sample from the boundary
and also to weaken the global convection effect that can occur along smooth
boundaries. Similar to an experimental setup for sheared systems [21], the
bottom wall (pink) moves sinusoidally in the horizontal direction with a shear
amplitude of $\gamma_{0}=0.05$ and a strain rate of $\dot{\gamma}=\text{0.1
s}^{-1}$, while the left and the right walls (blue) move accordingly. The top
wall (green) applies a pressure of $P=98100\text{g}\cdot\text{s}^{-2}$ on the
pack and it is constrained to move only vertically. This gives an inertia
number $I=\dot{\gamma}\bar{d}/\sqrt{P/\rho}$ is $3.6\cdot 10^{-4}$, where
$\bar{d}$ is the mean diameter of the particles, i.e., we are in the quasi-
static shear regime [3]. Typically we used for each values of $\mu$ 1-4
independent samples and did up to $9\cdot 10^{5}$ cycles to reach the steady
state. See Extended Data Table I for details.

Defining locked contacts: In practice we define a contact to be non-slipping,
i.e., “locked”, if the tangential force is inside the Coulomb-cone within a
relative accuracy $\varepsilon$, i.e., $|{\bf F}_{t}|<\mu(1-\varepsilon)|{\bf
F}_{n}|$, where we have chosen $\varepsilon=10^{-6}$. To improve the
statistics for obtaining the results for the probabilities shown in Fig. 3(b)
we averaged over 700 cycles. To calculate the integral in Eq. 1 we divided the
cycle into 100 equal intervals.

## References

  * [1] Bruno Andreotti, Yoel Forterre and Olivier Pouliquen  “Granular Media: Between Fluid and Solid”  Cambridge: Cambridge University Press, 2013  DOI: 10.1017/CBO9781139541008
  * [2] Scott V. Franklin and Mark D. Shattuck  “Handbook of Granular Materials” Google-Books-ID: EdOYCgAAQBAJ  CRC Press, 2016
  * [3] G… MiDi  “On dense granular flows”  In _European physical journal. E_ 14.4, 2004, pp. 341–65  DOI: 10.1140/epje/i2003-10153-0
  * [4] Gernot Michlmayr, Denis Cohen and Dani Or  “Sources and characteristics of acoustic emissions from mechanically stressed geologic granular media — A review”  In _Earth-Science Reviews_ 112.3, 2012, pp. 97–114  DOI: 10.1016/j.earscirev.2012.02.009
  * [5] Robert P. Behringer and Bulbul Chakraborty  “The physics of jamming for granular materials: A review” Publisher: IOP Publishing  In _Reports on Progress in Physics_ 82.1, 2018, pp. 012601  DOI: 10.1088/1361-6633/aadc3c
  * [6] J. Maxwell  “L. On the calculation of the equilibrium and stiffness of frames” Publisher: Taylor & Francis _eprint: https://doi.org/10.1080/14786446408643668  In _The London, Edinburgh, and Dublin Philosophical Magazine and Journal of Science_ 27.182, 1864, pp. 294–299  DOI: 10.1080/14786446408643668
  * [7] J.C. Phillips and M.F. Thorpe  “Constraint theory, vector percolation and glass formation”  In _Solid State Communications_ 53.8, 1985, pp. 699–702  DOI: 10.1016/0038-1098(85)90381-3
  * [8] M. Wyart  “On the rigidity of amorphous solids” Number: 3 Publisher: EDP Sciences  In _Annales de Physique_ 30.3, 2005, pp. 1–96  DOI: 10.1051/anphys:2006003
  * [9] Alessio Zaccone and Enzo Scossa-Romano  “Approximate analytical description of the nonaffine response of amorphous solids” Publisher: American Physical Society  In _Physical Review B_ 83.18, 2011, pp. 184205  DOI: 10.1103/PhysRevB.83.184205
  * [10] Kostya Shundyak, Martin Hecke and Wim Saarloos  “Force mobilization and generalized isostaticity in jammed packings of frictional grains”  In _Physical Review E_ 75.1, 2007, pp. 010301  DOI: 10.1103/PhysRevE.75.010301
  * [11] Jacco H. Snoeijer, Thijs J.. Vlugt, Martin Hecke and Wim Saarloos  “Force Network Ensemble: A New Approach to Static Granular Matter” Publisher: American Physical Society  In _Physical Review Letters_ 92.5, 2004, pp. 054302  DOI: 10.1103/PhysRevLett.92.054302
  * [12] L. Kondic et al.  “Topology of force networks in compressed granular media”  In _Europhysics Letters_ 97.5, 2012, pp. 54001  DOI: 10.1209/0295-5075/97/54001
  * [13] Aghil Abed Zadeh et al.  “Enlightening force chains: A review of photoelasticimetry in granular matter”  In _Granular Matter_ 21.4, 2019, pp. 83  DOI: 10/gng6zq
  * [14] Dapeng Bi, Jie Zhang, Bulbul Chakraborty and R.. Behringer  “Jamming by shear” Number: 7377 Publisher: Nature Publishing Group  In _Nature_ 480.7377, 2011, pp. 355–358  DOI: 10.1038/nature10667
  * [15] S.. Edwards and R… Oakeshott  “Theory of powders”  In _Physica A Stat. Mech. Appl._ 157.3, 1989, pp. 1080–1090  DOI: 10.1016/0378-4371(89)90034-4
  * [16] Adrian Baule, Flaviano Morone, Hans J. Herrmann and Hernán A. Makse  “Edwards statistical mechanics for jammed granular matter” Publisher: American Physical Society  In _Reviews of Modern Physics_ 90.1, 2018, pp. 015006  DOI: 10/gc5skk
  * [17] Ye Yuan et al.  “Experimental Test of the Edwards Volume Ensemble for Tapped Granular Packings” Publisher: American Physical Society  In _Physical Review Letters_ 127.1, 2021, pp. 018002  DOI: 10.1103/PhysRevLett.127.018002
  * [18] Yi Xing et al.  “Origin of the critical state in sheared granular materials”  In _Nature Physics_ , 2024  DOI: 10.1038/s41567-023-02353-4
  * [19] Leonardo E. Silbert et al.  “Geometry of frictionless and frictional sphere packings” Publisher: American Physical Society  In _Physical Review E_ 65.3, 2002, pp. 031304  DOI: 10.1103/PhysRevE.65.031304
  * [20] Antoinette Tordesillas et al.  “Transition dynamics and magic-number-like behavior of frictional granular clusters” Publisher: American Physical Society  In _Physical Review E_ 86.1, 2012, pp. 011306  DOI: 10.1103/PhysRevE.86.011306
  * [21] Binquan Kou et al.  “Granular materials flow like complex fluids”  In _Nature_ 551.7680, 2017, pp. 360–363  DOI: 10.1038/nature24062
  * [22] Yiqiu Zhao et al.  “Ultrastable Shear-Jammed Granular Material” Publisher: American Physical Society  In _Physical Review X_ 12.3, 2022, pp. 031021  DOI: 10.1103/PhysRevX.12.031021
  * [23] John R. Royer and Paul M. Chaikin  “Precisely cyclic sand: Self-organization of periodically sheared frictional grains” Publisher: National Academy of Sciences Section: Physical Sciences  In _Proceedings of the National Academy of Sciences_ 112.1, 2015, pp. 49–53  DOI: 10.1073/pnas.1413468112
  * [24] E. DeGiuli, J.. McElwaine and M. Wyart  “Phase diagram for inertial granular flows”  In _Physical Review E_ 94.1, 2016, pp. 012904  DOI: 10.1103/PhysRevE.94.012904
  * [25] Michio Otsuki and Hisao Hayakawa  “Shear modulus and reversible particle trajectories of frictional granular materials under oscillatory shear” arXiv: 2103.14457 version: 2  In _The European Physical Journal E_ 44.5, 2021, pp. 70  DOI: 10.1140/epje/s10189-021-00075-0
  * [26] Michael Sakellariou  “Granular Materials”  InTechOpen, 2017  DOI: 10.5772/66007
  * [27] Anaël Lemaître et al.  “Frictional Granular Matter: Protocol Dependence of Mechanical Properties” Publisher: American Physical Society  In _Physical Review Letters_ 126.7, 2021, pp. 075501  DOI: 10.1103/PhysRevLett.126.075501
  * [28] Yoël Forterre and Olivier Pouliquen  “Flows of Dense Granular Media” _eprint: https://doi.org/10.1146/annurev.fluid.40.111406.102142  In _Annual Review of Fluid Mechanics_ 40.1, 2008, pp. 1–24  DOI: 10.1146/annurev.fluid.40.111406.102142
  * [29] Kurt Binder and Walter Kob  “Glassy Materials and Disordered Solids: An Introduction to Their Statistical Mechanics” Google-Books-ID: lEOWIZeFSP4C  World Scientific, 2011
  * [30] Wolfgang Götze  “Complex Dynamics of Glass-Forming Liquids: A Mode-Coupling Theory”  Oxford University Press, 2008  DOI: 10.1093/acprof:oso/9780199235346.001.0001
  * [31] Vassiliy Lubchenko and Peter G. Wolynes  “Theory of Structural Glasses and Supercooled Liquids” _eprint: https://doi.org/10.1146/annurev.physchem.58.032806.104653  In _Annual Review of Physical Chemistry_ 58.1, 2007, pp. 235–266  DOI: 10.1146/annurev.physchem.58.032806.104653
  * [32] Jeppe C. Dyre  “Colloquium: The glass transition and elastic models of glass-forming liquids” Publisher: American Physical Society  In _Reviews of Modern Physics_ 78.3, 2006, pp. 953–972  DOI: 10.1103/RevModPhys.78.953
  * [33] Dietrich E. Wolf, Farhang Radjai and Sabine Dippel  “Dissipation in granular materials” Publisher: Taylor & Francis _eprint: https://doi.org/10.1080/13642819808205033  In _Philosophical Magazine B_ 77.5, 1998, pp. 1413–1425  DOI: 10.1080/13642819808205033
  * [34] Anton Peshkov, Michelle Girvan, Derek C. Richardson and Wolfgang Losert  “Reversibility of granular rotations and translations”  In _Physical Review E_ 100.4, 2019, pp. 042905  DOI: 10/gn8t38
  * [35] Zackery A. Benson et al.  “Experimentally Measuring Rolling and Sliding in Three-Dimensional Dense Granular Packings”  In _Physical Review Letters_ 129.4, 2022, pp. 048001  DOI: 10.1103/PhysRevLett.129.048001
  * [36] Glenn H. Fredrickson and Hans C. Andersen  “Kinetic Ising Model of the Glass Transition” Publisher: American Physical Society  In _Physical Review Letters_ 53.13, 1984, pp. 1244–1247  DOI: 10.1103/PhysRevLett.53.1244
  * [37] Walter Kob and Hans C. Andersen  “Kinetic lattice-gas model of cage effects in high-density liquids and a test of mode-coupling theory of the ideal-glass transition” Publisher: American Physical Society  In _Physical Review E_ 48.6, 1993, pp. 4364–4377  DOI: 10.1103/PhysRevE.48.4364
  * [38] F. Ritort and P. Sollich  “Glassy dynamics of kinetically constrained models” Publisher: Taylor & Francis _eprint: https://doi.org/10.1080/0001873031000093582  In _Advances in Physics_ 52.4, 2003, pp. 219–342  DOI: 10.1080/0001873031000093582
  * [39] Kuang Liu et al.  “Spongelike Rigid Structures in Frictional Granular Packings” Publisher: American Physical Society  In _Physical Review Letters_ 126.8, 2021, pp. 088002  DOI: 10.1103/PhysRevLett.126.088002
  * [40] Varghese Babu, H.. Vinutha, Dapeng Bi and Srikanth Sastry  “Discontinuous rigidity transition associated with shear jamming in granular simulations” Publisher: The Royal Society of Chemistry  In _Soft Matter_ 19.48, 2023, pp. 9399–9404  DOI: 10.1039/D3SM00725A
  * [41] Michael Naald et al.  “Minimally rigid clusters in dense suspension flow” Publisher: Nature Publishing Group  In _Nature Physics_ , 2024, pp. 1–7  DOI: 10.1038/s41567-023-02354-3
  * [42] Patrick Charbonneau et al.  “Spin Glass Theory And Far Beyond: Replica Symmetry Breaking After 40 Years” Google-Books-ID: gIrVEAAAQBAJ  World Scientific, 2023
  * [43] Edouard Hannezo and Carl-Philipp Heisenberg  “Rigidity transitions in development and disease”  In _Trends in Cell Biology_ 32.5, 2022, pp. 433–444  DOI: 10.1016/j.tcb.2021.12.006
  * [44] Lior Atia et al.  “Geometric constraints during epithelial jamming” Publisher: Nature Publishing Group  In _Nature Physics_ 14.6, 2018, pp. 613–620  DOI: 10.1038/s41567-018-0089-9
  * [45] Louis-Philippe Bigras, Michel Gamache and Gilles Savard  “The time-dependent traveling salesman problem and single machine scheduling problems with sequence dependent setup times”  In _Discrete Optimization_ 5.4, 2008, pp. 685–699  DOI: 10.1016/j.disopt.2008.04.001
  * [46] Yiqiu Zhao et al.  “Microscopic reversibility and emergent elasticity in ultrastable granular systems”  In _Frontiers in Physics_ 10, 2022  DOI: 10.3389/fphy.2022.1048683
  * [47] Elijah Flenner and Grzegorz Szamel  “Fundamental differences between glassy dynamics in two and three dimensions” Publisher: Nature Publishing Group  In _Nature Communications_ 6.1, 2015, pp. 7392  DOI: 10.1038/ncomms8392
  * [48] L. Berthier et al.  “Dynamical Heterogeneities in Glasses, Colloids, and Granular Media”, International Series of Monographs on Physics  OUP Oxford, 2011  URL: https://books.google.com/books?id=KWOS_qvgOFsC
  * [49] M.. Falk and J.. Langer  “Dynamics of viscoplastic deformation in amorphous solids” Publisher: American Physical Society  In _Physical Review E_ 57.6, 1998, pp. 7192–7205  DOI: 10.1103/PhysRevE.57.7192
  * [50] Hojin Kim et al.  “Dynamic-bond-induced sticky friction tailors non-Newtonian rheology”  In _Soft Matter_ 19 The Royal Society of Chemistry, 2023, pp. 6797–6804  DOI: 10.1039/D3SM00479A
  * [51] Corey S. O’Hern, Stephen A. Langer, Andrea J. Liu and Sidney R. Nagel  “Random Packings of Frictionless Particles”  In _Physical Review Letters_ 88.7, 2002, pp. 075507  DOI: 10.1103/PhysRevLett.88.075507
  * [52] Jonathan Barés et al.  “Local and global avalanches in a two-dimensional sheared granular medium” Publisher: American Physical Society  In _Physical Review E_ 96.5, 2017, pp. 052902  DOI: 10.1103/PhysRevE.96.052902
  * [53] A.. Thompson et al.  “LAMMPS - a flexible simulation tool for particle-based materials modeling at the atomic, meso, and continuum scales”  In _Comp. Phys. Comm._ 271, 2022, pp. 108171  DOI: 10.1016/j.cpc.2021.108171

Figure 1:  Translational and rotational dynamics: Time dependence of the
translational mean squared displacement (TMSD), (a), and of the rotational one
(RMSD), (b), for different friction constants $\mu$ (given in the legend). The
dashed and dotted lines are power-laws with exponent 1.0 and 2.0,
respectively. The RMSD is obtained by time-integrating the angular
displacements during the propagation of the system. Note that the value of
these MSD’s at $t=1$ is non-monotonic as a function of $\mu$, i.e., the
dynamics is re-entrant. Depending on the value of $\mu$, the long time
dynamics can be diffusive, super-diffusive, or caged. The Inset in (b) shows
the setup of the simulation. See Methods for details.  Figure 2:  Relaxation
times for the translational and rotational degrees of freedom as a function of
friction coefficient: (a) Slope of the TMSD and RMSD at long times, obtained
from the last three decades of the curves in Fig. 1. Error-bars indicate 95%
confidence interval. (b) Value of TMSD and RMSD at $t=1$. The straight line
indicates a power-law with exponent -2. Error-bars indicate the standard
deviation. (c) $\mu$-dependence of various relaxation times:
$\tau_{\text{cont}}$, the time for a particle to disconnect from a neighbor
with which it was in contact at $t=0$; $\tau_{\text{vor}}$, the characteristic
time for a particle to have a change in its (radical) Voronoi neighborhood;
$\tau_{\text{qs}}$, the decay-time of the self-overlap function (threshold
0.3$d_{s}$), and $\tau_{\text{TMSD}}$, the time for the TMSD to reach $0.1$.
The $t-$dependence of the corresponding time correlation functions is shown in
Extended Data Fig. 1. Error-bars indicate 95% confidence interval. Left Inset:
Same data for small $\mu$ in a log-lin plot. The gray straight lines are
exponentials with a slope 200. Right Inset: Same data for large $\mu$ as a
function of $\mu-\mu_{\rm max}$. The straight lines have slope -3.0. The
vertical red lines in the three panels indicate $\mu_{\rm min}$ and $\mu_{\rm
max}$.

Figure 3: Distributions of averaged number of constraints as a function of
friction coefficient: (a) Effective number of constraints per contact,
$\overline{C}$, number of constraints per particle, $\overline{C^{pp}}$, and
the averaged number of contacts, $\overline{Z}$. Solid lines show the time-
averaged results, while dashed lines and dotted lines show results at
$\gamma=0$ and $\gamma=0.05$, respectively. Error-bars indicate the standard
deviation. (b) Cartoons show three different states of a contact: Sliding,
breaking, locked. The curves show the probability that during a cycle a
contact is in a given state. (c) Distribution of the time-averaged number of
constraints per particle, $\overline{C^{i}}$. Inset: Zoom on the dotted
rectangular region. Cartoon: Chain-like motion corresponding to the peak at
$\overline{C^{i}}=4$. Figure 4:  Spatial structure of non-affine displacement
and of dynamic constraints on contacts for different $\mu$. Upper left part of
the panels show the network of cycle-averaged contacts, $\overline{C^{ij}}$,
with the color code given at the bottom left. Lower right part of the panels
show the non-affine displacement after one cycle, grey scale code at bottom
right. See main text for a discussion of these panels. Insets show
schematically the motion of the system in configuration space. Dark blue:
Mechanically stable regions where the particles have typically a coordination
number of 4. Light blue: regions with typically coordination number 3.
Orange/green arrows correspond to the forward/backward motion of the cycle.
Red lines are barriers that are incompatible with the dynamic constraints and
hence the system avoids them. Insets: (a) Many stable regions have
coordination number 4 and the absence of constraints allows the system to move
easily in configuration space. (b) The constraints make that the motion
becomes hindered resulting to a slowed down relaxation dynamics. (c) The
number of constraints is so high that motion is reversible and relaxation
suppressed. Formation of first zones with coordination number 3. (d) High
friction allows the system to access mechanically stable state in which
particles have a typical coordination number of 3. These states allow the
system to explore new pathways and hence to relax.

Figure 5: Percolation probability and length scales of the contact network.
(a): The number of contact in the largest connected cluster of
$\overline{C^{ij}}$ (normalized by the total number of contacts), full lines,
the extension of this cluster in the $x$-direction, dashed lines, and (Inset)
the fractal dimension of this cluster, dotted lines. The three colors
correspond to the three ranges of contacts given in the legend. Error-bars
indicate the standard deviation. (b): Network of contacts and non-affine
displacement field at maximal strain; $\mu=1.4\cdot 10^{-1}$. The micro-shear
bands where $D^{2}_{\rm min}$ is high after a quarter cycle (yellow boxes) are
basically reversible once the shear cycle is completed, see Fig. 4(c). Inset:
Zoom on the white square area.

Extended Data:

$\mu$ | cycles($10^{4}$) | $\mu$ | cycles($10^{4}$) | $\mu$ | cycles($10^{4}$) | $\mu$ | cycles($10^{4}$)
---|---|---|---|---|---|---|---
0 | 6 | 0.01 | 14, 48 | 0.1 | 6, 6 | 0.26 | 48, 72
$10^{-8}$ | 6 | 0.012 | 30 | 0.12 | 6 | 0.28 | 12, 48
$10^{-7}$ | 6 | 0.013 | 36, 84 | 0.14 | 6 | 0.3 | 12, 72, 84
$10^{-6}$ | 6 | 0.014 | 24, 60 | 0.16 | 24,24 | 0.32 | 6
$10^{-5}$ | 6 | 0.015 | 66, 72 | 0.18 | 18, 24 | 0.35 | 6
0.0001 | 6 | 0.016 | 24, 30 | 0.2 | 23, 24, 24 | 0.4 | 6
0.001 | 6, 6 | 0.018 | 18 | 0.21 | 24, 60 | 0.5 | 6
0.002 | 18, 30 | 0.02 | 2, 12, 18, 54 | 0.22 | 24, 72 | 0.7 | 6
0.003 | 6, 12 | 0.03 | 2 | 0.23 | 48, 90 | 1 | 6
0.005 | 16, 24 | 0.05 | 2 | 0.24 | 30, 60, 66 | 3 | 6
0.008 | 6 | 0.07 | 2, 6 | 0.25 | 54 | 10 | 6
Table 1: Table of simulation cycle numbers for each $\mu$: Each number
represent a run with independent random initial condition.

Extended Data Fig. 1: Time dependence of persistence probabilities: (a) Time
dependence of the self-overlap (defined using a threshold of 0.3$d_{s}$). (b)
Probability that a particle that at $t=0$ was in a given Voronoi cell, is
still in the same cell at time $t$. (c) The probability that at time $t$ a
contact that was present at time $t=0$ is still intact.  Extended Data Fig. 2:
Packing fraction, concentration of rattlers, relaxation times: (a) Packing
fraction as a function of $\mu$. $\phi$ decreases continuously with increasing
$\mu$, making it unlikely that the re-entrant dynamics is only an effect of
the packing fraction. At the upper bound of the glass-regime, i.e., $\mu_{\rm
max}$, $\phi$ drops quickly because the increased friction allows to stabilize
states in which the particles have a low coordination number. (b) Probability
that a particle is a rattler, i.e, particles with fewer than 2 contacts. This
probability is low ($\approx 1$%) at small and intermediate fraction, i.e., as
long as the typical coordination number is larger than 4. At high $\mu$ this
probability grows significantly due to the formation of the open network
structure discussed in the main text. Error-bars in (a) and (b) indicate the
standard deviation. (c) $\mu$-dependence of the relaxation times, presented in
Fig. 1 of the main text, normalized by $\tau_{\rm cont}$, the time it takes
for a particle to break a contact. Error-bars indicate 95% confidence
interval. Extended Data Fig. 3: Time dependence of $\chi_{4}$ for different
values of $\mu$: The dynamic susceptibility (self part) is used to quantify
the strength of the dynamical heterogeneities in the system and it has been
calculated from the variance of the overlap function. $\chi_{4}$ is relatively
small for small friction and increases rapidly if $\mu$ approaches the glass-
regime. Inside the glass-regime $\chi_{4}$ is basically zero, since the system
hardly relaxes. For $\mu\gtrsim\mu_{\rm max}$ it is again large and decreases
steadily with increasing friction. At the same time also the time scale at
which $\chi_{4}(t)$ peaks is non-monotonic in $\mu$. These results demonstrate
that the system has pronounced dynamical heterogeneities that peak at the
boundary of the glass-zone.  Extended Data Fig. 4: Probability distribution
function of the forces: Panels (a) and (b) show the PDF of forces at strain
$\gamma=0$ and panels (c) and (d) shows the ones at $\gamma=0.05$. The PDF’s
for ${\bf F}_{n}$ show a Gaussian decay for all values of $\mu$. The width of
the Gaussian increases in a monotonic manner as a function of $\mu$ at
$\gamma=0.05$, and fall on a master curve for $\gamma=0$. The PDFs of ${\bf
F}_{t}$ follow an exponential distribution with an increasing width. These
results indicate that the forces are not able to rationalize the re-entrant
dynamics of the system.

Extended Data Fig. 5: Network of the normal forces at the contacts: The upper
left part of the panels shows the magnitude of the normal component of the
force at the contacts while the lower right part shows the non-affine
displacement of the particles. Scales are given at the bottom of the figure.
Panels (a)-(d) are at strain $\gamma=0$ for different $\mu$ and panels (e)-(f)
at strain $\gamma=0.05$. No evident non-monotonic evolution is seen in the
force chains when $\mu$ is increased.

Supplementary Information

Mean field approximation Here we present the comparison between the mean field
result and the one of the simulation.

The number of constraints related to a contact between particle $i$ and $j$ at
time $t$ can be written as

$C^{ij}(t)=1\cdot\vmathbb{1}^{ij}(t)+2\cdot\vmathbb{2}^{ij}(t)\quad,$

where $\vmathbb{1}^{ij}(t)$ and $\vmathbb{2}^{ij}(t)$ are functions describing
the state of this contact at time $t$: $\vmathbb{1}^{ij}(t)=1$ and
$\vmathbb{2}^{ij}(t)=0$ if the contact is sliding; $\vmathbb{1}^{ij}(t)=0$ and
$\vmathbb{2}^{ij}(t)=1$ if the contact is locked; $\vmathbb{1}^{ij}(t)=0$ and
$\vmathbb{2}^{ij}(t)=0$ if the contact is broken.

The time-averaged constraint of this contact is

$\overline{C^{ij}}=\frac{1}{T}\int^{T}_{0}C^{ij}(t)dt\quad.$

We define the time-averaged constraint per particle as

$\overline{C^{pp}}=\frac{1}{2N}\sum_{i,j}\overline{C^{ij}}=\frac{1}{2NT}\sum_{i,j}\int^{T}_{0}(\vmathbb{1}^{ij}(t)+2\cdot\vmathbb{2}^{ij}(t))dt\quad,$

where $N$ is the particle number.

The system averaged contact number per particle at time $t$, i.e., the average
coordination number of a particle, is given by

$Z(t)=\frac{1}{N}\sum_{i,j}\vmathbb{1}^{ij}(t)+\vmathbb{2}^{ij}(t)\quad,$

and thus the time-averaged number of contacts is

$\overline{Z}=\frac{1}{NT}\sum_{i,j}\int^{T}_{0}(\vmathbb{1}^{ij}(t)+\vmathbb{2}^{ij}(t))dt\quad.$

Now we define the effective number of constraints per contact at time $t$ as

$C(t)=\frac{\sum_{i,j}\vmathbb{1}^{ij}(t)+2\cdot\vmathbb{2}^{ij}(t)}{\sum_{i,j}\vmathbb{1}^{ij}(t)+\vmathbb{2}^{ij}(t)}\quad,$

where the numerator is two times the total number of contacts at time $t$. The
$\mu-$dependence of the ensemble averaged value of $C(t)$ at $t=0$ and
$t=T/4$, which corresponds to strain $\gamma=0$ and 0.05, respectively, are
shown in Fig. 2a as a dotted and a dashed line, respectively.

The time-averaged effective number of constraints per contact is thus

$\overline{C}=\frac{1}{T}\int^{T}_{0}\frac{\sum_{i,j}\vmathbb{1}^{ij}(t)+2\cdot\vmathbb{2}^{ij}(t)}{\sum_{i,j}\vmathbb{1}^{ij}(t)+\vmathbb{2}^{ij}(t)}dt\quad.$

We now consider different way of averaging, by first calculating the averaged
number of total number of constraints in a cycle, and the averaged number of
total contact time in a cycle, and subsequently their ratio:

$\overline{C^{\prime}}=\frac{\frac{1}{T}\sum_{i,j}\int^{T}_{0}(\vmathbb{1}^{ij}(t)+2\cdot\vmathbb{2}^{ij}(t))dt}{\frac{1}{T}\sum_{i,j}\int^{T}_{0}(\vmathbb{1}^{ij}(t)+\vmathbb{2}^{ij}(t))dt}\quad.$

Note that

$\frac{1}{2}\overline{C^{\prime}}\cdot\overline{Z}=\overline{C^{pp}}\quad.$

Now we compare the difference between $\overline{C}$ and
$\overline{C^{\prime}}$. The $\mu-$dependence of the ensemble average of
$Z(t)$ at $=0$ and $t=T/4$, corresponding, respectively, to $\gamma=0$ and
$0.05$, are shown in Fig. 2a, dotted and dashed line, respectively. Note that
these two curves are basically identical, i.e., the $t$-dependence of $Z(t)$
is weak. This suggests that one can make the approximation

$\overline{Z}\approx Z(t)\quad.$

We have

$\overline{C}=\frac{1}{T}\int^{T}_{0}\frac{\sum_{i,j}\vmathbb{1}^{ij}(t)+2\cdot\vmathbb{2}^{ij}(t)}{NZ(t)}dt\approx\frac{1}{NT\overline{Z}}\int^{T}_{0}\sum_{i,j}(\vmathbb{1}^{ij}(t)+2\cdot\vmathbb{2}^{ij}(t))dt$

and the last expression is just $\overline{C^{\prime}}$, i.e., one has

$\overline{C}\approx\overline{C^{\prime}}$

Thus we conclude that

$\overline{C^{pp}}\approx\frac{1}{2}\overline{C}\cdot\overline{Z}.$