# On the critical nature of plastic flow: one and two dimensional models

## Abstract

Steady state plastic flows have been compared to developed turbulence because the two phenomena share the inherent complexity of particle trajectories, the scale free spatial patterns and the power law statistics of fluctuations. The origin of the apparently chaotic and at the same time highly correlated microscopic response in plasticity remains hidden behind conventional engineering models which are based on smooth fitting functions. To regain access to fluctuations, we study in this paper a minimal mesoscopic model whose goal is to elucidate the origin of scale free behavior in plasticity. We limit our description to fcc type crystals and leave out both temperature and rate effects. We provide simple illustrations of the fact that complexity in rate independent athermal plastic flows is due to marginal stability of the underlying elastic system. Our conclusions are based on a reduction of an over-damped visco-elasticity problem for a system with a rugged elastic energy landscape to an integer valued automaton. We start with an overdamped one dimensional model and show that it reproduces the main macroscopic phenomenology of rate independent plastic behavior but falls short of generating self similar structure of fluctuations. We then provide evidence that a two dimensional model is already adequate for describing power law statistics of avalanches and fractal character of dislocation patterning. In addition to capturing experimentally measured critical exponents, the proposed minimal model shows finite size scaling collapse and generates realistic shape functions in the scaling laws.

###### keywords:

plasticity, dislocations, self-organized criticality, Frenkel-Kontorova model, statistics of avalanches, intermittency, power laws^{1}

## 1 Introduction

It is known that crystalline and amorphous solids flow plastically when macroscopic stresses exceed certain thresholds. In plasticity, thermal relaxation may often be neglected while the driving can still be viewed as quasi-static. In this case the flow is an athermal process with rate independent dissipation. The singularity of this dissipative mechanism distinguishes plasticity from the regular dissipative phenomena such as viscosity or heat conductivity, and places it instead in the class of dry friction and fracture Cottrell (1953); Nadgorny (1988); Lubliner (1990); Kachanov (2004); Falk and Langer (2011).

At the macroscale rate independent plasticity can be perceived as a smooth process, whereas at microscale the plastic flow is known to exhibit fluctuations revealing isolated dissipative events Brechet (1994); Puglisi and Truskinovsky (2005). In crystals these fluctuations can be linked to nucleation and collective depinning of dislocations interacting through long range elastic stress fields Nabarro (1967); Hirth and Lothe (1982). In amorphous and random granular systems the nature of the ’quantum’ of plastic strain is still debated Argon (1979); Maloney and Lemaitre (2006); Liu and Nagel (2010); Barrat and A.Lemaitre (2011); Talamali et al. (2011) and therefore in what follows we mostly limit our discussion to crystals.

The smooth description of crystal plasticity adopted in engineering theories presumes spatial homogenization and time averaging. In some cases (bcc metals, tetrahedral covalent crystals, etc.) the obstacles are strong, the dislocation interaction is weak and the plastic flow can be viewed as a sum of uncorrelated events. The collective effects can then be neglected making classical homogenization an appropriate tool. Here we are interested in an alternative case (fcc metals, hcp crystals with basal glide, etc.) when hardening can be neglected, dislocational mobility is high and the elastic interaction among distant dislocations is significant. Then a specific collective behavior at the macroscale emerges as a manifestation of many correlated events at the microscale Zaiser et al. (2008); Seeger (1957); Kubin et al. (2002); Ananthakrishna (2007). In particular, such plastic flows exhibit in the steady state irregular isolated bursts and reveal apparently random localized active slip volumes, with both spatial and temporal fluctuations spanning all scales. The temporal intermittency manifests itself through acoustic emission with power law statistics of avalanches Weiss and Grasso (1997); Weiss et al. (2000, 2001); Richeton et al. (2006) while spatial fluctuations take the form of fractal dislocational cell structures Hähner et al. (1998); Zaiser and Hähner (1999); Mughrabi et al. (1986).

The ensuing macroscopic dynamics is usually characterized as critical because of the reasons that will become more clear in what follows. In physical terms the resulting behavior can be described as glassy. Since dislocational motion proceeds in a self similar manner in both space and time the importance of such notions as average dislocation density and average velocity in such regimes become questionable. Those notions, however, remain relevant in non-critical plastic regimes where dislocational microstructures are regular and are characterized by particular scales which depend on the level of loading, e.g. Kubin. (1993).

Despite the over-schematic description of dislocational activity in classical engineering theory, it has been remarkably successful in capturing the most important plasticity phenomenology such as yield stress, work hardening and shakedown Hill (1998); Rice (1975). However, the fractal patterning and the peculiar scale free structure of the temporal fluctuations remain invisible at the macroscale and are missed by the approaches based on ’intuitive’ spatial and temporal homogenization. As a result, a quantitative link between the fitting functions used in phenomenological plasticity and the microscopic picture of a crystalline lattice with moving line defects remains elusive. This gap has been a major obstacle on the way of quantitative estimate of plastic response for artificially designed materials.

The emergence of power laws suggests that the relation between microscopic and macroscopic pictures of plastic flow is rather complex and more akin to turbulence Cottrell (2002); Kubin et al. (2002); Choi et al. (2011) than to elasticity where classical homogenization of lattice models is usually sufficient to predict macroscopic response Born et al. (1998); Blanc et al. (2006). It is also clear that powerful methods of equilibrium statistical mechanics implying ergodicity and absence of correlations are not applicable for the description of small scale plasticity which appear to be highly nonequilibrium and strongly correlated phenomenon. The main problem with applying the classical continuum approach is that spatial renormalization symmetry (self-similarity) excludes separation of scales while non gaussian statistics prevents local description. More generally, we still lack analytic tools allowing quantitative description of evolutionary phenomena that depend on rare events spanning the whole domain and have convergence problems handling the uncertainty which grows algebraically with time which qualifies such systems as located at the ’edge of chaos’ Bak et al. (1991).

The fine structure of long range correlations in plasticity has been extensively studied in recent experiments which unambiguously established the power law statistics of avalanches and produced numerous examples of the scale-free nature of dislocation patterns Weiss and Grasso (1997); Weiss et al. (2000, 2001); Richeton et al. (2006); Dimiduk et al. (2006); Richeton et al. (2005a); Bharathi et al. (2001); Bharathi and Ananthakrishna (2002); Zaiser et al. (2008); Richeton et al. (2005b); Weiss et al. (2007); Ng and Ngan (2007); Brinckmann et al. (2008). Similar phenomena of self organization towards criticality have been also observed in other mechanical systems operating in nonequilibrium steady regimes including friction, fracture, porous and granular flows and martensitic transformations Turcotte (1999); Sornette (2000); Jensen (1998); Vives (1994, 1998).

While the general mechanisms for the formation of scale free correlations remain a major subject of research Henkel et al. (2008), several important ingredients of dynamic criticality, such as marginal stability, quasi-static driving, threshold type nonlinearity and nonlocal feedback, have been identified Papanikolaou et al. (2011); Dickman et al. (2000). All these components are present in crystal plasticity. For instance, marginal stability reveals itself through the fact that infinitesimal perturbations not only trigger small scale rearrangements of the ensemble of locked dislocations but can also induce a global transition between the jammed and the flowing regimes. The long range interactions responsible for the feedback have elastic nature, suggesting that, somewhat paradoxically, elasticity plays a fundamental role even in ’ideal’ plasticity.

The experimental findings of plastic criticality in crystals have been supported by a number of numerical simulations reviewed in Zaiser (2006); Weiss (2003). The two main approaches are: the discrete dislocation dynamics (DDD) accounting for dislocation interaction on different slip planes Miguel et al. (2001, 2002, 2001); Csikor et al. (2007); Beato et al. (2011) and the pinning-depinning models considering dislocations dynamics on a single slip plane Moretti et al. (2004); Koslowski et al. (2004a, b). Discrete mesoscopic models of crystal and amorphous plasticity implying some averaging and providing an effective description of the pinning-depinning and jamming-unjamming processes have also been shown to generate power law statistics of avalanches with realistic exponents Michael and Aifantis (2006); Zaiser and Moretti (2005); Dahmen et al. (2009); Baret et al. (2002); Chan et al. (2010); Talamali et al. (2011). Among those we particular emphasize the models of earthquakes aimed at reproducing the Gutenber-Richter law because the rapture-healing phenomena on a preexisting fault are very similar to 2D plasticity showing almost the same critical exponents Chen et al. (1991); Ben-Zion and Rice (1993); Fisher et al. (1997). In view of the implied analogy with turbulence, of considerable interest are also the models dealing with dynamics of continuously distributed dislocations and capturing the scale free effects already in a PDE setting Acharya (2001); Roy and Acharya (2005); Limkumnerd and Sethna (2008, 2006); Fressengeas et al. (2009); Choi et al. (2011); Chen et al. (2010).

