Machine learning depinning of dislocation pileups

One of the key ideas of many disciplines including materials science in particular is that the structure and properties of materials tend to be related.^1,2 For instance, much of metallurgy is about tuning the materials’ microstructure to design materials with desired mechanical properties.² Typically, the relation between, say, the precipitate content of a precipitation-hardened alloy³ and its yield strength is understood as an average property, and indeed, in the case of macroscopic bulk samples, the sample-to-sample variations in the yield stress of identically prepared samples tend to be small. However, this changes when dealing with samples with sizes down to the micrometer range and below: Recent micropillar compression experiments have revealed the fluctuating, irregular nature of small-scale crystal plasticity,^4,5 originating from the collective, critical-like dynamics of interacting dislocations mediating the deformation process. The fact that the plasticity of small crystals proceeds via an irregular sequence of strain bursts with a broad size distribution implies also that sample-to-sample variations of the plastic response or the stress–strain curves might be considerable⁶ even if the samples have been prepared and the deformation experiments performed using the same protocol.

This raises the important general question of what can be said about the relation between the initial structure and mechanical properties for individual samples of small crystals, when random variations in the initial (micro)structure become important. One way of framing the issue is to consider what we refer to as deformation predictability: Given a small crystalline sample with a specific arrangement of pre-existing dislocations and the possibility of some other defects interfering with dislocation motion, how well can one use such information to predict the plastic response of that sample. A key challenge here is the high dimensionality of any reasonable description of the disordered initial state, combined with the possibility of non-linearities in the mapping from the initial state properties to the ensuing response. This likely implies that simple empirical laws relating sample properties to its plastic response may be difficult to formulate. Challenges of this type are a key factor behind the recent emergence of machine learning (ML) algorithms as an important part of the toolbox of scientists in a wide range of fields including also physics and materials science.^7,8 ML has been demonstrated to be an efficient approach to address a wide spectrum of problems including materials’ property prediction, discovery of novel materials, etc.^9–15 These developments have led to the emergence of a novel research field referred to as materials informatics,¹⁶ where informatics methods—including ML—are used to determine material properties that are hard to measure or compute using traditional methods.

In general, many supervised ML algorithms including neural networks are capable of learning non-linear mappings from a high-dimensional feature vector to a desired output, and hence, these methods should be useful tools when addressing the question of deformation predictability. An important opening in this direction was recently achieved by Salmenjoki et al.,¹⁷ who applied ML to study the deformation predictability in the case of simple two-dimensional discrete dislocation dynamics (DDD) simulations. The key idea of Ref. 17 was that exploiting ML algorithms provides a useful set of methods to quantify the predictability of complex systems such as plastically deforming crystals. This was achieved by training ML models including neural networks and support vector machines to infer the mapping from the initial dislocation microstructure to the response of the system to applied stresses, characterized by the stress–strain curve. Following Ref. 17, ML has been further applied to predicting stress–strain curves of plastically deforming crystals¹⁸ and learning the interaction kernel of dislocations.¹⁹ Recently, ML has also been applied to the closely related problems of, e.g., predicting the local yielding dynamics of dry foams,²⁰ the creep failure time of disordered materials,²¹ and the occurrence times of “laboratory earthquakes.”²² 

By its nature, the problem of predicting the plastic deformation process of crystalline samples depends on details such as whether the crystal only contains pre-existing glissile dislocations (this was the case in Ref. 17) or if other defect populations interacting with the dislocations are present. The latter could include solute atoms, precipitates, voids, or even grain boundaries in the case of polycrystals. If present, a description of these additional defects needs to be included in the initial state of the system used as input for the predictive ML models. The presence of such static defects within the crystal may also change the response of the crystal to applied stresses and, thus, the process to be predicted in a fundamental manner. It has been shown that the deformation dynamics of “pure” dislocation systems are governed by dislocation jamming,^23,24 resulting in glassy dislocation dynamics characterized by “extended criticality,” with the cutoff scale of the size distribution of dislocation avalanches diverging with the system size at any stress level.^25–27 This can be contrasted with systems where significant pinning of dislocations due to other defects within the crystal may instead induce a depinning phase transition of the dislocation assembly, resulting in critical-like dislocation dynamics only in the immediate proximity of the critical point (stress) of the depinning transition.^28,29

