Correlation of plastic events with local structure in jammed packings across spatial dimensions.

In frictionless jammed packings, existing evidence suggests a picture in which localized physics dominates in low spatial dimensions, d = 2, 3, but quickly loses relevance as d rises, replaced by spatially extended mean-field behavior. For example, quasilocalized low-energy vibrational modes and low-coordination particles associated with deviation from mean-field behavior (rattlers and bucklers) all vanish rapidly with increasing d. These results suggest that localized rearrangements, which are associated with low-energy vibrational modes, correlated with local structure, and dominant in low dimensions, should give way in higher d to extended rearrangements uncorrelated with local structure. Here, we use machine learning to analyze simulations of jammed packings under athermal, quasistatic shear, identifying a local structural variable, softness, that correlates with rearrangements in dimensions d = 2 to d = 5. We find that softness—and even just the local coordination number Z—is essentially equally predictive of rearrangements in all d studied. This result provides direct evidence that local structure plays an important role in higher d, suggesting a modified picture for the dimensional cross-over to mean-field theory.

Mean-field calculations have emerged as a promising first-principles approach to understanding glasses and amorphous solids. Although these calculations are only exact in the limit of infinite spatial dimension (), they successfully capture many qualitative phenomena (1, 2) in d = 2, 3, including the time dependence of glassy relaxation, signatures of a Gardner transition to a marginal glass phase, and aging behavior within this phase. Furthermore, mean-field theory captures some aspects of jamming criticality in d = 2, 3 quantitatively, including critical exponents and the d dependence of the prefactors of scaling laws (1, 3).

Despite these successes, many questions remain regarding the cross-over from low-dimensional behavior to mean-field theory with increasing d. The existing literature collectively supports a dimensional cross-over picture in which localized effects thwart some mean-field predictions in low d but recede with increasing d, leading to smooth convergence to the mean-field limit (4–8). For example, near the zero-temperature (T = 0) jamming transition, rattlers and bucklers (particles with too few or just enough contacts to be locally stable) cause the scaling of the low tail of the contact force distribution to differ from the mean-field prediction, but these particles become exponentially rare with increasing d (4, 5).

An unexplored aspect of the dimensional convergence to mean-field behavior is the nature of rearrangement events. This aspect is key to glassy dynamics because such events, in which particles experience large, sudden displacements and change their relative positions, are the mechanism by which supercooled liquids relax. Particle rearrangements are also responsible for plasticity in athermal jammed packings under mechanical load, such as shear strain. In low d, rearrangements are localized, and local structure plays an important role (9–27). In particular, machine learning has identified a linear combination of local structural quantities named “softness,” which is highly predictive of rearrangements and provides insight into the underlying physics (18, 28–38).

In particular, local structure is predictive of localized rearrangements in systems under athermal, quasistatic shear in d = 2, 3 (18, 28, 37). Indirect evidence supports the prevailing dimensional cross-over picture, suggesting that this low-dimensional localized physics should diminish with increasing d, which leads to a decreasing correlation between local structure and rearrangements as d rises. First, at T = 0 in low d, each localized rearrangement corresponds to a quasilocalized vibrational normal mode whose frequency vanishes; as d increases, however, quasilocalized low-frequency modes are increasingly outnumbered by extended ones (6, 7). Second, structure becomes more homogeneous in higher d. Rattlers and bucklers, associated with soft deformations and deviations from mean-field behavior, disappear as d rises, and the SD of the relative excess coordination, , decays asymptotically as , suggesting that in high d, all particles have very similar local environments (5). Based on these observations and the existing dimensional cross-over picture, one would expect a cross-over in d from localized rearrangements strongly correlated with local structure to extended rearrangements that are necessarily much less correlated with local structure.*