Despite their ability to reproduce experimental findings, the above models cannot provide theoretical insight into the origin of criticality because they lack analytical transparency and depend on large scale numerical experiments. On the other hand, the physical adequacy of the proposed numerical schemes may still be disputed in some basic details. For instance, in case of DDD models, the fast microscopic topological changes (nucleation, annihilation, interaction with obstacles, etc) are treated by auxiliary hypotheses formulated phenomenologically in terms of local stresses without addressing the mechanism of barrier crossing Maloney and Lemaître (2004); Tanguy et al. (2006). Likewise, the theory of continuously distributed dislocations have difficulty representing kinks, locks, jogs and other individual entanglements with nontrivial topology. The numerical models can be made more adequate by involving various quasi-continuum schemes, however, it is clear that in order to understand plastic criticality mathematically, the existing models have to be simplified rather than further complicated.

The anticipation of a simple and transparent theory is based on the idea that truly scale free behavior is independent of either microscopic or macroscopic details of the system. The statistics of dislocation avalanches should then be an intrinsic feature of a particular crystal class described in terms of symmetry and dimensionality and not affected by specific material parameters and details of the loading geometry. It is then of interest to formulate simple prototypical models representing particular universality classes of plastic flows. The simplicity is understood here in the sense that the model is amenable to rigorous mathematical study while still capturing not only the observed exponents but also the characteristic shape functions in scaling relations Papanikolaou et al. (2011); Chen et al. (2011a).

At present, the only analytically tractable approaches to criticality in driven distributed systems are the models based on branching processes Sornette (2000), the mean field model Dhar and Majumdar (1990), the elastic depinning model amenable to FRG methods Le Doussal and Wiese (2009) and the Abelian sand pile type automata with integer valued fields Dhar (1999a); Redig (2005). The challenge is to reduce a realistic model of plasticity to one of these analytically transparent types.

Since the relation between the prototypical branching processes and the dynamical phenomena in distributed systems remains rather imprecise, we will not pursue this line of modeling here. In contrast, the observed critical exponents for crystal plasticity are so close to the mean field values that claims have been made that 3D plasticity is in the mean field universality class. This hypothesis is supported by the estimates of the critical dimension in various systems with long range elastic interactions Papanikolaou et al. (2011); Zapperi et al. (1998); Colaiori (2008); Dahmen et al. (2011). However, the question of validity of this approximation remains open because the match of the measured exponents is not perfect and the observed avalanche shapes appear to be less symmetric than the theoretical predictions Sethna et al. (2001). Also, the mean field values of exponents have not been supported by numerical experiments with colloidal crystals Salman and Truskinovsky (2011) and within phase field crystal model Chan et al. (2007).

The depinning models have a considerable relevance for plasticity because they study quasi-statically driven elastic objects interacting with a set of randomly distributed obstacles. The corresponding theoretical critical exponents for the case of dislocations have been computed with high precision Fisher (1998); Kardar (1998); Zaiser (2006); Pierre and Doussal (2010) and confirmed by molecular dynamics studies Patinet et al. (2011). These models, however, may not be in the same universality class as crystal plasticity because the critical exponents are different from those observed in plasticity experiments. Most probably, the origin of divergence is the neglect in depinning models of important interactions among dislocations on different slip planes.

Left with the only remaining option, we study in this paper the possibility to reduce the conventional dislocation mediated plasticity model to an integer valued Abelian automaton. The Abelian symmetry means that the order, in which the eligible ’spin variables’ are updated during the avalanche, is irrelevant for the avalanche outcome. Some Abelian automata are amenable to rigorous analysis due to partial ergodicity on a set of recurrent sequences which were fully characterized in several cases Dhar (1999a); Redig (2005); Jarai (2005); Dhar (1999b); Fey et al. (2010). The aim of this paper is to show that a reduction of a multi-slip-plane dislocational dynamics to an integer automaton with some form of Abelian symmetry allows one to reproduce the experimentally observed statistics of avalanches. An analytical study of the ensuing automaton will be reported elsewhere.

## 2 Summary of the main results

In search for the simplest representation of plasticity universality class(es) we begin with a series of one dimensional models allowing one to reproduce rate independent dissipation and hysteresis. Then we move to two dimensional models capturing also intermittency and power law statistics.

Our starting point is an assumption that the micro-scale dynamics is overdamped and that the rate of loading is much slower than the rate of viscous relaxation but much faster than the rate of thermal relaxation. Such models belong to the class of AQS (athermal quasi-static) systems that have been found relevant for the description of many rate independent phenomena from wetting to magnetism Maeda and Takeuchi (1978); Kobayashi et al. (1980); Maeda and Takeuchi (1981); Malandro and Lacks (1998); Yamamoto and Onuki (1998); Vives (2001).

The simplest 1D lattice model of this type, representing a chain of interacting particles (or a deck of interacting rigid cards), can be developed for transformational plasticity of shape memory alloys. To obtain in this setting a plastic response at the macro-level we must assume that the interaction potential has a double well structure which fully characterizes the local ’mechanism’ of the transformation. The presence of local bi-stability, which is also relevant for the description of the ’easy glide’ stage in crystal plasticity, places such mechanical system into the class of snap-spring lattices. The non-convexity of local interactions may be equally relevant for amorphous plasticity operating with a concept of STZ (shear transformation zones), which implies a multi-well structure of an effective potential Falk and Langer (2011); Maloney and Lemaitre (2006). We emphasize that the discreteness at this stage has to be understand broadly as describing a generic structural inhomogeneity, for instance, a grain structure.

In the simplest chain with nearest neighbor (NN) interactions the elastic elements are subjected to common stress field. The resulting mean field model captures the basic mechanism of rate independent dissipation and can be used to illustrate the idea of marginal stability of the underlying elastic system which is ultimately responsible for the plastic yield Puglisi and Truskinovsky (2005); Mielke and Truskinovsky (2011). A more realistic model must incorporate short range interactions and the simplest extension in this direction is a chain with next to nearest neighbors (NNN). The latter may be of ferromagnetic or anti-ferromagnetic type Truskinovsky and Vainchtein (2004a); Braides and Truskinovsky (2008) and in both cases the resulting mechanical system can be viewed as a soft spin generalization of the classical Ising model.

In the case of ferromagnetic interactions penalizing gradients we show that that the NNN model captures the important difference between the nucleation and the propagation thresholds and describes the characteristic nucleation peak Onuki (2003). In the case of anti-ferromagnetic interactions, favoring lattice scale oscillations, such model can reproduce periodically distributed shear bands. We find, however, that even in the presence of gaussian quenched disorder, the internal dynamics in all these 1D models is still too simple and none of them can generate scale free correlations in the overdamped setting. The situation does not improve if the double-well NN potential is replaced by a periodic potential.

Despite their failure to capture criticality, the 1D models are very important because they allow one to show that in the limit of quasi-static driving a viscous evolution in a rugged energy landscape takes the form of a stick slip dynamics adequately described by a discrete automaton. The replacement of the fast stages of dynamics by jumps leads to considerable simplification of the computational problem. It also shows that in quasi-static setting all dissipative mechanisms producing the same set of jumps are equivalent.

To simplify the system even further we replace a smooth multi-well potential by its piece-wise parabolic analog allowing one to solve the elastic problem analytically when the ’phase configuration’ is known. Since the phase configuration is prescribed by an integer valued field, such ’condensation’ of elastic problem allows one to reduce the original smooth dynamical system to an integer valued automaton of a sand pile (threshold) type. Because of the long range nature of elastic interactions, the discharge rules in this automaton are nonlocal which makes its rigorous analysis challenging.

While minimizing out the elastic fields in plasticity problems is a rather common approach leading to a reduced description in terms of either plastic distortion or dislocation density Ortiz et al. (2000); Koslowski et al. (2002); Berdichevskyître (2006); Limkumnerd and Sethna (2008); Choi et al. (2011); Chen et al. (2011b), the proposed integer automaton representation appear to be new. Its salient feature is that in place of straightforward time discretization of continuum dynamics Narayan and Fisher (1993) we utilize inherent temporal discreteness hidden behind the conventional gradient flow dynamics (see also Pérez-Reche et al. (2007); Puglisi and Truskinovsky (2005); Mielke and Truskinovsky (2011)).

The next level of complexity is exemplified by 2D lattice models. Here we abandon the double-well framework with its prototypical defects and consider the simplest setting allowing one to represent actual dislocations which interact through realistic long range elastic fields and are capable of generating collective effects. As it has been repeatedly proposed in the literature, the two dimensional models of plasticity may be sufficient to capture the corresponding universality class at least for FCC and hexagonal crystals Miguel et al. (2001); Zaiser and Moretti (2005); Ispánovity et al. (2010). For similar models of criticality associated with transformational plasticity, see Pérez-Reche et al. (2007, 2008, 2009)).

