# Diffusion in mesoscopic lattice models of amorphous plasticity

###### Abstract

We present results on tagged particle diffusion in a meso-scale lattice model for sheared amorphous material in athermal quasi-static conditions. We find a short time diffusive regime and a long time diffusive regime whose diffusion coefficients depend on system size in dramatically different ways. At short time, we find that the diffusion coefficient, , scales roughly linearly with system length, . This short time behavior is consistent with particle-based simulations. The long-time diffusion coefficient scales like , close to previous studies which found . Furthermore, we show that the near-field details of the interaction kernel do not affect the short time behavior, but qualitatively and dramatically affect the long time behavior, potentially causing a saturation of the mean-squared displacement at long times. Our finding of a short time scaling resolves a long standing puzzle about the disagreement between the diffusion coefficient measured in particle-based models and meso-scale lattice models of amorphous plasticity.

###### pacs:

–Over the past few decades, the notion of local shear transformations has been used to describe and explain the plastic flow of amorphous solids Bulatov and Argon (1994); Falk and Langer (1998). A class of mesocopic lattice models is built on this picture Baret et al. (2002); Talamali et al. (2012, 2011); Tyukodi et al. (2016); Picard et al. (2004); Martens et al. (2011); Nicolas and Barrat (2013); Nicolas et al. (2014a, b); Lin et al. (2014a); Puosi et al. (2014); Budrikis and Zapperi (2013); Budrikis et al. (2017a), (see Nicolas et al. Nicolas et al. (2017) for a recent review). In these lattice models, the system is partitioned into local regions, and any one of them may undergo a yielding event if loaded beyond some threshold. These models are designed to operate at a mesoscopic scale; slightly coarser than the particles, but not at a macroscopic scale where continuum thermodynamical models describe phenomena such as persistent shear localization Falk and Langer (1998); Manning et al. (2007); Bouchbinder and Langer (2009).

Avalanches of local shear transformations are observed in both particle-scaleMaloney and Lemaître (2006); Hentschel et al. (2010); Salerno et al. (2012); Salerno and Robbins (2013) and meso-scale models Talamali et al. (2011); Budrikis and Zapperi (2013); Lin et al. (2014b); Tyukodi et al. (2016); Karimi et al. (2017); Budrikis et al. (2017a) during slow steady shear. The cascades are caused by the elastically mediated redistribution of stress after a local yielding event Bulatov and Argon (1994); Picard et al. (2004); Eshelby (1957). The result is a broad spectrum of bursts of plastic activity Antonaglia et al. (2014) and fractal patterns of accumulated plasticity Sun and Wang (2011). Similar avalanching behavior is observed in many different dynamically critical systems Bak et al. (1989, 1987); Tanguy et al. (1998); Pindra and Ponson (2017); Benassi and Zapperi (2011); Fisher (1998).

Despite the quantitative agreement in the spectrum of avalanche sizes and the qualitative agreement in the spatial correlations in the plastic strain Talamali et al. (2012), one major discrepancy between particulate and mesoscale models has remained. It involves the diffusion coefficient, , of the motion of tagged particles. Lemaître and Caroli Lemaître and Caroli (2009) argued that the spatial correlations in the plastic strain field should give rise to a dependence of on the system length, . In quasi-static simulations of a Lennard-Jones glass, Maloney and Robbins Maloney and Robbins (2008) showed that out to particles. In a lattice model, Martens et. al. Martens et al. (2011), found a very different scaling with system length for the diffusion coefficient, . Nicolas and co-workers then Nicolas et al. (2014a) showed that including the effects of advection changes the scaling and suggested including advection was necessary to obtain agreement with particulate models. However, very little quantitative reconciliation has been done between the meso-scale and particulate models even for this case of advection.

To shed light on these inconsistencies, we have performed an extensive set of simulations of a simple athermal quasi-static mesoscopic lattice model. We find a short time regime where , and a long time regime where . The short time diffusive plateau ends after a characteristic time , characterizing the strain released in a system spanning event, which shrinks with system size. This reconciles the results of reference Maloney and Robbins (2008) with the results of reference Martens et al. (2011); Nicolas et al. (2014a). The crossover to the long time regime occurs at a size independent strain of order unity. At the same time, we also find consistency with recent results of Tyukodi et. al.Tyukodi et al. (2016): at the very longest times, the variance of the plastic strain field, and, correspondingly, the mean squared particle displacement, either continues to grow diffusively or saturates depending on whether or not the load redistribution kernel possesses null modes. We further show that the kurtosis of the displacement distribution decays with the size of the time window in the same way as in atomistic simulations and argue that this is a generic consequence of the fact that the displacement field is built from temporally uncorrelated shot noise with the spatial structure of each shot being a characteristic system spanning avalanche.

