Probing the transition from dislocation jamming to pinning by machine learning

Collective motion of dislocations is governed by the obstacles they encounter. In pure crystals, dislocations form complex structures as they become jammed by their anisotropic shear stress fields. On the other hand, introducing disorder to the crystal causes dislocations to pin to these impeding elements and, thus, leads to a competition between dislocation-dislocation and dislocation-disorder interactions. Previous studies have shown that, depending on the dominating interaction, the mechanical response and the way the crystal yields change. Here we employ three-dimensional discrete dislocation dynamics simulations with varying density of fully coherent precipitates to study this phase transition $-$ from jamming to pinning $-$ using unsupervised machine learning. By constructing descriptors characterizing the evolving dislocation configurations during constant loading, a confusion algorithm is shown to be able to distinguish the systems into two separate phases. These phases agree well with the observed changes in the relaxation rate during the loading. Our results also give insights on the structure of the dislocation networks in the two phases.


Introduction
While deforming, crystalline materials change irreversibly through discrete plastic events, i.e. avalanches, originating from the collective motion of dislocations -the topological defects of the crystal lattice [1]. These dislocation avalanches exhibit scale invariance with their distributions of sizes and durations encompassing several orders of magnitude [2]. This has lead to the discussion of plastic deformation as a non-equilibrium phase transition: below critical loading, the dislocations merely jump from one configuration to another, and the actual yielding of the crystal occurs at the critical point of diverging avalanches and uninhibited flow of dislocations. However, the dislocation movement has a highly complex nature arising from the interplay of the evolving, anisotropic interaction field produced by other dislocations and possible pinning field caused by disorder inside the crystal [3,4]. Thus, the collective dislocation behaviour is dictated by two competing phenomena − dislocation-dislocation interaction induced jamming and dislocation-obstacle induced pinning − that can be hard to distinguish from each other although they have fundamental differences [5,6,7].
Indeed in the case of dislocation jamming of pure dislocation systems, the interacting dislocations enter a state of 'extended criticality' where the system shows no distinct critical point but seems to recede in the constant vicinity of the transition independent of the loading force [5,8]. However, crystals are rarely completely pure, and introducing some disorder − such as precipitates − to the crystal to impede dislocation motion can increase the crystal's mechanical strength, and alter the criticality and ensuing avalanche behaviour of the system [9,10]. The key point here is that obstacles to dislocation motion may change the system behaviour by inducing dislocation pinning, which, if strong enough, results in a well-defined critical point of a depinning transition of the dislocation assembly [11].
Our recent study of 3D discrete dislocation dynamics (DDD) simulations of FCC aluminium with the inclusion of stationary fully coherent precipitates (see Fig. 1) showed that, by systematically increasing the strength or density of the precipitates, the system goes from the phase of dislocation-interaction dominated jamming to disorder-dominated pinning, and this transition can be observed in both constant load simulations as well as when quasistatically ramping up the external stress [7]. The related phenomenology depends on the loading protocol employed. For the quasistatic stress ramp simulations, one observes in general a sequence of strain bursts with a broad size distribution. In the jamming-dominated regime, the average strain burst size grows exponentially with the applied stress, while in the pinning phase we found a critical stress value where the average strain burst size exhibits a power-law divergence. Here, we focus on the creep-like constant loading simulations with varying precipitate density ρ p . There, the general behaviour in both of the phases, i.e. jamming and pinning, is on the one hand similar: In both phases the systems appear to possess a critical stress σ c (ρ p ) where one observes a power-law relaxation of the shear rate,ε ∼ t −θ . On the other hand, the relaxation becomes more rapid (larger θ) as the systems move further into the pinning phase [3,7]. This is illustrated in Fig. 2.
An open question regarding the phase transition from jamming to pinning is how exactly does it alter the dislocation structures in the systems? Furthermore, despite the apparent similarities in the response (i.e. power law relaxation) of the dislocation systems in the different phases, could one be able to distinguish them by their dislocation structures without specific a priori knowledge of the transition? To address this problem, we use machine learning (ML). ML is proving to be a flexible and useful tool for physics and materials science [12,13,14,15,16,17,18]. Using ML for the detection of phase transitions in statistical physics has given fruitful results [19,20,21] and here we applied the unsupervised 'confusion' scheme introduced in [22]. With the confusion algorithm, the only assumption one needs to make is that the system exhibits a phase transition in some control parameter range − in our case, the control parameter being the precipitate density ρ p − and the algorithm should be able to find the critical value ρ c p by using the states of the systems as input.
Here we followed the evolving systems by concentrating on both the fine details and the long-ranged structures of the dislocation network. To accomplish this we computed the dislocation junction lengths, geometrically necessary dislocation (GND) density and dislocation correlation, and used these separately to describe the microstructure for the ML algorithm. Our results show that the algorithm was able to find the phase transition from all of the used descriptors and the discovered values of ρ c p are in perfect agreement. Therefore, as the dislocation structures in the two phases evolve in notably different ways, we were able to quantify some of the changes in the systems by analyzing the used structural descriptors. The rest of this paper is structured as follows: the implementation of the ML method, along with the details of the DDD simulations and our approaches to characterize the dislocation structures, are presented in the next section. After the methodology, we proceed to show results obtained with the ML algorithm and we finish with some discussion.