Here, we test this hypothesis directly by quantifying the importance of local structure as a function of d via the prediction accuracy of machine-learned softness in athermal quasistatically sheared (AQS) jammed packings in dimensions d = 2 to d = 5. Although local structure is known to be predictive for AQS in low dimensions (18, 28), mean-field theory predicts plasticity with qualitative features matching low-d behavior (43, 44), making this a natural system to explore the dimensional cross-over to mean-field behavior. In contrast to the expectation that local structure only matters in low d, we find that softness is no less predictive of rearrangements as d rises. Moreover, softness increasingly coincides with the number of interacting neighbors or coordination number, Z, for each particle. Altogether, our results suggest a different picture of dimensional cross-over for plasticity, in which the distribution of Z plays an important role and a particle’s coordination number remains a valuable indicator of mobility even in the mean-field limit.

Approach

We first prepare bidisperse jammed packings of Hertzian particles at many pressures. Each system is then sheared athermally and quasistatically by applying small strain steps and minimizing the energy at each step. While past work considered static packings in spatial dimensions up to d = 10 (5–7), here we focus on 2 to 5 due to the much larger computational cost of quasistatic shear trajectories.

Next, we seek to identify structures that correlate with local fluctuations in the displacement field according to the framework established in ref. 38. In this scheme, we quantify local structure by counting the number of contacts and gaps at each distance in a triangulation of the packing. To describe the dynamics, we follow what has become standard practice (11, 18, 45) and train to predict the first rearrangement in each avalanche, rather than seeking to describe entire avalanches (37, 46). To do this, we compute the lowest eigenvector of the dynamical matrix immediately before the stress drop. Using this eigenvector, we compute the quantity for each particle, measuring its nonaffine displacement relative to its neighbors (45, 47) (a precise definition is given in Materials and Methods). Because the system is solid, a rearranging particle exerts a long-ranged strain field on the system, causing power law decay of with distance. Following previous work, we define a quantity , where the average is carried out over neighbors. Thus, describes local fluctuations in .

Finally, we use linear regression to find a linear combination of our structural descriptors S, which correlates with . The correlation coefficient R2 for this regression is not high, but our goal is only to predict the rearrangement itself, not the full field; we have shown previously that S trained this way is just as predictive of rearrangements as previous classification-based approaches but requires fewer training examples because all particles are used in the training rather than only a small subset (38). Further details on our approach can be found in Materials and Methods.

Results

If a system has an extensive field of structural defects, under athermal quasistatic shear only one is triggered at the beginning of the avalanche, while the vast majority does not rearrange. Thus, to evaluate the predictive power of S, we must check not whether every particle with high S moves but whether the particle at the center of each rearrangement (the particle with the maximum value of in the critical mode) has a high S (18, 48). We calculate the average percentile rank of S for rearranging particles, C, as our measure of correlation between structure and dynamics. A value of C = 1 corresponds to perfect prediction, while C = 0.5 corresponds to random guessing. We find that C calculated from linear regression is just as high as for previous classification-based approaches (38).

This quantification of prediction accuracy is superior to comparing classification accuracies on a training set or regression accuracies because those numbers depend on the definition of the training set, even if the resulting softness S does not change. An alternative approach is to set a threshold on so that particles above this threshold are “rearranging,” and then, to ask which fraction of these particles has high softness (37, 49). However, this means of quantifying prediction accuracy is problematic for us because different threshold values must be selected for different dimensions and pressures, making meaningful comparisons across d and P impossible. In thermal systems, the isoconfigurational ensemble may be used to average over noise, allowing for completely threshold-free measures of correlation between structure and dynamics (20, 27), but the facts that a rearrangement may localize at a single defect and that it is necessarily coupled to an elastic response, which dictates much of most particles’ motion (38, 50), preclude such an approach for athermal quasistatic shear.

In Fig. 1, we report C for packings in all d as a function of pressure P. The error bars show an estimate of the uncertainty in C due to sampling error (38). We find that C is high, ranging from 0.85 to 0.9, and does not depend strongly on either P or d within uncertainty.

