Yielding and Irreversible Deformation below the Microscale: Surface Effects and Non-Mean-Field Plastic Avalanches

Nanoindentation techniques recently developed to measure the mechanical response of crystals under external loading conditions reveal new phenomena upon decreasing sample size below the microscale. At small length scales, material resistance to irreversible deformation depends on sample morphology. Here we study the mechanisms of yield and plastic flow in inherently small crystals under uniaxial compression. Discrete structural rearrangements emerge as a series of abrupt discontinuities in stress-strain curves. We obtain the theoretical dependence of the yield stress on system size and geometry and elucidate the statistical properties of plastic deformation at such scales. Our results show that the absence of dislocation storage leads to crucial effects on the statistics of plastic events, ultimately affecting the universal scaling behavior observed at larger scales.


Introduction
Over the past years, experimental investigations have gathered increasing evidence that plastic deformation of crystalline materials proceeds through intermittent bursts of activity [1][2][3][4][5][6][7][8]. Plastic flow advances through a sequence of strain avalanches of broadly distributed sizes. Avalanches have been observed experimentally under stress control conditions at various scales. Although such plastic fluctuations may be hard to detect at macroscopic scales, they dramatically affect mechanical properties of crystalline materials at smaller scales, ultimately producing detrimental effects on material formability [9]. The introduction of microcrystal compression testing [4] has allowed access to the microscale and introduced sample size as a crucial variable in this scenario. Size effects have dramatic consequences on yield, making smaller samples harder to deform and more unpredictable [6]. It is then natural to wonder if such behavior would hold unchanged below the micrometer scale.
While several nanoindentation techniques have been developed to measure material resistance to irreversible deformation and plastic flow, the experimental observation of plasticity at the nanoscale still represents an enormous challenge [10]. Colloidal crystals, however, were proven to deform plastically by activating dislocation motion, in remarkable analogy with crystalline materials [11][12][13][14]. The unmatched advantage of working with such systems is brought in by their size. Micrometer-sized colloidal crystals may contain few thousands, or even hundreds, of particles. This aspect allows one to project the problem of nanomechanics into a length scale that is easily accessible on experimental grounds. In this work, we address some of the challenges posed by nanomechanics, by providing a theoretical study of the mechanisms of yield and plastic flow in inherently small systems. We propose a simple geometry, which can be reproduced in experiments on two-dimensional micrometer colloidal crystals, and provide robust input for mechanical testing of crystalline thin films below the micrometer scale. Our aim is to show that at very small scales the irreversible deformation of materials proceeds in a novel way, deviating from the allegedly universal behavior observed in larger systems. To this end, we investigate the dependence of the yield stress on system size and geometry and the statistical properties of plastic deformation and energy dissipation in uniaxially compressed two-dimensional small crystals, by means of atomistic simulations and analytical modeling.

