On local intrinsic dimensionality of deformation in complex materials

We propose a new metric called s-LID based on the concept of Local Intrinsic Dimensionality to identify and quantify hierarchies of kinematic patterns in heterogeneous media. s-LID measures how outlying a grain’s motion is relative to its s nearest neighbors in displacement state space. To demonstrate the merits of s-LID over the conventional measure of strain, we apply it to data on individual grain motions in a set of deforming granular materials. Several new insights into the evolution of failure are uncovered. First, s-LID reveals a hierarchy of concurrent deformation bands that prevails throughout loading history. These structures vary not only in relative dominance but also spatial and kinematic scales. Second, in the nascent stages of the pre-failure regime, s-LID uncovers a set of system-spanning, criss-crossing bands: microbands for small s and embryonic-shearbands at large s, with the former being dominant. At the opposite extreme, in the failure regime, fully formed shearbands at large s dominate over the microbands. The novel patterns uncovered from s-LID contradict the common belief of a causal sequence where a subset of microbands coalesce and/or grow to form shearbands. Instead, s-LID suggests that the deformation of the sample in the lead-up to failure is governed by a complex symbiosis among these different coexisting structures, which amplifies and promotes the progressive dominance of the embryonic-shearbands over microbands. Third, we probed this transition from the microband-dominated regime to the shearband-dominated regime by systematically suppressing grain rotations. We found particle rotation to be an essential enabler of the transition to the shearband-dominated regime. When grain rotations are completely suppressed, this transition is prevented: microbands and shearbands coexist in relative parity.

sequential and progressive process in which microbands form first in the early stages of loading and, depending on the initial density, shearbands may then develop as a result of cumulative local softening around a subgroup of microbands. However, they caution that confirming such causal links require a detailed analysis of these emergent structures. Accordingly, the first steps in such an analysis are attempted here.
On the whole, prior efforts to characterize deformation patterns in granular, and more broadly complex media, have mainly focused on the physical state space, both with respect to the space in which these patterns are analyzed and the metrics used to detect and quantify them 5,[9][10][11][12] . Most studies of deformation bands in the pre-failure regime adopt the conventional spatial gradient of deformation and other strain-based measures, as defined in Euclidean space. A perennial challenge confronting these approaches is the need to define a priori the size of the "window of observation" or spatial scale for the analysis. Patterns are inherently multiscale and what manifests at one spatial scale may be different at another scale. A further complication, supported by experimental evidence, is that the mechanisms underlying such patterns interact and give way to complex dynamics that change with loading history 7 . These complications may not only lead to novel patterns being missed but may also mask crucial aspects of coevolution among the known patterns of microbands and shearbands.
It is instructive to recall the known differences among these deformation bands. First, microbands have proved difficult to identify and characterize quantitatively due to their fine scale (i.e., a thickness of only a few particles wide) and dense criss-cross topology 9 . By comparison, system-spanning shearbands are easier to discern, since they are wider both in band thickness ( ≈ 8-30 grain diameters) and distance between bands [13][14][15] . Second, microbands prominently appear in the pre-failure or strain-hardening regime, albeit they exhibit highly transient dynamics therein 7,9 . Shearbands, by contrast, become persistent in the failure regime 16 . Third, microbands and shearbands differ in their orientation, as observed experimentally 7 and studied theoretically 7,12 . Lastly, the underlying mechanisms that govern microbands and shearbands also differ: the former is characterized by intense shear or slip events among particles 9 , while the latter is governed by the buckling of force chains [13][14][15] .
The question that then arises is: Is there a suitable state space and metric that can uncover and quantify kinematic patterns-across different spatial scales-while obviating the need to pre-define features of patterns, such as their spatial scale and geometry? To answer this question, we propose a new strategy for pattern discovery and characterization. Our approach shifts the analysis of kinematic patterns from physical space to the abstract state space of kinematics. To identify and quantify the relative dominance of distinct patterns in this state space, we introduce a novel metric which we call s-LID. This metric is based on the Local Intrinsic Dimensionality concept that was originally proposed by Houle 17,18 . In essence, s-LID accesses the intrinsic dimensionality of a reference point to its s nearest neighbors (sNN), applied here as a quantitative measure of how outlying a grain's motion is. By design, s-LID is relatively easy to use and yet has the unique capability of distinguishing concurrent deformation structures at multiple spatial scales without need to specify a priori a spatial scale to a target pattern.
We test the effectiveness of s-LID for a range of dense granular assemblies of spherical particles submitted to planar biaxial compression under constant confining pressure. Given the key role that grain rotations play in the development of instabilities in granular media, we suppress grain rotations to varying degrees to probe their influence on the uncovered patterns. The rest of the paper is arranged as follows. Methods are described in the next section followed by the results, discussion and a brief conclusion. Details of the studied data are also provided at the end of the paper.