In this paper, we study deformation predictability in a one-dimensional periodic model of a dislocation pileup, interacting with a quenched (frozen) random pinning landscape. The model is perhaps the simplest possible system including both interacting mobile dislocations and a quenched random pinning field interfering with dislocation motion, and therefore, it serves as a useful playground to explore the ideas of ML-based deformation predictability discussed above. It is known to exhibit a depinning phase transition at a critical flow stress value σ_flow, separating pinned and moving phases of the dislocation system.³⁰ While the general topic of applying ML in the context of disordered materials has recently gained some momentum,^{17,18,20–22,31} to our knowledge, no study addressing the question of deformation predictability in a dislocation system interacting with quenched disorder and exhibiting a depinning phase transition has been performed before. To this end, we generate a large database of stress–strain curves, each corresponding to a unique randomly generated initial microstructure, by simulating the model with a quasistatically increasing applied stress for different system sizes (or different numbers of dislocations, N). First, the statistical properties of the stress–strain curves and the strain bursts are analyzed, followed by training and testing of various predictive models ranging from linear regression to convolutional neural networks (CNNs; considering these different ML models allows us to assess the possible model dependence of deformation predictability) to establish mappings from the initial random microstructure (defined by both the static pinning landscape and the initial configuration of the dislocations) to the response of the system to applied stresses (i.e., the stress–strain curve). Notice that the deformation predictability as defined above is in general dependent on the ML model considered, but by studying the problem using multiple different ML models allows one to use the results from the one with the best predictive power as measures of the lower limit of the more general, ML model independent deformation predictability.

Our results show that the different predictive models employed are capable of learning the relation between the input and the stress–strain curves quite well, as measured by the correlation coefficient between the predicted and simulated stress values at a given strain. In particular, the sample flow stress (i.e., the sample-dependent finite size critical point of the depinning transition) can be predicted surprisingly well with the correlation coefficient reaching values as high as 0.89 for a regularized CNN. We also explore the predictability of individual strain bursts taking place during the deformation process. Critical avalanches are expected to be inherently unpredictable, and indeed, we find approximately zero predictability for most of the strain bursts belonging to the power-law scaling regime of their size distribution. Interestingly, we also find that strain bursts belonging to the cutoffs of the distributions as well as those taking place very close to the flow stress exhibit non-vanishing predictability.

The dislocation pileup model we study here is similar to the one in Refs. 30 and 32. It describes a system of N straight, parallel edge dislocations subject to an external shear stress σ gliding along the direction set by their Burgers vector (here, the x direction) within a given glide plane. By neglecting the roughness of the dislocation lines, the system may be described by a set of point dislocations with coordinates x_i representing the cross sections of the dislocation lines moving in one dimension. The dislocations interact repulsively with each other via long-range stress fields with the stress field magnitude inversely proportional to their mutual distance.³⁰ In addition, the dislocations are taken to interact via short-range forces with a set of N_p randomly positioned Gaussian pinning centers, playing here the role of quenched disorder. In real crystals, these pinning centers could be solute atoms or immobile forest dislocations threading the plane of the pileup. The overdamped equations of motion of the N point dislocations thus read

where χ is the effective viscosity, x_i is the position of the ith dislocation, μ is the shear modulus, b is the Burgers vector magnitude, and F is the pinning force landscape. We choose a pinning force consisting of Gaussian potential wells centered at pinning sites $x p$⁠, given by

where x_p is the position of the pth pinning point, E_p is a pinning energy scale, and s_p is the standard deviation of the Gaussian. For simplicity, we set χ = μ = b = 1 and $E p = s p 2 =0.5$⁠. The periodic boundary conditions within a system of size L are employed by performing an infinite sum over the stress fields of the periodic images of each dislocation,³⁰ 

The average spacing between dislocations is set to L/N = 16 and that between pinning sites is set to L/N_p = 2. The sum over pinning sites is truncated so that only pinning sites that are closer than a cutoff distance of 8 are included in the sum.

