Local yield stress statistics in model amorphous solids
We develop and extend a method presented in [S. Patinet, D. Vandembroucq, and M. L. Falk, Phys. Rev. Lett., 117, 045501 (2016)] to compute the local yield stresses at the atomic scale in model two-dimensional Lennard-Jones glasses produced via differing quench protocols. This technique allows us to sample the plastic rearrangements in a non-perturbative manner for different loading directions on a well-controlled length scale. Plastic activity upon shearing correlates strongly with the locations of low yield stresses in the quenched states. This correlation is higher in more structurally relaxed systems. The distribution of local yield stresses is also shown to strongly depend on the quench protocol: the more relaxed the glass, the higher the local plastic thresholds. Analysis of the magnitude of local plastic relaxations reveals that stress drops follow exponential distributions, justifying the hypothesis of an average characteristic amplitude often conjectured in mesoscopic or continuum models. The amplitude of the local plastic rearrangements increases on average with the yield stress, regardless of the system preparation. The local yield stress varies with the shear orientation tested and strongly correlates with the plastic rearrangement locations when the system is sheared correspondingly. It is thus argued that plastic rearrangements are the consequence of shear transformation zones encoded in the glass structure that possess weak slip planes along different orientations. Finally, we justify the length scale employed in this work and extract the yield threshold statistics as a function of the size of the probing zones. This method makes it possible to derive physically grounded models of plasticity for amorphous materials by directly revealing the relevant details of the shear transformation zones that mediate this process.
Despite numerous advances during the last two decades, a physical description of plasticity in amorphous materials, known to be quantitatively tied to well-characterized atomistic processes, remains a grand challenge Procaccia (2009); Barrat and Lemaître (2011); Rodney et al. (2011). All the constitutive laws describing the plastic flow of this large class of materials, such as glasses, amorphous polymers or gels, remain based on phenomenological assumptions. This fact is due to the lack of systematic characterization of elementary mechanisms of plasticity at the atomic scale. For the amorphous solids, the absence of crystalline order prevents, by definition, any identification of crystallographic defects such as dislocations. These defects are, however, the quanta of plastic deformation from which it has been possible to derive constitutive equations of crystal plasticity on robust physical grounds. In amorphous materials, plastic deformation manifests as local rearrangements Argon (1979) exhibiting a broad distribution of sizes and shapes Lerner and Procaccia (2009), nonaffine displacements Zaccone et al. (2014); Laurati et al. (2017) and connectivity changes between particles van Doorn et al. (2017) that lead to a redistribution of elastic stresses in the system Chattoraj and Lemaître (2013); Lemaître (2014). By analogy with dislocations, it therefore appears natural to try to describe the plastic flow from the dynamics of localized âdefectsâ commonly referred to as Shear Transformation Zones (STZs) Falk and Langer (1998).
A variety of metrics have been proposed to locate and characterize the defects that control plastic activity including structural properties (free volume Spaepen (1977), packing Jack et al. (2014), short-range order Shi and Falk (2007) and internal stress Tsamados et al. (2008)) and linear responses measures (elastic moduli Tsamados et al. (2009) and localized soft vibrational modes Widmer-Cooper et al. (2008, 2009); Tanguy et al. (2010); Ghosh et al. (2011); Manning and Liu (2011); Chen et al. (2011); Mosayebi et al. (2014); Ding et al. (2014); Schoenholz et al. (2014)). Unfortunately, the definition of these structural properties are often system-dependent and have shown a relatively low predictive power with respect to the plastic activity Tsamados et al. (2009); Jack et al. (2014). The approaches based on linear response measures (e.g. soft vibrational mode analysis) have shown that the correlation between these local properties and the location of plastic rearrangements decreases rather quickly as the system is deformed plastically since they are derived from perturbative calculations Patinet et al. (2016).
To address these problems, new methods have recently been proposed. They are based either on combinations of static and dynamic properties (atomic volume and vibrations) Ding et al. (2016), on the non-linear plastic modes Lerner (2016); Zylberg et al. (2017) or on machine learning methods Cubuk et al. (2015, 2016); Schoenholz et al. (2016). These approaches allow the calculation of local fields (respectively named by their authors flexibility volume, local thermal energy and softness field) which show abilities to detect plastic defects far superior to previous attempts. Independently of their high degree of correlation, they nevertheless have the disadvantage of giving access only to quantities that are not directly related to local yield criteria that are more commonly used in models of plasticity. Moreover, with the exception of works based on the vibrational soft modes Rottler et al. (2014); Smessaert and Rottler (2014), the vast majority of these previous approaches only gives access to scalar quantities, which by definition also neglect the important orientational aspect inherent in STZ activity.
It is precisely in this context that we have recently developed a numerical technique to systematically measure the local yield stress field of an amorphous solid on the atomic scale Patinet et al. (2016). This new approach responds to some issues raised previously by providing access to a natural quantity in the context of plasticity, i.e. a local yield stress, and shows an extremely strong correlation with the plastic activity. In addition, this method is non-perturbative and can investigate large strains. It also gives access to a tensorial quantity and is thus able to describe several possible directions of flow. It is therefore an ideal candidate to quantitatively characterize the relationship between structure and plasticity.
In this paper, we present in detail the principles of this method. The statistics of local yield stress are calculated in a model glass synthesized from different quench protocols. The correlation between local slip threshold and plastic activity is investigated as a function of the degree of relaxation of the system. The method is subsequently extended to the study of the amplitudes of the plastic relaxations. Additionally, the consequences of the orientation of the mechanical loading are examined. Finally, we address the effect of the length scale over which the local yield stress field is computed.
Ii Simulation methods
ii.1 Sample preparation
We performed molecular dynamics and statics simulations with the LAMMPS open software Plimpton (1995). The object of study is a two-dimensional binary glass which is known for its good glass formability Lançon et al. (1984); Widom et al. (1987). It has previously been used to study the plasticity of amorphous materials Falk and Langer (1998); Shi and Falk (2005, 2005); Shi et al. (2007). One hundred samples each containing atoms were obtained by quenching liquids at constant volume. The density of the system is kept constant and equals . As in Shi and Falk (2005), we choose our composition such that the number ratio of large (L) and small (S) particles equals . The two-types of atoms interact via standard Lennard-Jones interatomic potentials. In the following, all units will therefore be expressed in terms of the mass and the two parameters describing the energy and length scales of interspecies interaction, and , respectively. Accordingly, time will be measured in units of . In the present study, these potentials have been slightly modified to be twice continuously differentiable functions. This is done by replacing the Lennard-Jones expression for interatomic distances greater than by a smooth quartic function vanishing at a cutoff distance . For two atoms and separated by a distance :
Periodic boundary conditions are imposed on square boxes of linear dimensions . For the non-smoothed version of the interatomic potential, the glass transition temperature of this system is known to be approximately , where is the Boltzmann constant Shi and Falk (2005). This temperature corresponds to the mode coupling temperature which is an upper bound of the glass transition temperature Shi and Falk (2005). In order to highlight the links between the microstructure, the stability of glasses and their mechanical properties, three different quench protocols are considered. The first two kinds of glass are obtained after instantaneous quenches from High Temperature Liquid (HTL) and Equilibrated Supercooled Liquid (ESL) states at and , respectively. The last protocol consists in a Gradual Quench (GQ) in which temperature is continuously decreased from a liquid state, equilibrated at , to a low-temperature solid state at , over a period of using a Nose-Hoover thermostat Nosé (1984); Hoover (1985). Afterwards, the system is quenched instantaneously as well. All quench protocols are followed by a static relaxation via a conjugate gradient method to equilibrate the system mechanically at zero temperature. The forces on each atom are minimized up to machine precision. The same relaxation algorithm is used hereafter to study the response to mechanical loading.
This approach produces three highly contrasting types of amorphous solids. The greater the temperature from which the system has fallen out of equilibrium, the less relaxed the system Sastry et al. (1998); Debenedetti and Stillinger (2001). This fact is clearly reflected in the values of the average potential energies per atom of the generated inherent states, equal to , and for the HTL, ESL and GQ protocols, respectively.
ii.2 Mechanical loading: generation of plastic events
Beginning from a quenched unstrained configuration, the glasses are deformed in simple shear imposing Lees-Edwards boundary conditions with an Athermal Quasi Static method (AQS) Malandro and Lacks (1998, 1999); Maloney and Lemaître (2004, 2006) . We apply a series of deformation increments to the material by moving the atom positions following an affine displacement field such that and . After each deformation increment, we relax the system to its mechanical equilibrium.
In order not to miss plastic events, a sufficiently small strain increment equal to is chosen. Plastic events are detected when the computed stress decreases, a signature of mechanical instability. A reverse step is systematically applied after each stress drop to confirm that the strains generated in the solid are irreversible when the stress criterion is satisfied. The observed response is typical for amorphous materials and is characterized by reversible elastic branches interspersed by plastic events as illustrated in Fig. 1a. The more relaxed the system, the stiffer and harder the glass. In agreement with Shi and Falk (2005), we observe that the localization of plastic strain, under load, increases with the degree of relaxation of the initial state. Note that in the case of the GQ system a strong localization of the strain is observed around due to shear banding.
The average shear moduli are obtained from the ratio between the stress response of the entire system following a deformation increment , i.e. the slope at the origin of the stress-strain curves reported in Fig. 1a. equals , and for HTL, ESL and GQ systems, respectively. As expected, the stiffness of the system increases with its level of relaxation.
ii.3 Strain field computation
In order to quantify the correlation between deformation thresholds and plastic rearrangements two types of deformation fields are calculated. The first corresponds to the plastic deformation induced by a single plastic event, the second describes the total cumulative deformation. The displacement field of the former is calculated using the difference between the position of atoms after and just before that instability occurs minus the applied affine displacement increment. The displacement field of the latter is merely computed as the difference between the position of atoms in configurations at a given strain and the as-quenched state, that is to say the state that has not yet been deformed mechanically. The Green-Lagrange strain tensor is then evaluated from displacement fields of each atom following the coarse-graining method developed in Hinkle et al. (2017) based on the atomic gradient tensor evaluation Zimmerman et al. (2009). An octic polynomial coarse-graining function is employed Lemaître (2014). This function has a single maximum and continuously vanishes at where is the distance between the strain evaluation location and the atom positions. It is expressed as:
It is desirable to consider a large enough coarse-graining length scale so that continuum mechanics quantities make sense while keeping it as small as possible in order to account for heterogeneity and to preserve spatial resolution. To this aim, we choose . On this scale, a continuous description makes sense (Hooke’s law holds) but the solid is still anisotropic and heterogeneous Goldhirsch and Goldenberg (2002); Goldenberg et al. (2007); Tsamados et al. (2009).
To simplify the analysis, we choose to work with a scalar quantity by computing the maximum of the shear deformation . The positions of a plastic rearrangement are then defined as the position of the atom having undergone the maximum during a plastic event. This approach allows us to obtain the successive positions of localized plastic rearrangements during deformation from the quenched state as exemplified in Fig. 1b.
Iii Local yield stress measurement method
We used a method developed in Patinet et al. (2016) which allows us to sample the local flow stresses of glassy solids for different loading directions. Similar techniques have been employed to sample the local elastic moduli Mizuno et al. (2013) or the yield stresses along a single direction in model glasses Sollich (2011); Puosi et al. (2015). The principle of the numerical method is illustrated in Fig. 2a. It consists in locally probing the mechanical response within an embedded region of size (named region I) by constraining the atoms outside of it (named region II) to deform in a purely affine manner. Only the atoms within the region I are relaxed and can deform nonaffinely. Plastic rearrangements are, thus, forced to occur within this region and the local yield stress can be identified.
In Patinet et al. (2016), the embedded region I was centered on every atom of the system to test the reliability of the method. Here, to lay the groundwork for an up-scaling strategy, we rather sample the local yield stress on a regular square grid of size . Furthermore, the non-interacting atoms, located farther than a distance from the center of the probed region, are deleted during the local loading simulations, thereby speeding up the computation. Unless mentioned explicitly, we chose and , consistently with the coarse-graining scale used for the strain computation. The yield criterion and the AQS incremental method are the same as those described in section II.2 to shear the system remotely. At this scale the amorphous system is highly heterogeneous and the yield stress may not be the same for all orientations of the imposed shear. We thus sample the mechanical response using pure shear loading conditions in different loading directions . Eighteen directions, uniformly distributed between and , are investigated. corresponds to the remote simple shear direction .
To save computational time, the strain increment is this time equal to . The shear stress in the direction is computed over the atoms that belong to region I using the Irving and Kirkwood formula Irving and Kirkwood (1950) as a function of the applied strain. The local region is sheared up to the first mechanical instability occurring at a critical stress . Even at rest, the glasses feature non-zero internal stress due to the frustration inherent to amorphous solids (see inset of Fig. 2b). A more relevant quantity to link the local properties with plastic activity is thus the amount of stress needed to trigger a plastic rearrangement. The initial local shear stress state within the as-quenched glass is thus subtracted from the critical stress to get the local stress increase that would trigger an instability . Local shear stress-strain curves are exemplified for four different directions in Fig. 2b. It can be observed that the mechanical response depends on the loading orientation. As expected from elasticity theory, we verify that and that the local shear moduli as shown in the inset of Fig. 2b
We now want to consider the implications of the field of local for a particular direction of remote loading . For this purpose, we make the simplifying assumption of homogeneous elasticity within the system or, equivalently, of localization tensor equal to the identity tensor. Of course, the elasticity in this system at this length scale is heterogeneous and leads to non-affine displacements under remote loading as shown in Tsamados et al. (2009). This assumption is simply used to be able to estimate the stress felt by a local zone due to a remote loading. This assumption will be discussed further in section V. If the applied shear stress is homogeneous in the glass, plastic rearrangement that would be activated for a given site is the minimum (positive) projected along the remote loading direction. This may be expressed as:
Maps of local are shown in the Fig. 3 for the three quench protocols. One distinguishes a correlation length corresponding to the size of region I. Indeed, the same shear transformation zone can be activated for several grid locations if its threshold is smaller than others in its vicinity which leads to the assignment of close values on a scale . The influence of the size of region I will be addressed in section VI.
Iv Local rearrangement statistics
iv.1 Distributions of local yield stress
The effect of the quench protocol on the yield stress maps shown in Fig. 3 is remarkable. It is readily apparent that the lower the temperature at which the system falls out of equilibrium during its synthesis, the more mechanically stable the glass. An advantage of our method is that it allows one to assess stability not from the global scale, as in Fig. 1a, but locally. The HTL system shows an overabundance of small energy barriers characteristic of systems far from equilibrium. In contrast, the GQ system has a low proportion of soft zones embedded in a hard skeleton Shi and Falk (2005, 2006). As expected, ESL presents an intermediate situation. More quantitatively, the distributions of are computed for the three quenching protocols as shown in Fig. 4. The probability densities are noticeably shifted toward higher values with increasing system stability, weak areas being depopulated.
Through this method we are able to analyse the statistics of the sites that are about to rearrange plastically. Previous work (based on mean-field theoretical approaches Lin and Wyart (2016), atomistic simulations Karmakar et al. (2010); Hentschel et al. (2015) and mesoscopic simulations Lin et al. (2014, 2015)) proposed a scaling for these soft areas such as where is a non-trivial exponent. For systems at rest, i.e. after the quench, it was shown that Karmakar et al. (2010); Hentschel et al. (2015). From detailed inspections of our results, as shown as inset of Fig. 4, it seems difficult to extract this exponent with the exception of ESL. In any case, it appears that the behavior of in the limit of zero varies with the preparation of the system. In the case of the GQ system, a smaller exponent can even be observed while for HTL a larger one can be fitted for small threshold values. Several reasons may account for this disagreement, including the lack of statistics or the size of the strain increment . The fixed boundary conditions applied during local probing may also prevent some relaxation of the system. Also, notably, the previous approaches have considered the distribution of critical strains applied to the whole system which is strictly equivalent to our approach for an elastically homogeneous glass, a strong assumption at this length scale Tsamados et al. (2009). The answer to this question deserves more investigation, which is outside the scope of the present study.
iv.2 Correlation with plastic activity
The position of the first ten plastic rearrangements during remote loading are illustrated in Fig. 3. Plastic rearrangements clearly tend to occur in the soft zones, i.e. for areas in which are small. To quantify this correlation, we apply the same method as in Patinet et al. (2016).
We propose a correlation coefficient that allows us to relate the order of appearance of zones in which the plastic arrangements appear and a local scalar field, the local yield stress. The aim here is to compute the predictive power of a structural indicator for the location of successive rearrangements from the sole knowledge of the initial state of the system, i.e. before deformation. To achieve this, the correlation coefficient is computed from the value of the cumulative distribution function of corresponding to the point of the grid , i.e. the closest to the location of the plastic rearrangement (determined according to the method described in the Sec. II.2). The correlation coefficient is defined as:
where is the disorder average of the cumulative distribution function. indicates a perfect correlation, i.e. a localized plastic rearrangement on the lowest yield threshold grid point (), while means an absence of correlation. is calculated for all plastic events as a function of the deformation applied as shown in Fig. 5a. Note that relation (5) neglects the stress redistribution due to successive rearrangements. Moreover, it only takes into account the rearrangements producing the maximum of local shear strain located at and therefore ignores the possibility that a plastic event may be composed of several localized rearrangements. In agreement with Patinet et al. (2016), an excellent correlation is observed. Above all, this correlation shows a slow decrease indicating a persistence of weak sites.
We observe that the level of correlation depends on the preparation of the system. The more relaxed the system, the more robust the observed correlation. The first ten plastic rearrangements occur in areas belonging to the softest , and sites for the HTL, ESL and GQ systems, respectively. For slowly quenched glasses GQ, it is interesting to note that the correlation decreases sharply for a deformation corresponding to the softening due to the localization of the deformation. It can be argued that the origin of the best correlation observed in the most relaxed system comes directly from the distribution of local thresholds. The relaxed systems have a much smaller population of low yield threshold zones. They therefore exhibit larger shear susceptibility when compared to other zones or to mechanical noise, making it easier to predict the onset of plastic activity.
The problem of the preceding method is that it assumes the existence of individual and well-localized events. However, Refs. Maloney and Lemaître (2004); Dasgupta et al. (2012) have shown that if this hypothesis is relatively well satisfied for small deformation levels, it does not hold with the increase in deformation during which avalanches, through system spanning plastic events, are observed. To circumvent this problem, we deal directly with the correlation between the entire local yield stress field of the as-quenched state and the cumulative deformation field in the same spirit as Ref. Smessaert and Rottler (2014). The cross-correlation, or Pearson’s correlation, is calculated as a function of the applied strain as:
where is the number of points on the grid on which the thresholds are calculated, and are the standard deviations of and , respectively. The minus sign is added here to obtain a positive value since large are expected for locations where are small (i.e. anti-correlation). denotes the ensemble average of the quantity . Note that explicit dependence on of is omitted in the r.h.s for the sake of simplicity. Fig. 5b shows the evolution of as a function of the imposed deformation. The general trend is qualitatively similar to that of Fig. 5a. The correlation between local thresholds and plastic activity is greater for the more relaxed systems. There are, however, some differences. It can be observed that the correlation begins to increase as plastic rearrangements start to accumulate on weak sites. On the other hand, as expected, the decay of the GQ system is more marked upon the formation of shear bands, the latter concentrating the deformation.
The correlation appears smaller than the one expected from the computation based on local rearrangements in Eq. 5, however, we remark that this calculation is based on crude assumptions. For example, we have not dissociated the elastic part from the plastic part when calculating . Moreover, this approach does not take into account the distribution of amplitudes of plastic rearrangements. Still, we have verified that in both cases - correlation based on the maxima of the strain field (Eq. 5) and the cross-correlation based on the cumulative deformation (Eq. 6) - the correlations between and plastic activity are significantly better, and more persistent with deformation, than those obtained for the local classical structural indicators reviewed in Patinet et al. (2016).
iv.3 Distributions of local relaxation
In this section, we extend our method to study the amplitude of the local plastic relaxations that follow plastic rearrangements. The loading of region-I described in Sec. III is continued after the instability until the local stress increases again, signaling the end of the plastic rearrangement and the return to mechanical stability of the sheared zone. This final stress is also computed for all directions. Plastic relaxation in the direction is then deduced by simply subtracting from the stress just before the instability . This amplitude of relaxation is exemplified in Fig. 2b. Like the thresholds, it depends on the direction of shear. Note that this method can only give an estimate of the relaxation amplitude, as the frozen boundary conditions constrain some degrees of freedom during relaxation.
The plastic rearrangements actually observed during the shearing of the system correspond to the thresholds calculated in Eq. 4, and occur when the patch is loaded in the weakest direction . To place ourselves in a coarse-graining perspective, we want to derive a scalar indicator that corresponds to mechanical response in the remote loading direction and disregard for now the tensorial aspect of the problem. The amplitude of plastic relaxation is therefore calculated, in turn, by projecting in the direction according to the relation:
Note that this estimator of the stress relaxation slightly underestimates the plastic relaxations because of the projection. Nevertheless, it gives access to a sufficiently simple scalar indicator. We verified that the absence of projection does not qualitatively change our results (not shown here).
The distributions of stress relaxation amplitudes reported in Fig. 6a in lin-log scale for the three quench protocols show an exponential decay. The average stress drops increases with the relaxation of the system. The mean plastic relaxations calculated from exponential regressions are , and for HTL, ESL and GQ systems, respectively. The amplitude of plastic deformation, or slip increment, can also be estimated by computing the eigen-deformations of plastic rearrangements as where is the average shear modulus of the glass. Remarkably, the distributions collapse on a master curve of mean independently of the quench protocol as reported in the inset of Fig. 6a. These results justify the assumption of a characteristic relaxation commonly used in mesoscopic simulations or in mean-field models. It is also in agreement with previous atomistic computations based on different methods such as the mapping between elastic field and Eshelby inclusion model Albaret et al. (2016); Boioli et al. (2017) and automatic saddle point search techniques Fan et al. (2017).
We also take advantage of this analysis to study the dependence of the relaxation amplitude with the distance to thresholds as shown in Fig. 6b. The former increases on average according to the latter. Remarkably, the relationship observed does not seem to depend too much on the quench protocol. If it seems reasonable that the amplitude of relaxation increases with the increase of the local yield stress, the stored elastic energy being larger, we have no explanation to derive this relationship at the moment. An exponential dependence is adjusted empirically and gives .
Let us note finally that if we are sufficiently confident in the capacity of this local method to quantify the thresholds, the measurements of the relaxation amplitudes are more questionable insofar as the frozen boundary conditions prevent some relaxations. A more adequate treatment of this question would require developments that are beyond the scope of this work. One may for instance think about the implementation of quasicontinuum simulation techniques that, by relaxing elastically the surrounding matrix, will provide flexible boundary conditions to the atomistic region. The picture that emerges from Fig. 6 is nevertheless interesting and sheds new light on the plastic deformation as it greatly simplifies representation of relaxation in glassy systems.
V Orientation effects
v.1 Loading direction
So far, relatively few studies have addressed the issue of the variation of the local yield stress as a function of loading orientation. This is due to the fact that most of the proposed local indicators are scalar quantities. Only work based on soft modes attempted to explore susceptibility to loading orientation by taking advantage of the vectorial aspect of vibrational eigenmodes Rottler et al. (2014); Smessaert and Rottler (2014). In order to address this issue, we use here another asset of our local method which naturally gives us access to this directional information.
To test the dependence on the direction of the mechanical loading, the quenched glasses are deformed by the same AQS protocol but following different orientations. In addition to the simple shear described in Sec. II.2, the systems are deformed in pure shear by applying strain increments and simple shear in the negative direction by applying deformation increments . These remote loadings correspond in the infinitesimal strain limit to shearing along and directions, respectively. Pure shear thus produces a diagonally-oriented shear. The negative simple shear corresponds to a laterally-oriented shear but in the opposite direction with respect to the positive simple shear remote loading employed so far. The positions of the first ten plastic rearrangements for the different loading directions are exemplified for a GQ glass in the top row of Fig. 7. In agreement with Ref. Gendelman et al. (2015), the location of plastic rearrangements show a strong dependence on the loading protocol. Most plastic events occur in different areas for different loading protocols. Only occasionally rearrangements will appear in the same location.
At the same time, local yield stresses are also calculated using the formula 4 with the corresponding loading directions . is thus equal to , and for positive simple shear, pure shear and negative simple shear, receptively. The maps of are shown in the top row of Fig. 7. A strong dependence on the loading orientation is observed. The rotation of the shear results in the appearance (disappearance) of soft (hard) zones. For example, the areas close to the 1st and 3rd plastic rearrangements for positive simple shear in Fig. 7a disappear in the case of negative simple shear in Fig. 7c. Conversely, soft areas appear as those close to the 5th and 7th rearrangements in Fig. 7c. As with the simple shear detailed above, pure shear and negative simple shear show an excellent correlation between the soft zones and the zones where the plastic rearrangements occur. The quantification of correlations through Eq. (5) and (6) as described in Sec. IV.2 is quantitatively similar (not shown here).
In order to highlight the discrete aspect of the variation of the local yield stress field, we compute the threshold contrast existing between two loading directions. This contrast is defined locally as the ratio between their difference and their averages as
where and are two remote loading directions. Contrast maps are shown in the bottom row of Fig. 7. These maps feature the trends described qualitatively above. The change in loading angle clearly shows areas of marked contrasts as a function of the loading orientations considered. The change of loading direction “turns on” or “turns off” the soft areas which gives rise to large local contrasts.
We observe that the greater the difference between the angles , the greater the number and intensity of the contrasts. In order to quantify this trend, we compute the cross-correlation of the yield stress field as a function of the difference of the loading angles . The result is shown in Fig. 8. The trend observed confirms that the correlation of the yield stress field decreases rapidly with the loading angle, regardless of the quench protocol. We note, however, that the correlation is never completely zero, and is still significant even for the largest that corresponds to the correlation between a shearing direction and its opposite direction. We attribute this effect to the small correlation existing between stable (unstable) zones and their tendency to have large (small) slip barriers Rodney and Schuh (2009); Patinet et al. (2016). The decorrelation due to the directional aspect is nevertheless clearly the dominant effect.
This result shows that local stress thresholds are a very sensitive probe of the loading protocol insofar as it is possible to acurately predict the plastic activity as a function of the orientation of the load. Unlike Gendelman et al. (2015), we thus believe that these results are consistent with a plasticity-based view of shear transformation zones. Indeed, from our point of view, the dependence of the plastic activity upon the loading protocol does not rule out the existence of plastic deformation via discrete units encoded in the structure, and that these discrete units clearly preexist within the material prior to loading. Our results show rather that the plastic deformations of an amorphous solid, at least for the transient regime at small deformations, can be seen as a sequence of activation of discrete shear transformation zones having weak slip orientations.
The local information given by our method allows us to also study the fluctuations in the direction of plastic rearrangements around the loading direction . At first, we have verified that the distributions of thresholds do not depend on the angle . As expected, the glasses are isotropic on average in the as-quenched state. We then consider the angle minimizing Eq. (4) for a given , i.e. the weakest local direction for a given loading direction. The distributions of around , shown in figure 9a, are well described by a Gaussian function of standard deviation . The latter decreases slightly with the relaxation of the system.
We also examine the consequences of taking into account the different possibilities of rearrangement directions on the correlation between local yield stress and plastic activity. Three types of local indicators can be considered: the minimum over all directions, the threshold only along the loading direction and as previously defined in Eq. 4. The correlations calculated for these three quantities from the relation 5 are plotted in Fig. 10 as a function of the imposed deformation.
We observe that the minimum over all angles gives the best correlations only for the very first plastic rearrangements and then decreases rapidly with deformation. This indicator corresponds to isotropic excitation, such as fluctuations in thermal energy, and is therefore sensitive to small barriers. Conversely, shows a poorer correlation with the location of plastic activity for the first rearrangements as it misses the low thresholds which are slightly disoriented with respect to the loading direction of the system. On the other hand, the correlation is better for larger deformations. Finally, shows the best correlation for both small and large deformations. Due to the projection, it is sensitive to small thresholds while retaining information specialized for macroscopic loading direction for larger yield stresses. We see here the importance of having access to a directional quantity. The local yield stresses defined in this article are therefore a good compromise between simplicity (purely local and scalar) and performance that justifies our approach. Note that qualitatively similar results have been obtained as a function of the relaxation of the system or when the correlations are computed from Eq. (6) (not shown here).
Vi Length scale of the local probing zone
vi.1 Optimal size
We are interested here in the effects of the patch size (Region I) on which the local yield stresses are computed. is varied from to . The procedure is the same as described in Sec. III. The size of the grid on which is sampled is kept constant and equal to . We first investigate the correlations of the thresholds with the plastic activity using relation (5). To quantify the degree of correlation, three kinds of indicators are considered: the correlation of the first plastic rearrangements , the characteristic deformation on which the correlation decreases with imposed deformation and the average correlation over the investigated strain window . The variations of these three indicators as a function of are shown for the three quench protocols in Fig. 11.
The correlation of the first plastic rearrangement increases with the size . Indeed, the increase of the probing zone makes it possible to progressively integrate the elastic loading heterogeneities. The loading felt by the sheared zones converges with toward the effective loading produced by a remote loading, which makes it easier to identify the weak zones. For , we observe a marked drop of the correlation. For this small size, in addition to larger elastic loading heterogeneities, the frozen boundary conditions over-constrain the measurement of local shear stress thresholds.
The decorrelation strain is extracted from a fit of the curves with the expression: as shown in Fig. 5a. This deformation, corresponding to the characteristic deformation on which the glasses lose their memory of the quenched state, decreases with . Indeed, after the first plastic events, the use of a large loses information about the hard zones surrounding the softest ones. A small allows us, while still having good spatial resolution, to maintain a significant correlation for higher deformations, the hard zones being simply advected during plastic flow.
The last kind of correlation indicator is the average of computed as:
The upper bound of the interval of integration is chosen equal to the largest decorrelation strain , i.e. computed for the smallest . This is a global indicator that gathers information on the degree of correlation at the origin and during deformation as the glass loses its memory from the quench state. The results reported in Fig. 11 show an overall decrease of the average correlation with . This decrease is less marked between and . The maximum of average correlation is even found for for the quench protocols HTL and GQ.
These results show empirically that a patch size of is a good compromise in terms of correlation between and plastic activity. Calculating the stress thresholds over this scale allows one to precisely locate the first plastic events while preserving the spatial resolution and keeping the memory of the initial quenched state. The effect of quench protocols is qualitatively similar to our previous observations. A greater relaxation of the system results in both a greater correlation for the first plastic events as well as a larger characteristic decorrelation strain, resulting in a larger average correlation.
vi.2 Statistical size effects
We are interested here in the effect of the patch size on the slip barrier statistic. Several mechanisms such as mechanical and statistical size effects can be anticipated. Mechanical size effects correspond to elastic heterogeneities as well as the influence of frozen boundary conditions. Frozen boundaries affect the simulation in the following way: the closer an atom to the boundary, the more affine its displacement, thus deviating its trajectory with respect to non constrained simulations. Statistical size effects play a role insofar as the local yield stress is primarily controlled by the weakest zones in the patch since its amplitude is given by the smallest threshold contained in the patch. Maps of local computed for different are shown in Fig. 12 (top row) for the quench protocol GQ. We observe that the variation of modifies the global statistic of the thresholds. The distribution functions presented in Fig. 13 for the three quench protocols show that the increase of induces a significant shift of the distributions toward smallest values of .
Obviously, the maps obtained for the large can be explained by the spatial increase of the zones centered on weak sites. The softest areas tend to “invade” the glass as the radius of the area on which the threshold is computed increases. Hence, the statistical effect seems to be dominant. On the basis of this observation, we try to understand the variations of the distributions of the local yield stresses with . We choose to work from the observed distributions for a size . We make the simplifying hypothesis that all the thresholds of the grid points take the value of the smallest local minima located inside a disk of radius . For comparison, maps deduced by this procedure are given in Fig. 12 (bottom row). This purely geometric approach shows a remarkable agreement compared to the local yield stress maps calculated by actually varying .
This approach allows us to deduce the distribution of the yield stresses as a function of a given patch size from the distribution obtained for . The comparisons between the distributions computed for the three quench protocols for different and those estimated from our procedure reported in Fig. 13 show a satisfactory agreement. The variation of the distributions of is therefore dominated by statistical effects. The increase of plays the role of a low-pass filter for the thresholds, shifting their distributions toward smaller yield stress values. The agreement between the measured distributions and the deduced distributions is nevertheless slightly lower for the less relaxed systems and for the large values. We attribute this discrepancy to the larger elastic disorder and to the lower sensitivity of the soft zones due to narrower threshold distributions in these systems.
In this article, we describe a method for sampling local slip thresholds in model amorphous solids. A robust correlation is observed between the zones with small yield stresses and the locations where the plastic rearrangements occur. As expected, the more the state of the glass is relaxed, the more the barrier distributions shift towards the larger values, explaining the strengthening of glasses from their local stability.
This local method has been extended to measure the amplitude of the plastic relaxations. We show that the assumption of a characteristic mean plastic relaxation is reasonable, the relaxation amplitudes following exponential distributions. Interestingly, we have shown that the amplitude of the plastic relaxations increases on average with the yield stresses.
The effects of loading orientation have shown that the variation of the plastic activity with the direction of loading is well captured by the variation of the local yield stress field calculated using our method. Finally, the variation of the threshold statistics with the size of the probing zones can be reproduced with reasonable agreement on the basis of simple geometric arguments. These results reinforce the coherence of the amorphous plasticity modeling based on discrete flow defects that possess weak slip directions and which are encoded in the structure of the material.
The advantages of the method presented in this work are numerous. It allows to probe the local slip thresholds in a non-perturbative way over a well-defined length scale. Moreover, its generalization to other atomic systems does not seem to pose any particular difficulty since it is, in principle, transposable to all glassy solids. Finally, a last advantage of this method is its computational cost. While, for instance, normal mode analysis based-methods scale as the cube of the number of atoms, our method scales linearly with it. Furthermore, as it treats the different parts of the solid independently it is by construction suited for massively parallel simulations. It is therefore possible to handle extended systems.
The implementation of this local method opens up several promising perspectives. It will be interesting to compare quantitatively the predictive power of the plastic activity of this method with other recent works also providing robust indicators of plasticity Rodney and Schuh (2009, 2009); Fan et al. (2014, 2015, 2017); Ding et al. (2016); Cubuk et al. (2015); Schoenholz et al. (2016); Cubuk et al. (2016).
Future research could focus on the measurement of quantities on an atomic scale needed for coarse-grained approaches Hinkle et al. (2017). For instance, our method can provide the threshold statistics necessary to take into account the disorder in the mesoscopic models Bulatov and Argon (1994); Baret et al. (2002); Picard et al. (2002); Vandembroucq and Roux (2011); Talamali et al. (2012); Nicolas et al. (2014); Tyukodi et al. (2016); Nicolas et al. (2017) and could explicitly deal with the tensorial nature of the problem and the effect of the loading geometry Budrikis et al. (2017).
This work paves the way, for example, to study the correlation between local energy barriers and frequencies of thermally activated rearrangements simulated by molecular dynamics. Hence, an important question left for future work is to study the effect of thermomechanical history on the statistics of local yield stresses. Our method will allow us to test some of the many phenomenological hypotheses upon which continuum models are still based Rottler and Robbins (2005); Falk and Langer (1998); Sollich et al. (1997); Hébraud and Lequeux (1998) and thus significantly improve the multi-scale modeling of plasticity of amorphous solids.
Acknowledgements.M.L. and S.P. acknowledge the support of French National Research Agency through the JCJC project PAMPAS under grant ANR-17-CE30-0019-01. R.G.G. and A.H.G acknowledge the support of CNRS and French National Research Agency under grant ANR-16-CE30-0022-03. M.L.F acknowledges the support of the U.S. National Science Foundation under Grant Nos. DMR1408685/1409560.
- These symmetries are related to the rotation of the stress tensor and the elasticity tensor in the -plane through an angle about the center of the patch. Noting the new coordinates of a point after rotation, the shear stress along is equal to and thus . The shear modulus along is equal to and thus .
- I. Procaccia, The European Physical Journal Special Topics, 178, 81 (2009).
- J.-L. Barrat and A. Lemaître, Dynamical heterogeneities in glasses, colloids, and granular media, 150, 264 (2011).
- D. Rodney, A. Tanguy, and D. Vandembroucq, Modelling Simul. Mater. Sci. Eng., 19, 083001 (2011).
- A. Argon, Acta Metall., 27, 47 (1979).
- E. Lerner and I. Procaccia, Phys. Rev. E, 79, 066109 (2009).
- A. Zaccone, P. Schall, and E. M. Terentjev, Phys. Rev. B, 90, 140203(R) (2014).
- M. Laurati, P. Maßhoff, K. Mutch, S. Egelhaaf, and A. Zaccone, Phys. Rev. Lett., 118, 018002 (2017).
- J. M. van Doorn, J. Bronkhorst, R. Higler, T. van de Laar, and J. Sprakel, Phys. Rev. Lett., 118, 188001 (2017).
- J. Chattoraj and A. Lemaître, Phys. Rev. Lett., 111, 066001 (2013).
- A. Lemaître, Phys. Rev. Lett., 113, 245702 (2014).
- M. L. Falk and J. S. Langer, Phys. Rev. E, 57, 7192 (1998).
- F. Spaepen, Acta Metall., 25, 407 (1977).
- R. L. Jack, A. J. Dunleavy, and C. P. Royall, Phys. Rev. Lett., 113, 095703 (2014).
- Y. Shi and M. L. Falk, Acta Mater., 55, 4317 (2007).
- M. Tsamados, A. Tanguy, F. Léonforte, and J.-L. Barrat, The European Physical Journal E, 26, 283 (2008).
- M. Tsamados, A. Tanguy, C. Goldenberg, and J.-L. Barrat, Phys. Rev. E, 80, 026112 (2009).
- A. Widmer-Cooper, H. Perry, P. Harrowell, and D. R. Reichman, Nat. Phys., 4, 711 (2008).
- A. Widmer-Cooper, H. Perry, P. Harrowell, and D. R. Reichman, J. Chem. Phys., 131, 194508 (2009).
- A. Tanguy, B. Mantisi, and M. Tsamados, EPL (Europhysics Letters), 90, 16004 (2010).
- A. Ghosh, V. Chikkadi, P. Schall, and D. Bonn, Phys. Rev. Lett., 107, 188303 (2011).
- M. L. Manning and A. J. Liu, Phys. Rev. Lett., 107, 108302 (2011).
- K. Chen, M. L. Manning, P. J. Yunker, W. G. Ellenbroek, Z. Zhang, A. J. Liu, and A. G. Yodh, Phys. Rev. Lett., 107, 108301 (2011).
- M. Mosayebi, P. Ilg, A. Widmer-Cooper, and E. Del Gado, Phys. Rev. Lett., 112, 105503 (2014).
- J. Ding, S. Patinet, M. L. Falk, Y. Cheng, and E. Ma, Proc. Natl. Acad. Sci. U.S.A., 111, 14052 (2014).
- S. Schoenholz, A. Liu, R. Riggleman, and J. Rottler, Phys. Rev. X, 4, 031014 (2014).
- S. Patinet, D. Vandembroucq, and M. L. Falk, Phys. Rev. Lett., 117, 045501 (2016).
- J. Ding, Y.-Q. Cheng, H. Sheng, M. Asta, R. O. Ritchie, and E. Ma, Nat. Comm., 7, 13733 (2016).
- E. Lerner, Phys. Rev. E, 93, 053004 (2016).
- J. Zylberg, E. Lerner, Y. Bar-Sinai, and E. Bouchbinder, PNAS, 114, 7289 (2017).
- E. Cubuk, S. Schoenholz, J. Rieser, B. Malone, J. Rottler, D. Durian, E. Kaxiras, and A. Liu, Phys. Rev. Lett., 114, 108001 (2015).
- E. D. Cubuk, S. S. Schoenholz, E. Kaxiras, and A. J. Liu, The Journal of Physical Chemistry B, 120, 6139 (2016).
- S. S. Schoenholz, E. D. Cubuk, D. M. Sussman, E. Kaxiras, and A. J. Liu, Nature Physics, 12, 469 (2016).
- J. Rottler, S. S. Schoenholz, and A. J. Liu, Phys. Rev. E, 89, 042304 (2014).
- A. Smessaert and J. Rottler, Soft Matter, 10, 8533 (2014).
- S. Plimpton, J. Comp. Phys., 117, 1 (1995).
- F. Lançon, L. Billard, and A. Chamberod, J. Phys. F: Met. Phys., 14, 579 (1984).
- M. Widom, K. J. Strandburg, and R. H. Swendsen, Phys. Rev. Lett., 58, 706 (1987).
- Y. Shi and M. L. Falk, Appl. Phys. Lett., 86, 011914 (2005a).
- Y. Shi and M. Falk, Phys. Rev. Lett., 95, 095502 (2005b).
- Y. Shi, M. Katz, H. Li, and M. Falk, Phys. Rev. Lett., 98, 185505 (2007).
- S. Nosé, J. Chem. Phys., 81, 511 (1984).
- W. G. Hoover, Phys. Rev. A, 31, 1695 (1985).
- S. Sastry, P. G. Debenedetti, and F. H. Stillinger, Nature, 393, 554 (1998).
- P. G. Debenedetti and F. H. Stillinger, Nature, 410, 259 (2001).
- D. L. Malandro and D. J. Lacks, Phys. Rev. Lett., 81, 5576 (1998).
- D. L. Malandro and D. J. Lacks, J. Chem. Phys., 110, 4593 (1999).
- C. E. Maloney and A. Lemaître, Phys. Rev. Lett., 93, 195501 (2004a).
- C. E. Maloney and A. Lemaître, Phys. Rev. E, 74, 016118 (2006).
- A. R. Hinkle, C. H. Rycroft, M. D. Shields, and M. L. Falk, Phys. Rev. E, 95 (2017).
- J. A. Zimmerman, D. J. Bammann, and H. Gao, Inter. J. Sol. Struct., 46, 238 (2009).
- I. Goldhirsch and C. Goldenberg, Euro. Phys. J. E, 9, 245 (2002).
- C. Goldenberg, A. Tanguy, and J.-L. Barrat, EPL (Europhysics Letters), 80, 16003 (2007).
- H. Mizuno, S. Mossa, and J.-L. Barrat, Phys. Rev. E, 87, 042306 (2013).
- P. Sollich, in CECAM Workshop (ACAM, Dublin, Ireland, 2011).
- F. Puosi, J. Olivier, and K. Martens, Soft Matter, 11, 7639 (2015).
- J. H. Irving and J. G. Kirkwood, J. Chem. Phys., 18, 817 (1950).
- These symmetries are related to the rotation of the stress tensor and the elasticity tensor in the -plane through an angle about the center of the patch. Noting the new coordinates of a point after rotation, the shear stress along is equal to and thus . The shear modulus along is equal to and thus .
- Y. Shi and M. Falk, Scrip. Mater., 54, 381 (2006).
- J. Lin and M. Wyart, Phys. Rev. X, 6, 011005 (2016).
- S. Karmakar, E. Lerner, and I. Procaccia, Phys. Rev. E, 82, 055103 (2010).
- H. G. E. Hentschel, P. K. Jaiswal, I. Procaccia, and S. Sastry, Phys. Rev. E, 92, 062302 (2015).
- J. Lin, A. Saade, E. Lerner, A. Rosso, and M. Wyart, Europhys. Lett., 105, 26003 (2014).
- J. Lin, T. Gueudré, A. Rosso, and M. Wyart, Phys. Rev. Lett., 115, 168001 (2015).
- C. E. Maloney and A. Lemaître, Phys. Rev. Lett., 93, 016001 (2004b).
- R. Dasgupta, H. G. E. Hentschel, and I. Procaccia, Phys. Rev. Lett., 109, 255502 (2012).
- T. Albaret, A. Tanguy, F. Boioli, and D. Rodney, Phys. Rev. E, 93, 053002 (2016).
- F. Boioli, T. Albaret, and D. Rodney, Phys. Rev. E, 95, 033005 (2017).
- Y. Fan, T. Iwashita, and T. Egami, Nat. Comm., 8, 15417 (2017).
- O. Gendelman, P. K. Jaiswal, I. Procaccia, B. Sen Gupta, and J. Zylberg, EPL (Europhysics Letters), 109, 16002 (2015).
- D. Rodney and C. A. Schuh, Phys. Rev. B, 80, 184203 (2009a).
- D. Rodney and C. A. Schuh, Phys. Rev. Lett., 102, 235503 (2009b).
- Y. Fan, T. Iwashita, and T. Egami, Nat. Comm., 5, 5083 (2014).
- Y. Fan, T. Iwashita, and T. Egami, Phys. Rev. Lett., 115, 045501 (2015).
- V. V. Bulatov and A. S. Argon, Modell. Simul. Mater. Sci. Eng., 2, 167 (1994).
- J.-C. Baret, D. Vandembroucq, and S. Roux, Phys. Rev. Lett., 89, 195506 (2002).
- G. Picard, A. Ajdari, L. Bocquet, and F. Lequeux, Phys. Rev. E, 66, 051501 (2002).
- D. Vandembroucq and S. Roux, Phys. Rev. B, 84, 134210 (2011).
- M. Talamali, V. Petäjä, D. Vandembroucq, and S. Roux, C.R. Mécanique, 340, 275 (2012).
- A. Nicolas, K. Martens, L. Bocquet, and J.-L. Barrat, Soft Matter, 10, 4648 (2014).
- B. Tyukodi, S. Patinet, S. Roux, and D. Vandembroucq, Phys. Rev. E, 93 (2016).
- A. Nicolas, E. E. Ferrero, K. Martens, and J.-L. Barrat, (2017).
- Z. Budrikis, D. F. Castellanos, S. Sandfeld, M. Zaiser, and S. Zapperi, Nat. Comm., 8, 15928 (2017).
- J. Rottler and M. O. Robbins, Phys. Rev. Lett., 95, 225504 (2005).
- P. Sollich, F. Lequeux, P. Hébraud, and M. E. Cates, Phys. Rev. Lett., 78, 2020 (1997).
- P. Hébraud and F. Lequeux, Phys. Rev. Lett., 81, 2934 (1998).