The basic approach of the lattice models goes back to Eshelby who showed that the linear elasticity problem in which a local region undergoes a shift, , in its reference, stress-free configuration, is given by an integral convolution: where is the so-called Eshelby kernel.

Lattice models of amorphous plasticity then add a dynamical rule for the evolution of . One of our goals in this study was to develop a simple discretization of the Eshelby problem on a lattice which gives realistic displacements and compatible strains near the lattice site undergoing plasticity. Our approach is detailed in the supplemental material SI (), but briefly: i) we define our lattice model by partitioning space into square domains; ii) we define the strain on each square via a finite difference of a displacement field defined on the vertices of the square, and iii) the response, , to an increment of on a single square – i.e. the Eshelby kernel – is expressed analytically as a Fourier series on the square lattice and tabulated in real space for each lattice size, . This discretization scheme is similar to that used in studies of Martensitic transformations Kartha (1994); Shenoy et al. (1999) and to a scheme used recently in an amorphous lattice model by Nicolas Nicolas et al. (2015). At distances far from the yielding square, the shear component of this kernel gives the far-field solution of the Eshelby inclusion problem Eshelby (1957) and its shear component features a quadrupolar symmetry, i.e. in polar coordinates .

In this paper, we focus on two modes of shear with respect to the underlying lattice, ; ; which we call mode 2 and ; ; which we call mode 3. We assume the elastic constants have the Lamé form so that, in either case, , where is the elastic strain. We work in units where so that we can speak interchangeably of or . Note, denoting spatial averages with , that and would be the stress and strain of the sample measured by a load cell. As we will show below, despite residual correlations at long time in mode 2, the short time behavior of mode 2 and mode 3 in terms of the displacement and strain statistics and the avalanche spectrum (not studied here) is essentially indistinguishable.

We have studied different flavors of the model, characterized by different ways of introducing disorder and advancing the simulation in time.
For the stochastic ingredients, we have studied: i) random local stress thresholds, , with uniform increments in local plastic strain,
, and ii) random increments in with uniform .
For the dynamical rule we have used: i) an extremal protocol, where the total strain is adjusted uniformly across the system
(it may increase or decrease) at each step so that precisely one site is at threshold Baret et al. (2002), and ii) a synchronous protocol where all unstable sites are updated simultaneously while
is held fixed and this procedure is
iterated at the same until all sites become
stable before is incremented again (also uniformly as in
the extremal protocol). ^{1}^{1}1Note that in the latter
synchronous protocol avalanches are naturally produced at a discrete
set of with
monotonically
increasing during the protocol, while in the former extremal
protocol, need not
increase monotonically, and avalanches must be identified implicitly
a posteriori by grouping sets of extremal
events Talamali et al. (2011).
We have checked that the scaling exponents we define below do not depend on either
the stochastic model or dynamical update rule; although non-universal properties may. The data we present here
are for the case of random increments with uniform
and for the synchronous update protocol. We choose each
local increment of from a uniform distribution from
to .
For this study, we take
.
Because of the underlying
linearity, a shift in will simply shift all loading curves
and result in precisely the same sequence of shear transformations, so
we conventionally set and note that is the only non-trivial adjustable parameter in the
model Talamali et al. (2012); Budrikis et al. (2017b).

In Figs. 1a and 1b, we show the steady-state variance, , of the plastic strain field scaled by the length of the time window, , as a function of for various system lengths, , in mode 2 (a) and mode 3 (b) loading. For short times (small ), both modes of loading show a consistent, size independent, diffusion constant. We can make an ab initio estimate for the height of the plateau by assuming the probability distribution of local plastic strains is simply the uniform distribution corresponding to sites which have yielded precisely once plus a residue at zero corresponding to sites which have not yet yielded. This ab initio estimate gives a value of which is in excellent agreement with the measured plateau height. At later times, sites will eventually undergo more than one yielding event, and this estimate will break down.

There is a fall off from the plateau, starting at a strain of order (regardless of ) at which point each site has yielded approximately once on average. Beyond this fall from the plateau, the two loading modes show dramatically different behavior. The Mode2 curves all drop sharply. Each curve has a shoulder feature beyond which saturates and . The shoulder extends to longer for larger . The Mode3 curves show dramatically different behavior. After a subdiffusive regime, the curves again approach a diffusive plateau (with a lower diffusion constant than at short time), with the larger systems having a lower long-time diffusion constant. This behavior for mode 2 and mode 3 is consistent with Tyukodi et. al.’s recent work Tyukodi et al. (2016) where it was argued that the presence of null modes in the convolution operator (and associated stress-free slip lines) was necessary for to remain diffusive.