Fig. 1.

The S percentile C (A) and the –Z percentile C (B) of the particle with the highest in a plastic event as a function of dimension d and pressure P (d = 2, black triangles; d = 3, blue squares; d = 4, green pentagons; d = 5, red hexagons). Except for a weak decrease of C at low P and a more dramatic decrease at very high P, predictiveness depends little on P and d within statistical uncertainties. Except at high pressure, Z performs almost as well as S. Error bars are uncertainty in the mean value. (A, Inset) R2 for fitting of S to the local displacement fluctuations . Error bars show the full range over all pressures; all data are shown in (). Although the rearranging particle is predicted equally well in higher d, the correlation to the rescaled mobility field used to fit S is much worse in higher d. B, Inset shows that the correlation between C and C is strong, except at high P, with them nearly coinciding in d = 4, 5.

Although Fig. 1 shows that the local structure predicts the maximum of equally well in all dimensions studied, the linear regression coefficient R2, shown in Fig. 1 , Inset, decays dramatically between d = 2 and d = 5. Thus, in higher dimensions, fluctuations in are more poorly correlated with local structure, but the little correlation that remains is enough to identify a structural variable that correlates strongly with the rearrangement itself, as evidenced by the continuing high value of C. In other words, while S loses correlation with overall, it remains correlated with the extreme high tail of corresponding to rearrangements. For us, training on is only a means of finding a quantity (softness) that correlates strongly with rearrangements; the fact that C remains high in all d studied clearly demonstrates that structure is no less predictive in higher d. Note that C provides a lower bound on what can be achieved by machine learning methods.

Most of the predictive power of S comes from the number of contacts at distance (i.e., the coordination number Z of each particle). For comparison, we demonstrate the predictive power of Z by computing the average percentile of –Z, which we denote as C. In Fig. 1, we report C for the same set of particles used to measure C. Except at high pressures, we see that C is comparable with C and roughly independent of d; in 4d and 5d, they are equal within uncertainty (Fig. 1 , Inset). Thus, Z contains most of the local structural information in S, and it contains more and more of the available information with increasing d.

Various length scales diverge at the jamming transition (51), including the scale of spatial correlations of Z (52) and the length scale below which linear elasticity fails (53). Diverging length scales are also observed in nonlinear response (e.g., in relaxation, finite strain rate, and granular experiments) (54–56). Thus, an approach based on a fixed number of local structural descriptors might be expected to perform more poorly as since the number of necessary descriptors should diverge. Surprisingly, neither C nor C show a significant decrease with pressure. We further find that C shows no finite-size scaling at low pressures (), further indicating that the slight decrease in accuracy at low pressures is unlikely to be due to a diverging length scale (57).

In retrospect, it is not surprising that S reduces simply to the coordination number Z at higher d. In infinite d, it suffices to truncate the virial expansion at second order (i.e., to only include the effects of the nearest interacting neighbors) (58). With increasing d, a particle’s nearest neighbors are increasingly unlikely to also be neighbors of one another, so descriptors other than particle type, Z, and the number of nearest neighbor gaps should contain no information as .

What is surprising is that rearrangements still localize around low-Z particles as d increases, as shown by C. Because the low-frequency vibrational modes at zero strain become extended in high d, one expects rearrangements to become less localized in higher dimensions, causing a decrease of C (6, 7).

To determine the extent of localization of the initial rearrangement, we calculate the spatial correlations of and , along with those of S. At large distances, the of each rearrangement displays power law decay, consistent with continuum elasticity, but at short distances, there is exponential decay in correlations on the scale of the particle diameter (). This decay length gives a measure of the size of the rearrangement core, and we define it as the correlation length.

