Correlation and shear bands in a plastically deformed granular medium

Recent experiments (Le Bouil et al., Phys. Rev. Lett., 2014, 112, 246001) have analyzed the statistics of local deformation in a granular solid undergoing plastic deformation. Experiments report strongly anisotropic correlation between events, with a characteristic angle that was interpreted using elasticity theory and the concept of Eshelby transformations with dilation; interestingly, the shear bands that characterize macroscopic failure occur at an angle that is different from the one observed in microscopic correlations. Here, we interpret this behavior using a mesoscale elastoplastic model of solid flow that incorporates a local Mohr-Coulomb failure criterion. This differs from the interpretation of Le Bouil et al., which is based on purely elastic considerations ignoring the potential role of local friction on deformation patterns. We show that the angle observed in the microscopic correlations can be understood by combining the elastic interactions associated with Eshelby transformation with the local failure criterion. At large strains, we also induce permanent shear bands at an angle that is different from the one observed in the correlation pattern. We interpret this angle as the one that leads to the maximal instability of slip lines.

event. This preferential direction (with respect to the principal axis of the local shear event) is also the one along which a shear band should form in a system deformed at constant volume (in particular incompressible), as it corresponds to the direction of maximum macroscopic shear stress.
However, the symmetry that Le Bouil et al. 14 (see also ref. 15 for the experimental setup) observed in correlation patterns of the plastic activity seems to be at odds with the morphology of the fully formed band they also observed in strongly deformed granular media. While the latter was formally described as reflecting the Mohr-Coulomb macroscopic failure angle, the former revealed a distinct orientation. Moreover, both angles differ from the "canonical" 45°. In order to rationalize the preferential direction observed in correlations, the authors proposed to take into account the possibility of a local dilation in the plastic event, a feature that was already noted to affect the preferential directions in collective bursts plastic activity 16 . The authors proposed a phase co-existence scenario in which recurring mini bursts persist up until failure, which may be interpreted as a signature of discontinuous first-order or spinodal transition. However, this picture does not lead to a specific prediction for the angle of global failure and emergence of macroscopic friction in the system.
In this work, our aim is to provide, based on minimal ingredients, a scenario that explains simultaneously the deviation from the "canonical" 45° direction for the different observables, and the fact that different directions are observed for the correlations in intermittent activity and in permanent shear localization. Our analysis will be based on the hypothesis that, while the interactions between plastic events are mediated by the universal laws of linear elasticity, the triggering of events depends on a local failure criterion that is completely independent of these laws, and can be quite arbitrary. In granular media, for instance, with frictional resistance being an important mechanical property, the yield criterion may be stated as a critical ratio between the resolved shear and normal stress at any arbitrary material point. As a result, the correlations in plastic activity and its propagation involve a compromise between directions favored by the elastic interactions and those favored by the local failure criteria. The results that emerge from this compromise are nontrivial, in agreement with the experimental observations. From Le Bouil et al. perspective, friction is a topological property that emerges at global scales leading to a Mohr-Coulomb type failure; it seems to play no role in the dynamics of transient slip bands. In our approach, friction has a mesoscopic notion as well and, along with elastic couplings, governs triggering dynamics at intermediate scales.
This results in a frictional flow which is collective at macroscopic scales and appears as a shear-banding phenomenon with extended length-and/or time-scales associated with it. We also find that the "bulk" friction coefficient is different than the local ingredient which suggests a separation of scales between the meso and macro view of the plastic flow.
The minimal ingredients used in our analysis are the description of local plastic events as transformations of Eshelby-type 17 , combined with the use of a local Mohr-Coulomb criterion that introduces a pressure sensitivity in the propagation of plasticity. It can be seen as an extension to the case of Mohr-Coulomb failure of works that involve permanent damage as a cause of strain localization [18][19][20] . Remarkably, however, no permanent damage is required for inducing coexistence between transient micro events and fully developed shearing bands.
The organization of the paper is the following. In Sec. 1 we combine the hypothesis of Mohr-Coulomb local failure with the stress redistribution prescribed by Eshelby's elastic theory, and discuss hypothetically the consequences on the correlation patterns and macroscopic failure angle. In Sec. 2 we describe the numerical model that incorporates these basic ingredients. The numerical results and conclusions are given in Sec. 3 and Sec. 4, respectively.