The simulations are carried out by first choosing the positions of the N_p pinning sites at random from a uniform distribution along the one-dimensional system of length L. N dislocations are then placed in the system with the constant spacing of L/N = 16; this equally spaced dislocation configuration would be the minimum energy state of the pileup in the absence of pinning sites. The dislocations are initially allowed to relax to a metastable configuration with σ = 0; see Fig. 1(a) for an example of a relaxed configuration with N = 4, showing also the pinning energy landscape. Then, σ is increased quasistatically from zero with the rate

where m_σ = 10⁻⁴ is the maximum stress rate, v_th = 2 · 10⁻⁴ is a threshold velocity value, $v ¯ ( t ) = ( 1 / N ) ∑ i = 1 N v i$ (with v_i = dx_i/dt) is the spatially averaged dislocation velocity, and r = 100 is a shape parameter. In practice, this amounts to a continuous, smooth approximation of the step function, with dσ/dt ≈ m_σ for $v ¯ < v t h$ and dσ/dt ≈ 0 for $v ¯ > v t h$⁠, such that σ is increased only in between strain bursts and kept constant during them. The strain ε is accumulated by dislocation motion,

The simulations are run until the dislocations start to flow indefinitely, i.e., the strain burst size diverges, which happens at a sample-dependent flow stress σ_flow. The stress–strain curves σ(ε) are recorded by storing σ(t) and ε(t) every 20 dimensionless time unit during the simulations. Example stress–strain curves are shown in Fig. 1(b) for three different system sizes. Notice the characteristic irregular staircase-like structure of the curves as well as the larger size of strain bursts for smaller systems.

The resulting stress–strain curves, each corresponding to a different realization of the random pinning landscape and relaxed initial dislocation positions, are stored in a database consisting of 10 000 stress–strain curves for each system size to be used as training and testing data for the predictive models. For the largest systems (N = 64, …, 512), we also generate a separate dataset of 1000 stress–strain curves stored at a finer time resolution of Δt = 2 in order to better detect individual strain bursts. This latter dataset is used to analyze the statistical properties of the strain bursts as well as to train and test ML models for predicting strain bursts.

In what follows, we will introduce the predictive models used in this study. These encompass both linear and non-linear models using hand-picked features as input as well as a CNN that requires less feature engineering as it uses a complete description of the system’s initial state as input. Using these different models allows us to assess the dependence of our results on any particular model.

Before training, each input feature (each input channel in the CNN case) is standardized by subtracting the mean and dividing by the standard deviation, and the dataset is divided randomly into a training set containing 80% of the data points and a testing set containing the remaining 20%. Predictability is measured by calculating the correlation coefficient between the predicted and target outputs of the testing set and averaging the result over five training instances. Predictive models are implemented using the Python libraries scikit-learn³³ for linear absolute shrinkage and selection operator (LASSO) and Keras³⁴ for simple neural networks and CNNs.

The prediction models work by minimizing the loss, which mainly consists of the mean squared error. In addition, most models are L¹-regularized, which adds a penalty term to the loss to prevent overfitting. L¹ regularization often improves the result and encourages the model to focus on the most useful/promising input features while discarding the unimportant, and as a byproduct, it allows us to determine which input features are used by the model.

A linear regression model performs an affine transformation, where input features are first linearly combined by multiplication with a (weight) matrix and then biased by adding a translation vector. If the model is optimized using L¹-regularization for the weight matrix, it is commonly called LASSO (linear absolute shrinkage and selection operator³⁵). In this study, this model is applied by using sklearn.linear_model.Lasso, which has a built-in optimizer and only requires the user to choose the regularization parameter α. We found that values of α in the range 10⁻³ ≲ α ≲ 10⁻² gave good results and chose α = 10⁻².

When two linear regression models are stacked in such a way that the intermediate result is activated by applying a non-linear activation function, such as the rectified linear unit (ReLU), the combined model becomes a simple neural network that can, to some extent, learn general non-linear mappings. An intermediate result of such a model is called a hidden layer, and elements of the intermediate result are referred to as hidden units. L¹-regularization is applied to each layer with a parameter α_n describing the strength of regularization at each layer.

Tested variants for stress–strain curve prediction differ by the number of hidden units, hidden layers, and regularization. Most have one hidden layer with 64 hidden units and use α₁ = 10⁻² for the first weight matrix and α_n = 10⁻⁵ for the rest. All models use ReLU-activations for the hidden layers. One variant has less regularization, meaning that α₁ = 10⁻³ instead of 10⁻². Theoretically, one hidden layer is enough for reproducing any continuous mapping given enough hidden units, but a model with three hidden layers is employed here as well for comparison to other studies.^17,36 The models are trained using the Adam optimizer³⁷ with learning rate 10⁻³ for 100 epochs.