Methods
An overview of the proposed framework can be found in Fig. 1. At each strain stage, we map the data from the physical state space (PSS) to the displacement state space (DSS) where each grain is represented by a point whose coordinates are the components of the grain's displacement vector. A s-LID score is then computed for each point (particle) in DSS. s-LID estimates the intrinsic dimensionality of a data point with respect to its sNN. That is, the minimal number of latent factors that are required to represent a given neighborhood. As demonstrated in Fig. 1, three example data points with various s-LID values are shown. p 3 has s-LID equal to 1.28. This indicates the local submanifold surrounding p 3 has intrinsic dimensionality close to 1, which can be seen from the almost linear structure of the neighbors surrounding p 3 (close to a diagonal line). For p 2 , the s-LID is 1.88, close to 2, as shown in the almost uniform spread of neighbors in two dimensions around p 2 . For p 1 , the s-LID is 3.6 (close to 4), indicating that the distances from p 1 to its neighbors are more characteristic of what one might expect to encounter in a 4-dimensional space. This can be seen by the large number of neighbors which are close to the edges of the sNN circle, i.e., far away from p 1 . In this work, s-LID is applied as a measure of outlying degree such that data points that are far way from their sNN are characterized by high s-LID scores, in order to capture deforming regions with abnormal motions than the other areas.
The s is a parameter to control the length scale of the neighborhood to investigate. To uncover distinct patterns in different spatial scales, we use various kinematic neighborhood sizes, starting from s = 25 , and doubled all the way up to the maximum number of neighbors a particle can have, s max (one less than the total number of particles in the system). The aim of using low neighborhood size is to identify local abnormal patterns presented in DSS, such that the microscale deformation pattern-microbands-can be highlighted. Whereas applying high s promotes the presence of global outliers that are showing abnormal motions relative to most of the particles in the system, in order to capture the shearband pattern. The whole procedure is repeated for all strain stages during the loading history. Based on these outputs, we then study and summarize the coevolution of the multiscale deformation patterns by showing the development of their relative dominance.
Local intrinsic dimensionality. The intrinsic dimensionality of a data set is the minimal number of latent factors that are necessary to represent its data points without significant information loss 17,18 . That is, we say a data set X ∈ R N×M , consisting of N samples described by M dimensions, has an intrinsic dimensionality of www.nature.com/scientificreports/ D if these samples can be effectively approximated using D latent dimensions, where D is usually significantly smaller than M. Intrinsic dimensionality is a global geometric property that represents the complexity of a data set as a whole. In contrast, Local Intrinsic Dimensionality (LID) proposed by Houle 17,18 measures the intrinsic dimensionality around an individual data sample by assessing the relative growth rate of its neighborhood size as the relative distance from the sample of interest increases. The formal definition of LID originates from the classical expansion model 19 that is used for estimating the global intrinsic dimensionality of a data set. Intuitively, the volume of a D-dimensional ball grows proportionally to r D when its size is scaled by a factor of r in Euclidean space. The expansion model measures the volume growth rate with expanding of ball size and the expansion dimension D can be deduced as: To learn the local dimensional structure around a particular data sample, Houle 17 restricts the estimation of expansion dimension to only the neighborhood around the sample of interest. In addition, since the underlying distance measure should not be limited to Euclidean, the probability mass is used as a proxy to the volume. Formally speaking, given a data sample x ∈ X , let r > 0 be a continuous random variable denoting the distances from x to other (neighboring) samples, F(r) be the cumulative distribution of r and is continuously differentiable over r, the LID of F at r is defined as: wherever the limits exists and the LID of x is then defined as the limit when r tends to zero: (2) LID F (r) lim �r→0 + ln(F(r + �r)/F(r)) ln(1 + �r r )  www.nature.com/scientificreports/ To gain more intuition about LID, consider a pure case in which points in the neighborhood of x are distributed uniformly within a submanifold. Under this setting, the dimension of the submanifold would be equal to LID(x) . In general, however, data distributions are not pure and so the manifold model of data does not perfectly apply, and LID(x) is not necessarily an integer. Practically speaking, by estimating the LID at x we obtain an indication of the dimension of the submanifold containing x that best fits the distribution.
The above formulation is based on the assumption of a continuous data distribution (i.e., a continuous distribution of neighbors around a chosen query point). Since F is typically an unknown function, several estimators have been proposed by Amsaleg et al. 20 to efficiently estimate the LID of a data point from data. In this work, the Maximum Likelihood Estimator (MLE) of LID is employed because of its ease of implementation and useful trade-off between statistical efficiency and complexity. Specifically, given the data sample x and a pre-defined parameter, s, to control the size of the neighborhood around x , MLE estimates the LID of x based on the distances between x and its sNN as: where s ≥ 2 , d i (x) is the Euclidean distance between x and its i-th nearest neighbor. We further assume that there is no coincidence in x and its neighbors in the feature state space such that the s-LID value is well-defined as a unitless score ranging from 0 to +∞ . The MLE estimator was also earlier proposed by Levina and Bickel 21 and is equivalent to the Hill estimator popular in extreme value theory 22 .
Intuitively, s-LID assesses the complexity of the local structure around a data sample based on the distribution of distances to its sNN and can be used as a measure of outlyingness. If the s-LID is high, the structure of the sNN around data sample x is more complex (higher number of latent variables needed to characterize) and the data sample x is more outlying in comparison to its sNN. On the other hand, if s-LID is lower, the structure of the sNN around data sample x requires less degrees of freedom to characterize, and data sample x is more inlying in comparison to its neighbors. Several recent works have demonstrated that s-LID can measure outlyingness of data samples in heterogeneous data manifolds. For example, Ma et al. 23 demonstrated adversarial samples that are generated to fool deep neural networks are usually characterized by relative higher s-LID than the normal samples. Similarly, for machine learning tasks with noisy labels, the presence of mislabelled data can be detected by properties of s-LID during the training process 24 . s-LID has also recently been used to characterize the complexity of image representations in neural networks 25,26 . Quantifying kinematic outlying degree via s-LID. Given the power of s-LID in discriminating data points in different local structures in the feature space, we propose to apply this metric to capture particles that are abnormal in their relative motions compared to their neighborhoods in complex granular systems. Specifically, for each strain stage during the loading, s-LID scores are estimated for all particles in the system. The analysis is carried out in DSS that is defined by the cumulative displacements of particles within an axis strain interval �ε yy = 0.002 for samples 5K and 5K-SR, and �ε yy = 0.005 for samples 20K and 20K-NR. That is, at strain stage ε , the i-th particle, p i , is described by is the cumulative displacement on the x(y)-axis from ε − �ε yy to ε . Such window-based cumulative displacements capture the intermediate motion information in the system, at the same time, provide more smooth signals compared to using instantaneous displacements that are more unstable. Displacements of particles are fundamental signals reflecting the underling spatiotemporal dynamics in the system and analyzing in DSS provides physical transparency such that factors like particle size, shape and local topology can be omitted. Meanwhile, essential deformation patterns like microbands and shearbands can still be effectively captured since particles within microbands/shearbands are expected to have distinct relative motions compared to others. Based on this, we quantify how outlying the motion of particle p i is relative to its kinematic neighbors by using s-LID in equation (4), where the distance between p i and another neighboring particle p j is defined as the Euclidean distance between their displacements, i.e., d( The higher the s-LID, the more abnormal the motion is relative to its s nearest neighbors.

Identifying outlying structures across length scales in DSS and PSS.
A wide range of kinematic neighborhood sizes, s, are applied for detecting distinct patterns in different length scales in DSS and PSS at each strain stage. To study how well the pattern presented in s-LID coincides with microbands/shearbands, we compare the s-LID pattern associated with each s to the ground truth. The challenge is that there is no explicit information on the belongingness of particles to the actual microband/shearband patterns can be directly used. To the best of our knowledge, no metric can perfectly elaborate all constituting particles of microbands and shearbands, neither individually nor jointly. For example, findings from Kunh 9 suggested rapid rotations occurring within and near microbands in granular materials at low strains, and buckling of force chains (BFC) has been shown to be a distinctive mechanism that is unique to regions inside shearbands [13][14][15] . However, such metrics are only the necessary, but not sufficient, conditions to microbands/shearbands. As a result, we use the label extracted from the rotation (for microbands) and BFC (for shearbands) patterns as ground truth, but only compute the proportion of member particles labelled (imperfectly) as microbands/shearbands that are also characterized by high s-LID values as the effectiveness of s-LID. Note the connection between rotation and microbands is only revealed www.nature.com/scientificreports/ under low loading condition and BFCs only start to appear in the early precursory failure regime, we confine the effectiveness check between s-LID and rotation patterns at strain stages from the start of loading ( ε 0 ) to the end of the compression phase ( ε b ), whereas comparison to BFC patterns is limited to strain stages since the beginning of dilating process till the final failure stage, i.e., ε ∈ (ε b , ε f ] . Specifically, at a given strain stage ε , let P (ε) s be the s-LID pattern that is defined as the set of particles whose s-LID values are higher than the average s-LID score of all particles (see Supplementary Information for why average s-LID is used as the cutoff value), P (ε) r be the rotation pattern, which consists of particles having higher than the average accumulated rotation till strain stage ε , and P (ε) f be the BFC pattern that includes all particles presented in any BFC till stage ε , the effectiveness of P (ε) s in terms of detecting microbands is defined as: r | is the number of grains characterized by both the s-LID and rotation patterns. Likewise, we quantify the effectiveness of P (ε) The quantitative analysis enables a direct comparison among different s such that the optimal neighborhood size at the low ( s * l ) and high ( s * h ) ends can be learned.
As shown later in our results, P (ε) s can effectively highlight particles presented in microbands/shearbands by adaptively changing the neighborhood size. In order to understand how clear a s-LID pattern is and to learn the dominant pattern among a range of spatial scales, we measure the contrast of a specific P (ε) s associated with neighborhood size s at strain stage ε as: where γ s,ε is the average s-LID score of all particles, γ + s,ε is the average score of all particles in P (ε) s , and γ − s,ε is the average s-LID of all other particles. δ (ε) s indicates the strength of a given s-LID pattern and quantifies the separation between two clusters: the negative cluster that is occupied by regular grains that move in near-rigid body motion, and the positive cluster of particles that are located in deformation regions with highly abnormal relative motions (as marked by relative high s-LID score).

Results and discussion
In this section, we present the key findings of our study, starting with the distinct deformation patterns revealed by s-LID, followed by a comparison to strain and other known metrics for detecting strain localization structures. Furthermore, we present new insights obtained from our analysis, including the transition of relative dominance of microbands and shearbands throughout the loading history, and the critical role that particle rotation plays in such a transition. Note the details of the studied samples can be found at the end of this paper. Among the four studied samples, 5K and 5K-SR are well studied systems that have been reported in various publications, with comparisons made against those found in various real materials: sand (e.g., Ottawa, Hostun, Caicos Ooid, in triaxial tests), photoelastic disk assemblies 27,28 , and more recently, field scale radar data on landslides 29 .

Deformation band patterns across length scales in DSS and PSS.
We visualize the location of particles with higher than average s-LID scores for the four studied samples in Figs. 2 and 3, together with ground truth microbands and shearbands characterized by the rotation and BFC patterns. These cover the system's status at different stages: from the initial phase of compaction ( ε a , ε b ), followed by the onset of dilatation in the strain-hardening regime until peak stress ( ε c , ε d ), then to the strain-softening regime ( ε e ), and then finally the near-steady 'critical state' regime ( ε f ).
For all the samples, similar patterns from s-LID across different neighborhood sizes s can be observed. The patterns comprise system-spanning, criss-crossing bands and rhombus-shape areas encapsulated by these bands. As demonstrated in Fig. 4, particles located in bands highlighted by s-LID are outliers that are moving abnormally compared to their neighbors in DSS, whereas grains in the rhombuses show similar displacement to their neighbors, suggesting near-rigid body motions. Among these concurrent multiscale band patterns, two patterns draw attention: microband-like and shearband-like patterns. The former effectively identifies the presence of microbands with low s value ( ∼ 100), especially at the early stages of loading history. The latter detects mesoscale deformation structures, i.e., shearbands, by applying s-LID with a neighborhood size of a few thousands, in the order of the size of the system, during the final failure phase ( Table 1). The bands presented in these two patterns are different in their thickness and spatial separation. Generally, for the small samples (5K and 5K-SR), the bands in the microband-like structure are a few grain diameters wide, and the gap between two parallel bands is around 10 grain diameters ( ε a , ε b , s = 200 in Fig. 2a for sample 5K, and s = 100 in Fig. 2b for sample 5K-SR). However, with large s applied to the later failure stages, thicker bands (varying from 5 to 20 grain diameters wide) can be found in the s-LID patterns and the spacing between successive bands grows to as high as ∼ 50 grain diameters long (sample 5K-SR, ε e , ε f , s = 1600 in Fig. 2b). Similarly, the average band thickness increases from ∼ 6 to ∼ 25 grain diameters wide for sample 20K, with the distance to the next parallel band jumping from ∼ 15 to ∼ 50 grain diameters long ( ε a , ε b , s = 100 for microbands, ε e , ε f , s = 3200 for shearbands, Fig. 3a). Band thickness and spacing in the microband-like s-LID patterns in sample 20K-NR are similar to sample 20K. However, less growth can be found when shifting from microbands to shearbands (increased to ∼ 12 grain diameters wide in band thickness, and ∼ 30 grain diameters long in spacing), suggesting a less concentrated but more diffused shearband pattern given the particle rotations are completely blocked ( ε e , ε f , s = 1600 for shearbands, Fig. 3b). www.nature.com/scientificreports/ In the early stages of low to moderate strains, the microband-like s-LID patterns produce well-defined but comparably thin criss-crossing patterns which resemble well the microbands characterized by particle rotation. This is particularly evident in sample 20K (Fig. 3a), where similar, albeit less distinct, criss-crossing patterns can be seen in the rotation pattern (comparing the red pattern at s = 100 to the gray rotation pattern). With the growth of strain, shearbands start to dominate the kinematics, as suggested by the increasing number of buckling force chains. The shearband-like s-LID patterns become more and more pronounced and are steadily concentrated and localized in the region where the ultimate shearbands form at the failure stage. On the other hand, for the s-LID patterns with lower s, the clear regularized criss-crossing patterns in the early stages are gradually distorted and in some places smeared away. More importantly, by applying a larger neighborhood size s in the early stages prior to the stress peak ( ε a , ε b , ε c ), early indications of the impending shearbands can be observed in the s-LID patterns. Take sample 20K as an example. There are four main shearbands presented at the final failure stage, forming a giant rhombus in the middle of the sample (see Fig. 7c for the ground truth shearbands). We denote, respectively, the upper and lower short forward-incline bands by SB 1 and SB 2 , and the upper and lower long backward-incline bands by SB 3 and SB 4 . The following three key observations can be made. First, at ε c , out of the four main shearbands, both SB 1 and SB 4 are fully captured by the s-LID pattern with s = 3200 , and half of the SB 2 and SB 3 are also highlighted with higher than average s-LID values (Fig. 3a). Second, particles located in SB 1 − SB 4 are also frequently presented in the s-LID patterns when relatively smaller s are applied (e.g., s = 400 to s = 1600 ), with small rhombuses surviving in the middle of the four main shearbands. Third, with the growth of s, more and more of the small rhombuses disappear, giving way to the bigger and ultimate rhombus that is enclosed by the shearband. Thus, s-LID is capable of not only capturing the embryonic-shearbands formed in the early stages of loading prior to the stress peak, which highlights the impending shearband area, but also the convolution of deformation patterns at different scales. Similar results can also be found in the other samples. Criss-crossing microbands remain apparent at low s, while embryonic-shearbands are evident at high s, thus yielding early clues to the location of f , and the larger the θ , the more strict for a force chain to be considered as buckling.  Fig. 2a) and form a similar but wider band to the ultimate forward-inclined shearband (Fig. 7a). The ultimate shearband pattern in sample 5K-SR forms a V/-shape region (Fig. 7b). This consists of three major bands: a backward-inclined band and two parallel forward-inclined ones: the upper band is located along the main diagonal of the sample, while the bottom one distances itself from this band by up to around 20 grain diameters wide. As shown in Fig. 2b ( ε c , s = 1600 ), both forward-incline and the bottom half of the backward-incline bands are identified by the s-LID pattern. For sample 20K-NR (Fig. 3b), the s-LID pattern at ε c , s = 1600 identifies a collection of significant bands that are present in the ground truth (e.g., the bottom-left to top-center, bottom-center to top-right, top-left to bottom-center, and top-center to center-right ones). Note that the final shearband pattern in this sample is spatially diffused, and comprises a large number thinner bands that are interacting with each other (Fig. 7d). Such embryonicshearbands are invisible to the conventional visualization methods which rely on the underlying mechanisms of shearbands like buckling force chains (the blue patterns in Figs. 2, 3), which has understandably led to the common belief that such mesoscale failure patterns only appear at or close to the peak stress 5,7,30 . The evidence provided here clearly show that the embryonic-shearbands exist even earlier-early in the pre-failure regime-but is consistent with past work 31 , where embryonic-shearbands in sand were first detected.
Comparison to strain and other metrics used to detect strain localization. We compare s-LID to other metrics commonly used to visualize microbands and/or shearbands. The most common metrics are grain rotations 9,29,32 and various measures of strain 7,10,11,33-39 . Rotations appear to identify only the dominant pattern for the given stage of loading, and cannot distinguish the different concurrent patterns across multiple scales. In addition, the quality of the microband pattern identified by the particle rotation is questionable, especially for the two small samples (gray patterns at ε a , ε b in Fig. 2 for 5K and 5K-SR) where the criss-crossing bands are   Table 1. Effectiveness of s-LID in detecting microbands and shearbands. At each strain stage, the effectiveness of all neighborhood sizes are quantified and ranked, then the average ranking of them over stages in the corresponding regime are reported. For each sample, the lowest ranking value (the lower the value, the higher the effectiveness) is highlighted in bold, which gives the optimal neighborhood size for detecting the corresponding strain localization pattern ( s * l for microbands and s * h for shearbands). Note s max = 5097 for samples 5K and 5K-SR, and no rotation pattern P (ε) r can be used for sample 20K-NR since rotation is completely suspended. www.nature.com/scientificreports/ clearly evident in s-LID. To further support this, we examined samples where particle rotations are completely blocked (sample 20K-NR, gray patterns from ε a to ε c in Fig. 3b). s-LID clearly identifies the microband patterns, whereas the rotation field trivially fails to detect the microbands in this sample. The spatial distribution of buckling force chains only identifies shearbands. The spatial gradient of deformation, strain, is a well-established metric to identify strain localization patterns. Typical strain-based measures are backed up by classical or micropolar continuum theories 7,10,11,33-39 . The classical continuum strain focuses on the spatial gradient of displacements, which is averaged across an assumed representative volume element (RVE) in PSS. Micropolar strain and curvature are similarly based on an RVE but incorporates local rotation which introduces an intrinsic spatial length scale to the analysis of strain localization. On the whole, the main drawback of all these RVE-based metrics lies in their limitation in identifying concurrent strain localization patterns. Unlike strain-based measures, s-LID is a data-driven measure derived entirely in DSS, and thus bears no assumed spatial length scale. This means it can be easily applied to find concurrent strain localization patterns of different spatial scales by adaptively adjusting the kinematic neighborhood size s. In addition, s-LID does not depend on the tenets of continuum mechanics, which constrains and reduces the dimensionality of the motions under consideration. Instead, s-LID depends neither on the type of material (including its microstructure), nor on the underlying mechanisms of deformation (e.g., dislocation, fracturing, slip). Furthermore, by extending the analysis to a higher dimensional state space, additional kinematic features such as rotation and displacement rate can be easily incorporated. Consequently, s-LID holds great potential for studies of deformation patterns in complex materials where pronounced outlying motions can develop due to interactions between concurrent mechanisms at different spatial scales.