The Mohr-Coulomb Failure Criterion
Previous studies within the framework of mesoscopic elasto-plasticity 2,21 emphasized the role of elastic kernels in the yielding transition. The local yielding rule however, coupled with long-range elasticity, must have a strong relevance on the structure of correlations and macroscopic failure. In this section, we propose a simple theory taking into account these effects at both local and macroscopic scales. At the mesoscopic scales, we made a common assumption that, on average, the far-field stress patterns follow the classical Eshelby's solution in response to a localized shear. The actual response, however, displays significant fluctuations due to the disordered and discrete nature of amorphous media. Interestingly, statistical quantities at large scales, i.e. plastic-activity correlations, appear to be almost insensitive to the quenched disorder and can be, therefore, replicated within the context of elasto-plastic models and continuum framework.

Meso-scale Failure Criterion and Elastic Stress Redistribution.
We consider an infinite elastic matrix characterized by its bulk modulus K and shear modulus μ and an embedded inclusion going through a shear with an amplitude ε * and microscopic volume a d . Here δ αβ is a Kronecker delta and we focus on the two dimensional case d = 2.
The stress tensor at a material point can be expressed as σ αβ = −pδ αβ + σ(δ αx δ βx − δ αy δ βy ) + σ xy (δ αx δ βy + δ αy δ βx ) in terms of the pressure p and area preserving axial shear σ and diagonal shear σ xy . The far-field perturbation fields given by Eshelby The power-law decay with distance r marks the non-local long-range nature of the disturbance that also has distinct angular symmetries in terms of θ. Note the four-fold symmetry in σ which is contrasted by the bi-polar structure in p as illustrated in Fig. 1(a) and (b). The negative (blue) and positive (red) lobes represent regions with decreasing and increasing stresses in space, respectively. It should be noted that the stress patterns in a disordered solid become scattered due to the spatial heterogeneity and/or anisotropy in local elastic properties. It is only beyond a characteristic length-scale that fluctuations average out spatially and the amorphous system, at macroscopic scales, has a smooth elastic response.
If at a point within the bulk the shear stress on any plane becomes equal to the shear strength, failure will occur at that point. In the following, we will assume that the local shear strength can be expressed as a linear function of the confining pressure 22 , in accordance with the Mohr-Coulomb concept. Therefore, the distance to failure at any point can be expressed by the yield function f y : Assuming σ xy = 0, the far-field perturbations in f y to leading order following a shear transformation taking place at the origin is given by Obviously, failure is likely to localize into the overlapped sectors between positive lobes of δσ in Fig. 1(a) and negative lobes of δp in Fig. 1(b). Figure 1(c) displays perturbations in f y for a value φ = 65° of the local friction. Near the shear transformation site, the pattern retains the quadrupolar shape with its positive lobes slightly tilted upwards toward regions with decreasing pressure. Note that a positive δf y means that material points are pushed toward failure threshold. The direction of the maximum in the positive lobes (denoted by θ max hereafter) will depend on the internal friction φ and is given by for low values of φ. In the absence of friction, i.e. φ = 0, we recover θ max = 45° as reported in many glassy systems 23 . This implies that the change in f y takes its maximal value along the direction of maximum change in shear stress.
In addition to distortion, we now suppose that transformation zones may undergo dilation as well, i.e.
This will make an anisotropic contribution with the two-fold symmetry to the shear stress δσ Combining the effects of the local dilation and shear yields In comparison with Eq. 4, θ max in the above equation has an extra ingredient that includes ε ε ⁎ ⁎ 2 v the ratio between volumetric and shear strains. In the limit → ε ε , the term involving friction becomes dominant and we will retrieve Eq. 4.
Macroscopic Failure: Slip-line Instability Analysis. In a highly simplified setting, a fully operating band may be idealized as a quasi-linear object made up of evenly distributed Eshelby events 16,24 . Using an eigenmode based strategy, we now propose a simple theoretical argument that makes predictions on the band alignment. Point-like events with volume a d and amplitude ε * are evenly positioned at points → r i on a line at an angle α with the x axis: Taking a Fourier transform, the q-space representation reads q is the unit vector along → q , δ is the Kronecker symbol, and N is the total number of events. In order to keep the strain finite in the continuum limit, Na d is constant with N → ∞ and a d → 0. Here ε → αβ q ( ) takes a localized form as well but at a direction perpendicular to the band orientation in real space. Assuming pure shear distortion for the transformations, It is clear that the solutions in Fourier space will only depend on the direction θ (not on → q ) and take the exact same form as the effective source in Eq. 7, up to a normalization factor. In other words, the slip line does not induce any stress redistribution outside of the line itself, and the real-space response function is localized along the slip line at angle α As a result, the change in the yield function δf y is zero everywhere except on the slip line. We now argue that the most unstable lines will be those that maximize δf y , which is positive inside the band. This corresponds to a maximum amplification of the local "damage" that is caused by the operating shear band for an infinitesimal strain. The shear band angle θ sh is thus obtained by maximizing δf y with respect to α: This angle lies between 45° ≤ θ sh ≤ 60°, and for small φ is approximately θ sh ≈ 45° + φ 4 . The combination of the pure shear deformation and volumetric change in the transformation zones gives Note that a similar result in the case without friction was obtained in ref. 16 on the basis of an energy minimization argument. Inserting Eq. 10 into Eq. 9 and integrating δp over the active domain, it becomes evident that the total change in pressure is negative. This will effectively reduce the shear strength given by the Mohr-Coulomb yield surface inside the localization zone. With such a mechanism, yielded zones become likely places where next events tend to localize, hence making possible a permanent band-like formation. This is compatible with the numerical observations presented in the following sections.
From our discussions above it follows that, theoretically, θ max , the direction resulting from expected elastic-type correlations, will differ from the slip-line angle θ sh . This discrepancy appears to be consistent with the observations made in the granular experiment of Le Bouil et al. 14 . Note that the interpretation proposed by these authors was initially based on purely elastic considerations, with no allusion to the local friction. In this analysis, a large local dilatancy had to be assumed in order to make a sensible prediction of θ max . Our approach, by introducing the local friction angle, allows one to obtain a similar order of magnitude without invoking a large dilation. In a follow-up work 25 , it was suggested that the material anisotropy (or its combination with the volumetric strain) may explain the characteristic angle. The latter assumption was largely attributed to the presence of force chains that build up upon loading granular solids 26 . Our interpretation is somewhat different in that it is entirely based upon the dominant role of the friction angle on correlation patterns, and in addition makes a specific prediction for the orientation of the permanent shear bands.