The predictability of strain bursts is studied using a simple neural network having one hidden layer with 64 hidden units. The model output is in this case very big, consisting of a vectorized (one-dimensional) version of an avalanche map that is originally two-dimensional (see Sec. IV for details). We found that all parameters collapsed to 0 if α is too high, so a value of 10⁻⁵ is used for all weight matrices. A lower α parameter makes the model susceptible to overfitting. Therefore, models are trained only for 10 epochs but with a higher learning rate of 10⁻².

The LASSO and simple neural network models require hand-picked features for predicting stress–strain curves and strain bursts. It turns out that good, information-rich features for flow stress prediction are given by quantiles of distributions derived from the sample-specific pinning landscapes. Figure 2(a) shows the correlation of the flow stress with quantiles of the pinning energy E, the pinning force F, and $DF= d F d x$ and $D 2 F= d 2 F d x 2$⁠. In a non-interacting system, the flow stress would be controlled by the most negative pinning force F, which provides the strongest obstacle for dislocation motion. In an interacting system, it is possible for the dislocations to push on each other and overcome the strongest pinning force at a smaller applied stress. As a result, the correlation between the most negative value of F [given by the 0% quantile of F in Fig. 2(a)] and the flow stress is quite small. Interestingly, the 5% quantile of F exhibits a much larger negative correlation of ∼−0.8 with σ_flow, independent of the system size [Fig. 2(b)]. Other quantiles of the pinning force as well as of the other quantities listed in Fig. 2(a) also appear to contain useful information for flow stress prediction. We, therefore, give the entire quantile curve (sampled on 64 uniformly chosen points) of each field as input features.

In addition to these input features describing the pinning landscape, we also give input features that describe the initial relaxed configuration of dislocations ${ x j } j = 1 N$⁠. In particular, we use the values E(x_j), F(x_j), DF(x_j), and D²F(x_j) as well as the distances between successive dislocations, Dx_j = x_j+1 − x_j (with x_N+1 = x₁ + L). These features are given as sorted lists, so they can be interpreted as quantiles.

Finally, we choose a set of features given by local extrema of E and F, as well as the successive differences in energy extrema D(E extrema), also given as 64 uniformly spaced quantiles.

These models take discretized fields as input and have the ability to find out if the target output depends on the spatial variations within the input fields. The CNN structure used in this study is adapted for periodic input signals by applying a convolution-activation-pooling scheme that continues to down sample the input array down to a single spatial element. Initially, an input field has 1024 = 2¹⁰ spatial elements, so there are ten convolutional layers (each layer halves the number of spatial elements). This makes the result effectively independent of the particular choice of the origin of the input fields. Therefore, by contrast to more standard CNNs using fewer convolutional layers, this model’s final spatial element’s receptive field covers the whole system, and translational invariance is complete because the ten pooling operations leave no information left about the absolute locations of where the previous convolutional layers have activated. This is desirable in the case of periodic inputs, where circular shifting does not change the properties (such as the plastic response) of the system.

All convolutional layers employed here use a window size of 3, periodic boundary conditions, ReLU-activation, and non-overlapping max-pooling of size 2 (taking the maximum value within non-overlapping bins) and have 16 hidden units. The result from the final pooling layer is fed into a linear regression layer to obtain the model output. When predicting a single output feature from two input fields (channels), the CNN has 7185 learnable parameters.

The CNN is trained using the same Adam optimizer as for other neural network models, with a learning rate of 10⁻³, a learning rate decay of 10⁻⁵, and 50 training epochs. We found that even weak L¹-regularization caused the parameters to collapse to 0 in most cases, so we do not regularize these models. However, in the case of predicting flow stress only, a kernel regularization value of 5 · 10⁻⁴ for the convolutional layers is found to work.

The input to the CNNs consists of the pinning landscape fields E, F, DF, and D²F (although we found that giving one pinning field, such as E or F, suffices because the CNN can compute the higher derivatives by differentiation) as well as a field J defined as a sum of Gaussians centered on the relaxed dislocation positions x_i,