In Fig. 2, we show the diffusion coefficient, , of the displacement field scaled by vs . This rescaling collapses the data onto a short-time master curve which has the same shape for both loading modes. For short times, there is a diffusive plateau. As increases, the curve departs upward, superdiffusively, from the plateau. This superdiffusive regime sets in at a characteristic time scale when: .

To explain the scaling with , we recall and generalize the arguments of reference Maloney and Robbins (2008) which were motivated by reference Lemaître and Caroli (2009). In Fig. 3, we plot the incremental stress field in mode 3 loading for several consecutive non-overlapping time windows of a size corresponding to the end of the lower plateau in figure 2 at the initial stages of the superdiffusive regime. Similar features are observed for the other loading mode but rotated by degrees. The plasticity is organized into line-like features (either vertical or horizontal) which correspond to the directions where is large and positive.

Suppose the field for a typical time window at
short time is either zero if there has been no plasticity or composed
of a perfect line spanning the simulation cell if there has been
plasticity. On average, each site on the line has
^{2}^{2}2The factor of comes
from the uniform distribution from to .. Since
there are such sites in the line, the whole line will globally
relieve a strain precisely equal to:
(where is the lateral size of a
square element of the lattice). The displacement field,
, associated with that slip line is a linear profile
with a strain equal to so that (assuming, for
the sake of argument, a horizontal slip line centered at )
.

The variance of this displacement field is which is independent of . The rate at which these slip lines occur per unit strain, , has to be precisely enough so that, on average, so . If we are in a short time regime so that at most one of these slip lines has formed, then we have: So the simpleminded picture of elementary lines predicts a short-time characteristic strain, and a short-time .

In Fig. 2, we see that the scaling is only approximately correct and that scaling by and by gives a better quality data collapse for both mode 2 and mode 3 loading. We can explain this by slightly generalizing the argument above. If we imagine the short-time windows contain either no plasticity or a characteristic elementary event, then we still have that where is still the rate of events and is still the variance of a characteristic event. And we still must have balance between applied strain and plastic strain so that ; but now is the characteristic strain associated with an arbitrary characteristic event more general than a straight line: where is the number of sites involved in one of the elementary events. For lines, , while we generalize and let for fractal objects. So for the characteristic strain associated with an elementary line-like object, we have . And, finally, for the diffusion, we have . We must assume that the elementary events produce displacement fields whose variance is independent of , but given that assumption, we see that and . From our scaling collapse, we conclude that which would correspond to a fractal dimension of .

In Fig. 4 we again plot the diffusion coefficient vs the strain , but now with scaled by to collapse the upper plateau at long time. We see a crossover to the upper plateau at a strain of order unity regardless of . This occurs after the departure of from its short time plateau. The mode 3 case displayed in Fig 4 (b) remains perfectly diffusive for as long as we can simulate and we have no reason to believe it will do otherwise.

The mode 2 case shown in Fig 4 (a) shows a strikingly different behavior. For any finite , the mean square displacement eventually saturates and decays like at long enough . Despite this decay, we observe the emergence of an apparent upper plateau before the decay even for mode 2 at sufficiently large . Furthermore, the height of the plateau seems to obey the scaling as well.

The emergence of the upper plateau in mode 2 appears to be related to the spectral gap in the Eshelby convolution operator disappearing in the limit. For sufficiently large systems, there will be little difference between mode2 and mode3 loading, despite the lack of zero modes of the Eshelby convolution operator in the former, as long as remains below the onset of decay.

In figure 5, we plot the kurtosis, , of the
distribution of the (a) and (b) Cartesian
component of the displacement versus for
mode 3 ^{3}^{3}3We find virtually identical behavior for mode2 under
rotation of the Cartesian components by degrees. Both plots
show a striking initial behavior as observed
earlier in Lennard-Jones glasses Tsamados (2009); Tsamados et al. (2009).
The behavior can be explained as follows.
Suppose the displacement field is built up from a succession of
characteristic events which are spatially uncorrelated with each
other. Furthermore, suppose we are interested in a timescale
for which it is unlikely to observe more than one
event in a given window. Then a typical window of duration
contains either one event (with probability
) or no event (with probability
) where is the
characteristic strain release in the event. So any particular moment
of the distribution (in particular the second and fourth) should scale
like . So in
particular, for the kurtosis, which is precisely
what we see and explains the much earlier atomistic results from
Tsamados et. al.Tsamados (2009); Tsamados et al. (2009): at long times we recover , an indication of a Gaussian-like distribution. We would
expect our data for different system sizes to collapse when rescaled
by the characteristic strain, , which was found above, in
the analysis diffusion coefficient, to scale like .
However, we find the best collapse when is scaled by
and the discrepancy between the characteristic strain
inferred from the diffusion coefficient and the kurtosis remains an
outstanding puzzle.