Simulation Details
In order to test the predictions of the theoretical discussion above, we have made use of the Finite Elements based version of elasto-plastic models that we originally established in 6,27 . In our previous work, this model was used to study the deformation of amorphous media in which the local failure is governed by a standard, maximum shear stress criterion. Here, in order to follow the hypothesis of the previous section, the microscopic failure criterion for each element will be a Mohr-Coulomb condition. Each material point is therefore assigned with the shear strength parameters c and φ; while the former is randomly chosen from an exponential distribution with the mean value c , we allocate no disorder to the latter.
Below the failure limit at the material point in question, set by f y = 0 in the pressure-shear stress plane, the local stress trajectory is regulated by the imposed drift and non-local elastic interactions. We further presume a linear isotropic elastic response for the pre-failure dynamics. Figure 2 plots τ m = 1 2 (σ 3 − σ 1 ) against = − p 1 2 (σ 3 + σ 1 ) and represents any stress state by a stress point. Here σ 1 and σ 3 denote major and minor principle stresses, respectively. Upon yielding, the stress point will take on a path which is perpendicular to the p-axis and relax visco-elastically 6 toward a point representing the final state of stress. The released shear stress Δσ will create a localized net force that perturbs the force equilibrium in the medium. The perturbation will be of the generic form given in Eq. 1 plus a scattering term caused by currently yielded elements with decaying shear moduli 28 . The homogeneity approximation and subsequent theoretical derivations in the preceding section appear to be legitimate at least with regard to the fluctuation patterns and shear band morphology evaluated numerically in the following sections.
Simulations of shear deformations were performed by applying an area preserving axial shear rate ε ∼ −  10 4 to an L × L periodic cell. An irregular set of triangular elements with average size By setting a high value of the damping rate (in comparison with the vibrational frequency), we moreover impose an overdamped behavior during the relaxation phases, thus ensuring that inertial effects are negligible 6 . Prior to shearing, samples were prepared with random stresses assigned to each block followed by an equilibration within a purely elastic framework (no plastic events allowed) that resulted in a state of mechanical equilibrium.