This field J serves to describe the dislocation positions in terms of an input field with a similar structure as the pinning landscape, which is independent of the spatial resolution [see Fig. 7(a) for examples of these input fields]. Notice that by contrast to the case of hand-picked features, these fields are given in order rather than sorted.

In this section, we first consider the statistical properties of our dataset of stress–strain curves. Then, we report our results for training the predictive models and employ them to explore the predictability of the deformation dynamics of the pileup model. Finally, we also discuss our results on the problem of predicting individual strain bursts.

We start by quantifying the system size dependent sample-to-sample variations in the plastic response of the system. Figures 3(a)–3(c) show the bivariate histograms of the stress–strain curves for three system sizes (with N = 16, 64, and 256, respectively), illustrating that the sample-to-sample variability in the plastic response is more pronounced in smaller systems. The same applies to the probability densities of the sample flow stresses σ_flow, as shown in Fig. 3(d): The distributions are quite broad for small systems and get narrower upon increasing the system size, presumably approaching a delta peak in the thermodynamic limit. The average magnitude of the σ_flow-values is in agreement with previous results.³⁰ 

As the stress–strain curves consist of a sequence of discrete strain bursts or dislocation avalanches Δε, we also consider their stress-dependent size distributions P(Δε, σ). Quite generally, for systems exhibiting a depinning phase transition, one expects them to obey the scaling form

where τ is the critical exponent of the avalanche size distribution, f(x) is a scaling function that obeys f(x) = const. for x ≪ 1 and f(x) → 0 for x ≫ 1, and Δε^*(σ) is a stress-dependent cutoff avalanche size. The left panels of Figs. 4(a)–4(d) show the bivariate histograms of the avalanche count as a function of σ and avalanche size Δε for different system sizes (with N = 64, 128, 256, and 512, respectively). These illustrate that, indeed, the maximum strain burst size increases upon approaching the flow stress from below, and the increase close to the flow stress is more clear-cut for larger systems. Notice also that these histograms show an excess of strain bursts around σ ≈ 0.4. This feature is likely related to the characteristic pinning force magnitude due to an isolated pinning point given by Eq. (2). The right panels of Figs. 4(a)–4(d) show the corresponding stress-resolved strain burst size distributions P(Δε, σ), considering three different stress bins below the flow stress for each system size. For the largest system sizes considered, the data appear to be consistent with a power-law part characterized by an exponent τ ≈ 1.0 [see the right panels of Figs. 4(c) and 4(d)]. This value is somewhat lower than the expectation of τ ≈ 1.25 based on the scaling of the elastic energy of the pileup in the Fourier space,³⁰ suggesting that the pileup model would be in the same universality class with depinning of contact lines³⁸ and planar cracks.^39,40 As we focus here on studying rather small systems, this difference is likely due to finite size effects.

Overall, the statistical properties of the plastic response of the pileup model show two key features relevant for deformation predictability: (i) there are sample-to-sample variations in the response (that we will try to predict using machine learning) and (ii) the distribution of strain bursts Δε displays the characteristics of a critical system, suggesting that individual strain bursts may be inherently unpredictable. In what follows, we will use our trained predictive models to explore how well (i) the sample stress–strain curves and (ii) individual strain bursts can be predicted.

Figure 5 shows the mean correlation coefficient between the predicted and actual stress values as a function of strain for N = 64, using linear regression (LASSO) and different variants of simple neural networks. First, we note that the flow stress prediction is only slightly better than the magnitude of the correlation with the 5% quantile of the pinning force landscape found above (∼0.8, see Fig. 2), with the correlation coefficient assuming values in the range of 0.825–0.84 depending on the model used (the large-strain plateau in Fig. 5). Second, the correlation clearly depends on strain, exhibiting a rather pronounced dip reaching values as low as 0.65 or so for intermediate strains, along with two additional, lesser dips for smaller strains, before reaching a value close to 1 as the strain goes to zero. The latter is true especially for LASSO, while the neural network models seem to struggle a bit in capturing the deformation dynamics for very small strains. The reduced predictability for intermediate strains could be related to the importance of inherently unpredictable dislocation avalanches in analogy with the results of Ref. 17.