In general terms, our 2D model represents a cross section of a crystal with a single slip geometry. More precisely, we consider a set of parallel edge dislocations which can move only in horizontal direction. This implies that the crystal is highly anisotropic with the deformation in one direction strongly constrained. As a partial justification we mention that intermittency has been mostly observed under single slip conditions Zaiser and Moretti (2005); Ispánovity et al. (2010).

The ensuing discrete model with continuous dynamics can be viewed as a parallel set of coupled overdamped FK chains Miyazaki et al. (1982); Kovalev et al. (1993); Landau et al. (1993a); Landau (1994); Carpio and Bonilla (2003). It has been shown before that this model allows for generation and annihilation of dislocations and that it describes adequately their long range elastic interactions Carpio and Bonilla (2005); Plans et al. (2008). Most importantly, both the propagation and the nucleation/anihillation mechanisms in this model are associated with reaching the same strain thresholds. In particular, this model correctly accounts for Peierls stress associated with underlying periodicity, describes adequately sufficiently slow dislocational kinetics and does not require special rules governing the interaction of the dislocations with quenched disorder Carpio and Bonilla (2003, 2005). At the same time, this model still neglects tensorial nature of elasticity and ignores such important physical effects as cross-slip and climbing which may affect critical exponents.

In our numerical simulations we assume periodic boundary conditions and submit the 2D lattice to quasi-static cycling in the hard device by imposing periodic shear deformation. We show that after several loading-unloading cycles the system reaches a non-equilibrium steady state (plastic shakedown) where it exhibits critical behavior characterized by the power law statistics of avalanches. The associated spatial fractality is revealed through the study of spatial correlations. A study of the power spectrum of the energy fluctuations shows the characteristic noise which is another signature of scale free behavior. Overall, the critical manifestations of the plastic flow revealed by our simplified modeling are in full agreement with experimental observations and with previous more elaborate numerical studies in 3D.

In an attempt to understand the origin of criticality in analytical terms we performed a reduction of the system of ordinary differential equations describing our 2D visco-elastic lattice model to a threshold automaton with conservative dynamics. As in 1D case, we minimized out elastic fields and obtained an automaton model for a discrete variable which serves as an indicator of the number of dislocations passing a given point. Despite partial linearization implied by the use of piece-wise quadratic potential and the replacement of time-continuous dynamics by a sequence of jumps, the automaton model exhibits the same critical behavior as the original ODE model. Our numerical study suggests that the discrete automaton behind the FK model has a statistical (weak) Abelian property which makes it amenable in principle to rigorous mathematical treatment despite the nonlocal nature of the discharge rules.

A short announcement of some of the results of this paper can be found in Salman and Truskinovsky (2011).

## 3 The model

Classical phenomenological theory of plasticity is applicable in both crystal and amorphous settings because it circumvents an explicit reference to the nature of the defects carrying inelastic deformation. Instead, it operates with a concept of an averaged non affine deformation represented locally by a tensorial measure of plastic strain. The evolution of the plastic strain is governed by a dissipative potential which is assumed to be a homogeneous function of degree one making the dissipation rate independent. Other internal variables characterizing hardening, softening, additional inelastic rearrangements, etc. are often added as well, each equipped with its own dynamics. For instance, in crystal plasticity the effective plastic strain rate is parameterized by slip measures associated with different slip directions and evolved by coupled rate independent mechanisms. Likewise, in transformational plasticity of shape memory alloys the evolution of the phase/twin fractions is usually assumed to be governed by (coupled) rate independent dynamic equations.

Despite their success in fitting the observed memory behavior, the phenomenological models of plasticity do not describe adequately the internal dynamics revealed by endogenous fluctuations because they do not deal explicitly with microscopic time and length scales. In particular they miss temporal intermittency and spatial fractality.

An alternative program of describing plasticity directly in terms of defects responsible for the flow is far from being complete despite many important recent successes Motz et al. (2009); Csikor et al. (2007); Zhou et al. (2010); Chan et al. (2010); Needleman et al. (2006); Limkumnerd et al. (2008); Shishvan et al. (2011). First of all, in order to handle microscopic dynamics of dislocations one needs to go back to the elastic framework. Since the strain field of a dislocation decays rapidly away from the core region, the corresponding long-range distortions can be adequately described within the framework of linear elasticity. Therefore dislocations are modeled as (line) singularities of linear displacement fields bounding surfaces that carry the slip. The latter is represented in the elastic models as a quantized displacement discontinuity and this is the only place where the inherent discreteness enters the picture.

The localization of dislocations in space and time requires constitutive assumptions which are outside the realm of linear elasticity. Such assumptions are usually presented in the form of additional kinetic relations of phenomenological nature describing the associated mesoscopic, rate dependent dissipation O. and L. (2003). An additional problem in this approach is that independent phenomenological assumptions are also required to describe nucleation, annihilation and other topological transitions. Notwithstanding all this uncertainty, the simulation of realistic 3D flows requires tracking of a huge number of interacting defects which remains prohibitively expensive in terms of computational time.

In search for an alternative description we first recall that a simple hybrid discrete-continuous dislocational model can be formulated by using the Peierls-Nabarro framework Hirth and Lothe (1982); Nabarro (1951); Koslowski et al. (2002). Suppose that the bulk elastic energy is defined by

where is the displacement field and the function is quadratic. The field is allowed to have discontinuities which are penalized by surface energy of the form

In Peierls-Nabarro model the function is assumed to be periodic and the evolution of the field is governed by an appropriate gradient flow dynamics. The main technical difficulty in this approach is to track the a priori unknown surfaces of displacement discontinuity with arbitrary complex geometry.

In order to avoid the explicit ’tracking’ of the slip surfaces one can consider their ’capturing’ by incorporating discreteness already into the bulk energy term. To this end one can introduce an internal length scale and consider a mesoscopic energy which we schematically represent as

(1) |

Now is the lattice field and the function is periodic with infinite symmetry group Conti and Zanzotto (2004); Pérez-Reche et al. (2009) . The discrete energy contains information about both the bulk energy and the surface energy ; the mechanism of such splitting in the continuum limit was discussed in the context of fracture in Truskinovsky (1996) and in the context of plasticity in Pellegrini et al. (2010). It is important to keep in mind that the discrete energy describes multi-particle elastic units and is therefore different from the conventional molecular dynamics potentials dealing directly with interatomic interactions.

After the energy density in (1) is specified, one can formulate the rate dependent dynamics for the elastic fields carrying dislocations. Under the assumption that such dynamics is viscous and overdamped we obtain the following system of ODEs

(2) |

where is the ratio of the internal time scale and the time scale of the driving Pérez-Reche et al. (2007); Puglisi and Truskinovsky (2005). Notice that instead of classical visco-elasticity we consider here an ’environmental’ viscosity which allows us to avoid certain dynamic degeneracies in the system without quenched disorder Mielke and Truskinovsky (2011). Notice also that the system (2) is not autonomous because time enters through the equations for boundary units as, for instance, in the case of cyclic driving in a hard device.

The numerical solution of (2) may be rather challenging under generic loading. However, in the limit of quasistatic driving which can also be interpreted as the limit of zero viscosity, the dynamic system (2) can be replaced almost everywhere (in time) by the equilibrium system Puglisi and Truskinovsky (2005)

(3) |

The dynamics is then projected onto the branches of metastable equilibria which are defined as sequences of local minima (of the elastic energy) parameterized by the loading parameter. The crucial observation is that similar to the case of relaxation oscillations, the dissipation-free phases of such dynamics are necessarily interrupted by the branch switching events which are instantaneous at the time scale of the quasi-static driving. Such reduction of an original dynamic problem opens the way towards minimizing out the elastic variables and reformulating the problem in terms of a discrete integer valued variable labeling individual metastable branches Pérez-Reche et al. (2007); Puglisi and Truskinovsky (2005); Mielke and Truskinovsky (2011).

In the next section we illustrate this general approach by using the simplest one dimensional setting where our goal will be modest: to make the idea of marginal stability fully transparent and to obtain the simplest examples of fluctuating plastic flow. Then we formally expand the model to the second dimension and show that the ensuing minimal geometrical complexity is already sufficient to generate realistic fluctuational behavior.

## 4 One dimensional setting

The very first model of plasticity accounting for fluctuations was probably Prandtl’s zero-dimensional model describing a particle driven through a spring in a periodic landscape Prandtl (1928). This model involves marginal stability, captures hysteresis, and reproduces rate independent dissipation Puglisi and Truskinovsky (2005) but misses the crucial interactions among the defects. To obtain criticality in this model one needs to considerably complicate the effective landscape, replacing, for instance, a periodic potential by a Brownian motion type potential Colaiori (2008). Such fine ’tuning’ of the interaction potential places this model in the class of phenomenological models representing a natural development of the classical rheological models used in engineering plasticity.

