Quantitative Measure of Memory Loss in Complex Spatio-Temporal Systems

TL;DR

This paper introduces a practical measure for quantifying memory loss in complex systems, applies it to granular shear systems, and reveals irregular memory loss behavior linked to force network transitions and bifurcations.

Contribution

It proposes a new method to assess mixing rates in dynamical systems and demonstrates its application to granular shear systems, uncovering irregular memory loss patterns.

Findings

Memory loss rate varies erratically with shear rate.

Memory loss is linked to force network transitions.

Bifurcations influence the rate of memory loss.

Abstract

To make progress in understanding the issue of memory loss and history dependence in evolving complex systems, we consider the mixing rate that specifies how fast the future states become independent of the initial condition. We propose a simple measure for assessing the mixing rate that can be directly applied to experimental data observed in any metric space $X$. For a compact phase space $X \subset R^M$, we prove the following statement. If the underlying dynamical system has a unique physical measure and its dynamics is strongly mixing with respect to this measure, then our method provides an upper bound of the mixing rate. We employ our method to analyze memory loss for the system of slowly sheared granular particles with a small inertial number $I$. The shear is induced by the moving walls as well as by the linear motion of the support surface that ensures approximately linear shear throughout the sample. We show that even if $I$ is kept fixed, the rate of memory loss (considered at the time scale given by the inverse shear rate) depends erratically on the shear rate. Our study suggests a presence of bifurcations at which the rate of memory loss increases with the shear rate while it decreases away from these points. We also find that the memory loss is not a smooth process. Its rate is closely related to frequency of the sudden transitions of the force network. The loss of memory, quantified by observing evolution of force networks, is found to be correlated with the loss of correlation of shear stress measured on the system scale. Thus, we have established a direct link between the evolution of force networks and macroscopic properties of the considered system.

Equations

$$||F_m-F|{|}_{\infty }:=\underset{s\in \mathbb{R}}{\sup }|F_m(s)-F(s)|<C{2}^{-m},$$

$$|\mu (E)\mu (T)-\mu (f^{-m}(E)\cap T)|<{\epsilon }_m\mu (T),$$

$$||F_m-F|{|}_{\infty }<C{\epsilon }_m.$$

$${\mathbf{F}}_{i,j}^n=k_n{x}_{i,j}\mathbf{n}-{\eta }_n\overline{m}{\mathbf{v}}_{i,j}^n$$

$${\mathbf{F}}^t=\min (\mu |{\mathbf{F}}^n|,|{\mathbf{F}}^{t\ast }|){\mathbf{F}}^{t\ast }/|{\mathbf{F}}^{t\ast }|,$$

$${\mathbf{F}}^{t\ast }=-k_t\mathbf{\xi }-{\eta }_{t}m{\mathbf{v}}_{i,j}^t.$$

$$-{\mu }_b|\mathbf{g}|\frac{{\mathbf{a}}_i(t)}{|{\mathbf{a}}_i(t)|}$$

$${\mathbf{a}}_i=\frac{1}{m_i}\sum _{c_i}{\mathbf{F}}_{i,c_i}$$

$$\mathcal{A}(x,s,\epsilon ):=\left\{y\in X:s-\epsilon \le d(x,y)<s+\epsilon \right\},$$

$$\mu (\mathcal{A}(x,s,\epsilon ))\le K\epsilon .$$

$$V(R,\epsilon )=K_M((R+\epsilon {)}^M-(R-\epsilon {)}^M)=2K_M\epsilon \left(\left(\genfrac{}{}{0ex}{}{N}{1}\right)R^{N-1}+\left(\genfrac{}{}{0ex}{}{N}{3}\right)R^{N-3}{\epsilon }^2+\dots +{\epsilon }^{N-1}\right),$$

$$|\mu (E)\mu (T)-\mu (f^{-m}(E)\cap T)|<{\epsilon }_m\mu (T),$$

$$||F_m-F|{|}_{\infty }<C{\epsilon }_m.$$

$$F_m(s)=\sum _{T\in {\mathcal{T}}_m}\mu (\left\{x\in T:d(x,f^m(x))<s\right\})$$

$$F(s)=\sum _{T\in {\mathcal{T}}_m}\mu \times \mu (\left\{(x,y)\in T\times X:d(x,y)<s\right\}),$$

$$|F_m(s)-F(s)|$$

$$+\sum _{T\in {\mathcal{T}}_m}\left|\mu (\left\{x\in T:d(x_T,f^m(x))<s\right\})-\mu (B(x_T,s))\mu (T)\right|$$

$$+\sum _{T\in {\mathcal{T}}_m}\left|\mu (B(x_T,s))\mu (T)-\mu \times \mu (\left\{(x,y)\in T\times X:d(x,y)<s\right\})\right|.$$

$$\left\{x\in T:d(x,f^m(x))<s\right\}\Delta \left\{x\in T:d(x_T,f^m(x))<s\right\}$$

$$=f^{-m}(\mathcal{A}(x_T,s,{\epsilon }_m))\cap T.$$