It is interesting to analyze which features the predictive models end up using for prediction and how that might depend on strain. Since the models utilize heavy L¹-regularization, the analysis of the learned weights becomes straightforward. Features that contribute to predictability are allowed to have a high regularization penalty in exchange for decreasing loss, while unimportant features will not decrease loss sufficiently to justify a high regularization penalty. Therefore, if we were to omit features with a small L¹ penalty, the loss would not be affected and the model would presumably do just as well, so this analysis is closely related to studying feature ablations. Figure 6 shows how much penalty is generated by the weights corresponding to each feature category for the LASSO model. At the start, the pinning force derivative at the relaxed dislocation positions, DF(x_j), is the most active feature, but also the pinning energy at the initial dislocation positions, E(x_j), is used. This is reasonable as such quantities could be used to estimate the initial linear response of the dislocations. Interestingly, and perhaps somewhat counter-intuitively, the pinning energy at the initial dislocation positions becomes even more important for intermediate strains in addition to a larger diversity of features used. The flow stress prediction uses mainly the quantiles of the pinning force F as well as its local extrema points. This dependence of the important features on strain indicates that predicting the response of the system to infinitesimal stresses is of a different nature compared to predicting the flow stress.

We then proceed to employ CNNs as predictive models. Figure 7(b) shows the mean correlation as a function of strain as obtained using three different types of CNNs along with the LASSO result for reference. Here, the system size N = 16 is considered. First, we consider a CNN that predicts the entire stress–strain curve at once [a “generalist” CNN, shown as a red line in Fig. 7(b)]. This model actually underperforms the simple linear regression model for the majority of the strain values. This could be due to the fundamentally different nature of the problem for different strain values, as illustrated by Fig. 6, showing that different features are important for different strains. Hence, the “generalist” CNN might get exhausted and not be able to optimally learn the entire stress–strain curve at once.

Thus, we also consider a “specialist” CNN that is trained separately for each strain value to learn the corresponding value of stress. The result is shown as yellow circles in Fig. 7(b). Indeed, as compared to the generalist model, the specialist model is able to learn the deformation dynamics of the pileup model much better, especially for intermediate strains, where the generalist CNN exhibits a deep dip in the correlation coefficient, presumably due to the importance of hard to predict dislocation avalanches for those strains. It is interesting that the specialist model is able to reach a correlation coefficient exceeding 0.8 also in this regime.

Notably, neither of the models discussed above is able to predict the flow stress much better than by just considering the correlation with the 5% quantile of the pinning force F (see Fig. 2). We, therefore, try to further improve the specialist CNN by switching on L¹-regularization and using only the pinning force F as the input field (this was shown above to be the most informative quantity for flow stress prediction). The individual purple point at a large strain value in Fig. 7(b) is the resulting flow stress prediction with a correlation coefficient of 0.89, i.e., significantly better correlation than using any of the other models considered.

We also explore the possibility of predicting individual strain bursts along the stress–strain curves. First, we need to formulate a problem that is suitable for analysis using the predictive models at hand. To this end, we generalize the bivariate avalanche histograms (left panels of Fig. 4) to the case of individual samples. As such, these would be rather sparse sets of points in the space spanned by the avalanche size and stress. We, therefore, smooth out these maps by convolving with a Gaussian with a standard deviation of 0.03 units of stress in the vertical direction and 0.15 units of log₁₀(avalanche size) in the horizontal direction; both correspond to three times the bin size used in Figs. 4 and 8. The studied region is limited to bins where at least one avalanche appears within the dataset so that the Gaussian smoothing does not cause distribution boundaries to extend to impossible areas. Finally, the maps are standardized by removing the mean and dividing by the standard deviation at each bin. We then employ a simple neural network with 64 hidden units to find a mapping from the hand-picked input features to the smoothed sample-specific bivariate avalanche histograms.