To summarize, we have shown that lattice models for athermal quasistatic amorphous plasticity show good agreement with particle-based simulations for the system size dependence of the short-time diffusion coefficient. This is for two different stochastic prescriptions for the local energy landscape: random threshold or random plastic strain increment; two different dynamical update rules: synchronous or extremal; and two different orientations of the loading with respect to the lattice. Our results are also in agreement with Maloney and Robbins Maloney and Robbins (2008) who showed that the variance of the local strain field shows little size dependence, while the variance of the displacements shows dramatic size dependence.

At longer times, the diffusion coefficient shows a size dependence, which is similar to the observed by Martens et. al. Martens et al. (2011) In this long time regime, the behavior is different for the two different modes of loading. When loading along the axes of the lattice, the discretized Eshelby kernel has no null modes, so the variance of the plastic strain, and thus the variance of the displacements, saturates. When loading degrees away, the Eshelby kernel has proper null modes – perfect slip lines along the lattice axes which leave the stress field uniform, so the variance can continue to grow and the system can achieve a proper diffusive limit in agreement with earlier arguments by Tyukodi et. al.. Tyukodi et al. (2016) We note that even in the axial-load case where there are no perfect null modes of the kernel, a pseudo-diffusive-plateau develops at the very latest times. The extent of the pseudo-diffusive-plateau depends on system size with larger sizes maintaining a quasi-diffusive regime for a longer period of time, but a precise study of the long-time diffusive behavior is left for future work.

The present picture we put forward here of a separate early time diffusive regime crossing over to a distinct late time diffusive regime clarifies the apparent discrepancy between particle-based and lattice-based models. In particular, it appears that the introduction of convection discussed by Nicolas et. al. Nicolas et al. (2014a) is not necessary to recover the linear size scaling of the correlations observed in atomistic simulations. In light of our present work, it seems likely that a regime was already present in former advection-free lattice models based studies but that this early difusive regime was simply not analyzed. Of course, at very late times, advection should be important for a detailed comparison with particle-based simulations.

###### Acknowledgements.

This material is based upon the work supported by the National Science Foundation under Grant No. DMR-1056564 and in part by Grant No. PHY-1748958. CEM would like to acknowledge discussions with Kirsten Martens which motivated this study.## References