The simplest nontrivial ’first principle type’ model with the energy (1) describes the transformational plasticity of shape memory alloys Mueller and Villaggio (1977); Fedelich and Zanzotto (1992); Puglisi and Truskinovsky (2000, 2002, 2005). In this model each elastic element interacts with its nearest neighbors (NN) only through the total stress which then plays the role of a mean field. Having the advantage of utmost simplicity, this model possesses a nonphysical permutational invariance which introduces undesirable degeneracy. Therefore we subsequently augment this model by introducing linear interactions between next to nearest neighbors (NNN) and allowing them to be either attractive and repulsive Truskinovsky and Vainchtein (2004a); Braides and Truskinovsky (2008). Finally, we replace a bi-stable NN potential by a periodic one and introduce quenched disorder.

### 4.1 NN model

Consider a chain of particles forming in the reference configuration a regular 1D lattice. Denote the reference particle positions by where we imply that the atomic spacing is taken as the length scale (). Introduce discrete displacement field and the discrete strain field Our Fig.1 shows the simplest representation of this mechanical system. In plasticity framework one should rather view the chain as a deck of rigid cards interacting through nonlinear shear springs and view the displacement field as being perpendicular to the direction of the coordinate axis.

For the NN model the equations of motion (2) take the form

(4) |

where is the elastic force resulting from stretching of a spring. We assume that the system is placed in a hard device with

where is the dimensionless time playing the role of the loading parameter.

One can see that the total energy in the NN model is the sum of the energies of individual springs. This means that the springs interact only through the constraint imposed by the hard device. Indeed, in equilibrium where is the stress common to all springs; such infinite range of interactions is characteristic of paramagnetism.

To facilitate the subsequent elimination of elastic variables we assume that the double well potential describing individual springs is piece-wise quadratic:

(5) |

Here is the elastic modulus which we assume to be the same in both phases. The interval will corresponds to Phase 1 with equilibrium strain and the interval , to Phase 2 with equilibrium strain .

Our Fig. 2 shows the macroscopic strain-stress curve obtained from the numerical solution of (4) with sufficiently small. As we load the system from a homogeneous state in Phase 1, the phase transition starts at the critical strain , where only one spring changes phase. The transformation in one spring corresponds to a stress drop as the other springs relax. Each subsequent stress drop is also associated with only one spring changing phase during each jump. One can see that the resulting macroscopic stress strain curve forms a wiggly plateau with a rudimentary intermittency manifesting itself through equally spaced oscillations. This behavior was studied in Puglisi and Truskinovsky (2002, 2000) where it was shown that the individual springs can change phase only consequently.

As we mentioned in the Introduction, our model can be simplified further in the case of quasi-static loading which for the the dynamical system (4) corresponds to a zero viscosity limit . We begin with rewriting the piece-wise quadratic potential (5) as

(6) |

Here, we introduced a double-valued spin variable describing the phase state of a given spring with for and for . At a given spin field, the displacement field must satisfy the following set of linear equations

(7) |

where the right hand side plays the role of ’dislocation density’Berdichevskyître (2006). If the phase field is known, the tri-diagonal matrix on the left may be inverted and the strain field can be obtained explicitly. Moreover, due to permutational invariance of the problem, one can always assume that there is only one defect and therefore the equilibrium branches can be parameterized by a single variable which changes in the interval and describes the ’phase fraction’ in the mixture (see Puglisi and Truskinovsky (2000)).

Several metastable branches parameterized by are shown schematically in Fig. 3. The first branch with all springs in the first energy well (Phase 1) corresponds to . It ends at point where this homogeneous state becomes unstable. Suppose that we drive the system by quasi-statically changing the total strain. Then at point one spring changes phase and the system undergoes a dynamic transition to the new local minimum identified in Fig. 3 by point . This fast process is accompanied by a jump in stress; as a consequence of this jump a finite amount of energy dissipates into heat.

From point the driven chain evolves along the second equilibrium branch corresponding to which in turn becomes unstable at point . Here again one spring switches phase, the stress drops and the system finds itself in a new local minimum identified as point on the branch with . The process continues till all springs change phase and the system reaches the second homogeneous branch corresponding to . Interestingly, the way we introduced viscosity ensures that the transformation proceeds as a single front which is not the case if we directly deal with the degenerate problem at or if we consider classical visco-elasticity Mielke and Truskinovsky (2011).

A full analysis of the thermodynamics of the NN model including the computation of the ’heat to work ratio’ can be found in Puglisi and Truskinovsky (2005).

### 4.2 Automaton model

After the elastic fields are minimized out (either by the inversion of the matrix in (7) or by the Fourier transform as we show below), the ensuing problem in terms of the variables can be reformulated as a spin-valued automaton.

To be more precise suppose again that initially all the units are in the Phase 1, meaning that the phase variable vanishes everywhere. Since the field is known, one can compute by solving (7). Then one can find and locate the system on the equilibrium branch with . The stability condition in this case reduces to an inequality where is the upper spinodal stress Puglisi and Truskinovsky (2005). Suppose that the stability condition is violated at time . Then we need to update the phase configuration which leads to a stress drop from to . Notice that the choice of the flipping element remains arbitrary in the automaton why in the ODE framework it is uniquely prescribed by the dynamics. To overcome this (artificial) non-uniqueness in the automaton setting we later consider two regularizing mechanisms: NNN interactions and/or quenched disorder.

Along the new metastable branch one can continue computing till the spinodal stress is reached again and then it is necessary to make the next update. The stress oscillations continue till we reach the homogeneous Phase 2 with . During unloading the stability condition changes into where is the lower spinodal stress and the update rule becomes . The homogeneous branch corresponding to Phase 1 is reached again when . We do not present here the strain-stress curve obtained from the automaton-based simulations because it is practically indistinguishable from the results obtained by solving directly the dynamic equations (4) with shown in Fig. 2.

In order to show that on the yielding plateau the system evolves through a sequence of marginally stable states, we plot in Fig. 2 the energy of the system against the applied strain for one loading-unloading cycle. When the system is in one of the pure phases (affine or Cauchy-Born state) it resides in a global minimum of the energy until the Maxwell threshold is crossed (during either loading or unloading). Then the global minimum becomes a local minimum and eventually the homogeneous state becomes unstable as the system reaches either upper or lower spinodal stress. After that the system evolves through the sequence of non-affine states where different springs occupy different phases. Notice, however, that in this non-affine stage the stress in the chain never deviates considerably from one of the spinodal stresses (see points , , , , etc. in Fig. 3). It is not hard to see that in the continuum/thermodynamic limit the jumps become infinitesimal and the ensuing yielding plateau collapses the spinodal Puglisi and Truskinovsky (2005). Therefore, in the continumm limit the yielding system is marginally stable on the whole yielding plateau.

In summary, plasticity in this simplest model is associated with elastic states of minimal stability. This peculiar feature of plastic flow may be considered as a first manifestation of criticality. In addition, the model exhibits stress drops that can be interpreted as avalanches. The problem is that these avalanches are all of the same size and the corresponding microscopic transformation events are either un-correlated (zero viscosity limit) or over-correlated (finite viscosity case).

### 4.3 NNN model with ’ferromagnetic’ interactions

As we have seen, the NN model lacks the short range interactions capable of generating extended dislocation cores. To bring the missing internal length scale into the theory one can add to the NN model an additional ’ferromagnetic’ interactions favoring local homogeneity of the strain field Rogers and Truskinovsky (1997).

The simplest way to introduce short range interactions is to consider an NNN model with the energy density Truskinovsky and Vainchtein (2004a)

(8) |

Here is the spin variable introduced in the previous subsection and is the strength of short range interactions (we consider only the automaton version of the model). The ’ferromagnetic’ interactions ensuring the localization of phase boundaries corresponds to ; the case of ’anti-ferromagnetic’ interactions , favoring the formation of fine mixtures, will be considered in the next subsection (see also Truskinovsky and Vainchtein (2005)).

In the quasi-static limit discussed above, the equilibrium equations for the NNN model can be written in the form

(9) |

where again the ’dislocation density’ in the right hand side plays the role of a source of elastic strain. Due to non locality of the NNN model, we need to complement the hard device boundary conditions used in the NN model by two additional conditions. One way to avoid surface boundary layers is to add two fictitious springs just outside the boundaries of the chain ensuring the periodicity conditions and . This leads to additional boundary terms in the energy, however, the elastic problem can be again solved explicitly Truskinovsky and Vainchtein (2004a). The remaining problem for the spin field can be formulated as an automaton which is similar to the one described in the previous section and we omit the details.