Three different examples of smoothed avalanche maps are shown in the left column of Fig. 8 with the corresponding predicted avalanche maps in the right column. While the predictions are clearly not perfect, the algorithm is able to capture some key features of these avalanche maps. Figure 9 shows the correlation maps between the predicted and actual avalanche densities at four different system sizes, averaged over five training instances. These show a number of interesting features: (i) Large parts of the sample-specific avalanche histograms appear to be completely unpredictable, with the correlation coefficient assuming a value very close to zero (shown in black in Fig. 9). This is especially true for most of the regions corresponding to the power-law scaling regime of the avalanche size distribution in agreement with the idea that critical avalanches should be intrinsically hard to predict.⁴¹ (ii) Avalanches with the largest size for each stress value (i.e., those belonging to the cutoff of the stress-resolved avalanche size distribution) appear to exhibit some degree of predictability, although the related correlation coefficient values are somewhat lower than those found for the predictability of the stress–strain curves. This is in agreement with the idea that avalanches in the cutoff are not critical and hence somewhat predictable. (iii) Surprisingly, we also find the non-vanishing predictability of avalanches of any size occurring in the immediate proximity of the flow stress (bright horizontal segments in the upper parts of the panels of Fig. 9). This is surprising because those events are expected to be critical and hence should be unpredictable.

In order to shed some light on the origin of the key features of avalanche predictability discussed above, we finish by considering the correlation between the avalanche count and flow stress. The logic here is that since we have demonstrated above that the sample-specific flow stress can be predicted quite well, if the avalanche maps are correlated with the flow stress, these should be predictable as well. Figure 10 reveals the two stand-out features: (i) The avalanches taking place in the immediate proximity of the flow stress exhibit clear positive correlation with the flow stress and (ii) the largest avalanches for stresses smaller than the flow stress (i.e., avalanches belonging to the stress-dependent cutoff of the avalanche size distribution) are negatively correlated with the flow stress. This means that both a higher than average number of avalanches very close to the flow stress and a lower than average number of large avalanches at smaller stresses imply a higher than average flow stress value. The good level of predictability of the flow stress demonstrated above thus translates into the reasonable predictability of these features of the sample-specific bivariate avalanche histograms. Notice that the widths of the predictable bands in the avalanche maps (Fig. 9) appear to become thinner with increasing system size, so we cannot exclude the possibility that this avalanche predictability would be a finite size effect.

We finish by considering the effect of small perturbations of the pinning point positions on the flow stress value. Notice that these perturbations also affect important features such as the extreme quantiles of the pinning force landscape as large pinning force values originate from overlapping force fields of more than one pinning point. This, therefore, serves as a measure of the robustness of our results on deformation predictability in situations where the initial microstructure can be characterized with only a finite precision. To that end, we perturb the initial pinning point positions by random amounts drawn from a normal distribution with mean 0 and standard deviation δx. Figure 11 compares the resulting flow stress with the unperturbed one as a function of δx, by measuring the relative difference [Fig. 11(a)] and the correlation between the perturbed and unperturbed flow stress [Fig. 11(b)]. The results are calculated for the system of 64 dislocations using a set of 100 reference cases, which are compared with the corresponding perturbed systems with varying δx. Although slightly noisy, the overall shape of the curves is clearly visible.

Figure 11 shows that perturbing the pinning point positions has essentially no effect below a standard error of δx ∼ 10⁻², and a significant effect of perturbing the pinning point positions becomes visible only when the perturbation magnitude clearly exceeds that threshold value. This suggests that the dislocation pileup model and hence our results concerning deformation predictability in that model are essentially unchanged even if the pinning point positions (or, more generally, the features describing the pinning potential) are perturbed a little.

We have established that predictive models ranging from linear regression to CNNs can be used to predict the stress–strain curves of the pileup model with a good accuracy. While the different models employed give rise to somewhat different results, a practical conclusion is that the model giving the highest correlation coefficient for a given strain provides a lower limit for the deformation predictability of the pileup system. Notably the “specialist” CNN discussed above results in correlation coefficients exceeding 0.8 over the entire range of strains considered, and adding regularization further improves the correlation coefficient of the flow stress prediction to 0.89. These should be interpreted as lower limits of deformation predictability as we obviously cannot exclude the possibility that a hypothetical predictive model not considered here might outperform our models. For flow stress prediction, the 5% quantile of the pinning force landscape turned out to be a surprisingly information-rich feature, highlighting the fact that our method of using quantiles as features might be useful also in other contexts.

