Learning local, quenched disorder in plasticity and other crackling noise phenomena

When far from equilibrium, many-body systems display behavior that strongly depends on the initial conditions. A characteristic such example is the phenomenon of plasticity of crystalline and amorphous materials that strongly depends on the material history. In plasticity modeling, the history is captured by a quenched, local and disordered flow stress distribution. While it is this disorder that causes avalanches that are commonly observed during nanoscale plastic deformation, the functional form and scaling properties have remained elusive. In this paper, a generic formalism is developed for deriving local disorder distributions from field-response (e.g. stress/strain) timeseries in models of crackling noise. We demonstrate the efficiency of the method in the hysteretic random-field Ising model and also, models of elastic interface depinning that have been used to model crystalline and amorphous plasticity. We show that the capacity to resolve the quenched disorder distribution improves with the temporal resolution and number of samples.


INTRODUCTION
When Gibbs proposed the famous measure for quantifying the probabilities of individual gas microstates, 1 the primary interest was on figuring out the microscopic principles that gave rise to well defined thermodynamic state variables, even though actual microstates kept changing through Brownian motion. Gibbs' original inquiry has re-emerged in the context of material deformation of crystals but also other mechanically deformable systems: the far-from-equilibrium response to mechanical deformation depends critically, especially at non-self-averaging small scales, on the initial configuration of defects. Here, we provide a way that can be used, in principle, to estimate the probability of quenched microstates as a function of material parameters (plastic strain, size, and hardening coefficients). 2 The evidence for quenched disorder in initial defect microstructures has been accumulated through observations of abrupt plastic events or material-crackling noise in a large variety of materials, such as FCC and BCC crystals, 3-6 amorphous solids 7 and also earthquake geological faults. 8,9 This crackling noise 10 has been commonly explained by random field models 11,12 or interface depinning ones, [13][14][15][16][17][18] where the major component is homogeneous solid elasticity, but also a spatially inhomogeneous and random distribution of local, quenched disorder (typically entering local flow stress information) 4,17,[19][20][21] and the allowed microstates are characterized by its stress and strain and minimize the elastic energy functional. The evidence of crackling noise has led to an extensive study of the local, statistical properties of abrupt events and their properties, such as their sizes, durations, average shapes, and their critical exponents. However, the major concern has been the fact that while homogeneous elastic properties are relatively straightforward to measure and test at virtually any scale, 22 the model distribution of local, quenched disorder is elusive, despite its commonly observed signature response of stochastic plastic bursts. [3][4][5]7 In this paper, we propose a feasible approach to "learn" the quenched disorder distributions directly from load-response timeseries: We argue that the full characteristics of the timeseries may unveil the information on the form of the quenched disorder distribution which is not available through typical temporally local observables (such as abrupt event size/duration). 23 While the major motivation of this work stems from plastic deformation, this method is generally applicable across crackling noise phenomena, defined through timeseries of an applied field (magnetic field, force, stress) and the associated response variable (magnetization, displacement, and strain). We label this method as time series-machine learning (TS-ML). TS-ML is clearly limited by the physical applicability and completeness of the utilized model. In its current form, for clarity purposes, TS-ML makes use of an unsupervised machine learning approach through principal component analysis (PCA) and k-means clustering. However, other unsupervised ML approaches could be successfully applied (see Supplementary Information for some  examples). For demonstration purposes, TS-ML is applied to a typical model of plasticity, 14 as well as to a general example of crackling noise. 24 Materials at small scales display size effects in material properties such as the strength and hardening coefficients. 21,[25][26][27][28][29][30] Stochastic yielding has been known to vanish with increasing system size. 31 However, its mechanism has been debated. 4,26 A plausible phenomenological scenario that explains the overall behavior is that the stochastic yield distribution at a representative volume has non-trivially large tails 32 at the nanoscale (that could either be a fat-tailed/large-kurtosis distribution or special distributions such as Weibull or Gumbel) and it gradually self-averages into the central limit behavior as the volume increases (cf. Fig. 1). While still a phenomenologically-driven hypothesis, the validity of this scenario would explain the self-consistent emergence of size effects and stochasticity in small-scale plastic deformation, in a material and sample independent manner. 4,5,33 In contrast to this picture, it is common in multiscale modeling of material deformation 22,34 to assume that the statistics of yield parameters in a representative material volume are dominated solely by their averages. Such a severe assumption is clearly not valid in a scenario where the yielding distribution has a large variance, or it is non-Gaussian in nature. Given that all current multiscale modeling approaches are fundamentally based on estimates of atomistic mechanisms that are active at the nanoscale, it is natural to consider and pursue the validation of such a scenario. For this reason, robust and swift methods are required for the derivation of local yield distributions as function of sample size, prior loading history and time.
In this paper, we propose a systematic way to solve the inverse problem of deriving the stochastic yield distributions from timeseries that may be the outcome of mechanical deformation or other crackling noise experiments. We use typical methods of unsupervised machine learning that naturally depend on physical modeling's completeness and the possible universality class. 4,23 The method does not aim to identify novel physical mechanisms of crackling noise; instead, it may provide a classification of parent quenched disorder distributions despite the existence of coexisting universal behavior. In the following, we perform simulations in two characteristic models of crackling noise, the T = 0 mean-field Random Field Ising model 24 (RFIM) and the elastic long-range dipolar interface depinning model(LRDIDM), 14,35 designed for crystalline or amorphous solids. In all these cases, we observe a capacity of the method to cluster model data despite the similarity of avalanche distributions, and then classify them according to the quenched yield distribution at a testing volume. Finally, the method may be also directly applied to actual experimental data on uniaxially compressed Ni micropillars, following prior work 6 (see Section 4 and associated Fig. S6 of the Supplementary  Information).