DDD simulations
We study the effect of varying the precipitate density ρ p on the nature of the collective dislocation dynamics within 3D DDD simulations using our modified version of the ParaDiS code [23,24]. ParaDiS implements the dislocation interactions by approximating the continuous dislocation lines by a set of straight dislocation segments. The segments interact through the stress fields arising from the linear elasticity theory, while the diverging fields at dislocation cores are replaced by the use of results from molecular dynamics simulations. To cope with the long-range elastic forces, ParaDiS uses multipole expansion. Our version of ParaDiS also enables including disorder to the system, in the form of spherical precipitates [24]. The precipitates are frozen pinning sites for the dislocations that produce a short-range radial force where A is a constant, r is the distance from the precipitate to the dislocation and r p is the radius of the precipitate. In the context of transition from dislocationdominated jamming to disorder-dominated pinning, the relevant parameters are the precipitate density ρ p and the precipitate strength A [7]. For our simulations, we set parameters to approximate those of FCC aluminium with precipitates of fixed strength and size in a simulation box with periodic boundaries − A was especially chosen so that the system exhibits both jamming and pinning-dominated response depending on ρ p [7]. The parameters are presented in Table 1.
The simulations started with two relaxation periods, the first with only the dislocation networks and the second with also the precipitates present, to ensure the systems reached meta-stable states. After the initial relaxation, the systems were driven by applying a constant external stress σ. Depending on the magnitude of driving force, the systems tend to either get stuck (exponential decay of strain ratė ε with small σ) or reach linear creep-like conditions (constantε with large σ). However independent of precipitate density, all of the systems possess also a critical value σ c (dependent on ρ p , see Table 2) that leads to a power-law relaxation ofε, as seen in Fig. 2 [3,7]. The effect of precipitate density is seen in the rate and starting time of the power-law decay: there is a transition between the behaviour of less disordered systems withε ∼ t −0.3 and more disordered systems with the more rapid decay starting earlier. To see how this transition affects the system, we characterized the dislocation structure and observed its evolution during the constant stress loading with σ = σ c .