While the predictability of the stress–strain curves turns out to be quite good, it is not perfect. This is to be expected because the stress–strain curves consist of critical-like strain bursts that are expected to be inherently hard to predict. Indeed, our attempts to predict the sample-level bivariate avalanche count resulted in essentially zero predictability for most of the critical dislocation avalanches belonging to the power-law distributed part of the avalanche size distribution. On the other hand, these sample-level avalanche distributions also showed an interesting regularity: A larger than average flow stress value indicates a lower than average number of the largest avalanches for stresses below the flow stress and a larger than average number of avalanches in the immediate vicinity of the flow stress. It would be interesting to extend such an analysis from our stress-controlled simulations to strain-controlled ones. One should point out that because of the unpredictable nature of critical dislocation avalanches, the ML models cannot fully replace the integration of the equations of motion as they are unable to fully resolve such strain bursts. However, the fact that the overall shape of the stress–strain curve is still predictable to some degree suggests that the trained ML models could be used for simple classification tasks like determining if a given sample is “weak” or “strong” compared to the average, using significantly less computational resources than what is required for integrating the equations of motion.

An interesting issue that is also of practical relevance has to do with the accuracy of the input data used for prediction. Here, we have access to the precise microscopic information of the pinning force landscape, and indeed, the features we use have been constructed from such precise information that might not be available in real experimental systems. However, our robustness analysis of the pileup model suggests that extreme precision is not needed for good predictions, and the same conclusion can be made by remembering that our hand-picked features are constructed in a way that does not preserve the information of the precise locations of the pinning points. Thus, we expect our general methodology to be applicable also in situations where the initial microstructure can be characterized with a finite precision only.

To put our work in the context of the broader literature on applying ML to physics problems, let us point out some details: First, incorporating the inductive biases from different physical systems into machine learning architectures has recently received attention,^42–44 the key idea being to use ML to learn missing physical relations from data by constructing models that mimic the equations of motion and then applying them to interpolate or extrapolate the original data. One should note that our aim here is simply to find a mapping from features characterizing the initial microstructure to the plastic response of the sample, without explicitly incorporating information about the laws of physics in the predictive models. When constructing the CNN, we make use of the fact that the system we study is periodic, but otherwise it is up to the ML model(s) to learn the above-mentioned mapping without explicit input related to the physics of the system. A second relevant development in ML has to do with incorporating the domain knowledge in ML models.⁴⁵ In the context of our study, the feature engineering process that resulted in the features we use for LASSO and simple neural networks can be argued to incorporate some “domain knowledge,” as the features used reflect, at least to some extent, our physical intuition as to what the important features might be. Third, our study dealing with a rather simple dislocation pileup model allowing us to calculate an extensive training database with relative ease could perhaps serve as a “sister system” for related systems such as more realistic/complicated models as well as experiments where data are scarce, in analogy to the approach of Hoffmann et al.³¹ For instance, in systems with a higher dimensionality than the one considered here, such an approach might work well especially for small strains, where the dislocation dynamics is governed by the linear response independent of the system’s spatial dimensionality. For larger strains, the formation of tangled dislocation structures resulting in strain hardening might, however, complicate the situation, as such features are not captured by our one-dimensional pileup model.

We finish by noting that the pileup model studied here has the advantage that the dislocation system will sample the entire one-dimensional random energy landscape, and hence, extracting relevant descriptors of the sample-specific quenched disorder is straightforward. This is in contrast to other systems exhibiting a depinning phase transition such as magnetic domain walls in disordered ferromagnetic thin films,⁴⁶ planar crack fronts in disordered solids,^39,40,47 or individual dislocation lines interacting with point defects.⁴⁸ In such systems, driven by an external force in a direction perpendicular to the average elastic line direction, the disorder sampled by a given realization is not known a priori. Nevertheless, it would be interesting to extend the present study to consider the predictability of depinning dynamics in such systems. Finally, we point out two crucial outstanding questions related to deformation predictability: The problem remains to be addressed in realistic three-dimensional DDD simulations, with^29,49 or without²⁶ additional defects interfering with dislocation motion, as well as in experiments where descriptors of the initial sample microstructure could be extracted, e.g., by x-ray measurement techniques.^50–52 

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

The authors wish to thank Henri Salmenjoki for technical help on machine learning. We acknowledge the financial support of the Academy of Finland via the Academy Project COPLAST (Project No. 322405).