Fig. 2 shows the correlation lengths for , and S as a function of P in all d. At low P in low d, the correlation length grows but not as a power law and with no sign of diverging at unjamming. At high P and higher d, it appears to approach the S and correlation lengths, which are approximately equal and pressure independent. The correlation lengths all decrease with increasing d, indicating that rearrangements become more localized spatially. In , we also provide the inverse participation ratio of the critical modes and compare them with previous reports at zero strain () (6, 7, 59). These results show that the delocalization of the lowest-frequency modes as at zero strain does not necessarily imply delocalization of rearrangements at nonzero strain in higher d. Fewer and fewer quasilocalized modes exist (6, 7) with increasing d, but it appears that shearing picks one of the few remaining quasilocalized modes to go unstable, suggesting that the quasilocalized modes are more sensitive to shear strain than the extended disordered modes. This may not seem surprising when one considers that even in low dimensions, quasilocalized modes, which exist at arbitrarily low frequencies in the thermodynamic limit , are outnumbered by acoustic phonons but nevertheless, control plasticity. Our results suggest that the limits and do not commute; taking the former limit first yields localized rearrangements, while taking first eliminates localized modes that could produce such rearrangements, leading to extended ones.

Fig. 2.

Correlation lengths extracted from the fits to initial exponential decay of correlation functions as a function of pressure. In 2d and 3d, the correlation length is much longer than that of S or ; in d = 4, 5, they appear almost equal. All correlation lengths appear to decrease with increasing spatial dimension. Black, 2d; blue, 3d; green, 4d; red, 5d. The x’s indicate , the open squares indicate , and the circles indicate S.

It is possible to construct disordered mean-field models with localized low-energy modes or modes that are either localized or extended that correlate with local structure (39–42). However, existing numerical work on our system (7) suggests convergence to the spectrum of the perceptron, with zero density of low-frequency localized modes, as at γ = 0, N finite. Thus, mean-field models with a finite density of localized or quasilocalized low-frequency modes may not apply to these systems, where the density of states at the zero-strain ground state appears to converge to that of the perceptron, but nonetheless, a single localized mode dominates as the critical strain is approached, suggesting the existence of a subextensive but relevant set of localized modes.

We also note recent work on the convergence to mean field in the (thermal) random Lorentz gas (60). This work found more complicated dimension dependence in dynamical properties due to an interplay between mean-field caging and percolation in finite-dimensional systems. In particular, however, they find that while mean-field theory correctly describes the dynamics at any finite time in large d, it fails to describe the percolation transition in any large finite dimension because the limit fails to commute with that of , a scenario analogous to ours.

Fig. 2 shows that the correlation lengths of and S are the same, indicating that structural correlations set the size of fluctuations (61). Earlier work concluded that the spatial correlations of are the same as those of S (61); here, we find that spatial correlations of are longer ranged in d = 2, 3. However, the systems studied in ref. 61 did not exhibit elastic correlations over as long a range as we find here due to friction or thermal fluctuations. The correlation lengths of and in those systems are, therefore, likely to be far more similar.

It has previously been noted that “bucklers,” particles with the minimum number of contacts for local stability, are associated with deviations from mean-field behavior and become vanishingly rare in higher-dimensional packings (5). This appears to contradict our finding that C does not decrease with increasing d. To reconcile these results, we examine the probability that the rearrangement is located at a given small particle as a function of Z in each d in Fig. 3. Since the pressure dependence is weak except at high P, we bin together data from all P before the drop in C; full data are in . Comparing systems at the same pressure is not necessarily meaningful since P has different units in different d. The prediction accuracy of local structure, however, appears to be roughly constant as long as the pressure is not extremely high, justifying the average over P.

Fig. 3.

Probability of rearranging for small particles as a function of in varying spatial dimension as well as the statistics and σ for small particles used to scale Z. (A) Probability of rearranging for small particles as a function of in varying spatial dimension (d = 2, black triangles; d = 3, blue squares; d = 4, green pentagons; d = 5, red hexagons). After rescaling, data from all pressures before C decreases are binned together. As d and P are varied, the dependence of the rearrangement probability on this rescaled Z is similar, indicating that particles with relatively low Z are always more likely to rearrange, even in higher d where they are not bucklers. Error bars indicate statistical uncertainty in the mean. (B and C) Alternate attempts at data collapse under different scenarios for the d dependence described in the text. These collapses are worse than that in A, and as shown in , they become still worse when data from all pressures are plotted separately.