Methods
We first consider the compression of perfect crystals of various sizes and aspect ratios. Crystals are simulated as two-dimensional aggregates of short-range interacting monodisperse particles in their lowest energy configuration, that is a triangular lattice in the xy plane. Boundaries parallel to the y direction are in contact with rigid walls, while those parallel to the x direction are free (see Fig. 1). Uniaxial compression is applied symmetrically along the x direction. We consider two different compression protocols: i) displacement-control, in which rigid walls are quasi-statically displaced at equal constant velocities in opposite directions; ii) force-control, in which the force exerted on the walls is slowly increased at a constant driving rate. We consider a set of N particles located at the nodes of a triangular lattice and confined between two rigid walls (as in Fig. 1). Particles interact pairwise through a short-range potential V i j . For simplicity, we use a Lennard-Jones potential, but any other short-range one would lead to qualitatively similar results. We implement overdamped dynamics simulations. This choice is inspired by experiments on colloids, in the presence of a viscous carried fluid. The equation of motion for each particle i at position r i where C is the viscous friction coefficient of the carrier fluid, f is the interparticle force, d the characteristic size of the particles, and f d is the corresponding driving force for each deformation protocol. We choose as units of space and time d and t 0~C d 2 =V 0 respectively, with V 0 the amplitude of the interparticle potential, and we measure the driving force f d in units of V 0 =d. As a consequence we can take V 0~1 and d~1 without loss of generality. In dimensionless units, the linear system sizes L 0 x considered range from L 0 x~1 6 to L 0 x~7 2. Particles are confined by two rigid walls, that we model as two extra columns of particles of the same size, commensurate with the crystalline planes (see Fig. 1). In displacement control simulations, rigid walls are driven quasi-statically with a constant velocity v~0:005 in simulation units, whereas in the force control simulations the force exerted on the walls is slowly increased at a constant driving rate equal to 0:005. In both cases, these values are sufficiently small to avoid the overlap of plastic events. Smaller driving rates do not yield significantly different results. Upon compression, the strain c is measured as the ratio between the deformation along x and the size of the undeformed sample, i.e. c~DL x =L 0 x (engineering strain), while the applied stress s is the force per unit surface (L 0 y d) exerted on the rigid walls.
The coupled Eqs. (1) for i~1,:::,N are integrated numerically with an adaptive step size fifth-order Runge-Kutta method with precision 10 {6 . Thermal and hydrodynamic effects on the particles have been neglected. Thermal fluctuations can be disregarded because the characteristic time for dislocation motion is much faster than the characteristic time for thermal diffusion. Moreover in the problem at hand, elastic interactions would prevail over hydrodynamic interactions between the particles, if in suspension.

Elastic loading and plastic yield
In both protocols, the response is initially elastic. In a perfect crystal, the elastic limit is reached as soon as the motion of a pair of opposite sign edge dislocations is activated, as in Fig. 2, marking the beginning of the irreversible or plastic flow regime (see Supporting Videos S1 and S2). We use this limit to define the yield stress s y and the yield strain c y of the crystal.
By moving, dislocations allow the system to slip plastically and emit/dissipate part of the stored elastic energy. The value of both the yield stress s y and strain c y are found to be independent of the deformation protocol. Figure 3 shows the dependence of the yield stress on the geometry of the sample under examination. Different systems sizes L and aspect ratios r~L 0 x =L 0 y are considered. Even below the microscale, we qualitatively recover the inverse size dependence of the yield stress that makes smaller samples stronger.
According to the literature, in larger systems the dependence of the yield stress on the system size follows a power law [6], which corroborates the view of plastic yield as a critical phenomenon, dominated by scale free rearrangements of dislocation patterns over finite-size effects [15]. Plasticity is dominated by bulk phenomena and surfaces are somehow peripheral to such mechanisms. At smaller scales, our results show instead that s y is strongly affected by both size and geometry of the specimen. Size effects at this length scale were experimentally observed in compression tests in gold nanopillars [16], which were able to relate the increase in hardness to the dislocation starvation mechanism. Figure 2 shows that in our systems yield is mediated by a single dislocation pair (or even by a dislocation alone in the case of an imperfect crystal) and the onset of plasticity does not require the collective motion of a complex dislocation network. Dislocation interplay with boundaries thus becomes crucial. Early studies showed that the interaction of dislocations with rigid substrates is inherent to the strengthening of sheared thin films [17,18]. Nevertheless, the role of rigid boundaries in uniaxial compression experiments is still a matter of investigation.
The connection between the yield stress and the boundary effects can be easily visualized in our simple model system. Up to values of the stress very close to s y , the behavior of the system is elastic and the strain energy per unit length is approximated by where S is the specimen surface and Y r is the effective Young modulus. Here and in the following we consider energies per unit length, given the quasi two-dimensional nature of the problem. At yield the onset of dislocation motion is almost instantaneous compared to the dynamics of steady loading. Dislocations must account for the stress distribution inside the compressed sample, giving rise to the elastic energy density E. However, before dislocation motion is initiated, the elastic energy stored in the dislocated system must be comparable to the energy E e right before yield, as approximately no dissipation has occurred yet. This leads to the simple relation Equation (3) establishes a connection between the yield stress and dislocation strain distribution. It also bears implicitly the information about the dependence of the yield stress on the system size and geometry. The essence of our problem then lies in the stress distribution that accounts for E. In particular, an anomalous stress concentration is required at rigid boundaries, in order to enforce the condition of vanishing displacements. It is evident that boundary conditions alter dramatically the energy landscape E in the specimen and ultimately affect the yield stress, as prescribed by Eq. (3).

An estimate of the yield stress
By means of elasticity theory, we can demonstrate that the stress fields produced inside the sample by an edge dislocation close to a rigid boundary are long ranged and decay as *1=r, thus showing no screening effects, unlike dislocations close to free boundaries [19]. For simplicity, we consider a positive straight edge dislocation with Burgers vector perpendicular to the rigid wall, a distance l apart from the wall, as in Fig. 4 (left). The problem has translational invariance along the z direction (plane strain conditions): the wall is at x~0, the dislocation is placed in (l,0,z). We can calculate the exact displacement field u as u~u inf zu img zw where u inf is the displacement field of the original positive dislocation in (l,0,z) in an infinite medium, u img is the displacement of the image dislocation of opposite sign in ({l,0,z) and w is and additional field, which is analytic in xw0 and satisfies the equilibrium elastic equations with boundary conditions such that the full u(0,y,z)~0, being m the shear modulus and n the Poisson ratio. Eq. (4) with the given boundary conditions is commonly known as the 1{plane problem in linear elasticity. Several techniques can be devised to find exact solutions [20] and the complexity resides solely on the complicated form of boundary conditions. Here we follow [21] as reviewed by [20] and observe that Eq. (4) can be rewritten as a Laplace equation in the form where q~+ : w is the dilatation. The problem can then be solved using the Green function method, being the Green function for Eq. (5) in the 1{plane geometry known from classical electrostatics. While further details will be given in a future  publication [22], for the problem at hand we obtain From the above result, the elastic strain tensor e ij , the stress tensor s ij and the elastic strain energy density E can be calculated using linear elasticity. Interestingly, while the stress components derived from u inf zu img decay like *1=r 2 (being that the stress of a dislocation dipole), the stress due to w will decrease like *1=r. Overall, the stress field of an edge dislocation close to a rigid boundary will hence be long-ranged, as opposed to the case a free boundary. As a consequence, stress fields are able to sample the entire system size and are responsible for the size-sensitivity of the yield stress.
In the case of our simulations, however, Eq. (6) may seem of limited help in calculating the strain energy E at yield for our model system, where we have two fixed boundaries, two free surfaces and two dislocations of opposite signs, as the exact solution of such a problem would be great complexity. Nevertheless, an approximation of the integral in Eq. (3) can be provided as follows. Fig. 4 (right) shows the result of the superposition of the solutions for the shear stress of two edge dislocations of opposite signs, as calculated from Eq. (6), arranged in a configuration which mimics the one observed at yield in our simulations for the perfect crystal. For simplicity we consider Burgers vectors along x, but the general solution would lead to similar results [22]. In such a configuration, the condition of zero displacement at the vertical walls is not met anymore, however one can prove that the deviations from zero affect the stress and the strain energy in a negligible way if the dislocations are far apart. At the same time, Fig. 4 (right) shows that the two dislocations behave like a dipole, in the sense that the stress goes rapidly to zero outside the region enclosed by them, approximating the stress field close to free boundaries. Such observation is verified analytically and is valid for all stress components.
We can conclude that the configuration in Fig. 4 (right) provides an acceptable approximation for the elastic problem in our simulations and an estimate for the strain energy E in Eq. (3). Finally, upon inverting Eq. (3) and computing numerically the integrals of E, we evaluate the yield stress and its dependence on the system size and geometry. Predictions from our theory are given in Fig. 3 (solid lines) and are in remarkably good agreement with our simulation results, presented below (discrete points) , suggesting that boundaries and shape effects are essential in our problem. Given the nature of the approximations involved, we stress that ours is just a dimensional estimate. However the agreement with simulations corroborates the view that the role played by boundaries is integral to nanoplasticity.

Plastic avalanches
As soon as the yield point is reached, the response of the system to further loading differentiates depending on the deformation protocol. Under conditions of displacement control, stress-strain curves are characterized by serrated yielding, while they assume a staircase shape under conditions of stress control. We emulate realistic realizations of compressed samples by introducing randomness at free boundaries as follows. The initial state of each realization is obtained from the perfect crystal by extracting a random number of particles from one free surface and relocating them at random positions on the opposing surface. In this way both the number of particles and the linear size of all simulated specimens are kept constant, while the morphologies of their free surfaces are allowed to vary stochastically. Strain plateaus Dc in force control or stress drops Ds in displacement control always correspond to plastic events in which dislocation motion is reactivated in order to reduce the elastic energy stored during prior loading phases (see Fig. 1). Remarkably, a statistical analysis of event sizes over several realizations of surface disorder ( Due to the limited system size, moving dislocations easily leave the sample through free boundaries. Pioneering studies have shown that in sub-micrometer Ni samples, pure mechanical loading can induce dislocation depletion within the sample [23]. In our case plastic flow proceeds through the activation and motion of few dislocations at a time, while no storage is observed (see Supporting Video S3). Such phenomenology differs from the widely accepted picture of plastic flow at larger scales. Yet, our statistical study leads to broadly distributed plastic events for several system sizes, suggesting that kinematic constraints and long range dislocation interactions still rule plastic flow, as shown in Fig. 5. The cumulative distribution is defined as p c (x)~ð  Plastic event sizes are commonly quantified in experiments by looking at the amount of energy W dis released during each event or, under stress control conditions, by recording the magnitude of platen displacements. By defining the stored potential energy as U~X ivj V ij , each variation of U in our simulations will be given by the energy balance relation Here the dissipated energy W dis can be measured as the time integral of the dissipated power, i.e. W dis : (force time velocity for all particles), whereas W ext is the work done by the external load, W ext~ð sdc. The strain is the sum of its elastic and plastic components c~c el zc p . As in our system no defect accumulation is encountered, we have Ds!Dc el (linear elasticity) and evidently DU*(Dc el ) 2 .
Let us first consider the case of force control. If focusing on a single platen displacement event, we have Ds^0 and Dc el^0 . This implies that the change in stored elastic energy is roughly zero, DU^0, and, according to Eq. 8, W dis^Wext , which means that the dissipated energy comes essentially from the work done during the platen displacement event of magnitude Dc^Dc p . Then evidently W dis *Dc p , as the amount of slip due to a moving dislocation is proportional to the energy dissipated in the process [24]. This relation indeed is verified in our simulations , as shown in Fig. 6(a). Experimental studies at larger scales show that both quantities follow a power-law distribution, decaying with an exponent close to the mean-field value t~3=2 [24]. The striking aspect of the exponent t is its universality [9]. In our model, which aims to reproduce nanoscale plastic flow, the exponent t c *0:9 deviates substantially from the universal value, suggesting that at such scales the microscopic dislocation dynamics giving rise to plastic deformation are qualitatively different. Signals of univer-  sality breakdown also come from experiments of compressed highpurity LiF micropillars [25].
Under displacement control conditions, instead, the statistical analysis of plastic flow and dissipated energy can be performed by looking at the distribution of stress drops (Fig. 5(b)). In Fig. 6(b) we show the evolution of the energy balance during compression. For quasi-static driving, stress drops occur almost instantaneously, i.e. at constant total strain (Dc^0, or W ext^0 ). This implies that, during the stress drop, a sudden decrease of the elastic strain Dc el of the crystal must occur at the expense of increasing its plastic strain, jDc el j^jDc p j. According to Eq. 8, being W ext^0 during the stress drop, the energy dissipated by plastic avalanches must balance the change of internal energy, DW dis^D U. Then, as in the the case of force control, since there is no accumulation of large numbers of defects, Ds!Dc el and DU*(Ds) 2 . Being stress drops distributed according to a power-law of exponent t s , we obtain the following exponent relation t E~( t s z1)=2 for the dissipated energy distribution, which for t s *0:7 yields t E *0:85, in good agreement with the results obtained for t c in the force control protocol. We emphasize that stress drops in displacementcontrolled deformation are in principle different from stress increments between bursts in force-controlled tests. Recent experimental studies have shown the stress increments are Weibull-distributed in Mo micro-and nanopillars [26].
Compared to numerical studies of size effects in dislocation dynamics at the microscale [9], which highlight a strong size dependence of the avalanche distribution cutoff, our small systems span a very narrow linear size range. Plasticity is thus always mediated by equally small numbers of dislocations and maximum avalanche sizes for systems with initial size L 0 x~1 6, 24 and 32 do not give rise to appreciable differences in the avalanche-size cutoff as plastic flow advances. It should also be noted that by definition of c, platen displacements are always Dcv1, and similarly under displacement control, large stress drops are bounded by the very nature of the short-range interatomic potential considered.
As for the origin of the anomalous avalanche exponents, we should remark that the novel behavior is related to the inability of the system to store large numbers of dislocations. In the absence of collective behavior, plastic flow departs from the traditional picture of cooperative dislocation organization. In fact our simulations show a behavior which approximately recalls a sequence of load-unload events, much in the spirit of stick-slip dynamics or fracture/failure mechanics. We notice that an energy release exponent very close to 1 has been recently encountered in simulations of yielding and failure of heterogeneous materials in the plastic steady state [27]. However the rationale behind the analogy with our finding remains to be ascertained. Finally, a statistical analysis of surface roughness of strained systems would provide further insights into the nature of plasticity at such scale. However, linear sizes of our free surfaces are only one order of magnitude larger than the lattice constant and such a study would not lead to significant results. Further studies in this direction are currently under way.

Conclusions
In conclusion, we have shown that the onset of plasticity at small scales is mediated by few dislocations. The number and arrangement of nucleated dislocations must account for the distribution of stress stored inside the crystal during the elasticloading regime, allowing one to estimate the dependence of the yield stress on sample size and geometry. Our results confirm that both size and shape are crucial factors in determining the strength of materials at these scales. We find that plastic flow occurs in an intermittent manner reminiscent of irreversible deformation at larger length scales. Plastic avalanches of broadly distributed sizes are still observed, however, the absence of dislocation storage has important effects on the scaling characteristics of viscoplastic dynamics, which ultimately violate the universal mean-field behavior observed at larger scales. Our results are thus a significant example of source-limited deformation [23] and arise naturally from the impossibility of such a small system to store dislocations . The new exponent values are obtained for a twodimensional crystal geometry and should be relevant for thin films of several self-assembled nanoparticles under external loading conditions and amenable of experimental analysis in colloidal systems.

Supporting Information
Video S1 The first appearance of a dislocation pair signals the onset of yield for the perfect system. As soon as each dislocation reaches the opposite rigid boundary, plastic activity stops and the the first event is over. The animation consists of 7 snapshots of the dynamics. Top: dislocations are represented as pairs of 5-and 7-