Numerical Results
The results of shearing tests can be used to determine the bulk shear strength characteristics together with the structure of interactions between transient slip events. A number of tests has been performed at the same (initial) pressure = p c 8 and φ = 65°, and the resulting average stress strain curve is displayed in Fig. 3(a). For the given p and φ, the material shear strength has the average value of σ ≈ c 8 y . We also report ε in units of the yield strain ε = σ μ y 2 y . A peak stress is typically reached around ε ≈ ε y in every sample followed by a reduction in strength as the loading continues. With increasing strain, the shear strength ultimately falls to a residual value at large deformations.
A strong plastic activity in the form of extended linear structures is present at all times after the initial yielding for ε > ε y . Spatial maps of active sites in Fig. 3(b-e) demonstrate the highly intermittent and non-local nature of bursts during plastic flow (see Supplementary movie). Following the stress peak, multiple shearing bands start to evolve within which most of plastic activities are taking place. The orientation of these bands are quite distinct from maximum shearing directions, 45° or 135° in our loading set-up. The quantification of the slip directions will be the subject of the next section, which presents an analysis similar in spirit to 29 . Structural Characterization. Here we focus on the number density of active sites and associated spatial correlations so as to quantify the spatial structure of the plastic activity. Let ρ δ . These variations can be characterized by a two-point density correlation function between two different positions → r and → ′ r in space: Here the angular brackets 〈...〉 correspond to an average over different realizations and a spatial average. As a result of translational invariance → → ′ ≡ → − → ′ S r r S r r ( , ) ( ). A Voronoi cell analysis was performed that enabled interpolations of the density field onto fine regular grids. We subsequently used a 2 d Fourier transform in order to compute the correlations.
We shall naturally expect that density fluctuations are highly correlated along shear band directions, inducing strong anisotropies in → − → ′ S r r ( ) . Figure 4 displays the evolution of density correlations at different loading stages averaged over 16 samples at = p c 8 and φ = 65°. Angular symmetries seen in the post-failure regime are fairly stable features showing only weak fluctuations with increasing strain. The banded regions contain correlations that are explicitly longer-ranged (as opposed to other orientations) indicating the system-spanning nature of the localization. Figure 5 displays the anisotropic part S(θ), an averaged → − → ′ S r r ( ) over different distances → − → ′ r r . There are two marked maxima in every data set at 45° < θ < 135° that delineate the degree of anisotropy. We quantify the positions of these peaks by fitting the data to a sum of two Gaussian peaks, as in Fig. 5(b-e). These fits enable us to obtain the locations of the peaks-denoted by θ sh -and to follow their evolution upon shear loading. This is shown in Fig. 5(a). Apart from the initial transient part prior to failure, the peak locations are essentially strain independent at larger strains. The dashed lines in Fig. 5(a) mark the directions of failure given by Eq. 4.
The basis for the permanent localization observed here can be understood in simple terms. An incident avalanche tends to depressurize currently-damaged blocks, with a pressure drop proportional to sin φ, making them vulnerable spots against further deformations. In the standard framework of elastoplastic models, it is possible to observe a similar behavior (however with the classical 45° orientation for the shear band by adding some permanent or transient weakening mechanism to the model). Examples include models based on the damage factor 18 , weakened stress thresholds 19 , or lingered restoration time 21 , just to name a few. In our simulations, the weakening . .  . .  As stipulated by the Mohr-Coulomb phenomenology, the theoretical angle between the major principal stress direction and the plane of failure must be θ = + φ  45 MC 2 r . Inserting φ r = 40°, it follows that θ MC = 65° which is off from θ sh by about 6°. Therefore, we find that applying the macroscopic Mohr-Coulomb theory using the observed values of the macroscopic failure envelope tends to overestimate the actual inclination of the yield surface, in agreement with experimental observations 30 .
Pre-failure Patterns of Plastic Activity. Beyond the yield strain, correlations in plastic activity are easily identified based on the analysis of a snapshot taken at a single time or strain, as described above using the function → − → ′ S r r ( ) . This is related to the existence of well correlated linear regions of plastic activity that are operating simultaneously, eventually giving rise to localized shear bands. In the pre-failure regime, however, plastic activity is much more scattered, and the study of a single configuration does not reveal any established pattern. A calculation of → − → ′ S r r ( ) does not allow one to identify any preferred direction. The structure of the plastic activity, however, can be revealed by the study of two time correlations between configurations separated by a fixed strain interval Δε. In practice, the statistics is improved by averaging over different simulation and performing an average over a strain window for the initial configuration. We define here the brackets denote the average over different realizations. The averaging interval [ε 1 ,ε 2 ] is taken before the strain peak, ε 1 = 0.25ε y and ε 2 = 0.7ε y , in order to avoid the contamination of the correlation functions by the formation of permanent shear bands in the post yield regime. Statistically speaking, ε → Δ C r ( , ) will provide spatial details about the most likely position of an event that is triggered following a local slip event at the origin. In Fig. 7(a), correlations are displayed for a small strain difference Δε = 5 × 10 −4 ε y . The results exhibit a four-fold structure, which persists up to a strain interval of the order 10 −3 before fluctuations become uncorrelated at higher strain differences. The four-fold structure is of course expected from standard elasticity theory, and has been observed in experiments and simulations of a number of glassy systems. However, we observe here a marked deviation from the 45° that would be predicted by pure elasticity, reflecting the influence of the friction angle in the yield criterion.
In Fig. 7(b), the angular dependence of ε → Δ C r ( , ) is shown together with a fit to a Gaussian function, similar to the one used for the analysis of S(r, θ). The fit function has a peak centered around 52° which is in close agreement with the theoretical prediction θ max . Note that in order to reduce the noise in the data, the integration over r was limited to small values of the distance < r L 16 . The structure of correlations preceding the ultimate failure reflects the hypothesis of our model. A localized event induces stress fluctuations in the surroundings triggering plastic rearrangements nearby in the medium. events are preferentially triggered in directions in which the changes to the yield function are maximal, as a result from non-local elastic couplings combined with a local frictional yielding rule. While one may have expected that these fluctuations would strongly influence the collective behavior that emerges upon macroscopic failure, the analysis of the shear band orientations shows that they represent a distinct phenomenon, and must be analyzed within a more collective perspective.