The macroscopic strain-stress curves for the NNN model with are shown in Fig.4. In addition to the horizontal plateaus, that again correspond to the regimes when an isolated phase boundary sweeps through the chain (propagation), one can also see two peaks at the onset of the direct and inverse transformation (nucleation) and two smaller anti-peak describing the disappearance the existing phase boundaries (annihilation). In the computational experiment shown in Fig.4 four springs out of 24 participate in the nucleation event while the propagation involves successive jumps of individual springs. The necessity of the peak in the NNN model and the nature of the nucleation event in the continuum limit were discussed in Truskinovsky and Vainchtein (2004a).

In Fig. 5) we show the macroscopic energy-strain relation for this model (loading only). In contrast to what we have seen in the NN model, here the energy increases in the systematic manner as the system evolves along the plastic plateau. This means that plastic dissipation is accompanied by an additional energy build up which can be vaguely qualified as the ’cold work’ Puglisi and Truskinovsky (2005). The origin of this effect in our model is the particular form of the NNN term in the energy (8) which, in addition to penalizing gradients, also contributes to homogeneous elasticity. Interestingly, similar continuing growth of the energy along the flat-stress yielding plateau have been observed in several microscopic (MD-type) models of amorphous plasticity Procaccia (2009).

In order to avoid such hidden contributions to bulk elasticity, one can introduce short range interactions in a different form. The simplest way to keep the ferromagnetic coupling intact while removing the NNN contribution to homogeneous elasticity is to introduce a three-particle interaction mimicking the strain gradient (SG) term in the discrete setting

(10) |

It is not hard to see that in the case this energy is equivalent to the NNN energy (8) with up to a ’local’ term depending on only. This means that the ’cold work’ is eliminated in the SG formulation which is confirmed by our Fig. 5(b). Notice that outside the nucleation (annihilation) events the behavior of the NN and SG models are quite similar with isolated phase boundaries propagating along the chain as the total strain increases. The main difference is that in the NN model a phase boundary (mimicking dislocation) is atomically sharp, while in the SG model its core has a finite size.

Despite all these new features, neither NNN nor SG model are successful in predicting the realistic statistical structure of avalanches. Indeed, both models exhibit one correlated ’snap’ event (nucleation) and a succession of equally sized ’pop’ events (propagation). It is clear that this fluctuation/avalanche structure is still too simple to be compared with realistic experiments.

### 4.4 NNN model with ’anti-ferromagnetic’ interactions

As it is well known, the nonlocal kernel describing linear elastic interactions in dimensions 2 and 3 has both ferromagnetic and anti-ferromagnetic contributions Kartha et al. (1995); Ren and Truskinovsky (2000); Lookman et al. (2003). Therefore it is of interest to see how plastic behavior in our toy model changes if ’ferromagnetic’ NNN interactions with are replaced by ’anti-ferromagnetic’ NNN interactions with . In the 1D setting, anti-ferromagnetic interactions favoring formation of small scale microstructures, should be understood as a poor man’s attempt to capture the gradient nature of the strain (elastic compatibility) and the effect of the constraints placed by higher dimensional boundary conditions Ren and Truskinovsky (2000).

Our numerical results for the anti-ferromagnetic NNN model driven in a hard device (automaton version) are summarized in Fig. 6. The most important feature is that in addition to two homogeneous phases the model exhibits the third phase representing binary lattice scale mixture which is homogeneous only in average (see Fig. 7). This new phase shows binary oscillations with neighboring springs occupying different phases; the energetic preference of such non Cauchy Born (non-affine) phase arrangements follows from the studies of the global minimum of the elastic energy, e.g. Cahn and Vleck (1999); Pagano and Paroni (2003); Braides and Gelli (2006); Braides and Truskinovsky (2008); Charlotte and Truskinovsky (2002); Truskinovsky and Vainchtein (2005).

The salient feature of the macroscopic response in the anti-ferromagnetic NNN model is the disappearance of the nucleation peak. One can see that on the first half-plateau of the strain-stress curve (Fig. 6), the homogeneous phase is progressively replaced by the growing ’binary phase’ and the corresponding front is shown in Fig. 7. Once again, each drop in stress profile corresponds to one spring changing phase. Notice that the thresholds for the nucleation of the binary phase and for its propagation are similar. By the end of the first half-plateau the binary phase occupies the whole domain and the step in the middle of the strain-stress curve (Fig. 7) describes purely elastic deformation of the binary mixture. Along the second half-plateau the mixture phase is gradually replaced by the second phase and we again observe a propagating front. The dynamics of such fronts has been studied in Cahn et al. (1999); Vainchtein and Van Vleck (2009).

From Fig. 6 we see that ’ferromagnetic’ model, the bulk contribution to the energy in the ’anti-ferromagnetic’ model changes its sign. The main plastic phenomenology, however, remains the same even though the microscopic spatial arrangement of phases changes from phase separation to fine scale mixing.

For completeness, we also present here the results of the simulations for the anti-ferromagnetic SG model with , see Fig. 8. The strain profiles in this model are basically the same as in the NNN model with anti-ferromagnetic interactions (), however, the energy profile shows a new feature: a negative slope in the first segment and a positive slope in the second. This implies that the negative contribution to the cold work along the first half-plateau is compensated by the equal in magnitude positive contribution on the second half-plateau.

To summarize, the 1D anti-ferromagnetic NNN model and its SG analog describe rate independent dissipation and show a new feature which can be interpreted as a propagation of regular slip band microstructures. Still, both models predict highly coherent fluctuations structure and therefore none of them captures the complexity of the experimentally observed spatial and temporal correlations.

### 4.5 Quenched disorder

In modeling of realistic material response one cannot neglect disorder due to various chemical and mechanical imperfections at the level of a lattice. Such disorder can influence the macroscopic behavior of the system by facilitating inhomogeneous nucleation of defects on one side and by creating a variety of pinning obstacles for the propagating defects and defect microstructures, on the other side. In this subsection we show how the quenched disorder influences the macroscopic response in the basic NNN model with ferromagnetic interactions.

A disorder can be incorporated into the energy (10) in different ways, for instance, by randomizing the size of the barriers separating individual energy wells Puglisi and Truskinovsky (2002). In this paper, due to technical reasons, we follow a different approach originating from the RFIM model Sethna et al. (1993) where disorder is represented by a random field biasing different phases in different spatial points. The ensuing RSSM model (random soft/snap spring model Pérez-Reche et al. (2008)) is characterized by the energy

(11) |

where and are independent random numbers distributed according to

In all our numerical illustrations we put ; we have checked that the value of dispersion does not affect our statistical results even though it affects some macroscopic features, for instance, the degree of hardening (see Mielke and Truskinovsky (2011)).

A typical strain-stress relation for a particular realization of disorder is shown in Fig. 9. We observe that the nucleation peak persists, however the successive stress jumps due to phase boundary advances have random amplitudes and each avalanche now involves different number of springs. Between the jumps the system is trapped/pinned in metastable configurations where it evolves purely elastically. The instabilities can be again associated with a succession of minimally stable states where at least one element have reached the spinodal state. Notice that the disorder in this model does not evolve and therefore the shakedown state is reached already during the first cycle.

We observe that as in the case without disorder, each jump at time is associated with a finite dissipation which can be computed as the energy difference before and after the jump

(12) |

At a given the value of may be either zero (elastic deformation) or positive (irreversible jump). In a system with disorder a single jump, originating from a phase change in one unstable element, does not necessary lead to a stable state and therefore an avalanche may contain several elementary jumps. We have already encountered this phenomenon while dealing with nucleation requiring simultaneous transition in several elastic elements. In the case of nucleation, however, we had to deal with massive transformation at one value of the loading parameter. In a disordered system typically only one element becomes unstable at a given value of the load, however, several such jumps may take place before the next elastic stage is reached.

A typical time series for the dissipated energy is shown in Fig. 10. We observe that jumps are separated by quiescent time intervals where . The shape of several typical avalanches is shown in Fig. 10 where we do not see any significant variation of time and length scales. Notice also that avalanches are very short, corresponding to 10-20 time steps. All this indicates that despite marginal stability of the system we record a succession of local events which fail to trigger global rearrangements.

The total dissipated energy during an avalanche separated by two silent intervals can be written as

where summation is over the avalanche duration. The cumulative probability distribution Clauset et al. (2009) of the variable is shown in Fig.11. One can see that the distribution is close to Gaussian if we exclude few large scale events corresponding to macroscopic nucleation/annihilation. The fact that we do not record a power law signal implies that this 1D system does not exhibit criticality and that outside highly synchronized nucleation events the avalanches remain largely uncorrelated.