RESULTS
TS-ML assumes N samples that generally would include the discrete timeseries of the applied field f ij and the material response m ij with time index i ∈ [1, T] and sample index j ∈ [1, N]. The complete data matrix is built through a multiscale description of the timeseries that may be performed through sequential estimation of the sliding window time-varying moments of order n, at scale p and for sample j 36 where X could be either m or/and f. Through the construction of these moments, all fluctuations are captured (that could be equivalently unraveled through histogram distributions (cf. Fig.  1b)). For each sample j, a vector may be constructed that contains all moments up to a maximum resolution scale p max and max order n max , with total length M = n max p max . The parameters n max and p max are controlled by the resolution of the timeseries and can be estimated through: p max ' log 2 ðTÞ=2 and n max ∈ , 3,8 given that identified moments remain non-zero for the given resolution. Then, the effective data matrix D eff is built through these vectors and has dimensions N × M. D eff is used towards unsupervised machine learning through principal component analysis (PCA) and k-means clustering. Machine learning (ML) has had important success in various fields of science and engineering, but also materials science, for example accurate predictions of phase diagrams, crystal structures, and materials properties. 37 Unsupervised learning through using PCA may work independentently of the input data types, as soon as the data cluster into distinct, distinguishable spaces that can be tracked by a clustering technique such as k-means. 38 PCA reduces the dimensionality of a dataset 39 by identifying subspaces that demonstrate characteristic data variation. The success of the method depends on the capacity of data that originate in different quenched disorder distributions, to be spatially separated into distinct clusters in PCA space. The identification and quantification of the clusters is done through the k-means method, which finds k clusters that minimize the pairwise squared deviations of points in the same cluster. It is worth noting that our utilized ML is just a general way to statistically distinguish different signals, where the differences are caused by quenched disorder. However, other ML methods, may be similarly applicable (see Supplementary Information, Section 1 and Fig. S1), or may be expected to perform much better, such as deep learning methods. 40,41 These approaches will be studied in detail elsewhere.
The models TS-ML is applied on two characteristic models of crackling noise: a the elastic long-range dipolar interface-depinning model 13 (LRDIDM) and b) the hysteretic Random-Field Ising Model (RFIM). 24 The interactions in both models can be short or long ranged, altering the universality class. The purpose of this paper is to demonstrate examples, so some limiting behaviors are considered: mean-field interactions in the LRDIDM 42 and a complex, anisotropic and long-range kernel for the RFIM. 14 The method is identically applicable for quite general interaction kernels in both models. In the LRDIDM, the Talamali et al.'s approach 14 in two dimensions and assume that plastic deformation in disordered solids is modeled by the xx-component of the strain tensor ϵ ðpÞ ϵ ðpÞ xx ¼ Àϵ ðpÞ yy . The interaction due to local slip is the stress generated by local deformations of a random medium, 43 which takes the formF int ðk; ωÞ = ðÀcosð4ωÞ À 1Þεðk; ωÞ where k, ω are the polar coordinates in Fourier space 44 andε;F int are the transforms of strain and interaction respectively. We initialize the L × L system with ϵ(x) = 0 ∀ x and stress thresholds f p (x) are taken from a uniform distribution [0, 1). The external field F is increased  Fig. 1 Possible scenario of size effects and stochasticity in uniaxial compression of crystalline nano and micro pillars. a Uniaxial stress-strain curves have been shown to be increasingly stochastic -with clear bursts-and with higher apparent strength, as the sample/probe volume decreases, in the range where the probed volume has effective diameter 0.5, 10, or 100 μm. [3][4][5]7 The yield stress in such samples ranges from 50 to 500 MPa. (inset): Avalanche bursts, quantified through their strain magnitude S, have been shown to follow power-law distributions with a cutoff that decreases with sample volume. 3,5 b Increased strength and stochasticity at nanopillars could possibly originate into a nanoscale quenched yield distribution with non-trivial wide form, where the "most probable" yield is displaced from the bulk yield point. This distribution should evolve into a normal distribution as sample volume increases, according to the central limit theorem. The very existence of quenched disorder manifests in stochastic events that are common to describe through avalanches. 4 The green arrows display the direction of increasing representative volume that is being deformed Learning local, quenched disorder in plasticity S Papanikolaou adiabatically, and at each time-step, for all x that the total local force f(x) = F int (x) + F − f p (x) becomes positive, there is a strain increase dϵ(x) randomly picked from a uniform [0, 1) distribution. The external stress is decreased by k/L 2 at each time-step to cutoff avalanches-this cutoff mechanism resembles typical machine response at slow nominal straining. 14,45 In the RFIM, the approach follows the basic algorithm, 46,47 for which the zero temperature H the applied field, J the homogeneneous mean-field coupling, and h i the quenched random field that follows the distribution ρ(h) (to be specified). The local spin s i gives rise to the response m ¼ P N i¼1 s i . Starting with s i = −1∀i, H is increased until s i = +1∀i. In both models, it is possible to define abrupt events in terms of the size of response S (strain in the LRDIDM, magnetization per site in RFIM) during an event, as well as other quantities such as duration and energyrelease. 4 The distribution of S has been studied extensively 13,14,24 and it is known to be described, in the pinned regime, by P(S)S −τ P(S/S 0 ) where τ is~1.3 44 for the LRDIDM and 1.5 for RFIM. 24 The cut-off function P(x) typically resembles an exponential function and S 0 is a nominal maximum event size. In the depinned regime, there are system-spanning events, and therefore P(S) may include an additive component P inf that scales with the system size. In this paper, we consider various behavioral regimes.
Quenched disorder distributions Various possible random distributions are considered in both models. We consider: (a) normal distributions with mean μ and variance σ, (b) flat in a specified range, c)Weibull with probability density function g(f) = α λ ðf =λÞ αÀ1 e Àðf =λÞ α where δf~N −1/α defines the range of the distribution. Here, α is the Weibull exponent 48 , d) Gumbel distributions with g(f) = e Àðf ÀμÞ=β e Àe Àðf ÀμÞ=β =β with mean μ + 0.57721β and variance π 2 β 2 /6. All these distributions are very typical in materials science, depending on the character and origin of randomness. 49 First, in order to demonstrate the applicability of TS-ML method in realistic plastic yield problems, we consider the behavior of the LRDIDM. As shown in Fig. 2, the input random distributions that range from Weibull, normal, and flat (cf. Fig. 2c), give a range of nonlinear stress-strain curves (cf. Fig. 2a), consistent with the commonly observed behavior in crystalline and amorphous systems. 4 However, abrupt event distributions in terms of the strain burst sizes, are characteristically independent of the imposed distribution (cf. Fig. 2b). This result is consistent with previous studies in various avalanche models. 50 The presence of the abrupt event distribution validates TS-ML and facilitates a self-consistent physical picture. In order for TS-ML to be applicable (since high-order temporal moments need to be non-trivial), the presence of noise is required in the behavior. If such uncertainty was absent, then a simple constitutive formula Learning local, quenched disorder in plasticity S Papanikolaou would suffice to perform parameter fitting of the behavior. The very presence of abrupt events in the behavior validates the application of TS-ML, since it implies that the overall behavior is driven by noise-thus, a noise distribution should be identified. In crystal plasticity, the presence of such abrupt events further validate their origin into a quench disorder distribution, therefore TS-ML may be attempted on reasonable physical grounds. In such a noisy case (like, for example, in the experimentally testable case of uniaxially compressed Ni micropillars, 6 studied in Section 4 of the Supplementary Information), it makes no physical or practical sense to apply parameter fitting for parameters such as work hardening or yielding-these features are caused by a quenched disorder distribution that can be identified using TS-ML. In a natural way, work hardening is caused by the form of the quenched disorder distribution, which can be directly associated to the distribution of e.g., dislocation sources. If no such abrupt events are observed, TS-ML should be regarded as unphysical/ unnnatural, since there would be no motivation for a quenched disorder distribution in material modeling. In this sense, the selfconsistency of the physical picture may be tested. The application of TS-ML on 100 samples (20 each random distribution) is shown in Fig. 2, justifying the applicability of the method. Then, the RFIM is discussed, which does not include any spatial resolution. In this model, as shown in Fig. 3, 20 samples of seven different distributions are produced, (see Fig. 3c), totaling 140 samples. Then, the type of imposed randomness naturally influence the form of the M − H curves (cf. Fig. 3a), however the avalanche distributions (cf. Fig. 3b) demonstrate either a pre-depinning power-law behavior, or a spanning behavior, depending on the ratio of the disorder variance and the mean-field coupling strength. Then, all these samples are clustered and classified using TS-ML, and a very clear clustering effect that captures completely the original quenched threshold distribution is seen (cf. Fig. 3a).
The behavior of the RFIM is analyzed in further detail, by considering the limits of the applicability of TS-ML (for the method validation, see also Supplementary Information, Sections 2 and 3 and associated Figs. S3-S5). As shown in Fig. 4, the performance of the method is influenced as the number of samples S, the timeseries resolution R, and the number of degrees of freedom in the crackling system N. Three relevant observables are considered: 1. The "inertia" of k-means clustering D h i, which amounts to the sum of distances from the center of the assigned cluster; 2. The F β score with β = 2 which amounts to the success of the classification approach (k-means clustering) to identify correctly which samples are in the clusters; 51 The definition of the F 2 score amounts to a weighted average of precision and recall: can be seen, D h i increases with resolution R (cf. Fig. 4a) and number of samples S (cf. Fig. 4b), which amounts to naturally "larger" clusters in the PCA space. However, F 2 decreases with both R, S (cf. Fig. 4a,b), signifying that the larger clusters are also mutually well displaced, keeping the reliability of the method robust. Moreover, D h i and F 2 decrease drastically with the number of system degrees of freedom N. Finally, the variance is clearly distributed among the first three components, with more than 97% of the data variance.