Characterizing dislocation structures
In the characterization of the dislocation structures, we used three distinct descriptors. First, we exploited the fact that disorder causes the dislocations to stretch when parts of the dislocations get pinned. With this in mind, we measured the length of dislocation links between two junction nodes [25] along the dislocation segments l along , and compared this to the shortest possible length between the nodes l shortest . Thus, we define parameter J, which represents the roughness of a dislocation and by collecting its distribution inside a system provides information on the dislocation structure. As an example, Fig. 3a shows the distribution of J in the simulated systems. The second used descriptor was GND density [26,16]. We computed the local GND density (the total GND density is constant throughout the simulation [27]) by first evaluating the Nye tensor α in voxels by where V voxel is the voxel volume, b is the Burgers vector, l is the line direction giving also the segment length and the sum is over all dislocation segments i inside the voxel. Then, the GND density ρ GN D was calculated from the Nye tensor. The resulting GND density fields, for instance the one with 10×10×10 voxels illustrated in Fig. 3b, are quite system specific, and as we are interested especially in the changes in the dislocation structure, we focused on the evolution of GND density, i.e. ρ GN D (t) = ρ GN D (t) − ρ GN D (0). Moreover to remove the effect of periodic boundaries, we took the Fourier transform of ρ GN D (t) as we collected the data.
As the third and final descriptor, we calculated the dislocation spacing correlation according to [28] where ρ is the total dislocation density and L(r) is approximated by computing the mean line length in spheres of radius r centered at random points along the dislocation structure. We note that in the case of 2D DDD simulations, the average dislocation-dislocation correlation function changed drastically when mobile solutes (pinning points) were introduced to the system [29]. Here we focused on longerrange correlations to avoid possible effects caused by the assigned segment length restrictions of the computations. Fig. 3c shows the dislocation correlation in systems with varying ρ p . We proceed by collecting the descriptors listed above during the loading with σ = σ c (ρ p ) at intervals of t = 10 −9 s = 2.6 · 10 5 GM , where the times are given in the units of shear modulus G times the dislocation mobility M = M edge = M screw . Due to computational challenges of 3D DDD simulations, we simulated only 19 systems for every value of ρ p and σ c . The time reached in every simulation was at least 4.7 · 10 6 GM although some systems were able to run even longer in their allocated simulation time.

Unsupervised learning of the phase transition
To observe the transition from dislocation jamming to pinning in an unsupervised manner, we used the confusion method presented in [22]. The idea is that, assuming the studied system experiences a transition in a control parameter range (in our case [ρ 0 p , ρ 1 p ]) with some value ρ c p , one expects that the different systems below and above ρ c p are distinguishable from each other. Thus by appointing trial values ρ p in the range [ρ 0 p , ρ 1 p ], the sample systems are assigned to classes depending on whether ρ p is below or above ρ p . This way, a machine learning classifier trained on the trial samples in supervised fashion should perform best near the critical point ρ p ≈ ρ c p where the systems are truly distinguishable. Correspondingly further from ρ c p , the classifications should get worse as some of the samples are wrongly labeled. If for instance a system was in jamming state (with ρ p < ρ c p ), trial value ρ p < ρ p would lead to the system being mislabeled to the pinning state with samples that actually belong to the pinning state. Then this labeling would be especially challenging for the classifier to learn because some of the samples in jamming state should be classified as jamming but some as pinning -therefore the confusion. Observing the accuracy of the classification in the range [ρ 0 p , ρ 1 p ] should therefore be somewhat W -shaped, as the accuracy is good at the transition but also at the beginning and at the end of the range (as large majority of the samples are labeled to one class, the classifier gets high score by simply predicting always the majority class).
As we were dealing with a data set of 190 systems with more than one thousand of collected features, we applied some dimensionality reduction before teaching any classifier. The three distinct data sets (different descriptors) were cast to lower dimensions by principal component analysis (PCA). In PCA, every feature of the data is first scaled to zero mean and unit variance, and then the entire dataset is represented by n orthogonal linear combinations of the original data which maximize amount of explained variance. This happens in descending order, so that with the first principal component (PC), the explained variance is the largest. Fig. 4 already shows that by projecting the data to the space of the two first PCs, there is a rather smooth transition in dislocation structures from less to more disordered landscapes with all of the used descriptors.
Supervised classifier for the confusion method For a classifier, our choice was based on linear discriminant analysis (LDA) [30]. In this simple case of 2-class classification, LDA builds one linear decision boundary into the input space according to where x is the feature vector of a sample, w and w 0 are the weights and bias of the classifier and y(x) = 0 is the boundary. We used the implementation by scikit-learn [31], that computes the boundary parameters by assuming that the samples inside different classes are Gaussian distributed, i.e. probability of a sample with features x when belonging to class k is where d is the length of x, µ k is the class-specific mean of features and Σ k is the covariance matrix. Moreover to obtain a linear boundary, the different classes are assumed to have identical covariance matrices, so in our case of two classes, k = −1 or k = 1, Σ −1 = Σ 1 = Σ. The weights for the decision boundary are obtained by applying the Bayes theorem, as at the boundary the probabilities of different classes given the sample are equal P (y = −1|x) = P (y = 1|x) and, thus, the log-probability ratio is log P (y = 1|x) P (y = −1|x) = log P (x|y = 1)P (y = 1) P (x|y = −1)P (y = −1) = 0.
From this, the final weights, w and w 0 , are obtained by substituting the probability distribution of Eq. 6 and comparing to Eq. 5, The LDA classifiers were evaluated by the straightforward accuracy, i.e. score S = number of correctly predicted test samples / number of test samples, and trained by 2-fold cross-validation to provide some tentative confidence intervals.