### 4.6 Explicit elimination of elasticity

In the case when the elastic potential is piece-wise quadratic, one can solve the linear elastic problem analytically and present the automaton explicitly in terms of spin variables. Given our periodic boundary conditions, it is natural to invert the linear elasticity operator in the Fourier space. The discrete Fourier transform of a lattice field is defined by where is the 1D wave vector with values , . The inverse Fourier transform is given by

Consider first the NN model. By transforming (7) with added quenched disorder into Fourier space we obtain

(13) |

where

The field is the Fourier transform of the discrete spin field and the field is the transform of the random field ; both fields represent sources of elastic strains. The difference between these two sources is that the field is fixed (quenched disorder), while the field is evolving in accordance with the automaton rules which we formulate below.

First, we decompose our linear elastic field into a homogeneous part associated with prescribed time dependent affine deformation on the boundary and an inhomogeneous part which is a solution of a problem with periodic boundary conditions. More precisely, we assume that the component of the strain field is controlled by the loading while the component can be obtained by solving a periodic problem with given sources.

The expression for the non affine component of the displacement field can be obtained from (13)

(14) |

The corresponding strain field where

can be written as

(15) |

The affine component of the deformation field , is fully defined by the loading conditions

(16) |

The total strain is then equal to

The same method can be applied to the NNN model and the SG models with the replacement of by either or . We obtain

(17) |

and

(18) |

Such elimination of the elastic fields allows one to formulate the condensed automaton model for the spin variable in explicit form.

### 4.7 SG model with periodic potential

As an application of the reduction method presented in the previous section we consider here a discrete chain with periodic NN potential. Once again we circumvent dynamic simulations and move directly to the study of the automaton model corresponding to the inviscid limit .

We chose the periodic potential to be piece-wise quadratic (see Fig. LABEL:atomic111) and assume that in one period

(19) |

where is a new integer-valued spin variable which describes quantized slip; each value of defines a quadratic energy well with a period . At field given, the elastic strains must satisfy

(20) |

To regularize the NN model we add to the elastic energy the ferromagnetic SG term as in (10) and introduce a quenched disorder as in (11). The elastic field can now be decomposed into three components

(21) |

Here the affine component is given by (16); the Fourier transform of the residual strain due to quenched disorder can be written as

(22) |

where the Green’s function is

(23) |

and finally the Fourier representation of the strain field due to plastic slip is

(24) |

Now, according to (20), a metastable configuration with a given distribution is defined (is stable) within the following limits

(25) |

where

(26) |

When reaches one of the thresholds, the integer parameter must be updated

(27) |

where

(28) |

In a driven system an initially stable distribution of slip can get destabilized because the thresholds evolve with time. The update of the slip given by (28) leads to a redistribution of the elastic strains . In Fourier space the corresponding ’discharge rules’ take the form

In Fig. 13 we show the physical image of the discrete kernel for (ferromagnetic SG model) and (anti-ferromagnetic SG model). One can see that if an elementary positive slip takes place in an element, the elastic strain is increased in several neighboring elements and this ’influence domain’ increases with absolute value of . In NN model where the strain of only the unstable unit increases while the strain of all another units decreases which is a typical mean field behavior. In the regularized model with the local response is more complex as one can see in Fig. 13.

One of the crucial points in the automaton formulation is the choice of the update strategy when several elements become simultaneously unstable. To avoid this problem we introduce slip dependent disorder which signifies that random bias fields may be different in different energy wells. Physically this means that imperfections may affect unevenly the repeated slips along the same slip direction. Similar assumption is usually made in the meso-scopic models of friction/amorphous plasticity where new random thresholds are assigned to the newly formed bonds, e.g. Bak et al. (1991).

In Figs. 14 and 14 we show the strain-stress relations in the periodic problem with (SG ferromagnetic model) and = (NN model). In both cases, the macroscopic response is very similar: deformation is irreversible, springs always remain close to the marginal stability limit and the transformation proceeds intermittently through random avalanches. However, the strain profiles in the two models are different as shown in the inserts in Figs. 14 and 14. Thus, in the local model with , the strain field is chaotic at the lattice scale due to the presence of disorder. Instead, in the nonlocal model which penalizes interfaces we see a more orderly (but still random) arrangement of the domains with different level of slip. In both cases the statistical distribution of dissipated energy during avalanches remains Gaussian as in the case of transformational plasticity (a model with a double well potential).

It is also of interest to analyze the power spectrum of the signal Jensen (1998); Laurson et al. (2005)

(29) |

where is the total duration of the time series, is the frequency. The log-log plot of the power spectrum for the SG model with is shown in Fig. 15. We see the signature of a typical random process when the energy distribution is constant over a broad range in the low frequency domain. The rapidly decaying part of the spectrum at high frequencies is just a sign that the exerted power is finite. This behavior suggests a short range of temporal correlations which is incompatible with criticality. We recorded similar structure of the power spectrum in all other 1D models, with and without short range interactions.

## 5 Two dimensional setting

Even though our driven 1D models exhibited marginal stability, we could not detect in our numerical experiments any signs of scale free correlations. This may be associated with the simplicity of the geometrical structure of interaction which prevents local events from triggering global rearrangements. As we have seen, this conclusion remain valid even in the presence of strong Gaussian disorder that is not sufficient in the 1D setting to generate avalanches with power law statistics. It is then natural to try to augment the model by introducing an additional spatial dimension.

### 5.1 The model

The simplest 2D extension of the models presented in the previous sections can be viewed as a series of 1D NN chains linked transversally by linear elastic springs (see Fig.16 ); the ensuing 2D lattice model can also be interpreted as an array of coupled FK chains Suziki (1968); Landau et al. (1993b).

Notice that in 2D the short range NNN or SG terms are not necessary because the nonlocal elastic interactions take place already in the NN model by virtue of an additional geometrical dimension. We also emphasize that while in 1D we could discriminate only homogeneous slips (see Fig.17), in 2D one can also describe the structure of the dislocation cores (see Fig.18).

We assign to each particle with coordinates defined on a square lattice () a horizontal scalar displacement field . We associate with linear longitudinal springs the strain measure and with nonlinear shear springs - the strain measure The total energy of this highly anisotropic lattice can be written as

(30) |

where the potential is assumed to be quadratic in and periodic in . More specifically, we set

(31) |

where is a periodic function (see Fig. LABEL:atomic111), in particular, in our dynamic (ODE based) model we use

(32) |

The two lattice functions are independent Gaussian random variables introducing quenched disorder.

Observe that in the limit , we must have and our 2D model reduces to a one dimensional model with periodic potential which we have already discussed in the previous sections (see Fig. 17). When the variables and are no longer independent and the corresponding long-range interactions can be revealed by minimizing out a linear variable which is lacking in the 1D model Kartha et al. (1995); Rasmussen et al. (2001).

An elementary slip in 2D can be illustrated already in the constrained model with if the adjacent atomic planes are displaced by one or several lattice spacings (Burgers vector) without changing the energy, see Fig. 17. Within the classical interpretation, one usually relabels the nearest atoms after the slip takes place: notice that the bond in Fig. 18 between the atoms and is broken, therefore, a new bond has to be added between the atoms and , which are now nearest neighbors. This approach gives rise to discontinuous displacement fields and requires additional phenomenological hypotheses to deal with evolving discontinuities. In our model we do not perform the relabeling and consider compatible representation of the lattice before and after slip, see Fig. 17. Notice that in this case, one cannot use linear elasticity and has to introduce a periodic potential which ensures lattice invariance with respect to finite shears.

If the plastic flow does not occur by a slip of complete atomic planes. Instead the slipped atomic planes end inside the crystal forming dislocations. In Figs. 18 and 18 we illustrate the role of relabeling in the classical representation of dislocations and compare it with the compatible description of large shears adopted in this paper.

To complete the formulation of the problem we assume that the lattice is subjected to a time dependent shear in a hard device. In our framework driving with constant strain rate can be presented as displacement controlled boundary condition

for . This boundary condition can also be rewritten in terms of the overall shear strain as

where is again the slow time playing the role of loading parameter. In longitudinal direction we assume the periodic boundary conditions given by

for . The overall axial strain must then satisfy

Observe that our elastic energy depends on two material parameters: , representing the ratio of the bulk and shear moduli and , the ratio of the Burgers parameter and the interatomic distance. In addition, the behavior of the system may be affected by the variance of the quenched disorder . The natural question is where in the space of parameters one can expect to see the critical regimes. For large the model becomes one dimensional, while for large , the dislocations disappear. When is small the model again becomes one-dimensional and for small we lose lattice trapping. At very large disorder, dislocation mobility is very limited (POP regime, Pérez-Reche et al. (2008)) and the critical behavior can hardly be expected. On the other hand, at small disorder, one can expect synchronization (SNAP regime, Pérez-Reche et al. (2008)). The detailed reconstruction of the corresponding boundaries presents a formidable task which is outside the scope of this paper. A limited parametric study shows robustness of the scale free behavior, in particular, we observed that the power law statistics was disorder insensitive in the interval .