Avg. ranking in E(P
One merit of strain-based measures is that their tensor form provides more information on the anisotropic microstructures in solid systems like crystalline material. In comparison, s-LID quantifies the outlying degree of a grain's relative motion as a scalar, with less information on the development of anisotropy in the material. That said, to some extent, s-LID indirectly captures the presence of 'kinematic anisotropy' . One example is p 3 in Fig. 1, which resides in a cluster of points lying roughly along a linear band in DSS. Additionally, in the samples studied, we observe clear preferential motions of particles in distinct directions at all stages of loading. This results in clusters in linear formation in DSS which in turns leads to s-LID values close to 1. Last, it is possible to extend s-LID into tensorial form, as one can estimate the s-LID of a grain based on its displacement component in each direction separately. To better demonstrate the advantages of s-LID in visualizing deformation patterns than strain-based measures, we present in Fig. 5 the spatial field of deviatoric strain that is commonly used in the literature for visualizing the microbands and shearbands for samples 20K and 20K-NR at strain stage ε c , including both the accumulated and incremental deviatoric strain fields, which are commonly used for capturing both short-and long-lived patterns. Compared to the clear system-spanning criss-crossing microbands ( s = 100 , Fig. 3a) and early detected embryonic-shearbands ( s = 3200 , Fig. 3a) presented concurrently in s-LID, the strain field method captures partially the embryonic-shearbands, but is clearly incapable of fully realizing all the microbands, especially in sample 20K (Fig. 5a). Similarly, much sparser criss-crossing microbands pattern can be found in the strain field www.nature.com/scientificreports/ than the s-LID pattern with s = 100 (Fig. 3b vs. Fig. 5b). More importantly, similar to other measures, the strain field also fails to distinguish between the patterns at different scales, i.e., shearband and microbands. A more advanced method for microband detection is presented by Kunh 9 , which is based on the void-graph model [40][41][42] to detect microbands as a set of thin obliquely trending chains of void cells within which slip deformations are most intense. Compared to this algorithm, the proposed s-LID is considerably more straightforward to implement, since it relies directly on the displacement data. Moreover, the s-LID method obviates the need to pre-define a priori any geometric features of the pattern, like the angle of inclination of bands. Instead, our approach automatically identifies deformation bands of various slopes. Overall, the comparison above conclusively demonstrates the superiority of s-LID over the other common schemes in capturing concurrent and coevolving deformation structures at multiple scales.
Hierarchy and coevolution dynamics of microbands and shearbands. The discovery and characterization of deformation patterns in granular media have mainly focused on the evolution of either microbands or shearbands in the physical (spatial) state space. This limitation is likely due to the different spatial scales of the patterns. Distinguishing the patterns from each other-especially when they interact or coevolve with a changing relative dominance as loading advances-is difficult when relying on tools that only consider the contrast between magnitudes of Euclidean measures (e.g., strain). In such methods, only the dominant pattern will likely be picked up. In this context, it is not surprising that a sequential evolution is the picture that emerges, whereby first the microbands are believed to form before the stress peak. As the sample reaches the stress peak, the more dominant (and connected) microbands are believed to coalesce leading to the formation of the shearband that spans the entire sample. The other microbands are thought to disappear due to the release of elastic energy (unloading) outside the shearband region. Inside the shearband, microbands formed between and through vortices 10,43 , but these strictly apply to the failure regime when the shearband is fully formed and not the whole of loading history.
The proposed s-LID method, however, puts forward a more nuanced hierarchy of coexisting patterns across all of loading history by considering the outlierness of kinematic measures regardless of their magnitude. The microbands and the pattern of embryonic-shearbands are seen to form concurrently in the early stages of deviatoric loading prior to the stress peak. Evidence of this embryonic-shearbands in the nascent stages of loading in sand and virtual samples has been previously reported in a study 31 which used the metric of closeness centrality, a property that is based on non-Euclidean geodesic distances between nodes of a kinematic complex network. A subsequent study of the coevolution dynamics of force chains and their supporting 3-cycles confirmed this pattern and elucidated a possible cause 32 . A progressive spatial bias in the degradation of the most stable subgroup (3 grains in mutual contact which persist in time), initiated by symmetry-breaking vortices, formed at the onset of loading 44 . This bias then predisposed colocated force chains to collective buckling which in turn gave way to shearbands in situ. None of these studies, however, uncovered a connection between embryonic-shearbands and microbands.
Here we are able to show not only the coexistence of microband and shearband patterns, but also quantify their coevolution by measuring their relative dominance. Specifically, the criss-crossing microband patterns are dominant in the early stages, while as the stress peak is reached, the embryonic-shearband pattern starts to overtake the microbands, leading to the fully formed localization pattern spanning the entire system in later postpeak stages. To demonstrate this in a quantitative manner, we measure and compare the contrast (as defined in equation (6)) between the microband and shearband patterns in Fig. 6. Across all the samples, a higher contrast can be seen in the microband pattern (red curves) in the early stages of loading, compared to the contrast of shearband pattern (blue curves). For the two small samples, 5K and 5K-SR, shearbands start to take the lead over the microbands as early as ε b , at the onset of dilatancy (Fig. 6a,b). The transition from the microband-dominant regime to the shearband-dominant regime initiates a little later for sample 20K (since the peak stress stage, ε d , Fig. 6c) but more intensely, given the contrast of shearbands are consistently doubled than that of microbands in this sample. That said, the shearband pattern in sample 20K-NR is less pronounced in the later stages of loading, as both the contrast of microband and shearband patterns fluctuate around 0.7 (Fig. 6d). Additionally, it is interesting to observe that, despite the dominance of the shearband, the microband criss-crossing patterns are still discerned, albeit in a mangled and distorted form, in the later stages of loading. Note the presence of small rhombus-shaped structures in the s-LID ( s = s * l ) patterns (e.g., ε c , ε d for sample 5K, ε d , ε e for sample 5K-SR, and ε e , ε f for sample 20K, Figs. 2, 3), which suggest deformation patterns at smaller spatial scales. As a result, in contrast to the common belief of a sequential formation of microbands and shearbands, patterns from s-LID suggest that they are coexisting structures since the early stages of loading, while the dominance of the shearbands over microbands is progressively amplified and promoted, as governed by a complex symbiosis among them.
Finally, we note that across all studied samples, there is a concomitant burst to a peak in the contrast of the shearband pattern during unjamming events when stress ratio drops and dissipation rises (not shown). We refer readers to the evolution in the amount of dissipated energy in Figure 3 in Tordesillas 13 for sample 5K, Figure 3a in Tordesillas and Muthuswamy 45 for sample 5K-SR, and Figure S5 in Supplementary Information for samples 20K and 20-NR, respectively.

Influence of particle rotations.
To study the influence of particle rotations in the coevolution of microbands and shearbands, we pay special attention to the comparison between systems with free and suppressed rotations, i.e., 5K v.s. 5K-SR, and 20K v.s. 20K-NR. As shown before in Figs. 2 and 3, clear criss-crossing microbands can be observed in all samples in the early microband-dominated stages, regardless of the extent by which particle rotations are suppressed. This suggests that: (a) the formation of microbands is independent of the exist- www.nature.com/scientificreports/ ence of particle rotations, and (b) a causal relation may be possible whereby the particles rotations follow and are promoted by the microband patterns, and not the other way around. However, unlike the 5K and 20K samples where thick, concentrated and localized shearbands are formed in the failure phase, the fully formed shearbands in the samples with suppressed particle rotation (partially for 5K-SR, completely for 20K-NR) are thin and spatially dispersed (see Fig. 7 for the shearbands formed at the failure stage). Such phenomenon is consistent with previous observations 14,46 . More importantly, it can be seen that the shearbands in these two samples are less dominant, as noted by the relatively smaller difference in the contrast of microbands and shearbands (Fig. 6). Consequently, it seems that particle rotation, as an energetically cheaper deformation mechanism, plays a critical role in this regime transition from the microband-dominated to the shearband-dominated, to the extent that this transition is altogether prevented when grain rotations are completely suppressed. This observation calls for more detailed studies into the significance of grain rotations in the coevolution dynamics of plastic deformation bands.

Conclusion
We presented a new and powerful method for identifying and quantifying localization patterns from kinematic data called s-LID. Direct comparisons with strain-based measures demonstrate the superiority of s-LID in revealing localization patterns invisible to the conventional strain-based methods, which renders s-LID the perfect tool for a wide range of deformation patterns from data for heterogeneous, complex media.
Unlike the conventional strain-based methods, s-LID describes the correlation between the motions of each particle and its s nearest neighbors in the state space of kinematics. By varying s, we uncovered a hierarchy of concurrent kinematic localization patterns of different length scales in samples submitted to planar biaxial compression under constant confining pressure. Contrary to common belief, patterns corresponding to both microbands and precursory embryonic-shearbands coexist in the early stages of loading, albeit the former is dominant. As the deviatoric loading unfolds, this imbalance is progressively tipped towards the embryonicshearbands until the shearband is fully formed in situ and becomes the dominant structure at the onset of the critical state regime. Results reveal that the transition from the microband-dominated regime to the shearbanddominated regime is obstructed when grain rotations are completely suppressed.
Our findings reinforce an enduring hallmark of complexity, viz., a hierarchical organization of structures at different levels and scales which coevolve. A more detailed study of this coupled evolution between microbands and shearbands is now the subject of an ongoing study and will be reported in future work.

Input data and DEM simulations
We presented the details of the four studied samples in Fig. 7 and Table 2. All samples exceed the size needed for realistic microband and shearband patterning. The samples 5K and 5K-SR are from a well-studied set of 3D discrete element simulations of a granular assembly of 5098 polydisperse spherical particles subjected to planar biaxial compression, under constant confining pressure 13  www.nature.com/scientificreports/ a rolling resistance and a sliding resistance act at the contacts, both of which are governed by a spring up to a limiting Coulomb value of µ|f n | and µ r R min |f n | , respectively, where f n is the normal contact force, R min is the radius of the smaller of the two contacting particles, µ and µ r are the parameters to control the sliding and rolling friction, respectively. A summary of all the interaction parameters governing the contacts and other quantities relevant to this sample is presented in Tordesillas 13 . The samples 20K and 20K-NR are simulated using YADE software 47 . The samples consist of 20,000 spherical particles confined within rigid walls that apply a strain-controlled biaxial loading. The contact law includes normal and tangential linear springs with latter limited by a Coulomb friction limit. The two samples differ in that particles in the sample 20K can rotate freely while the rotation is prevented in the sample 20K-NR.  The inset in panel c shows the same information for the early stages of loading, with strain in log scale for ease of presentation. ε a , . . . , ε f indicate different stages during the loading: ε b is the stage where volumetric strain hits its minimum, ε d is the peak stress stage, and ε f is the end of post-peak softening stage, as marked by the switching of volumetric strain to a flatten curve. The shearbands are identified as the accumulated buckling force chains exceeding the threshold buckling angle θ = 1 • till the final stage. Table 2. DEM simulation parameters and material properties for the four studied samples. Identical initial conditions are applied in samples 5K and 5K-SR, except for that 10 times larger rolling friction is used in sample 5K-SR, which results in suppression of rotation. Samples 20K and 20K-NR share the same initial conditions and differ only in particle rotation, which is blocked in 20K-NR. The contact stiffness for 20K and 20K-NR depend on r = r 1 r 2 /(r 1 + r 2 ) where r 1 and r 2 are the radii of the contacting particles.