$$|\mu ((\mathcal{A}(x_T,s,{\epsilon }_m))\mu (T)-\mu (T\cap f^{-m}(\mathcal{A}(x_T,s,{\epsilon }_m)))|\le {\epsilon }_m\mu (T).$$

$$\mu (T\cap f^{-m}(\mathcal{A}(x_T,s,{\epsilon }_m)))\le {\epsilon }_{m}K\mu (T)+{\epsilon }_m\mu (T).$$

$$\left\{x\in T:d(x_T,f^m(x))<s\right\}=f^{-m}(B(x_T,s))\cap T.$$

$$\left|\mu \times \mu (T\times B(x_T,s))-\mu \times \mu (\left\{(x,y)\in T\times X:d(x,y)<s\right\})\right|\le K\epsilon \mu (T)$$

$$\left\{T\times B(x_T,s)\right\}\Delta \left\{(x,y)\in T\times X:d(x,y)<s\right\}\subset T\times \mathcal{A}(x_T,s,{\epsilon }_M)$$

$$|F_m(s)-F(s)|\le 2(K+1){\epsilon }_m.$$

$$F_m^N(s)=\frac{1}{N-m}\sum _{i=0}^{N-m}{\chi }_s(d(x_i,x_{i+m})).$$

$$F_m(s):={\int }_{x\in X}{\chi }_s(d(x,f^m(x)))d\mu .$$

$$\underset{N\to \infty }{\lim }\frac{1}{N-m}\sum _{i=0}^{N-m}{\chi }_s(d(x_i,x_{i+m}))={\int }_{x\in X}{\chi }_s(d(x,f^m(x)))d\mu ,$$

$$f(x):=2\min \left\{x,1-x\right\}:[0,1]\to [0,1].$$

$$||F_m-F|{|}_{\infty }\le C{2}^{-m}$$

$$F(s)=\left\{\begin{array}[]{ll}0\text{ if\;}s<0,\\ 2s-s^{2}\text{ if\;}0\leq s\leq 1,\\ 1\text{ if\;}s>1.\end{array}\right.\\$$

Full text

I Introduction.

Understanding the global dynamics of nonlinear systems is typically challenging, especially if the governing equations are not known and only experimental data are available. A significant challenge arises from chaotic behavior as this makes precise long term prediction of future states impossible. While low dimensional chaotic dynamics is relatively well understood via the geometric theory of differential equations Guckenheimer and Holmes (2013); Robinson (1998), higher dimensional systems are usually treated using the ergodic theory Walters (2000); Petersen (1989).

In this Letter we propose a quantitative measure to assess the rate at which an ergodic system looses its memory, i.e. how fast the future states become independent of the initial condition. The main idea is to analyze distributions of distances between consecutive states of the system sampled at different sampling rates. We stress that knowledge of the governing equations of the system is not required and the method can be applied to experimental data. However, to obtain a rigorous upper bound of the mixing rate for the underlying dynamical system $f\colon X\to X$ we have to require that: (i) $X\subset R^{M}$ is compact, (ii) the system has a unique physical invariant measure $\mu$ and is strongly mixing with respect to this measure, i.e. $\mu(f^{-m}(A)\cap B)$ converges to $\mu(A)\mu(B)$ as $m\to\infty$ for any two measurable sets $A,B\subset X$. There exist a few slightly different definitions of the mixing rate. We consider it to be the rate at which the strong mixing coefficient $\alpha$ converges to zero  Bradley et al. (2005). Our result provides an upper bound on convergence rate of $\alpha$ and we employ it to quantify the memory loss of the system.

First we familiarize the reader with our method by considering discrete dynamical systems generated by the tent map and the logistic map. However, our main motivation comes from the dynamics of a large number of interacting particles modeling granular systems. Granular systems have attracted a lot of attention over the last centuries because of their importance to everyday life. However, many of their properties remain obscure even today. In particular, any description of a granular system in terms of particles’ positions and moments is necessarily incomplete, since it sheds limited light on the properties of the particle interactions; even for the simplest case of a static system of (frictional) granular particles, the interaction field is not uniquely determined by the particle positions. Both physical experiments and simulations have shown that the interaction field consists of complex interaction networks that are known to be crucial for understanding mechanical properties of the system Papadopoulos et al. (2018). Properties of these networks can be well described using tools of algebraic topology Dijksman et al. (2018). In particular, persistence diagrams Carlsson (2009); Edelsbrunner and Harer (2010) provide quantitative succinct descriptions of the changes in topology of an interaction network as the force level changes from infinity to zero. For a detailed treatment of the subject we refer the reader to Kramár et al. (2014a) while a more compact presentation can be found in Kramár et al. (2013). The space of persistence diagrams is a metric space and one can study the dynamics of interaction networks Kramár et al. (2014b); Kondic et al. (2017) in this space.

This Letter presents an application of our measure of the rate of memory loss to a simple shear flow of granular particles in the regime where inertial effects are strong, and particles are stiff. We show that dependence of the rate of memory loss on the parameters of the system is nonlinear and not at all obvious. Our study shows that the global dynamics of the system changes erratically with the control parameter suggesting a presence of bifurcations even if the inertial number, $I$,  da Cruz et al. (2005) is kept fixed. We also show that the rate of memory loss is closely connected to the frequency of abrupt changes in the persistence diagrams of the interaction networks, triggered by sudden reorganization of the force network. Moreover, we show that these changes are correlated to system-wide average measures such as shear stress. We document this by demonstrating a relation between decay of autocorrelation of shear stress and our measure of the memory loss.

II Mixing rate.

We propose that the mixing rate can be inferred from the convergence rate of the cumulative distribution functions (CDFs), $F_{\tau}$, of the distances $d(x(t),x(t+\tau))$ along a trajectory $x$. To provide motivation and intuition for our result we consider the discrete system generated by the tent map $f(x):=2\min\left\{{x,1-x}\right\}$ on the interval $[0,1]$. The physical invariant measure $\mu$ of this system coincides with the Lebesgue measure on $[0,1]$. We begin by computing a finite sample $\left\{{x_{i}}\right\}_{i=0}^{N}$, $N=5\times 10^{6}$, of a trajectory starting from $x_{0}=0.1$, see the Section VII.4 for more details. From this sample, we calculate the CDFs, $F_{m}^{N}$, of the distances $\left\{{d(x_{n},x_{n+m})}\right\}_{0}^{N-m}$ for different values of $m$. Figure 1(a) suggests that $F_{m}^{N}$ converge as $m\to\infty$. In Section VII.3 we prove that for the CDFs $F_{m}$ of the random variable $d(x,f^{m}(x))$, with $x$ distributed according to $\mu$, there exits $C>0$ such that

[TABLE]

where $F$ is the CDF of the distance $d(x,y)$ between two random variables $x$ and $y$ which are i.i.d. according to $\mu$. Figure 1(b) shows that $||F^{N}_{m}-F||_{\infty}$ exhibits the decay predicted by (1) only for $m\leq 5$. For $m>5$ a larger sample is needed to properly approximate ’spatial averages’ $F_{m}$ by ’time averages’ $F_{m}^{N}$.

It might seem that the decay given by (1) is closely connected to the Lyapunov exponent of the system. The tent map locally stretches the space by factor of two and its Lyapunov exponent is $\ln 2$. However, the following theorem shows that the convergence rate reflects the mixing rate, which has a more global nature than the Lyapunov exponent.

Theorem 1**.**

Let $X\subset{\mathbb{R}}^{M}$ be compact and suppose that $f\colon X\to X$ has a unique invariant measure, $\mu$, whose Radon-Nikodym derivative is continuous with respect to the Lebesgue measure on ${\mathbb{R}}^{M}$. Let $\left\{{\varepsilon_{m}}\right\}_{m=1}^{\infty}$ be a sequence of positive numbers converging to zero. If there exits a sequence of partitions $\left\{{\mathcal{T}_{m}}\right\}_{m=1}^{\infty}$ of $X$ such that diameter of every set $T\in\mathcal{T}_{m}$ is less than $\varepsilon_{m}$ and for every measurable set $E\subset X$

[TABLE]

then there exits a constant $C>0$ such that

[TABLE]

The relevant consequence of this theorem is that the convergence rate of $F_{m}$ provides an upper bound on how fast the sets with diameter $\varepsilon_{m}$ are mixed. As the system is mixed the trajectories starting from similar initial conditions become independent. Thus our ability to forecast the future states decreases, and we can interpret the mixing rate as a rate of memory loss of the system. For further intuition the reader is referred to Section VII.5 where a study of the logistic equation is presented.

III Continuous dynamics.

Our model for the granular dynamics takes the form of a differential equation $\dot{x}=h(x)$ where $x\in{\mathbb{R}}^{M}$. For a short time scale $\tau$ the dynamics tends to be well approximated by the linearization of the system and $d(x(t),x(t+\tau))\approx\tau||h(x(t))||_{2}$ where $||\cdot||_{2}$ is Euclidian norm in ${\mathbb{R}}^{M}$. To detect the longest time scale for which the system is well approximated by its linearization, let us consider the CDFs $F_{m}$ of the distance $d(x(t),x(t+m\tau))$ along an orbit $x$. If the dynamics of $\dot{x}=h(x)$ is well approximated by its linearization at the time scale $m\tau$, then $d(x(t),x(t+m\tau))\approx md(x(t),x(t+\tau))$ and $F_{m}(s)\approx F_{1}(ms)$. This scaling is lost at the time scale at which the non-linear effects become important. Moreover, if the functions $F_{m}$ are essentially constant for $m>m_{0}$, then it is likely that the system has a bounded attractor and its subsets with small volume are well mixed at the time scale $m_{0}\tau$.

IV Linearly sheared systems.

Now we consider simulations of sheared granular particles. The details of the simulation protocol are given in Section VII.1; here we provide an overview. In order to avoid as much as possible the complications involved in inhomogeneous flows, but also governed by the goal to consider experimentally realizable configurations, we consider the following setup. Two-dimensional frictional bidisperse circular particles (elastic disks) are placed between solid walls that impose shear flow (shear rate $\dot{\gamma}$) by moving to the right (top) and left (bottom) with the same speed, $v$, (all the relevant quantities in what follows are non-dimensionalized using average particle diameter, mass, and binary collision time, $\tau_{c}$, as the length, mass, and time scale, respectively, see Section VII.1. To ensure uniform flow, the particles are placed on a solid substrate that moves with a linear velocity profile, similarly as in recent experiments Bi et al. (2011). The walls are subject to the applied pressure, $P$, chosen in such a way that the interaction forces between the particles are orders of magnitude stronger than the particle/substrate forces. In the present work, we vary both the shear rate and applied pressure while keeping the inertial number $I=\dot{\gamma}\sqrt{m/P}$ constant. The question is whether we can understand the memory loss and interaction network evolution in this simple setup.

There are three relevant time scales in the setup considered: $t_{s}={1/{\dot{\gamma}}}$ (shear time), $t_{I}=\sqrt{m/P}$ (inertial time) and $t_{c}=\sqrt{m/k_{n}}$ (binary collision time, also comparable to the time needed for a signal to propagate through a grain). One could also think of the typical contact time between the particles as a relevant time scale, however, for a simple shear flow considered here, this time is comparable to $t_{I}$. Out of these three time scales, one can produce two independent parameters; one possibility is to define $I={t_{I}/t_{s}}$ and $\kappa={t_{I}/t_{c}}$. Since $t_{c}\ll t_{I}\ll t_{s}$, $I\ll 1$ and $\kappa\gg 1$ (for the reference case, $I\approx 10^{-4}$ and $\kappa\approx 10^{5}$, see Section VII.1. For such large values of $\kappa$, it has been argued that the main features of the flow are $\kappa$-independent, and that the effect of finite elastic modulus of the particles can be ignored da Cruz et al. (2005). Note that two simulations with different values of $P$ and $\dot{\gamma},$ but with $I$ and $\kappa$ the same are identical after rescaling the time with $t_{s}$.

The reference case of our simulations corresponds to $P=1$ and $v=v_{ref}\approx 2.5\times 10^{-5}$. We refer to the system sheared with $v=uv_{ref}$ by $S_{u}$. For different values of $u$, we record the positions of the particles and forces acting between them at times $t_{i}=2i\tau_{c}$, where $i=0,\ldots,N=5\times 10^{6}$. By applying persistent homology to the force networks recorded at times $t_{i}$, we obtain a sequence $\left\{{x_{i}}\right\}_{i=0}^{N}$ that captures the evolution of the topological properties of the interaction network along the sampled trajectory. To assess the differences between the persistence diagrams, we use the Wasserstein $d_{W^{2}}$ distance which mitigates the influence of noise Kramár et al. (2014a).

V Memory loss of linearly sheared systems.

We use Theorem 2 to infer the rate of memory loss for the systems $S_{u}$. Using the sample $\left\{{x^{u}_{i}}\right\}_{i=0}^{N}$ of the system $S_{u}$, we compute the CDFs $F^{N}_{m}[u]$ of $\left\{{d_{W^{2}}(x^{u}_{n},x^{u}_{n+m})}\right\}_{i=0}^{N-m}$. Figure 2 suggests that $F^{N}_{m}[1]$ and $F^{N}_{m}[5]$ converge. For every $S_{u}$ we verified that the functions $F^{N}_{m}[u]$ do not depend on the choice of initial conditions and change very little if $m\geq 10^{4}$. Hence, we approximate the limiting distributions $F[u]$ by $F^{N}_{2^{15}}[u]$. As documented by Fig. 3(a), the limiting distributions $F[u]$ vary with $u$ and there is no obvious trend. Thus, the asymptotic dynamics of the systems $S_{u}$ changes in a rather complicated manner.

To study the memory loss we denote a CDF of $\left\{{d_{W^{2}}(x^{u}(t),x^{u}(t+{\Delta t}))}\right\}$ by $F^{N}_{\Delta t}[u]$. Figure 3(b) shows $||F^{N}_{\Delta t}[u]-F[u]||_{\infty}$ as a function of $\Delta t$ for different systems $S_{u}$. It indicates that the rate of memory loss increases with $u$, but in a non-smooth fashion. This suggests a presence of bifurcations along the parameter $u$. The rate of memory loss increases with $u$ as we pass over the consecutive bifurcation values while it stays essentially constant away from this values.

Based on the earlier discussion about expected minor influence of the parameter $\kappa$ measuring elasticity of particles, it is natural to analyze memory loss by considering the time scaled by $t_{s}$. To be specific, we rescale time by the wall speed in such a way that the wall moves by one particle diameter during one time unit. Using rescaled time $t_{u}=t/(uv_{ref})$, we consider $||F^{N}_{\Delta t_{u}}[u]-F[u]||_{\infty}$ as functions of $\Delta t_{u}$. They indicate that the decay rate of $||F^{N}_{\Delta t_{u}}[u]-F[u]||_{\infty}$ varies erratically with $u$.

This finding is surprising since, as mentioned earlier in the text, one would expect essentially identical behavior of the considered systems since $I$ is kept fixed, and the values of $\kappa$ are large. To better understand this behavior, we recall that at the physical time scale the rate of memory loss does not change much with $u$ unless the parameter $u$ passes through an expected bifurcation value. At those values the rate increases with $u$, see Fig. 3(b). On the other hand, Fig. 4(a), shows that at the time scale $t_{u}$ the rate of memory loss decreases with $u$ between the consecutive bifurcations which is consistent with da Cruz et al. (2005). However, the rate of memory loss still increases with $u$ as it crosses over a bifurcation value.

Finally, we connect the unexpected behavior shown in Fig. 4(a) directly to the evolution of interaction network. The animations, indicate that this evolution exhibits slow-fast dynamics: slow dynamics is dominated by build up of force ‘chains’ that buckle and lead to large and fast rearrangements. We proceed to analyze these transitions and their possible connection to the memory loss.

To quantify the transitions, we consider the evolution of the systems at the fastest timescale given by our sampling. Let $D_{u}(t_{u})$ be the $d_{W^{2}}$ distance between the state of the system $S_{u}$ sampled at time $t_{u}$ and the next sampled state. The peaks of $D_{u}$, visible in Fig. 5, indicate presence of abrupt transitions of the system which otherwise exhibits a relatively slow evolution. We define individual abrupt transitions as distinct peaks of $D_{u}$ that exceed a specified threshold; the number of detected transitions depends on the choice of threshold, however we find that the following results suggesting universal dependence of the memory loss on the number of abrupt transitions hold for thresholds in $[0.2,0.4]$. Note that if $\bar{D_{u}}$ is the mean of $D_{u}$ and $\sigma_{D_{u}}$ is its standard deviation, then the interval $[\bar{D_{u}}+\sigma_{D_{u}},\bar{D_{u}}+3\sigma_{D_{u}}]\subset[0.2,0.4]$ for each $D_{u}$. We use the threshold value of $0.3$.

Let $\delta t_{u}$ be an average time between two consecutive abrupt transitions of the system $S_{u}$ and consider the rescaled time $t^{*}=t_{u}/\delta t_{u}$. Note that $t^{*}$ does not depend explicitly on $u$, but only on the number of transitions. Figure 4(b) shows that values of $||F^{N}_{\Delta t^{*}}[u]-F[u]||_{\infty}$ are very similar for all systems $S_{u}$. This suggests that the rate of memory loss depends universally on the frequency of the transitions. Thus, we conjecture that the memory is predominately lost during the abrupt transitions caused by localized buckling of the interaction networks.

Now, we investigate influence of the abrupt transitions on global measures such as the stress tensor. Figure 6(a) shows the temporal autocorrelation of the shear stress. We see that the ordering of the autocorrelations for various values of $u$ is similar to the ordering of the curves $||F^{N}_{\Delta t_{u}}-F||_{\infty}$, see Fig. 4(a). Moreover, only a handful of transitions occur before the autocorrelation drops below $0.2$. Figure 6(b) shows the average number $E_{u}$ of transitions of the system $S_{u}$ during the time that the wall moves by two particle diameters. Figures 6(a)-(b) indicate that the systems with large $E_{u}$ decorrelate faster.

VI Conclusion.

In summary, we formulated a method for estimating the rate of memory loss in complex spatio-temporal systems. We employed it to study the time evolution of force networks of granular systems in the space of persistence diagrams. We showed that dependence of the dynamics on the control parameters is very complex even for the simple case of uniform flow in a planar geometry. The rate of the memory loss depends erratically on the stiffness of the particles even if the particles are stiff and the inertial number is fixed. While proper understanding of this dependence should be the subject of future research, we found that the rate of the memory loss is strongly correlated with the frequency of abrupt transitions of interaction network. The fact that even a simple planar flow of circular particles exhibits an extremely complex dependence on the control parameters suggests that significant new research is still needed in the field of dense granular matter. In particular, much more work will be needed to fully understand the systems exposed to more complex flows or built from more complex particles.

**Acknowledgements ** M.K. was supported by ERC Gudhi (ERC-2013-ADG-339025). L.K. and L.K. were partially supported by NSF Grant No. DMS-1521717, DARPA contract HR0011-15-2-0033, and ARO Grant No. W911NF1810184. K.M. was supported by NSF-DMS-1248071, DMS-1521771, DMS-1622401, DMS-1839294 and DARPA contracts HR0011-16-2-0033 and FA8750-17-C-0054.

VII Appendix

VII.1 Particle interactions and simulation protocol

The particles in the considered numerical system are modeled as 2D soft frictional inelastic disks that interact via normal and tangential forces, specified here in nondimensional form. We use the average particle diameter, $d_{ave}$, as the lengthscale, the binary particle collision time $\tau_{c}=2\pi\sqrt{d_{ave}/(2gk_{n})}$ as the time scale, and the average particle mass, $m$, as the mass scale. Force constant, $k_{n}$ (in units of ${mg/d_{ave}}$), is set to a value corresponding to photoelastic disks Geng et al. (2003). The parameters entering the linear force model can be connected to physical properties (Young modulus, Poisson ratio) as described e.g. in Kondic (1999).

Dimensionless normal force between $i$–th and $j$–th particle is

[TABLE]

where $\boldsymbol{v}_{i,j}^{n}$ is the relative normal velocity, $\overline{m}$ is reduced mass, $x_{i,j}=d_{ave}-r_{i,j}$ is the amount of compression, with $d_{ave}={(d_{i}+d_{j})/2}$ and $d_{i}$, $d_{j}$ diameters of the particles $i$ and $j$. The distance of the centers of $i$–th and $j$–th particle is denoted as $r_{i,j}$. Parameter $\eta_{n}$ is the damping coefficient in the normal direction, related to the coefficient of restitution set to $e=0.5$.

We implement the Cundall–Strack model for static friction Cundall and Strack (1979). The tangential spring $\boldsymbol{\xi}$ is introduced between particles for each new contact that forms at time $T=T_{0}$ and is used to determine the tangential force during the contact of particles. Due to the relative motion of particles, the spring length $\xi$ evolves as $\xi=\int_{T_{0}}^{T}\boldsymbol{v}_{i,j}^{t}(t)dt$ with $\boldsymbol{v}_{i,j}^{t}=\boldsymbol{v}_{i,j}-\boldsymbol{v}_{i,j}^{n}$ and $\boldsymbol{v}_{i,j}$ being the relative velocity of particles $i,~{}j$. The tangential direction is defined as $\boldsymbol{t}={\boldsymbol{v}_{i,j}^{t}/|\boldsymbol{v}_{i,j}^{t}|}$. The direction of $\boldsymbol{\xi}$ evolves over time and we thus correct the tangential spring as $\boldsymbol{\xi}=\boldsymbol{\xi}-\boldsymbol{n}(\boldsymbol{n.\xi})$. The tangential force is set to

[TABLE]

with

[TABLE]

Viscous damping in the tangential direction is included in the model via the damping coefficient $\eta_{t}=\eta_{n}$. The value of the normal spring constant is $k_{n}=4\times 10^{3}$ and parameters $\eta_{n}$ and $k_{t}$ are set to $\eta_{n}=1.4$ (consistent with the specified value of $e=0.5$) and $k_{t}=0.8k_{n}$. Friction coefficient is set to $\mu=0.7$.

In the simulations the particles are placed on a base that moves with the speed that varies linearly from [math] to the speed of the top wall. The purpose of the base is to ensure linear velocity profile across the domain. There is a dissipative effect from friction between the particles and the base; for the $i$-th particle

[TABLE]

where $\mu_{b}$ is the friction between particle and base and $\boldsymbol{a}_{i}$ is the acceleration of the $i$-th particle. The magnitude and direction of $\boldsymbol{a}_{i}$ is determined here from the interaction of the particle $i$ with all particles in contact

[TABLE]

where $c_{i}$ runs over all particles in contact with particle $i$.

In our simulations we integrate Newton’s equations of motion for both the translational and rotational degrees of freedom using a $4$th order predictor-corrector method with time step $\Delta t=0.02$. Initially, particles are placed on a grid and given random initial velocity; particles are bidisperse with the ratio of the large to small particle diameter $1.4$. Approximately $1/3$ of the particles have large and $2/3$ of the particles have small diameter. There is $\approx 1200$ particles and the rectangular domain is $54$ particle diameters wide and $23$ particle diameters high (in terms of $d_{ave}$).

Our simulations start by slowly compressing the domain with a specified pressure, $P=\rm k^{2}P_{0}$, applied on the top wall, until the top wall reaches a steady position. The system is then sheared with velocity $v=kv_{0}$, where $v_{0}=2.5\cdot 10^{-5}$ (in the units of $d_{ave}/\tau_{c}$). Walls are built of monodisperse particles with diameters of size $d_{ave}$ placed initially at equal distances, $d_{ave}$, from each other; in the horizontal direction the boundary conditions are periodic. The value of $P_{0}$ is found by compression of the top wall up to a packing fraction $\rho=0.80$ ensuring a dense packing above jamming point. We show results for different values of $k\in[1.0,3.0,4.0,\sqrt{20},5]$; note that the ratio of the applied pressure and shearing velocity guarantees fixed inertial number (as discussed in the main body of the paper). Specifically, the shear rate, $\dot{\gamma}=H/v$, where $H$ denotes the height of the granular system and inertial number $I=\dot{\gamma}\sqrt{m/P}\approx 10^{-4}$.

VII.2 Forces and Effective Friction

Interaction networks of the granular systems with the fixed inertial number, $I$, are expected to evolve in the same manner. In particular, the velocity profile during shear should stay the same regardless of the pressure or shearing velocity applied and the measures such as probability distribution function of forces, $\rm PDF$, and effective friction, $\mu^{*}$ found as a ratio of shear and normal stress, attain the same value da Cruz et al. (2005).

We show that, indeed, these global measures confirm the previous findings and that our conclusions do not depend on the subtle differences such as the existence of the frictional base. Specifically, the ratio of the typical interparticle force, $\langle F\rangle$ and the force between the particle and frictional base, $F_{b}$, is $\langle F\rangle/F_{b}<10^{-2}$. We also find that the ratio of the kinetic and elastic energy in each system considered is $\ll 1$, showing that the force dynamics is not influenced by the existence of the frictional base in the dense flow regime, characterized by $I\leq O(10^{-3})$ da Cruz et al. (2005). Figure 7 shows the probability distribution functions of forces, $\rm PDF$, averaged over complete simulation. As expected in the case of a fixed inertial number, the $\rm PDF$ curves collapse for all the systems considered.

Next we discuss the effective friction, $\mu^{*}$. It has been shown da Cruz et al. (2005) that the value of $\mu^{*}<\mu$ for dense flows and that it remains constant for a large range of inertial numbers, $I\leq(10^{-3})$. Since in our case $I\approx 10^{-4}$ we expect to see similar results (i.e. $\mu^{*}<\mu$). Figure 8 shows $\mu^{*}$ for different systems as a function of the wall displacement, $d$ (in the units of $d_{ave}$). We observe $\mu^{*}\approx 0.45$ which is smaller than the interparticle friction $\mu=0.7$. Moreover, the value of $\mu^{*}$ fluctuates around the same value regardless of the pressure applied.

We point out the fact that all the considered measures yield the same results for all the systems that we consider. In other words, the conclusions found in da Cruz et al. (2005) are confirmed here as well. However, our findings strongly suggest that not all essential information can be gathered by means of global measures and one needs to consider the topology of interaction networks as well.

VII.3 Strongly mixing systems

In this section we consider a dynamical system generated by $f\colon X\to X$, where $X\subset{\mathbb{R}}^{M}$ is compact. We suppose that the system has a unique invariant measure $\mu$, and its dynamics is strongly mixing with respect to this measure. We also require that the measure $\mu$ has a continuous Radon-Nikodym derivative with respect to the Lebesgues measure $dx$, i.e. there exists a continuous function $g\colon X\to X$ such $\mu(A)=\int_{A}gdx$ for every measurable set $A\subseteq X$. Under these assumptions we get the following upper bound on the measure of the intersection of an $\varepsilon$ neighborhood of a sphere and subsets of $X$.

Lemma 1**.**

Let $X\subset{\mathbb{R}}^{M}$ be a compact set and $\bar{\varepsilon}>0$. For $x\in X$, $s\in{\mathbb{R}}$ and $\varepsilon>0$ we define

[TABLE]

where $d$ is the Euclidean distance. If $\mu$ is a measure on $X$ whose Radon-Nikodym derivative is continuous with respect to the Lebesgues measure on ${\mathbb{R}}^{M}$, then there exists a constant $K$ such that

[TABLE]

for every $\varepsilon<\bar{\varepsilon}$

Proof.

The compact set $X$ is contained in some ball with a sufficiently large radius $R$. Thus, volume of every set $\mathcal{A}(x,s,\varepsilon)$ is smaller than the volume of an $\varepsilon$ neighborhood of an $M$-dimensional sphere with radius $R$ which is given by

[TABLE]

where $K_{m}=\frac{\pi^{\frac{M}{2}}}{\Gamma(\frac{M}{2}+1)}$. This implies that, for $\varepsilon<\bar{\varepsilon}$, the Lebesgues measure of the set $\mathcal{A}(x,s,\varepsilon)$ is bounded by $K^{\prime}\varepsilon$ where $K^{\prime}=2K_{M}\left({N\choose 1}R^{N-1}+{N\choose 3}R^{N-3}\bar{\varepsilon}^{2}+\ldots+\bar{\varepsilon}^{N-1}\right)$. Finally, by continuity of the Radon-Nikodym derivative $g$ of the measure $\mu$ we get that $\mu(\mathcal{A}(x,s,\varepsilon))\leq K\varepsilon$ for $K=K^{\prime}\max_{x\in X}g(x)$. ∎

We use this lemma to prove convergence of CDFs, $F_{m}$, of the random variable $d(x,f^{m}(x))$ with $x$ distributed according to $\mu$. The main theorem of our Letter (restated below) guarantees that the mixing rate of the system provides an upper bound on $||F_{m}-F||_{\infty}$ where $F$ is a cumulative distribution of the random variable $d(x,y)$ with $x,y$ independently distributed according to $\mu$. Thus, the theorem implies that for strongly mixing systems $||F_{m}-F||_{\infty}\to 0$ as $m\to\infty$.

Theorem 2**.**

Let $X\subset{\mathbb{R}}^{M}$ be compact and suppose that $f\colon X\to X$ has a unique invariant measure, $\mu$, whose Radon-Nikodym derivative is continuous with respect to the Lebesgues measure on ${\mathbb{R}}^{M}$. Let $\left\{{\varepsilon_{m}}\right\}_{m=1}^{\infty}$ be a sequence of positive numbers converging to zero. If there exits a sequence of partitions $\left\{{\mathcal{T}_{m}}\right\}_{m=1}^{\infty}$ of $X$ such that diameter of every set $T\in\mathcal{T}_{m}$ is less than $\varepsilon_{m}$ and for every measurable set $E\subset X$

[TABLE]

then there exits a constant $C>0$ such that

[TABLE]

Proof.

It follows from additive property of the measure that the CDFs $F_{m}$ and $F$ can be calculated as:

[TABLE]

and

[TABLE]

where $\mu\times\mu$ is the product measure. For every set $T\in\mathcal{T}_{m}$ we choose a random point $x_{T}\in T$ and define $B(x_{T},s):=\left\{{x\in X\colon d(x_{T},x)<s}\right\}$. The following inequality follows from Equations (13) , (14) and the triangle inequality:

[TABLE]

To show that Inequality (12) holds we need to properly bound the individual terms in the above inequality.

We start by estimating the terms in the first sum. We denote the symmetric difference of the sets $A$ and $B$ by $A\Delta B$ and recall that $|\mu(A)-\mu(B)|\leq\mu(A\Delta B)$ for all measurable sets $A$ and $B$. Now, we estimate the symmetric difference of the sets

[TABLE]

The inclusion follows from the fact that $B(x_{T},s-\varepsilon_{m})\subset B(x,s)\subset B(x_{T},s+\varepsilon_{m})$ for all $x\in T$. It follows from (11) that

[TABLE]

By applying Lemma 1 we get

[TABLE]

Therefore, the first sum in (VII.3) is bounded by $\varepsilon_{m}(K+1)$.

Now, we turn our attention to the second sum . Note that

[TABLE]

It follows from (11) that $\left|\mu(f^{-m}(B(x_{T},s))\cap T)-\mu(B(x_{T},s))\mu(T)\right|<\varepsilon_{m}\mu(T)$ and the second sum is bounded by $\varepsilon_{m}$. Finally, we consider the last sum. By definition $\mu(B(x_{T},s))\mu(T)=\mu\times\mu(T\times B(x_{T},s))$ and to bound the last sum by $K\varepsilon_{m}$ we just need to show that

[TABLE]

for every $T\in\mathcal{T}_{m}.$ By similar reasoning as above

[TABLE]

and Inequality (16) follows from the fact that $\mu(T\times\mathcal{A}(x_{T},s,\varepsilon_{M}))<\mu(T)K\varepsilon_{m}.$

By combining the estimates for individual sums we get that

[TABLE]

for all $s\in{\mathbb{R}}$. Hence Inequality (12) holds for $C=2(K+1)$.

∎

In practice we can only approximate the CDFs, $F_{m}$, from a finite sample $\left\{{x_{i}}\right\}_{i=0}^{N}$ of a trajectory starting from some initial condition $x_{0}$. The approximation based on this sample is defined by

[TABLE]

where $\chi_{s}$ is the characteristic function of the set $\left\{{x\in{\mathbb{R}}\colon x<s}\right\}$. We close this section by proving, that if $f$ is a Borel measurable function, then $\lim_{M\to\infty}F_{m}^{N}(s)=F_{m}(s)$, for every $s\in{\mathbb{R}}$ and $m\in{\mathbb{N}}$. By definition

[TABLE]

Hence to prove that $F_{m}^{N}$ converges to $F_{m}$ point-wise we just need to show that

[TABLE]

for $s\in{\mathbb{R}}$. This follows from Birkhoff ergodic theorem for almost every $x_{0}\in X$ under the assumption that $\chi_{s}(d(x,f^{m}(x)))$ is a Borel measurable function. If $f$ is Borel measurable, then the map $x\to(x,f(x))$ is Borel measurable. Moreover, the distance function $d$ is Borel measurable with respect to the product algebra and so the functions $\chi_{s}$ are Borel measurable as well.

VII.4 Tent map

In this section we consider the dynamical system generated by the ten map

[TABLE]

The invariant measure of this system coincides with the Lebesgues measure and we can use Theorem 2 to estimate the convergence rate of $||F_{m}-F||_{\infty}$. Lets us consider a family of partitions $\left\{{\mathcal{T}_{m}}\right\}$ of $[0,1]$ where $\mathcal{T}_{m}$ consists of the intervals $T_{m}^{i}=[i2^{-m},(i+1)2^{-m}]$ for $i\in\left\{{0,1,\ldots 2^{m}-1}\right\}$. The radius of each interval $T_{m}^{i}$ is equal to $2^{-m}$. Moreover, every interval $T_{m}^{i}$ is uniformly starched to $[0,1]$ by $f^{m}$. This implies that for a fixed $m$ and every measurable $E\subset[0,1]$ the measure of $\mu(f^{-m}(E)\cap T_{m}^{i})$ is the same for all $i$. Because $f$ preserves $\mu$ it has to be equal to $\mu(E)2^{-m}$ which is exactly $\mu(E)\mu(T_{i})$ and so $|\mu(f^{-m}(E)\cap T_{m}^{i})-\mu(E)\mu(T_{i})|=0$. Now it follows from Theorem 2 that

[TABLE]

which is the inequality presented in our Letter.

The fact that the invariant measure $\mu$ coincides with Lebesgue measure on $[0,1]$ makes it easy to calculate the limiting distribution

[TABLE]

This allows us to study the convergences rate $||F_{m}^{N}-F||_{\infty}$ of the approximate distributions obtained from a finite sample $\left\{{x_{i}}\right\}_{i=0}^{N}$, $N=5\times 10^{6}$, of a trajectory starting from $x_{0}=0.1$. Due to the binary floating-point representation of the numbers in a computer, the direct iteration of the tent map does not produce a good sample of the trajectory. Instead we use the fact that the tent map is homeomorhpic to the logistic map with parameter $r=4$. If we denote the logistically evolving variable by $y_{n}$, then the iterations of the tent map are given by $x_{n}={\tfrac{2}{\pi}}\sin^{{-1}}(y_{{n}}^{{1/2}})$. We showed in this Letter that the bound on the convergence rate given by (17) does not hold for $||F_{m}^{N}-F||_{\infty}$ if $m>5$. This is due to the fact that the sample is not large enough and $F_{m}^{N}$ does not approximate $F_{m}$ with necessary accuracy required to observe the predicted decay rate.

VII.5 Logistic map

Our main theorem shows that convergence rate of $||F_{m}^{M}-F||_{\infty}$ can be used to infer the mixing rate of the system. The mixing rate is known to influence decay of correlations between the states of the system Wiggins and Ottino (2004). For the logistic map we demonstrate that convergence rate of $||F_{m}^{M}-F||_{\infty}$ is closely related to the decay of correlations and thus well suited for measuring the memory loss of the system.

The logistic map

[TABLE]

is a classical example showing that complex chaotic behaviour can be produced by very simple non-linear systems. Dependence of the dynamics on the parameter $r$ is very complicated and exhibits a large number of different bifurcations. This is documented by erratic behaviour of the Lyapunov exponents, shown in Fig. 9(a), which measure the complexity of the dynamics. However, this measure is rather local and is not directly connected to the mixing rate.

Figure 9(b) depicts the decay of correlations for different values of $r$. Generally, the decay rate increases with $r$ and does not reflect the erratic behavior of the Lyapunov exponents. The decay of $||F_{m}^{M}-F||_{\infty}$ is shown in Fig. 9(c). As expected the systems with faster decay of correlations are mixing faster and thus $||F_{m}^{M}-F||_{\infty}$ also decays faster for these systems. This suggest that our measure is indeed well suited for assessing the rate of memory loss.