### 5.2 The dynamic problem

As we have seen in the 1D case, a quasistatic evolution of the driven system can be studied either in a dynamic setting with small but finite viscosity or in an automaton setting corresponding to the limit of infinitely slow driving. In the 2D setting we first formulate the dynamical model with small viscosity and smooth periodic potential (32). Later we compare the results with the automaton model implying zero viscosity and operating with a piece-wise quadratic approximation of the periodic elastic potential.

To facilitate the Fourier representation we reintroduce the finite difference operators and

(33) |

Then the axial and shear stresses can be written as

(34) |

and the equation of motion (2) takes the form

(35) |

where

(36) |

In order to avoid numerically challenging non-local spatial derivative terms in (35), we use the Fourier method. The discrete Fourier transform of a 2D lattice function is defined by

(37) |

where is the wave vector. By mapping our dynamic equation (35) into Fourier space we obtain

(38) |

where we defined

and added an inhomogeneity due to quenched disorder

To solve (38) numerically we use a semi-explicit Euler time discretization. The axial strain is a linear term with respect to the displacement and we can treat it implicitly in Fourier space. The non linear term is first calculated in real space and then transformed into Fourier space and treated explicitly. As a result, we obtain

(39) |

where

and is the time step. After the field is known, the strains can be computed in Fourier space by simple multiplication. To obtain the real space solution we need to invert the Fourier transform by using the formula

(40) |

Notice that the FFT method automatically enforces periodic boundary conditions for the displacement field . For instance, the conditions are automatically satisfied. Driving in shear can be again implemented by means of decoupling of the homogeneous shear strain component

### 5.3 Numerical results in the dynamic problem

In this subsection we present the results of our numerical experiments with implicit-explicit FFT method described above. In our tests we always take and while our initial state is always dislocation free.

As a first illustration we show in Fig. 19 the strain state due to internal prestress inside an ellipsoidal inclusion. In the driven system such inhomogeneously stressed regions play the role of nucleation centers for dislocational dipoles (mimicking dislocation loops in 2D).

In Fig. 20 we show a fragment of the deformed lattice with two dislocation dipoles nucleated around an ellipsoidal imperfection as a result of the application of homogeneous shear strain on the boundary of the domain. The local shear stress given by , is indicated by colors: red, for positive shear stress and blue for negative shear stress. One can see that the stress reaches its maximum inside the four distinct dislocation cores. As the applied strain is increased the dislocations of different sign move in opposite directions. Eventually they reach the boundaries which means that two full atomic layer have slipped.

In Fig. 21, we show a typical dislocational configuration in a crystal with random quenched disorder subjected to a finite shear. Once again the colors indicate the location of dislocations with different signs. Here we already see several atomic planes that have experienced either singular or multiple slip. Notice also that the dislocational dipoles have a tendency to concentrate around regions with strong inhomogeneity.

The corresponding stress field at larger scale is shown in Fig. 22. We observe that dislocations, which interact through long-range elastic fields, organize themselves into complex microstructures which appear random but turn out to be highly correlated. In our simulations we see that in the steady state regime plastic activity reduces to intermittent dislocational exchanges between stable clusters containing multiple dislocational dipoles. While these clusters remain mostly unchanged from one cycle to another, they have a finite lifetime as in experimental observations reported in Jakobsen et al. (2006).

In Fig. 23, we show the evolution of the dislocation density during the first six cycles of loading/unloading which we compute by localizing and identifying individual dislocation cores. We observe an initial overshoot and the subsequent stabilization which indicates that the system has reached the shakedown state. In Fig. 23 we present the corresponding macroscopic strain-stress curves exhibiting plastic hysteresis. The size of the hysteresis loops decreases with cycling, which indicates that the rate of dissipation diminishes as the system get stabilized in the shakedown steady state.

One can see that the macroscopic description of plasticity in our 2D model is quite realistic. At the microscale we see that the initial incipient inhomogeneity is augmented and modified during the cycling. This leads to the formation of statistically stable distribution of residual stress which is largely independent of the initial disorder.

### 5.4 Statistics of avalanches in the dynamic model

The fluctuations associated with intermittent dislocational dynamics are detected in experiment through the acoustic emission (AE) Weiss et al. (2000, 2007). The of intensity of the AE is accessed by measuring the square of the voltage which is expected to be proportional to the dissipated energy. Therefore by computing the rate of dissipation

(41) |

in our simulations we can make comparison with experimental data. The typical time series for is shown in Fig. 24 at two different scales. One can see that the signal has a characteristic spiky appearance and that the complexity of the time series is not reduced with magnification which is an indication of scale free behavior.

To separate individual avalanches we introduce an irrelevant threshold and define the avalanche energy by integrating the dissipation rate over the duration of the avalanche :

(42) |

The complexity of the intermittent signal is characterized by the probability distribution of avalanches which in our case stabilizes after few first cycles. In Fig. 25 we show the stationary distribution of avalanches exhibiting a robust power law pattern for over 4 decades with exponent obtained by the maximum likelihood method Clauset et al. (2009). This value of the exponent is in perfect agreement with experiments in ice crystals and fits the generally accepted range 1.4-1.6 Dimiduk et al. (2006); Zaiser et al. (2008); Brinckmann et al. (2008); Richeton et al. (2005a); Dimiduk et al. (2006); Csikor et al. (2007). It is also consistent with the value obtained for 2D colloidal crystals Pertsinidis and Ling (2005) and with the value obtained in phase field crystal simulations Chan et al. (2007). Curiously, the Gutenberg-Richter law for earthquakes potencies/moments, interpreted as statistics of dissipated energy, gives a very close value of the exponent Bak et al. (1991); Benzion (2008). This is not too surprising, however, in view of the essentially two dimensional nature of the fault friction/plasticity.

Due to very small scale of stress fluctuations the dissipated energy during an avalanche is approximately proportional to the size of the total plastic slip before the system stabilization. The probability distribution of the plastic strain increments

(43) |

is shown in Fig. 25 and, as expected, we again find the power law distribution with exponent . In this graph the abscissa is normalized by the system size, and therefore one can see that the biggest events traverse the whole system. On the other hand, the lower bound in our scaling region corresponds to the smallest plastic strain increment which scales with the unit cell size (viscous scale is even smaller). This suggests that the cut off have purely geometric nature which is a signature of a scale free regime.

In addition to energy it is also of interest to study the statistics of the avalanche durations Maslov (1996); Pérez-Reche et al. (2007); Dimiduk et al. (2010). The results shown in Fig. 25 suggest a power law scaling with the exponent , see Fig. 25. This value can be compared with measurements in micro-crystals where the duration exponents were found to be consistently bigger than Dimiduk et al. (2010).

We now turn to the study of the power spectrum for dissipated energy which carries information about temporal correlations not only between avalanches but also inside a single avalanche Hamon et al. (2002). The log-log plot of the computed power spectrum is shown in Fig. 26. It has a characteristic structure of the noise with . Such signals, also known as pink or flicker noise, are encountered in many self organized critical systems ranging from astronomy to physiology Bak et al. (1987), moreover the self organized criticality was originally proposed as an ’explanation of 1/f noise’ Bak (1987).

It is important to point out that the time series obtained in our dynamical model result from slightly overlapping avalanches. This finding is consistent with the results of the numerical simulations in ”running” sand pile models reported in Kardar (1992), where the appearance of noise was explicitly linked to the avalanches overlap at finite driving forces and where a very close RG estimate of the exponent in 2D case was obtained. Our results also agree with the scaling relation derived in Kertesz and Kiss (1990) under the assumption of random superposition of individual avalanches.

Another important signature of a universality class is the average shape of avalanches with a given duration . The avalanche shapes have been studied in different physical systems exhibiting critical behavior Houston (2001); Colaiori et al. (2004); Papanikolaou et al. (2011), in particular, plastic avalanches are known to be slightly asymmetric (skewed) Laurson and Alava (2006). The average shapes of avalanches in our dynamic model are shown in Fig. 27. We observe that the avalanche shapes in the 2D model show much more variability with respect to durations than in the 1D model. For instance, at longer durations the avalanche in our model are getting more symmetric, which is also a prediction of the mean field model Papanikolaou et al. (2011), however, the agreement is not complete because some asymmetry always remains. In our calculations, we observed that the number of units close to the marginal stability limit is increasing as the avalanche evolves towards its maximum which may be the origin of the tail responsible for the observed asymmetry. In the dynamic model we were not able to perform the scaling collapse because of the large computational cost associated with explicit resolution of the fast events. This will be done later on in the automaton model.