Results
The confusion curves with the different dislocation structure descriptors in Fig. 5a show the expected W -shape. What is striking, is that every curve shows a possible transition in the form of local maximum at the same spot, ρ c p ≈ 3 · 10 19 m −3 . Moreover, the classifying accuracy there is extremely good as every descriptor achieved score larger than 0.95 at the local maximum. Comparing the position of the transition to the relaxation curves with different ρ p of Fig. 2b and their power-law part represented by the exponents θ presented in Fig. 5b, we see that the relaxation behaviour is distributed nicely to the two phases so that θ is close to constant on one side (jamming) of the transition, while on the other side it starts to increase (pinning) [7]. Of course with ρ p = 5.1 · 10 19 m −3 θ seems to still be close to the constant value of the jamming-side of the transition, but there the error of θ, arising from the fact that the σ c used in simulations is impossible to get spot on to the one producing power-law relaxation, is notably higher than with other ρ p .
The used number of PCs for the best confusion curves (i.e. the curve with the highest maximum accuracy somewhere else than the ends of the range) was 5 for junctions and GND density, and 10 for correlations. Interestingly, the confusion curve obtained with junction lengthening data shows another distinct maximum near ρ p ≈ 3 · 10 20 m −3 , although there the accuracy is not as good as at ρ c p . Similar fluctuations from the pure W -shape are also observed in Fig. 6 which shows the confusion curves with different amount of PCs used for the classifying task. Basically all of the secondary maxima are positioned to the more disordered side with ρ p > ρ c p . Most likely this arises from the fact that in the pinning phase, the systems get more and more pinned with growing ρ p yielding faster relaxation with larger θ causing these systems to possess some distinguishability from each other despite being in the same phase. This also explains the tendency of slightly asymmetric W -shaped curves in Fig. 6, as the LDA score does not drop as much in the pinning phase as in the jamming phase. But as was ensured by the choice of the best confusion curves, the dominant maximum is indeed near ρ c p . We can also study how the ability of the confusion scheme to distinguish the two phases using different microstructure descriptors evolves in time by computing the confusion curves based on single snapshot structures, presented in Fig. 7. There the classifiers were trained by using two PCs of the dislocation structure at the specific times. Starting from the junction lengthening in Fig. 7a, there seems to be a short transient time until the single time step curves have converged to close to the shape of the best confusion curve in Fig. 5. This indicates that the junction lengthening shows early on the signs of distinct jamming and pinning phases. On the other hand GND density in Fig. 7b, which was measured as the difference to the initial density field, shows that the phases are separated well in the immediate beginning of the driving. However, the information about the transition is lost if looking at a momentary GND density field compared to one before loading. Trying the confusion scheme to GND density field without extracting the initial field or difference in the field of subsequent time steps yielded no observable phase transition (not shown here). Finally with the observed dislocation correlation functions in Fig.  7c, the behaviour is similar as with junction lengthening: There is now a longer time during which the transition is not observed, but after that the curves start to resemble the best confusion curve with maximum at ρ c p . Again, this is quite evident because the correlation functions focused on long-range structures, so it takes time until the systems have evolved structures that are noticeably different in the two phases. Notable here is also that the converged confusion curves are quite flat in the pinning phase.