- Bulatov and Argon (1994) V. V. Bulatov and A. S. Argon, Modelling and Simulation in Materials Science and Engineering 2, 167 (1994).
- Falk and Langer (1998) M. L. Falk and J. S. Langer, Phys. Rev. E 57, 7192 (1998).
- Baret et al. (2002) J. C. Baret, D. Vandembroucq, and S. Roux, Phys. Rev. Lett. 89, 195506 (2002).
- Talamali et al. (2012) M. Talamali, V. Petäjä, D. Vandembroucq, and S. Roux, C.R. Mécanique 340, 275 (2012).
- Talamali et al. (2011) M. Talamali, V. Petaja, D. Vandembroucq, and S. Roux, Phys. Rev. E 84, 016115 (2011).
- Tyukodi et al. (2016) B. Tyukodi, S. Patinet, S. Roux, and D. Vandembroucq, Phys. Rev. E 93, 63005 (2016).
- Picard et al. (2004) G. Picard, A. Ajdari, F. Lequeux, and L. Bocquet, Eur. Phys. J. E 15, 371 (2004).
- Martens et al. (2011) K. Martens, L. Bocquet, and J. L. Barrat, Phys. Rev. Lett. 106, 156001 (2011).
- Nicolas and Barrat (2013) A. Nicolas and J. L. Barrat, Faraday Discussions 167, 567 (2013).
- Nicolas et al. (2014a) A. Nicolas, K. Martens, L. Bocquet, and J. L. Barrat, Soft Matter 10, 4648 (2014a).
- Nicolas et al. (2014b) A. Nicolas, K. Martens, and J. L. Barrat, EPL (Europhys. Lett.) 107, 44003 (2014b).
- Lin et al. (2014a) J. Lin, E. Lerner, A. Rosso, and M. Wyart, Proc. Natl. Acad. Sci. 111, 14382 (2014a).
- Puosi et al. (2014) F. Puosi, J. Rottler, and J. L. Barrat, Phys. Rev. E 89, 042302 (2014).
- Budrikis and Zapperi (2013) Z. Budrikis and S. Zapperi, Phys. Rev. E 88, 062403 (2013).
- Budrikis et al. (2017a) Z. Budrikis, D. F. Castellanos, S. Sandfeld, M. Zaiser, and S. Zapperi, Nature Commun. 8, 1 (2017a).
- Nicolas et al. (2017) A. Nicolas, E. E. Ferrero, K. Martens, and J.-L. Barrat, (2017), arXiv:1708.09194 .
- Manning et al. (2007) M. L. Manning, J. S. Langer, and J. M. Carlson, Phys. Rev. E 76, 056106 (2007).
- Bouchbinder and Langer (2009) E. Bouchbinder and J. S. Langer, Phys. Rev. E 80, 031132 (2009).
- Maloney and Lemaître (2006) C. E. Maloney and A. Lemaître, Phys. Rev. E 74, 016118 (2006).
- Hentschel et al. (2010) H. G. E. Hentschel, S. Karmakar, E. Lerner, and I. Procaccia, Phys. Rev. Lett. 104, 1 (2010).
- Salerno et al. (2012) K. M. Salerno, C. E. Maloney, and M. O. Robbins, Phys. Rev. Lett. 109, 105703 (2012).
- Salerno and Robbins (2013) K. M. Salerno and M. O. Robbins, Phys. Rev. E 88, 062206 (2013).
- Lin et al. (2014b) J. Lin, A. Saade, E. Lerner, A. Rosso, and M. Wyart, EPL (Europhys. Lett.) 105, 26003 (2014b).
- Karimi et al. (2017) K. Karimi, E. E. Ferrero, and J. L. Barrat, Phys. Rev. E 95, 013003 (2017).
- Eshelby (1957) J. D. Eshelby, Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences 241, 376 (1957).
- Antonaglia et al. (2014) J. Antonaglia, W. J. Wright, X. Gu, R. R. Byer, T. C. Hufnagel, M. LeBlanc, J. T. Uhl, and K. A. Dahmen, Phys. Rev. Lett. 112, 155501 (2014).
- Sun and Wang (2011) B. A. Sun and W. H. Wang, Applied Physics Lett. 98, 201902 (2011).
- Bak et al. (1989) P. Bak, K. Chen, and M. Creutz, Nature 342, 780 (1989).
- Bak et al. (1987) P. Bak, C. Tang, and K. Wiesenfeld, Phys. Rev. Lett. 59, 381 (1987).
- Tanguy et al. (1998) A. Tanguy, M. Gounelle, and S. Roux, Phys. Rev. E 58, 1577 (1998).
- Pindra and Ponson (2017) N. Pindra and L. Ponson, Phys. Rev. E 053004, 1 (2017).
- Benassi and Zapperi (2011) A. Benassi and S. Zapperi, Phys. Rev. B 84, 1 (2011).
- Fisher (1998) D. S. Fisher, Physics Reports 301, 113 (1998).
- Lemaître and Caroli (2009) A. Lemaître and C. Caroli, Phys. Rev. Lett. 103, 1 (2009).
- Maloney and Robbins (2008) C. E. Maloney and M. O. Robbins, J. Phys. Cond. Matt. 20, 244128 (2008).
- (36) See Supplemental Material at http://www.SI-address.
- Kartha (1994) S. Kartha, World, Ph.D. thesis (1994).
- Shenoy et al. (1999) S. Shenoy, T. Lookman, A. Saxen, and A. Bishop, Phys. Rev. B 60, R12537 (1999).
- Nicolas et al. (2015) A. Nicolas, F. Puosi, H. Mizuno, and J. L. Barrat, J. Mech. Phys. Solids 78, 333 (2015).
- (40) Note that in the latter synchronous protocol avalanches are naturally produced at a discrete set of with monotonically increasing during the protocol, while in the former extremal protocol, need not increase monotonically, and avalanches must be identified implicitly a posteriori by grouping sets of extremal events Talamali et al. (2011).
- Budrikis et al. (2017b) Z. Budrikis, D. F. Castellanos, S. Sandfeld, M. Zaiser, and S. Zapperi, Nature Comm. 8, 15928 (2017b).
- (42) The factor of comes from the uniform distribution from to .
- (43) We find virtually identical behavior for mode2 under rotation of the Cartesian components by degrees.
- Tsamados (2009) M. Tsamados, PhD Thesis, Ph.D. thesis (2009).
- Tsamados et al. (2009) M. Tsamados, A. Tanguy, C. Goldenberg, and J. L. Barrat, Phys. Rev. E 80, 026112 (2009).