DISCUSSION
In this paper, we suggested a possible scenario for size effects in crystal plasticity and also, we proposed an explicit method for correlating the shape of noisy load-response curves in crackling noise phenomena to the shape of the quenched disorder distribution, which is the natural cause of the noise in the response timeseries. In the case of crystalline deformation of nanocrystals, we argued that this method and approach can investigate in detail the uncertainty in the work hardening behavior during statistical sampling of the deformation, in order to provide an estimate of the originating quenched disorder distribution that may relate to various physical quantities (e.g., precipitate defects or pre-existing immobile dislocation forests). This approach can be equivalently thought as an unbiased "fitting" approach for work-hardening curves using stochastic distributions (instead of constitutive formulas that would only utilize average behaviors). In this paper, we argued that in various physical phenomena (especially ones that display crackling noise effects), the very nature of the problem enforces the use of stochastic distributions in order to be able to pursue multiscale modeling using physically accurate (statistically) microscopic information.
The TS-ML method provides a unique description of the quenched disorder distribution in the limit of a large number of test samples, in the sense that the PCA clusters separate distinctly and can be easily distinguished by a classifier (e.g. K-means). In the context of our theoretical models, we provide such evidence (cf. Fig. S3) in both the RFIM and the LRDIDM plasticity models (cf. Fig. 4). Also, it is natural to check the cluster centers (identified through the centers of the K-Means clusters) and how they progress with the increase in the number of samples. As it is discussed in the Supplementary Information (cf. Fig. S4) the average cluster position for all clusters in the RFIM case (cf. Fig. S2) is quite stable with the increase of the number of samples. Analogous results are found for the LRDIDM model. Moreover, uniqueness is critically dependent on the physical completeness of the utilized theoretical model. In that respect, the uniqueness of the TS-ML method can be established by the generality of its applicability: the LRDIDM and RFIM models can be generalized by adding short or long range interactions, valid for various applications. We tested TS-ML for a number of additional possibilities, such as the 3D-RFIM hysteretic model, 47 the model of ref. 35 for continuum dislocation plasticity, and also discrete dislocation dynamics models of crystal plasticity. 52 TS-ML appears to be generally applicable towards identifying quenched disorder distributions, appropriately defined for each model of reference. Finally, in continuum plasticity models the threshold distribution Learning local, quenched disorder in plasticity S Papanikolaou edge is critical for the yield behavior during avalanches. Thus, if two distributions have large overlap near the edge, then the limit of number of samples needed to distinguish different quenched disorder behaviors is prohibitively large (see SI's Section 3 and Fig.  S5 for a characteristic example, where overlapping Gaussian distributions lead to corresponding distinguishable but highly overlapping PCA maps).
The validation of the method can be performed by separating data sets into "experiment" and "simulations", with the experimental data sets being the ones that will be assessed, while all the rest contributing to the classification scheme. As it is shown in Fig.  S3, we separated statistical samples into two equal subsets for both the RFIM and the LRDIDM plasticity model, ("experimental"/ testing and "simulations"/training), and then the experimental data sets are projected on the trained PCA components. The testing data set successfully correlates in a percentage of around 90% for this size of the data set (which includes 50-100 samples per quenched noise distribution). However, we acknowledge that the method requires further validation, especially in the experimental front: In the current work, we showed how the method may be practically applied to actual experimental data on uniaxially compressed Ni micropillars, following prior work 6 (see Supplementary Information, Section 4 and Fig. S6). Through this investigation, a quenched disorder distribution with Weibull statistics and α ' 2) was shown to be consistent with various theory-based suggestions, originating on Weibull statistics of dislocation sources (single-arm etc.) (see for example ref. 53 ). Overall, the method is a concrete, unsupervised classifier of noisy stress-strain curves (and possibly other timeseries) that takes into account all moments of the noise, beyond the average behavior.
Through this work, it has become clear that it is viable to distinguish random threshold distributions from stochastic fieldresponse timeseries in typical crackling noise models. TS-ML can be efficiently implemented in order to distinguish quenched stochastic yield distributions in plastic deformation, that may originate in nanoscale experimental data. For this reason, high throughput experiments are required in order to efficiently probe the uncertainty in well defined deformed volumes. These distributions help define the material properties, if a generalized definition of a "material" is used where processing and deformation history is considered as a defining characteristic. In the future, the target would be to classify these distributions for various experimental cases in the processing mechanics of various materials and also various prior loading histories, aiming at producing a library of stochastic yield distributions that can be implemented in multiscale mechanics models, 22,34 in a similar way that interatomic potentials libraries 54 exist for molecular dynamics models. 55 The crucial importance of such a physical picture is the transition to predictive models that are based on the selfconsistent and intrinsically out-of-equilibrium statistical mechanics of crystal and amorphous plasticity.

METHODS
For the plasticity model, a custom Python code is used, solving the model on a 256 × 256 grid. The samples were loaded on well equilibrated strain configurations, that were generated by repeated loading/unloading of the samples. k is selected as 0.01 for these simulations. Typical cellular automaton rules are used, as described in the text. For the random-field model, the coupling J is chosen to be unity and the number of spins is chosen to be 512. For applying PCA, the sklearn library of Python is used: Non-linear PCA is accomplished by the application of Singular Value Decomposition, where the data matrix may be decomposed to a diagonal matrix of singular values S and left/right singular vectors V/U, D eff = V T SU. The V vectors that correspond to the largest singular values, capture the most characteristic temporal behavior. For solving the k-means problem, Lloyds algorithm 56 is used; the average complexity is given by O(knT), were n is the number of samples and T is the number of iteration, while the worst case complexity is given by O(N k+2/M ). 57 For more details, see also the Supplementary Information. The computational codes for the application of TS-ML on generic stress-strain or other timeseries are available upon e-mail request at stefanos.papanikolaou@mail.wvu.edu.

Data availability
The data that support the findings of this study are available from the corresponding author upon reasonable request.