Discussion
As the results with the unsupervised ML scheme showed, the dislocation configurations can be separated into two phases with different relaxation rates even though the general response, i.e. the power-law relaxation, is similar in the two cases. The confusion scheme succeeded extremely well, as it was able to achieve accuracy > 0.95 at the observed transition indicating that the systems where the dislocation-dislocation interactions dominate are significantly different from the precipitate-dominated systems. This was further supported by the fact that all the three dislocation structure characterization metrics considered captured the transition happening at the same value of ρ c p where also the relaxation starts to turn more rapid.
The success of all three descriptors reveals some of the notable differences between dislocation structures in the two phases. Firstly, the distribution of the junction lengthening J captures the bowing of the dislocation lines and, clearly, the pinning points cause more stretching and bowing of junctions than the other possible obstacles, namely the jamming dislocation structures, as depicted already in Fig. 3a. Secondly, the spacing correlation of the dislocations, C(r) shows that even longrange structures are slightly affected, although there the differences seem to arise more from the magnitude (and scaling by the total dislocation density in Eq. 4) than the shape of the correlation functions which are plotted in Fig. 3c.
Thirdly, the evolution of the local GND density finds similar structural changes as the other two descriptors: On one hand, the bowing dislocations are seen as a 'spreading' density of GND, while on the other hand with only few precipitates dislocations tend to move more in their straight forms. This is illustrated in Fig.  8 which shows the probability of a computational voxel having a non-zero GND density as a function of simulation time for different ρ p . The systems in the jamming phase show more or less constant number of active voxels as dislocations keep their shape while in the pinning phase the number is clearly increasing as dislocations bow. This happens despite the fact that the total GND density stays constant during the simulations. Undoubtedly, the effectiveness of GND density as a descriptor of the phase transition is also enhanced by the fact that in the pinning phase σ c is larger (faster changes in the dislocation structure right in the start of the simulation) but relaxation is more rapid (more constant structures on longer time-scale). However as Fig. 9 shows, the confusion scheme seems to be quite robust with respect to the resolution of GND density computation: even sparse number of voxels reveals the changes in the evolving structures.
To conclude our findings, we have studied the transition between dislocationdislocation interaction dominated jamming and disorder dominated pinning. By tuning the disorder content through precipitate density and strength, the system changes the mechanical response and yielding which is also seen in the power-law relaxation rate during the plastic flow with constant loading. Here we have been able to distinguish the simulated systems to the two phases of jamming and pinning solely by their dislocation structures during the constant stress simulations and, thus, highlighted the changes in the microstructure caused by the phase transition. These results offer two obvious prospects for future study: first, to conduct further simulations of the borderline case system where neither dislocation-dislocation nor dislocation-precipitate interaction dominates over the other. The second one is that our results tell that the dislocation structures are different in the two phases. This means that one can correlate these with the most interesting engineering quantity, the yield strength, possibly on a sample-to-sample basis as well. One should thus use the dislocation structure -oriented approach in the experimental verification of the different phases of crystal plasticity and for strength prediction.

Availability of data and materials
The data that support the findings of this study are available from the corresponding authors on reasonable request.      1.40 1.0 · 10 19 1.70 2.0 · 10 19 2.25 5.1 · 10 19 3.50 1.0 · 10 20 4.50 2.0 · 10 20 6.25 4.1 · 10 20 9.00 1.0 · 10 21 14.0        Figure 9 Confusion curve dependency on the number of GND density voxels. Only slight decrease of the maximum in the curve is observed depending on the number of voxels (i.e. voxel size from 280b to 3500b) when measuring GND density.