{'fraction_non_alphanumeric': 0.028540149702123326, 'fraction_numerical': 0.0048882325984011404, 'mean_word_length': 4.512068481616615, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Shearing of loose granular materials: A statistical mesoscopic model

A two-dimensional lattice model for the formation and evolution of shear bands in granular media is proposed. Each lattice site is assigned a random variable which reflects the local density. At every time step, the strain is localized along a single shear-band which is a spanning path on the lattice chosen through an extremum condition. The dynamics consists of randomly changing the `density' of the sites only along the shear band, and then repeating the procedure of locating the extremal path and changing it. Starting from an initially uncorrelated density field, it is found that this dynamics leads to a slow compaction along with a non-trivial patterning of the system, with high density regions forming which shelter long-lived low-density valleys. Further, as a result of these large density fluctuations, the shear band which was initially equally likely to be found anywhere on the lattice, gets progressively trapped for longer and longer periods of time. This state is however meta-stable, and the system continues to evolve slowly in a manner reminiscent of glassy dynamics. Several quantities have been studied numerically which support this picture and elucidate the unusual system-size effects at play.


I. INTRODUCTION
Modeling the rheology of granular media using continuum solid mechanics, has reached a high degree of sophistication in terms of constitutive equations [1]. Whatever the complexity of load paths being studied, an accurate account of the experimental stress strain relationship can now be achieved provided enough parameters or internal variables are included in the constitutive laws. However, such approaches are descriptive and leave unanswered questions pertaining to the scale of grain sizes.
In parallel to such phenomenological descriptive theories, a lot of effort has been spent in recent years in developing powerful computer models able to simulate granular systems at the individual grain level. Molecular dynamics approaches [2], or other techniques such as "contact dynamics" [3,4], now offer the possibility of dealing with several thousand particles, and provide extremely realistic pictures of the detailed micro-mechanics.
Such numerical techniques can be used to accurately investigate displacement fields, resolved both spatially and temporally [5]. The latter reveal an intriguing feature: namely that even in the most simple tests, such as a simple steady shear imposed over large strains, the local displacement appears very unsteady, with short quiescent periods where the displacement field is spatially smooth, separated by sudden changes where the configuration of grains reaches a local instability and undergoes a rapid reorganization through significant displacements at the grain level. This temporal variability manifests itself in giant stress fluctuations observed experimentally when particles and walls are stiff, and when the high frequency part of the stress signal is not filtered out [6]. Numerical simulations indicate [7] that instantaneous strain fields consist essentially in localized strains occurring along one or a few shear bands. However, as the strain increases (over the moderate range accessible in the simulation), there seems to be little or no correlation between successive shear bands, so that the time average of the displacement field erases these discontinuities and produces smooth strain fields.
Such fluctuations are obviously ignored in continuum modeling. And indeed it may appear that the identification of these instabilities is relevant only for discussing fine details of microscopic and transient features. However, their relevance can be judged only at a mesoscopic level of modeling, since the microscopic numerical techniques are far from being able to reach the relevant time scales. In fact, in the following we will argue that these instabilities may have a significant impact, both on large scale heterogeneities of the medium itself, and on a systematic slow time evolution of the macroscopic friction angle. A short account of some of our results has appeared in a previous publication [8].
The paper is organized as follows: in Section II, we will recall some features observed experimentally or numerically that we consider essential, and, in Section III, we progressively introduce the rules of a model whose aim is to describe some statistical aspects of shearing of loose granular media over large strains. In Section IV, we present in details the different quantities studied numerically for this model. We conclude in section V with a summary of our results and a discussion on possible experimental checks. The shear band is parallel to the shear direction z due to the periodic boundary conditions in this direction. We sum up along this direction to get a two dimensional sample in the xy plane.

II. THE SHEAR PROCESS IN LOOSE GRANULAR MATERIAL
We will address here the question of the behavior of granular media subjected to a simple shear for large strains. We restrict ourself to the simplest granular medium one may consider, namely rigid (undeformable) grains with Coulomb friction. This refers experimentally to dry sand subjected to a low confining pressure. We are concerned here with large strains, and thus in order to avoid the problem of boundary conditions which would limit the maximum strain, we consider an annular shear cell. To simplify the problem further, we consider only the case where the problem is invariant along the shear direction. As shown in Figure 1, the displacement is a single function of the coordinate of a radial cross-section (x, y), and constant along the orthoradial direction z (traditionally this situation is termed "antiplane"). Moreover, we are interested only in the quasistatic regime, i.e., time as such is irrelevant, and only the total strain matters. Thus, in what follows, what is referred to as "time" is to be seen here as a practical means of parameterizing the total strain being imposed on the medium.
One of the important observations of soil mechanics concerning such media is the concept of a critical state [9,10]. Depending on the preparation of the sample, the behavior under shear may differ considerably. For loose sand, (low density), the deviatoric stress to be applied increases with the total shear strain and simultaneously, a densification is observed [11,12]. However, as the shear strain increases, the density and shear stress seem to reach a plateau independent of the initial density. This state is called the "critical state". On the contrary, if the initial density is large, a single shear band forms, while the rest of the medium remains frozen [13]. The formation of the shear band is preceded by a volume expansion of the medium [14], but after the band is formed, all further properties remain quasi constant. A detailed experimental investigation has revealed [15] that inside the shear band, the density tends to approach the critical state density. This concept of the critical state has received considerable experimental evidence over the years, and is implemented in a number of continuum constitutive laws. Experiments however mostly deal with a rather moderate total strain well below unity.
A simple picture which is consistent with the critical state concept is that both the friction and the dilation angle increase with the density, and that the critical state is the density for which the dilation angle is zero (no change in volume under shear). Retaining the density as the only internal variable is an approximation. Other characteristics of the texture of the medium such as the fabric tensor (which has information about the orientation of contact normals), certainly play a significant role. For the purpose of simplicity, we will in the following only retain one single scalar internal variable governing the friction angle. It could be either simply the density or a combination of density and texture. Nevertheless, in all cases we will refer to this internal variable as "density", irrespective of its precise meaning.
As mentioned above, numerical simulations seem to reveal [7] evidence for the existence of instantaneous shear bands even in loose granular media. On the other hand, in experiments the strain appears to be homogeneous and not localized. The resolution of this apparent paradox is that the shear bands change rapidly, and may visit the entire medium in the process. Thus during an increment of shear, which can be observed experimentally, only a time average over many such shear bands is seen. In our modeling, we introduce a basic time scale for each elementary procedure. This time step is then clearly much shorter than most experimentally accessible time scales. However, our model attempts to achieve a qualitative rather than a precise quantitative mapping. One of the main features of the model is to show that these two apparently unrelated facts: the existence of instantaneous shear bands at early times and its localization at late times, are actually related, with a slow transition between these two limiting cases. This slow dynamics is reminiscent of slow ageing properties encountered in glassy systems, and indeed, we will see that a breakdown of ergodicity does appear in this model.

A. Motivation and definition
At every instant the two-dimensional medium is characterized by a single, scalar internal variable, the density ̺(x, y). This represents an average of the density along the orthoradial direction z. From this density, we deduce a corresponding local friction coefficient µ(x, y). The latter is assumed to be a single monotonically increasing function of the density [16]. For simplicity, we may assume a linear relationship in the following although this is inessential.
The strain is imposed on the shear cell through pre-scribed displacements of the bottom and top planes. As particles are considered rigid, (no elastic deformation), the shear cell can only move if the shear force exceeds a threshold value proportional to the normal pressure. This limit stress is given by the "weakest internal surface". Indeed, in our anti-plane geometry, the shear strain will localize on the surface (i.e. path in the (x, y) plane) which will fail first. The latter is assumed to be given by the following algorithm. For each directed path P spanning the entire cross-section along the x-axis, we compute the maximum shear force it can support according to the local density. Assuming that the local slope of the path is always small, this maximum force F (P) is simply proportional to the sum of local friction coefficients, and thus making use of the assumed linear variation of the friction coefficient with the density, F (P) is proportional to the sum of local densities, where the sum runs over the sites along the path. Among all the possible paths, the weakest P * (for which S(P * ) = min) will fail first, and this fixes the value of the shear force F = F (P * ). In agreement with the previously mentioned observation, at every basic time step, the shear strain is realized along a single shear band. Away from this shear band, the strain rate is zero, and thus the density is kept constant in time. However, inside the shear band, there is a relative motion of grains, and thus the density is susceptible to evolve. The next step is now to determine how the density inside the shear band evolves with time. Even though it is observed that large strains are necessary in order for a system to reach its critical state, we argue that at the microscopic level, the evolution of the medium cannot depend on the total imposed strain. Thus the evolution rules for the density within the shear band should be designed in such a way that they do not depend on the past history, but only on the present state (density field). As the density ought to contain the basic information of the local characteristics, we propose that within the shear band, the local density ̺(x, y) is randomly modified. More precisely, in one elementary time step, corresponding to the "life-time" of the shear band in a very loose granular sample, we assume that the density along the shear band acquires random uncorrelated values picked from a statistical distribution p(̺). The uncorrelated character of the distribution is however justified only on a mesoscopic scale.
After the elementary strain event, we have a new density map ̺(x, y). We now simply reiterate this procedure as long as desired: Namely we identify the new path which minimizes S(P), and update the value of the densities along this path randomly. As the purpose of the present article is to illustrate some statistical aspects of this dynamics, we do not try to mimic any specific granular system by imposing a realistic density distribution or initial correlations in ̺. We will choose here a simple uniform distribution between 0 and 1 for p(̺). The mean value and variance of the distribution p can be chosen arbitrarily, since a translation and rescaling of ̺ does not affect the result.
A key assumption of our model which may appear as precluding the occurrence of a slow evolution toward a critical state is the selection of the density values within a shear band from uncorrelated, smooth distributions. In fact, we will show below that, on the contrary, a collective and purely statistical effect produces a slow increase of the mean density over large strains.
Our model is furthermore discretized on a regular square lattice. We have looked at two different kinds of square lattices to check the robustness of our results. In the first one the value of the density ̺ is carried by the bonds. The orientation of the lattice is chosen so that the principal directions lies at π/4 with respect to the (x, y) axis as shown in Fig. 2 a). In the other version, density values are assigned to the sites of a square lattice. In this case the minimal path can be connected through the next nearest neighbours too as shown on Fig.  2 b). Both square lattices give exactly the same results so in the following we just refer to them as square lattice realizations.
We mention here that a third type of lattice, the hierarchical diamond lattice [17] was also studied. The numerical results are surprisingly, essentially unchanged by the unusual topology of this recursively constructed lattice. The easier construction of this lattice however allows us to solve the model analytically thus giving us a quantitative picture of the behaviour of the system. These results are presented in [18] where the intimate relationship of our model to other models of statistical physics is also discussed.
The rules of our model, finding the extremal directed spanning path at every instant, is similar to finding the ground state of a directed polymer in a random potential [19]. However, in our case this potential is uncorrelated only at the beginning; it changes in time through the process described above, of ascribing new densities to all sites along the minimal path. It is clear from this re- lationship between the models, that the shear band is self-affine with a Hurst exponent ζ = 2/3 at the beginning, i.e. the transversal fluctuations of the band grow with the size (L x ) of the sample width as L ζ x . We will discuss the time evolution of the roughness later in this paper.

IV. NUMERICAL RESULTS
We first show the density map, Fig. 3, of the system at different times, t/L ranging from 4 to 4000. The grey scale chosen focuses on the vicinity of 1 so as to highlight the progressive densification. It may appear counterintuitive at first that the rest of the medium shows a densification at all, when the only dynamics consists of finding a minimal path and updating sites along it randomly. However the reason is simply that this update systematically hunts out the sites with the lowest density values and replaces them.
At early times, we observe an apparently uncorrelated field. However, as time proceeds, it is possible to distinguish preferential channels of high density aligned along the direction (x-axis) of the minimal paths. These channels however have a significant width which shows that though the minimal path has been confined to this zone, it has enough freedom to explore different neighboring configurations and achieve a significant local densification. We also see within these wide and dense channels, a single path with a lower density. This has been the last active minimal path in the channel. As time passes, the number of channels increases, and so does their width. They get partly interconnected, leaving always the same scars of low density paths. Finally, at the latest time shown on the figure, the average density is quite high, and traces of ancient minimal paths are still visible. Nevertheless what is striking is the occurrence of islands entrapped by these high density channels, where the density map looks like at the very early stage of the time evolution. This signals that these regions have basically not been visited by the minimal path during the entire history of the system. These features are quite generic, and they reveal that the spatial (and temporal) organization of the activity is rather complex. The rest of the study is devoted to a more quantitative account of this activity, of the resulting kinetics of compaction, and the unexpected finite size effects which appear in this problem.
In the following subsection, we will introduce the main measurements performed numerically on the model.

Average density
The most important quantity is the average density of the sample that we define as the mean of the density of the inactive sites, i.e. the sites not belonging to the shear band. We denote this by ̺ . This definition is convenient because ̺ increases monotonically by the rules of our model. In experiments one of the most frequently measured quantities is the volumetric strain which is just the change in the inverse of the average density. In our model, when p(̺) is chosen to be the uniform distribution in the interval [0, 1], the average density is bounded by ( ̺ ≤ 1). It is easy to see that in finite systems the steady state (the asymptotics) cannot be other than a system with maximal densities everywhere except for a path which will be always chosen as the minimal path. This state is equivalent to ̺ = 1.
We are interested in the approach to this asymptotics so we plot the quantity (1 − ̺ ) as a function of time. Within the granular medium context, this means that we mainly study a loose initial state and its convergence to the critical state [9,10]. However, we will also present results obtained when one starts from a high initial density later (Section IV G). Figure 4 shows that the difference of the average density from its asymptotic value first remains almost constant during a first stage t/L ≪ 1, and then it decreases steadily. This first increase of ̺ is well captured by a reduced time equal to t/L. However, as time progresses, the average density increases more and more slowly. Quite strikingly, the larger the system size, the slower the increase in density. Further down this will be interpreted as a breakdown of ergodicity.

Shear band density
It is also natural to define the density of the shear band that we denote by ̺ SB . This is just the average density of the sites along the minimal path (before updating them). As already mentioned in the introduction, we assume that the maximal static shear force is a single function of the density. Thus the density of the shear band can be related to the shear stress in experiments. Figure 5 shows the evolution in time of the difference (0.5 − ̺ SB ). As expected, as time proceeds, the density along the shear band will tend toward the average of the random densities which are used to refresh the sites or bonds along the shear band. Using a uniform distribution of densities between 0 and 1 implies that this average is 0.5. Note that in contrast to the previous case, the reduced time t/L accounts nicely for the time evolution of this quantity for all system sizes for t/L ≤ 10 4 . This is a second puzzle we will try to address further in the following sections as well as in [18].

Mean Hamming distance
We also calculate the Hamming distance of two successive shear bands which is defined as the number of sites (or bonds) by which two consecutive shear bands differ. We denote this distance by d. The natural normalization is to divide this distance by the total length of the path, L x . As we will see this quantity is very useful in characterizing the time evolution of the localization process.    6 shows that the mean Hamming distance is close to unity (i.e. two consecutive paths do not overlap at all) at early times, and decrease towards 0 for t/L ≫ 1. We recall that when the distance is equal to 0, then the two consecutive conformations of the shear band are identical, in spite of the total renewal of random densities along them. This indicates that the shear bands have a tendency to remain more and more persistent as the system "ages". We will analyze further the complete statistical distribution of the Hamming distance later in this paper.

Cumulative shear
An experimentally relevant quantity is the cumulative shearing denoted by σ cum . The numerical procedure we follow to obtain this quantity in our model is the following (see Fig. 7): We mark a line in the y-direction (see Fig.  2). We measure the total activity n a (y) along the line, i.e. the number of instances when the shear band went through a point y on the line. From this, we define σ cum (y) = y j=0 n a (j). ( By definition σ cum (L y ) = t, since at every instant, the shear band has necessarily to pass through one or the other site on a cut along the y axis. The fluctuations of σ cum (y) about its mean value (t/L y )y then reflect the inhomogeneity of the shear process within the sample. In Fig. 8, we track the time evolution of σ cum (y), after subtracting out the mean value. As indicated, a snapshot of this quantity encodes the history of the process of shearing in this system.

B. Early time regime
It is apparent from Figures 4, 5 and 6 that the initial behavior of the model is very different from the late stages. In the former regime, the average distance d (Fig.  6) is equal to the system size indicating that successive shear bands do not overlap at all. In other words there is an effective strong repulsive interaction between them.
The density on the shear band (Fig. 5) provides an explanation for this behavior. From the directed polymer [19] picture it is known that the first shear band has a mean density of ̺ SB (t = 0) = e * ≈ 0.22 on the tilted square lattice. The state of the shear band after the densities of the sites along it have been changed, is independent of its previous state and the mean density is now 0.5 (since it is just an average of L x independent random numbers taken from the uniform distribution [0,1]). This implies that at the next instant, the chosen shear band will have a very low probability of sharing bonds/sites with the previous one, since there will still be many spanning paths with density smaller than 0.5. Thus, until all the sites are visited at least once, successive shear bands have few sites/bonds in common (d = L x in Figs. 6). Also in this regime, since the path sweeps all sites through the cumulative shearing appears to be homogeneous (Fig. 8 first plot). As mentioned in the introduction, this uniform shear strain is consistent with the experimental observations that no well-defined shear-bands persist over observable time scales in loose samples.
The characteristic time to build up correlations is set by the sweeping through the sample, i.e., it is given by L y . This is the reason why in the early time regime the plots in Figs. 4-6 can be scaled together with L y .

C. Localization
The above described behavior is drastically changed as time goes on and the shear band gets localized for very long times at the same location.
As the average density increases, the probability of choosing a minimal path with density less than 0.5 starts decreasing, and thus for the new path, it becomes more favorable to overlap to a greater and greater extent with the previous one. As a result, the density of the region where the path is located, is further increased, hence trapping the minimal path in canyon-like structures surrounded by extremely high density regions (see Fig. 3 b-d).
The early repulsive interaction is thus now inverted to an attractive one as is shown by the rapid 1/t like decrease of the Hamming distance (Fig. 6), as well as the decrease of 0.5 − ̺ SB and of ̺ from its initial value (Figs. 5 and 4 respectively) in this regime.
There are a number of consequences of the localization. First, as time goes on, the shear band (which was not visible at all in the density map of the system (see Fig. 3 a)) becomes more and more apparent until finally it becomes localized at a given position for macroscopic times. This is in accordance with experiments and with the critical state concept [9,10]. Secondly, since the same path for the shear band is chosen most of the time, its density saturates to its asymptotic value 0.5 (Fig. 5). Again this is consistent with the experimental observation that the density within the shear band tends to achieve a well defined value somewhat smaller that the rest of the medium for dense granular media. Simultaneously, the shear stress saturates to a constant value since this is imposed by the shear band itself [1]. However, we note that in our model the global density of the system continues to increase in time, albeit very slowly. This extremely slow trend may well be out of reach experimentally. However, it is known [1] that the shear stress saturates much faster than the volumetric strain (the average density in our case) which is clearly justified by the numerical results.
On the cumulative shear which is a straight line with small statistical fluctuations in the early time regime, there appear step-like structures, indicators of progressively more persistent localization (Fig. 8). However, this localization is not everlasting since the shear band may perform big jumps to other local minima. This can be seen on the series of cumulative shear curves in the form of certain steps disappearing and others becoming more prominent. This prediction could easily be tested experimentally.

D. Systems with different aspect ratios
All changes take place along the shear band, which is aligned along the x direction, and thus we may anticipate that the x and y direction will play different roles. Therefore in this section, we study the influence of the width and length of the system. In what follows, we use the terms long for samples with L y < L x and wide in the opposite case (L y > L x ).
With very long and wide systems we are able to separate the two kinds of dynamics described in Section IV C. If one considers a long sample (L y /L x is small), we expect the time evolution to be independent of L x for all quantities of interest since the lattice can be split into subparts placed in series. Thus one may expect the large jumps to disappear and the average density to scale solely with L y (Fig. 9). A wide system, on the other hand, might be expected a priori to behave like a number of competing subsystems.

Long systems
On Fig. 9 we have plotted the time dependence of the density in long samples with L y = 4 (lower curves) and L y = 8 (upper curves). The t/L y scaling is excellent in both cases. However, the densification obeys a different time evolution for different L y . The rate at which the density increases is slower as the width increases.

Wide systems: Breakdown of ergodicity
Wide samples can be considered as subsystems placed next to each other and coupled in parallel. In contrast with the previous case, we will see that the evolution of the different subsystems cannot be accounted for by a simple average. In the case of the wide systems the same plot as before (Fig. 10 a) shows no t/L y scaling for small system sizes.
We could imagine the following construction: suppose we split a given wide system into two subsystems of size L x × (L y /2). Provided L y is large enough, we can ignore the interaction between the two subsystems, and thus we could study independently the time evolution of both subsystems. Now, if we merge them again, we realize that the only reason why the resulting densification could differ from the measurement on the separate subsystems is that the time t 1 the shear band has stayed in subpart 1 is far from being equal to t/2. In other words, the breakdown of the data collapse of 1 − ̺ vs. t/L y is a breakdown of ergodicity. This is naturally associated to what we termed "localization" earlier.
In order to make this concept more explicit, let us con-sider an extreme version of such a breakdown of ergodicity. During a first stage, up to t/L y of order 1, the activity is evenly spread over the system. Then we assume that after such a time, the activity remains confined in a subsystem of size L x × ℓ. All other subsystems are assumed not to be visited by the shear band, and thus their density is quenched at their value reached at the onset of localization, ̺ 0 , at time t 0 = θL y . The global density will thus obey where f ℓ (t) = 1 − ̺ (t) describes the densification of the representative cell of size L x × ℓ where the activity is confined, and thus 1 − ̺ 0 = f ℓ (θℓ).
In this crude scenario, we note that the global density does not converge to 1 as time goes to infinity, but rather remains stuck at a value such that (1 − ̺ (t)) → f ℓ (θℓ)(1 − ℓ/L y ). In more quantitative terms, we tried to carry out such a procedure, and indeed for a fixed L x , it is possible to account for the time evolution of systems of different width using ̺ 0 as free parameter and calculating ℓ from the L y dependence. It turns out that ℓ changes for small values of L x but becomes constant (ℓ ≃ 30) above the system size of L x ≃ 30. The test of this analysis can be seen on Fig. 11. However, the asymptotic density turns out to depend on L x .
The conclusion is that although, such an extreme modeling of the localization is able to capture part of the strong size effects observed numerically, it is too crude to provide a quantitative account of the densification. The hierarchical lattice provides us with a convenient case where an analytical investigation of this breakdown of ergodicity can be made. It is shown in Ref. [18], that the local "age" distribution assumes a multifractal distrbution whose spectrum can be computed exactly. This property can then be used to provide an expression of the density evolution in time.

E. Time evolution of the Hamming distance distribution
We study here the distribution of Hamming distances as a function of time, P (d, t) for both lattices. This quantity is the analogue of an "avalanche distribution", such as is usually studied in self-organized critical(SOC) systems. As we will see below, this quantity does indeed decay as a power-law like in SOC systems. However, since a steady state is never reached, the power-law decay is multiplied by a time-dependent prefactor.
The distribution P (d, t) is shown in Fig. 12. At early times, this quantity is peaked around the maximum value (d = L x ) while in the localized regime, it becomes peaked at the minimum value (d = 0). This corresponds to the transition from repulsive to attractive effective interaction between consecutive path conformations.  Fig. 10 with a different, effective density ̺ * = ̺0(Ly − ℓ). Note that the plot displays nice scaling as well as a clean power law decay over at least three decades. However if the simulation is continued further the average density increases over ̺ * .
At fixed (large) times, the distribution P (d, t) decays as a power-law of d, as can be seen in Figure 13. The measured exponent is in addition to which there exists a peak at d = 0 the amplitude of which varies significantly with time. The decay of average d as 1/t (see Fig. 6) found earlier, implies that the time dependence of the d = 0 part is p(d, t) ∝ 1/(td 3 ). Thus, including the different scalings with L x and L y , we obtain finally the asymptotic form Here again, the hierarchical lattice allows us to compute this distribution analytically (for L x = L y ), and we find that a similar behavior is obtained [18].

F. Roughness exponent of the shear band
From the directed polymer analogy we know that the shape of the shear band is self-affine with an exponent of ζ = 2/3 for infinitely large systems. It is an interesting question whether this property of self-affinity is conserved in the time evolution of our system. We have investigated this question and have found self-affine scaling for all times, albeit with a time dependent Hurst exponent (Fig.  14).
We have estimated the value of ζ by measuring the width of the shear band, w L (t), for different system sizes, L × L. The width is defined as the standard deviation of the y coordinate of the active path. The latter is expected to scale as w L (t) ∝ L ζ for a self-affine object. To estimate   ζ, the roughness exponent, we have computed the ratio of two such widths for lattice sizes differing by a factor of 2, and used the following estimate where w L (t) is the width of a shear band in a L×L lattice at time t. The results obtained for the tilted square lattice can be seen on Fig. 14. It starts from ζ(t = 0) = 2/3 as expected from the directed polymer result and has an asymptotic value of ζ ≃ 0.8.

G. Systems with high initial densities; relevance of initial conditions
We have so far studied the situation when a very loose granular medium compactifies under shear, while simultaneously a shear band gets quasi-localized in the system. It is also of interest to study samples with higher initial density where it is known experimentally that a shear band is localized from the very beginning. In this section we study the interesting crossover to that state from the previously described dynamics.
We choose initial densities from the interval [̺ init : 1] with a uniform distribution and varying ̺ init . In Fig. 15 the time evolution of the difference of the average density from its asymptotic value is plotted using different initial conditions. Since all previous arguments hold we assume that in finite systems the asymptotic value of the average density is 1 irrespective of the initial density.
The striking result of these simulations is that if ̺ init < 0.5, all curves coincide in the decreasing regime. However, if ̺ init > 0.5 a different time evolution is observed for large times. Thus we can assume that in the early time regime, when the shear band is swapping uncorrelatedly, it visits all sites that have a value less than the expectation value of the refreshing density distribution = 0.5. After the first regime, as all small values are eliminated, the system effaces the initial condition almost entirely.
On the other hand, starting from initial conditions with ̺ init > 0.5 we largely eliminate the possibility of big jumps. The shear band fluctuations are now very small, involving changes in a very few sites, and hence, the density change of the sample is extremely slow.

H. Summary of the numerical results
We have seen that the system densifies with time so as to approach a unit density, i.e. the maximum available density from the distribution used to refresh the sites. The kinetics of the densification is slow (slower than any power-law). Moreover, after a first transient where the reduced time t/L accounts for the L dependence, the compaction process depends on the system size in a nontrivial way. The width of the system is the parameter which really controls this anomalous behavior, signaling that the competition between parallel paths may somehow play a key role in this breakdown of system size rescaling. This competition is a subtle one however, lying somewhere in between complete localization of the path (which, as we saw in section IV B, is too crude to mirror the actual scenario) and complete delocalization (which, as mentioned earlier, accounts only for the early time behaviour).
The density maps display an interesting organization of "canyon-like" paths with density much lower than their immediate surroundings (where the density approaches 1 quite uniformly). Moreover, large regions are left quiescent, being systematically avoided by the minimal paths. This contrast of high and low activity within the same system is at the heart of the breakdown of ergodicity observed after an initial transient. As remarked however, the distribution of Hamming distances between consecutive paths displays a somewhat simpler behavior where the role of the width and the length of the system can be simply accounted for.

V. CONCLUSION AND DISCUSSION
Though very simply defined, our model seems to capture some essential features of granular shear and provides at the same time several predictions. The model demonstrates the self-organized mechanism of the localization of the shear band in loose granular materials. As the sample ages, very high fluctuations in density appear where we can observe some kind of screening effect: more resistant regions of higher density, protect the looser ones. This model also gives an insight into a dynamics that exhibits very non-trivial system size effects.
Our stochastic model makes an attempt to describe the large strain behavior of sheared loose granular matter on a mesoscopic level. The rules of the model do not include any dependence on the total amount of shear imposed on the medium, nevertheless, a constant friction angle and slow densification is observed -a property referred to as "ageing" -which reproduces the experimental results qualitatively [11]. By construction, the strain takes place through local shear bands which ini-tially travel throughout the medium homogeneously (and hence produce a uniform shear), but which progressively become more permanent giving rise to more steady shear bands, a feature also observed experimentally [13]. Our model reproduces further features seen in experiments and numerical simulations, including the high frequency fluctuations of the local shear [6].
In addition, we predict a complex self-organization of these shear bands, displayed in the inhomogeneities in the local density. This feature can be studied experimentally, in particular through the use of X-ray tomography, to access the local density of a sheared medium. The use of tracer particles could also be helpful in identifying the inhomogeneous ageing and localization of the shear bands as well as their sudden changes.
Most of the results presented here for the Euclidean lattice are closely mirrored by the results on the hierarchical lattice studied in [18]. The recursive topology of this lattice allows a quantitative analytical understanding of many of the quantities studied numerically in this paper. This includes elucidating the mechanism for the breakdown of ergodicity, the slow density evolution, as well as the behaviour of the of the Hamming distance at late times.