Conclusion and Discussion
In this work, we have taken a coarse-grained, mesoscopic view of the deformation in a disordered granular medium, by integrating the notion of a deformation taking place through local shear transformations with that of a local criterion of failure described by the Mohr-Coulomb friction condition. As a result, the mechanical response is described by two ingredients: (i) the elastic moduli of the medium which, according to the Eshelby's picture, describe the response to a localized shear transformation and (ii) the local friction angle that characterizes the Mohr-Coulomb condition. The latter feature is, in fact, material specific and supplements Le Bouil et al. 's line of reasoning which relies entirely on non-local couplings mediated by elasticity 14 .
Based on this picture, we proposed a scenario in which spatial correlations of the scattered plastic activity that takes place before the yield point, and the spatial structure of the permanent failure planes that can be observed beyond the yield point are explicitly calculated as a function of these simple ingredients, but occur at different directions. The discrepancy between the morphology of permanent bands and transient correlations is in full agreement with Le Bouil et al. 14 , who proposed a scenario based on the phase-coexistence between correlated mini bursts and persistent shearing bands. Within this picture, the flowing, post-yielding state is intrinsically different from the states explored at small strains, and the yield transition appears similar to a first order, discontinuous transition.
The failure mechanism we have proposed can be interpreted as the maximal instability of uniform lines of slip. Localization of deformation in this form is favorable, as no further stress is built up in the surrounding medium 21 . Furthermore, the intense shearing effectively lowers the yielding strength inside the localization band giving rise to a weakening process. The maximally unstable modes are dictated by both the elastic interactions 31 and the Mohr-Coulomb plasticity criterion (more specifically local friction) and may be viewed as a major source of mechanical instability 32 . It is noteworthy that the incurred damage is an ingredient of the model through the failure criterion, and does not enter as an extra material parameter.
The applicability of our analysis to a real granular medium may be questioned on two accounts: firstly, one may argue that a real granular medium is not described by an linear elastic continuum, due to the existence of force chains and micro plasticity in the form of contact breaking and formation. Second, the very existence of a Mohr-Coulomb criterion at the local scale is a rather arbitrary assumption, although in general a pressure sensitivity of the shear transformation is expected. Therefore, we have numerically tested our proposition by building a lattice-based model that incorporates rigorously, if perhaps artificially, these ingredients. The results of these simulations are in good agreement with the theoretical expectations, and, when compared to experiments, can lead to a prediction of an effective local friction coefficient.
In addition to confirming the theoretical analysis, the simulations offer the possibility to perform an analysis of the macroscopic failure envelope. This analysis yields an orientation, labeled as θ MC , that is incompatible with θ sh , the shear inclination. While the former is purely based on the measured bulk stresses, the later arises from kinematic considerations only. This discrepancy is in contradiction with the Mohr-Coulomb phenomenology at the global scale. Similar experimental observations were made (see 30 and the references herein) leading to ad-hoc remedies that incorporate extra material parameters, i.e. the dilatancy angle 33 , into the theoretical framework.