### 5.5 Spatial correlations in the dynamic model

The spatial counterpart of the observed time correlations is the fractal structure of dislocational patterns. Fig. 28 shows the energy density map of the system after the fifth cycle. We observe that the energy distribution exhibits spatially separated local peaks associated with dislocation-rich regions which are rather stable in time. The distribution of these peaks is only partially guided by the quenched disorder and is largely a result of self organization.

In order to characterize quantitatively the spatial complexity of this distribution we constructed a level set representation for the energy density and identified the boundaries of the spatially separated local peaks by selecting an irrelevant threshold (see the insert in Fig. 29). We then computed the energy density associated with these regions and obtained its probability distribution which is a power law with exponent , see Fig. 29. The obtained scaling is very robust spreading over five decades.

Another way to quantify the fractal clustering is to compute the correlation function of the actual dislocation distribution shown in Fig. 30. We define

(44) |

where is the total number of dislocations and is the number of pairs of dislocations with separation less than Hentschel and Proccacia (1983). If with non-integer exponent , this exponent is called the fractal dimension. In particular, for a randomly distributed (not correlated) set of points, the fractal dimension is equal to the dimension of space ( in our case).

In our numerical simulations we observed that during the first loading cycle , which is expected given the random nature of the quenched disorder. With cycling, however, the long range correlations developed progressively and in the shakedown regime we recorded independently of the initial disorder, see Fig. 31.

Interestingly, the dislocation patterns with have been observed experimentally in crystals with multiple slip systems; in simulations with a single slip system fractal patterning has been previously linked to the possibility of dislocation multiplication Groma and Bakó (2000) which is operative in our model.

Another well-known method of identifying the fractal structure of a set is computing its box-counting dimension Barnsley (1988) . It requires a covering of a set by boxes of size and counting the asymptotics of the number of boxes as . If then is the desired fractional dimension. Application of this method to the contour plot shown in the insert in Fig. 29 gives again , see Fig.31 which confirms the consistency of different measures of spatial correlations.

In summary, our 2D dynamical model shows a broad interval of self organized critical behavior characterized by a series of critical exponents revealing temporal and spatial correlations (, , and ). In the domain of parameters where we studied our system, these exponents do not depend on and and are not affected by the level of disorder . The system self-tunes to criticality by developing a highly correlated structure of persisting inhomogeneities which ensures effective energy dissipation.

### 5.6 The automaton model

Despite the conceptual transparency of the dynamical model, the mechanism of reaching the critical regime remains obscure in this largely computational formulation. In an attempt to obtain a mathematically more transparent version of the model we reduce in this subsection our dynamic equations to an automaton model and compare the results of conceptually similar numerical experiments. As in 1D case, we utilize the fact that in the quasi-static limit viscous relaxation is instantaneous and the driven system can be viewed as almost always residing in the state of mechanical equilibrium. This dramatically simplifies the tensorial problem because instead of dynamic equations of visco-elasticity one needs to deal only with equilibrium equations of elasto-statics

In order to solve the equilibrium equations analytically we again replace the smooth periodic potential shown in Fig. 12 by its piece-wise quadratic approximation defined in each period by

Here is the integer-valued spin variable describing a quantized slip. In the dynamic version of the model we we always implicitly assuming that .

Suppose first that the lattice field is given. Then the equilibrium equations can be written in the form

(45) | ||||

This linear problem for the displacement field can be mapped into the Fourier space giving an algebraic equation

(46) |

where and all other notations have been already introduced in 5.2. By solving (46) we obtain

(47) |

where . The Fourier image of shear strain can be then calculated as .

It is now straightforward to reformulate the dynamical problem as an integer valued automaton. We first define the residual shear strain due to quenched disorder as . Then we observe that the variable

representing shear strain associated with the slip must be confined between the limits

One can see that the stability of a particular slip configuration is controlled by two random thresholds which also evolve with the loading (with time).

When one of the thresholds is reached in one lattice point, the integer field is updated

where

After each increment of loading one has to check the stability of all units and continue the update until all units are stabilized.

The update of the ”slope” field can be represented in Fourier space as

Here the kernel

can be viewed as the analog of the toppling matrix in the sand-pile models. The inverse Fourier image of this kernel in the real space is highly anisotropic, long-range and conservative (see Fig. 32). The conservativeness is understood here in the sand pile terms and means that after each re-distribution of shear strain (discharge), the total deformation , which is controlled by the boundary conditions, remains unchanged.

In Fig. 32 we show the discrete map of the real space of the Green’s function . The discrete map indicates the distribution of the strain for each discrete unit (spring) in response to an elementary slip which can also be viewed as a localized force dipole. From this map one can see see that if an element reaches a threshold an strain increase takes place in this unit. Away from this element, in the direction, the strain is decreased in two neighboring elements but then it again increases reaching zero at infinity. In the direction, the elastic strain simply decays. Observe that while the structure of our is different from the quadrupolar structure of the 2D elastic kernel in isotropic elasticity Talamali et al. (2011), the main anti-ferromagnetic structure in the direction of the slip is preserved. In the limit the kernel vanishes everywhere except , where for . One can see that the 2D model reduces in this limit to the 1D model with the kernel (23) (NN version with ).

### 5.7 Numerical results for the automaton model

As we have already noticed, the automaton representation greatly reduces the complexity of numerical computations which, in particular, allows one to deal with much bigger domains. In the automaton model fast depinning events are replaced by jumps and outside the jumps the system evolves through a succession of equilibrium states. One can expect that such simplification of the model preserves the basic correlations between the avalanches but may affect the avalanche shapes. It may also affect the power spectrum particularly at high frequencies. In this section we begin a systematic comparison of the two models.

The stress-strain curves in the automaton model are shown in Fig. 33.

We observe that the macroscopic response remains basically the same as in the dynamical model. The small difference in the configuration of the hysteresis loops is due to the fact that in both models we compute the stress-strain relation only at some selected time steps. For instance, in the automaton model we know exactly the silent intervals and therefore can capture the strain/stress curve with higher resolution. Instead, in the dynamic model we randomly chose our test points which leads to stronger fluctuations.

A typical steady dislocation pattern in the automaton model is shown in Fig.34. One can see that the dislocations again form clusters without an obvious scale. A computation of the correlation function for this set of points gives a slightly bigger fractal dimension than in the dynamic model.

The study of the temporal correlations in the automaton model can be based on the time series representing the dissipated energy (12). The typical signal in the shakedown state is shown in Fig. 35.

Since now we have access to data for different domain sizes we can perform the finite-size scaling (FSS) analysis of the avalanche distribution Fisher (1971). In the critical state the power law probability distribution is expected to have a size dependent cut-off

(48) |

where is a universal shape function and is the cut-off energy with the scaling . According to the FSS hypothesis both, the exponents ( and ), and the shape function , define the universality class.

If the FSS hypothesis is valid we must obtain for large that , where . The exponent can be computed from the slope of the log-log plot of vs. and the slope of with respect to q is the cut-off exponent . As we see from Fig. 39, the linear behavior starts from with a slope . The exponent is close to the value obtained for plastic strain increments in Miguel et al. (2001); Zaiser and Nikitas (2007). Once is calculated, using , one finds the scaling relation that leads to . The data collapse in the coordinates and , presented in Fig. 36, shows an excellent agreement with the criticality hypothesis. The knowledge of the exponent allows one to introduce a correlation length and study the size effect associated with collective interaction of dislocations Roux (1996); Beato et al. (2011).

As we have already seen, one can also access temporal correlations through the analysis of the power spectrum for the dissipated energy.

The power spectrum for our automaton is shown in Fig. 37. Observe that the 2D power spectrum is different from the 1D spectrum shown in Fig. 15. It is, however, similar to power spectrum one obtained in our dynamic model although with a larger exponent (which is, interestingly, very close to the value obtained for the 2D BTW sand pile Laurson and Alava (2006)). The scaling part of the power spectrum is limited on the low frequency side by the duration of the longest avalanche. The absence of correlations in the low frequency region shows that the power spectrum reflects correlations within a single avalanches while missing correlations between different avalanches (see also Ref. Laurson et al. (2005)).

We now move to the study of the avalanche shapes. Recall that in a critical state one can expect to see avalanches of all sizes and the average avalanche shapes must scale with their durations and magnitudes in a universal way Kuntz and Sethna (2000). This scaling hypothesis can be thoroughly tested in the automaton model.

The first step is to average the avalanches of a given duration and plot the resulting function against the rescaled time . In the critical state one should see the following scaling Mehta et al. (2006)

(49) |

where is a universal scaling function and is the exponent in the power law relation linking the average avalanche size