The curves collapse within error bars when Z is standardized by subtracting its mean and dividing by its SD σ. This collapse demonstrates that it is a particle’s coordination number relative to the distribution that determines its propensity to rearrange, rather than solely its status as a buckler or a rattler. The lower the standardized Z is, the higher the probability is to rearrange. In the limit , the standardized Z of rattlers and bucklers approaches , and they are never observed. Most rearranging particles have a value of Z that lies low in the distribution, consistent with Fig. 1. Indeed, our data suggest that it is particles with values of Z lying one to two SDs below the mean that are most likely to be the locus of a rearrangement in all d.

Fig. 3 addresses two other possible asymptotic behaviors of . Defining , there are three different scenarios that could prevail in high d. 1) Our results favor the scenario summarized in Fig. 3, where . Since the dynamics are a function of rather than , x is insufficient to describe the physics; as , but goes to zero at the same rate as . So, rearranging particles are still atypical, and is constant. 2) One may have naïvely expected that only particles with extremely low , such as bucklers, act as structural defects. This scenario is made mathematically precise by assuming that . Since , the typical for rearranging particles goes to one, and (no correlation) as . This scenario is consistent with the existing picture that localized physics becomes irrelevant in high dimensions. We do not, however, find a decrease in C with spatial dimension, and Fig. 3 and especially show that the collapse of predicted by this scenario is not as good as the first. 3) It has been found in the past that the distribution of contacts in large d has a large-deviation form (5). One might expect an analogous form , where near the bulk of the distribution to reproduce the apparent exponential dependence of . If and are expanded near the typical x = 2, one predicts (i.e., ); rearranging particles in large dimension are not necessarily bucklers but are still extremely atypical, and since while is finite, as . Thus, in this third scenario, the coordination Z predicts rearrangements perfectly in the mean-field limit. We do not see evidence for this increase in accuracy with d, and Fig. 3 and show again that this form does not collapse the data as well as scenario 1.

In , we show the dependence of and σ on P and d. We see that as expected. However, we find that in seems to be crossing over from scaling like d to the expected scaling (5). There is thus a small caveat to our results that at least σ has not yet reached its asymptotic scaling in d = 5. Nonetheless, the collapse in Fig. 3 and suggests that has converged to a form that implies that Z remains an accurate predictor in all d.

Discussion

Our results suggest a revision of the picture relating low-dimensional to mean-field behavior, at least in jammed packings. In contrast to expectation, local structure remains important in determining particle rearrangements for all d. In high d, the contact number Z for each particle becomes a good structural predictor of rearrangements, implying that rearrangements are controlled only by nearest neighbors that interact directly with a particle. Although the contact number distribution narrows relative to the mean contact number , the SD diverges, scaling as (5). Thus, different particles can still have different Z in high dimensions as , and low Z particles may still have a higher propensity to rearrange. It should be possible to capture such behavior in a mean-field description of structural glasses. For example, recent work has shown how mean-field spin glass models can have low-energy modes that localize on sites with a low local magnetic field (which acts as a heterogeneous structural variable) (41, 42). In those models, this correlation is proven by studying the low-energy vibrational modes; as discussed above, however, our picture suggests that such a computation for structural glasses must be carried out at nonzero applied strain before the limit is taken. Another possible route suggested by our work is a computation of in mean field [e.g., by computing for some suitable identification of R with mobility and then using Bayes’ theorem].

Previous work has assumed or supported a scenario, such as scenarios 2 and 5–7, but our results do not. Rather, they suggest the possibilities of scenarios 1 and 3. Scenario 1, in which rearrangements localize at particles of low Z and C is independent of d, is most consistent with our data. Either scenario 1 or 3 would reconcile the softness picture of dynamical behavior in low d with mean-field theory; in both cases, a particle’s coordination number is the local structural variable that controls the particle’s propensity to rearrange in high dimensions.

Materials and Methods

Simulations.

We prepare bidisperse jammed packings of Hertzian particles at many pressures in each spatial dimension: that is, with the pair potentialwith σ the radius of particle i and the step function. The two particle radii are σ and .

We study N = 4, 096 particles in 2d, N = 8, 192 in 3d and 4d, and N = 16, 384 in 5d. Initially, particles are assigned random positions () before quenching to an energy minimum, producing a force-balanced, zero-temperature state. The system is sheared quasistatically by applying small strain steps () and reminimizing the energy E with respect to at each step. Commonly, an affine deformation field is used as the “guess” for the new positions before minimization at each step; during the linear elastic branches, we speed up the simulation by using the nonaffine displacements from the previous strain step as the initial guess.

When the pressure P is large, it is insensitive to shear (), allowing us to hold the volume V fixed. However, when P is small, it fluctuates with shear, so we fix P by minimizing the enthalpy, , with respect to both and V (62). When a plastic event is detected in the form of a drop in shear stress, the resulting state is discarded, and the stress drop is approached with a shear step (63, 64). We repeat this procedure until the stress drop is approached with a shear step of , and the main obstruction to realizing true quasistatic shear is the force tolerance in the energy minimization algorithm (65, 66).

Structural Descriptors.

We quantify local structure by finding a combination of local structural descriptors that captures each particle’s propensity to rearrange. In previous studies, populations of rearranging and nonrearranging particles were identified from a collection of configurations and used to train a classifier to sort particles into these two groups (28). Here, we use linear regression to identify structures that correlate with local fluctuations in the displacement field (38). This method needs far less data to construct a good training set as it can use all particles rather than a small subset.

We describe each particle’s structure by counting the number of contacts and gaps at each distance in a triangulation of the packing (38). Once the Delaunay triangulation of the packing is constructed, each particle j in the triangulation is assigned a discrete distance d to particle i such that d = 0, d = 1 for neighbors, d = 2 for particles sharing a neighbor, etc. Each edge in the triangulation is either a contact between particles or a gap and is assigned a distance from particle i, . The numbers of gaps and contacts at each distance up to , together with the particle radius, form a set of 17 descriptors that we use to encode the local structural environment around each particle.

Identifying the Initial Rearrangement and Computing .

One complication is that quasistatically sheared packings exhibit avalanches in which one localized rearrangement triggers a second and so on at each stress drop (46). To focus on localized rearrangements rather than avalanches, we follow what has become standard practice (11, 18, 45), training to predict the first rearrangement in each avalanche by computing the lowest eigenvector of the dynamical matrix immediately before the stress drop. This eigenvector is the normal mode whose frequency vanishes at the onset of the rearrangement, signaling the rearrangement. Using this eigenvector, we compute the quantity for each particle, measuring its nonaffine displacement relative to its neighbors (45, 47).

We define using the two nearest neighbor shells: that is, for a displacement eigenvector u,where and Λ is chosen to minimize .

Definition of and Linear Regression.

Because the system is solid, a rearranging particle exerts a long-ranged strain field on the system, causing power law decay of with distance. A particle near the rearrangement shifts more than one far from it, regardless of their respective structural environments. To cancel this effect, we rescale for each particle by dividing by the average of its neighbors (38); we denote this quantity as . Thus,

We find that correlates far better with local structure than (38).

Next, we perform linear regression to identify a linear combination of our structural descriptors, the softness S, which best correlates with . Although the correlation coefficient R2 is not high, our training variable does not directly correspond to rearranging and non-rearranging particles, so R2 does not quantify the prediction of rearrangements.