Elementary processes in dilatational plasticity of glasses

Materials typically fail under complex stress states, essentially involving dilatational (volumetric) components that eventually lead to material decohesion/separation. It is therefore important to understand dilatational irreversible deformation -- i.e., dilatational plasticity -- en route to failure. In the context of glasses, much focus has been given to shear (volume-preserving) plasticity, both in terms of the stress states considered and the corresponding material response. Here, using a recently-developed methodology and extensive computer simulations, we shed basic light on the elementary processes mediating dilatational plasticity in glasses. We show that plastic instabilities, corresponding to singularities of the glass Hessian, generically feature both dilatational and shear irreversible strain components. The relative magnitude and statistics of the strain components depend both on the symmetry of the driving stress (e.g., shear vs.~hydrostatic tension) and on the cohesive (attractive) part of the interatomic interaction. We further show that the tensorial shear component of the plastic strain is generally non-planar and also extract the characteristic volume of plastic instabilities. Elucidating the fundamental properties of the elementary micro-mechanical building blocks of plasticity in glasses sets the stage for addressing larger-scale, collective phenomena in dilatational plasticity such as topological changes in the form of cavitation and ductile-to-brittle transitions. As a first step in this direction, we show that the elastic moduli markedly soften during dilatational plastic deformation approaching cavitation.


I. INTRODUCTION
Glassy materials are ubiquitous in the natural and technological world around us, and include various noncrystalline solids such as oxide glasses, glassy polymers, organic glasses and metallic glasses.These intrinsically disordered materials possess notable properties and hence find an enormous range of engineering applications [1][2][3][4][5].Processing glassy materials [6] and more so their performance, durability and structural integrity in various applications require deep and fundamental understanding of their mechanical deformation and failure modes.Failure typically involves complex stress states, essentially involving dilatational (volumetric) components that eventually lead to material decohesion/separation.This is evident from extensive experimental observations regarding cavitation in glasses under a wide variety of failure conditions (e.g., [7][8][9][10][11][12][13][14][15][16]).
Our goal in this work is to study the fundamental micro-mechanics, geometric and statistical-mechanical properties of the elementary processes that mediate dilatational plasticity in glasses.In terms of driving forces, we invoke the hydrostatic tension test and compare it to its well-studied simple shear counterpart, both shown in Fig. 1.The hydrostatic tension test applied to computer glass samples physically represents a mesoscopic portion of a macroscopic glass, on which the surrounding material exerts predominantly hydrostatic forces.In Fig. 1a, an ensemble-averaged shear stress-strain curve under simple shear is presented, obtained through 3D computer simulations using athermal quasi-static (AQS) deformation [17][18][19][20].The resulting curve features a small-strain linear regime and a smooth, monotonic transition to a steady-state flow plateau upon shear plastic yielding.The curve does not feature a stress overshoot (and the accompanying stress drop), a situation that is typical for poorly annealed (i.e., rapidly quenched) glasses, which indeed corresponds to the employed glass ensemble [55].In Fig. 1b, an ensemble-averaged dilatational stressstrain curve under hydrostatic tension (pure dilation) is presented, obtained using the very same ensemble of computer glasses as in panel (a).The resulting curve features a small-strain linear regime, followed by an abrupt stress drop and continuous strain softening.It qualitatively differs from its simple shear counterpart in Fig. 1a, despite performing the deformation simulations on the very same ensemble of glass realizations.One manifestation of this qualitative difference is the emergence of a cavitation instability in the hydrostatic tension test, as illustrated in the inset (see figure caption for additional details).A major goal of this work is to study the elementary irreversible processes that contribute to the observed differences.
Our methodology, involving computer glasses of N particles and initial volume V 0 , is discussed in Appendix A, where additional technical details are provided in [55].As material cohesion is essential for understanding dilatational plasticity and failure, our computer simulations employ a class of recently introduced potentials, which feature both repulsion and cohesion/attraction, where the strength of the attractive part is continuously adjustable through a dimensionless parameter r c [55][56][57][58].The smaller r c is, the stronger the cohesive/attractive part of the interaction, as demonstrated in the inset of Fig. 2a.It was recently shown that reducing r c can lead to a ductile-to-brittle transition [34].
Simple shear AQS deformation in a given direction is controlled by a shear strain γ, where a representative simple shear stress-strain curve σ(γ) is presented in Fig. 1a.Hydrostatic tension (pure dilation) is controlled by a dilatational (volumetric) strain ϵ, representative hydrostatic tension stress-strain curve −p(ϵ) is presented in Fig. 1b (p is the hydrostatic pressure).Unstable plastic eigenmodes u(r) (where r is a position vector relative to the center of the eigenmode [55]), corresponding to a zero crossing of an eigenvalue of the Hessian M, are identified during AQS deformation, controlled either by γ or ϵ (see Appendix A and [55]).u(r) features a highly-nonlinear, disordered core of linear size a (i.e., of volume V ∝ a 3 ) and a power-law decay |r| −2 in the far field (e.g., [20]), |r| ≫ a, associated with a linear elastic continuum behavior.
The irreversible deformation inside the nonlinear core is quantified through an eigenstrain tensor E * in the framework of Eshelby's inclusions formalism [59,60].A recently developed method [61], based on a class of contour integrals evaluated in the far field |r| ≫ a, allows to extract V E * .Our primary goal is to study the properties of E * and V (or alternatively a) as a function of the symmetry and magnitude of the driving force, quantified by γ and ϵ, and as a function of the strength of the cohesive/attractive part of the interatomic interaction, quantified by r c .

II. PLASTIC EIGENSTRAIN TRIAXILITY
The eigenstrain tensor can be additively decomposed into its dilatational (volumetric) part, E * dil , and deviatoric (shear, volume-preserving) part, dev .In 3D, E * is characterized by 3 independent amplitudes, one characterizing the dilatational (isotropic) part, E * dil = ϵ * dil I and two characterizing the deviatoric part, which features tr(E * dev ) = 0. We denote by ϵ * dev,1 the largest (in absolute value) of the three eigenvalues of E * dev and by ϵ * dev,2 the second largest, which has a different sign compared to ϵ * dev,1 (the third eigenvalue is not independent, but is rather determined through tr(E * dev ) = 0).
The 3 independent amplitudes that characterize E * provide basic information about the geometry of plastic rearrangements/instabilities in glasses.Of particular interest in the present context is the amplitude of the dilatational eigenstrain component ϵ * dil , and more specifically its relative magnitude compared to the deviatoric component.Consequently, we aim at constructing a ratio of the two components in order to quantify their relative magnitude.This goal naturally fits the contour integrals method [61] that -as noted above -allows to extract only the product V E * , but not V and E * individually.In fact, as u(r) is a normalized mode (by construction), considering absolute eigenstrain values is not immediate.Yet, as will be shown later, the fully nonlinear u(r) does contain information that allows to estimate the core volume V.
In order to construct a ratio that quantifies the relative magnitude of the dilatational and deviatoric components of the eigenstrain tensor E * , we draw analogy with the well-known and widely-used stress triaxiality measure [62,63] -constructed at the macroscopic scale for a similar goal, but with respect to the stress tensorand define where We term the ratio R t in Eq. ( 1) the plastic eigenstrain triaxiality.Note that R t is a signed quantity, as will be further discussed below.
In Fig. 2, the probability distribution p(R t ) is plotted for both the simple shear and hydrostatic tension tests, and two values of r c .p(R t ) was obtained by calculating R t of Eq. (1) for 200 independent glass samples (made of N = 128K particles each) per r c value, where for each sample up to the first 50 plastic instabilities/events were detected [55].Under dilation, we focused on strains below the cavitation strain (corresponding to the peak stress under hydrostatic tension, cf.Fig. 1b), implying that not all of the detected plastic events were analyzed.Overall, the presented distributions were generated using between 3000 and 10000 events each.Statistical convergence and possible effects of the magnitude of the applied strain are discussed in [55], also for other observables to be considered below.
In Fig. 2a, we show p(R t ) for r c = 1.5, which corresponds to the canonical Lennard-Jones interatomic potential (see inset).It is observed that under simple shear (red circles) p(R t ) is symmetric and narrowly peaked around R t = 0 (i.e., a vanishing dilatational eigenstrain, ϵ * dil = 0).Note that ϵ * dil > 0 corresponds to isotropic core expansion/dilation and ϵ * dil < 0 to isotropic core contraction [55].It is observed that finite values of ϵ * dil exist, though R t is small, indicating that the deviatoric eigenstrain component is significantly larger than the dilatational one.Under hydrostatic tension (green circles), for the very same ensemble of glass realizations, the symmetry between core expansion and contraction, R t → −R t , is broken, and p(R t ) is biased toward positive R t values.p(R t ) is peaked at relatively small values, again indicating a dominance of the deviatoric component over the dilatational one, even when the global driving force is dilatational (hydrostatic tension).
The results presented in Fig. 2a demonstrate that plastic instabilities in glasses generically feature both dilatational and shear components.A corollary of this finding is that at least part of dilatational plasticity is carried by the same micro-mechanical objects that carry shear plasticity.Moreover, the relative magnitude of the dilatational and deviatoric components depends on the stress state, i.e., the microscopic plastic strain is not an intrinsic/geometric property of plastic instabilities (for a given glass composition and interatomic interaction).Finally, hydrostatic tension gives rise to larger dilatational plastic strains, though for the canonical Lennard-Jones potential the deviatoric plastic strain is dominant.
We then repeated the calculations leading to Fig. 2a for an ensemble of glasses generated with r c = 1.2.This value of r c corresponds to an interatomic potential featuring a stronger cohesive/attractive part, see inset of Fig. 2a.In the presence of stronger cohesion, one expects that in order to trigger an irreversible rearrangement, particles at the core of a plastic mode (saddle configuration) would feature a larger separation, i.e., p(R t ) is expected to feature larger values.The results in Fig. 2b support this expectation.
It is observed that under simple shear (red squares) p(R t ) is still symmetric around R t = 0, but becomes bimodal with peaks at R t values significantly larger than the typical R t values for r c = 1.5 (cf. the red circles in Fig. 2a).Moreover, this physical effect significantly intensifies under hydrostatic tension, as demonstrated by p(R t ) in Fig. 2b (green squares).It is observed that under these conditions, r c = 1.2 glasses feature almost only positive R t values, and p(R t ) is broad, where plastic events with a sizable dilatational component exist with a nonnegligible probability.The geometry of plastic cores with different R t values is illustrated in the visual insets of Fig. 2b, see figure caption for details.
Overall, the results of Fig. 2 show that a dilatational plastic strain generically exists in plastic instabilities in glasses, along with a deviatoric plastic strain, and that their relative magnitude depends on both the driving conditions and the cohesive/attractive part of interatomic interactions.Moreover, our findings show that at least part of dilatational plasticity in glasses is carried by the same micro-mechanical objects that carry shear plasticity, and consequently that these provide a coupling mechanism between the two.
The generic existence of a dilatational plastic strain in plastic events in glasses has also been demonstrated previously in 2D Lennard-Jones glasses studied under uniaxial tension/compression [64] and under simple shear [65], and in 3D computer models of amorphous silicon studied under simple shear [66].Amorphous silicon has been modelled in [66] using either the standard/original Stillinger-Weber (SW) potential [67] or a modified SW (termed SWM therein) potential [68,69].The latter potential, which features a three-body term twice as large as the original SW potential, has been developed in relation to the fracture of crystalline silicon, where it was found that the SW potential leads to a ductile behavior, while the SWM potential to a brittle one [68,69].This is reminiscent of the ductile-to-brittle transition induced by reducing r c , as discussed in [34].The ductile-to-brittle transition in both cases is also accompanied by a reduction in Poisson's ratio (see Table I in [66] for the SW and SWM values, and Table S2 in [55] for the different r c values).
Interestingly, it was reported in [66] (cf.Fig. 7 therein) that the relative magnitude of the dilatational and deviatoric components of plastic events under simple shear (termed "shear-tension coupling") increases from SW to SWM.This is precisely the trend observed in Fig. 2 (red symbols in the two panels) with decreasing r c .Consequently, while "ductility" and "brittleness" are clearly collective phenomena that are affected by both the history dependence of a glass (i.e., its initial state of disorder emerging from the quench self-organization, being in itself affected by the interatomic potential) and by spatiotemporal interactions of an extensive number of plastic events during deformation [16,34], they might also have some signature in the geometry of individual plastic events.

III. PLANARITY OF THE DEVIATORIC EIGENSTRAIN TENSOR
The plastic eigenstrain triaxility R t , studied above, constructs out of the 3 independent amplitudes that characterize the eigenstrain tensor E * a measure of the relative magnitude of the dilatational and deviatoric eigenstrain components.Next, we consider another geometric property of the core of plastic instabilities.As the dilatational component is isotropic, i.e., it features 3D spherical geometry, we focus our attention on the geometry of the deviatoric eigenstrain tensor E * dev , which is characterized by the 2 independent amplitudes ϵ * dev,1 and ϵ * dev,2 , defined above.
A lot of previous insight into glass plasticity has been gained using computer simulations in 2D, which -as highlighted above -focused mostly on shear plasticity.In 2D, the deviatoric eigenstrain tensor is characterized by a single independent amplitude.There were some preliminary indications in the literature that plastic instabilities in 3D feature such a planar deviatoric eigenstrain tensor as well (e.g., [70]).That is, in terms of the abovedefined quantities, it was suggested that ϵ * dev,1 ≃ −ϵ * dev,2 (such that the magnitude of the third eigenstrain is much smaller than |ϵ * dev,1 |) also characterizes plastic instabilities in 3D.Yet, to the best of our knowledge, this issue has not been systematically investigated in the past.
To address this basic issue, we define the planarity ratio of the deviatoric eigenstrain tensor E * dev of plastic instabilities as As explained above, in the purely planar limit one eigenvalue of E * dev vanishes and we have ϵ * dev,1 = −ϵ * dev,2 , i.e., R p = 0.The opposite limit, the least planar one, corresponds to ϵ * dev,2 /ϵ * dev,1 = −1/2, i.e., R p = 1.These limiting cases are illustrated in the visual insets in Fig. 3a (see figure caption for additional details).In these terms, the suggestion that the deviatoric eigenstrain tensor remains planar in 3D corresponds to p(R p ) → δ(R p ), i.e., to a probability distribution p(R p ) that is strongly concentrated near the planar limit R p ≃ 0, approaching a delta-function distribution.Our goal is to study p(R p ) as a function of the external driving force (simple shear vs. hydrostatic tension) and r c .The analysis was performed on the same set of plastic instabilities/events discussed in relation to Fig. 2, and the results are presented in Fig. 3.In Fig. 3a, p(R p ) under simple shear is plotted for two values of r c (those previously used in Fig. 2, see legend in Fig. 3b).It is observed that for r c = 1.5, p(R p ) features a broad peak around R p ≃ 0.3, significantly deviating from δ(R p ).Moreover, for r c = 1.2 (corresponding to a stronger attractive part of the interatomic potential), p(R p ) is also broad, but is peaked at significantly larger values of R p , away from the planar limit.The corresponding results under hydrostatic tension are presented in Fig. 3b, and are similar (note that p(R p ) is more sharply peaked around R p ≃ 0.7 compared to the corresponding result under simple shear, cf.panel (a)).Overall, the results indicate that the geometry of the deviatoric eigenstrain tensor E * dev is generally non-planar, i.e., that the dimensionality of E * dev is the same as space dimensionality.

IV. THE CORE SIZE OF PLASTIC INSTABILITIES
As pointed out above, the contour integrals method, which employs the far field linear elastic properties of unstable plastic modes u(r) to infer the core properties, does not allow to separately extract the core volume V and plastic eigenstrain tensor E * .Yet, the core volume V -or equivalently the linear core size a (with V ∝ a 3 ) -is an important glassy lengthscale.For example, it has been argued that V influences various physical properties of glasses [71][72][73].Consequently, it is important to extract V ∝ a 3 itself.In fact, the particle-level (atomistic), fully nonlinear u(r) can be used to estimate a.
The unstable plastic mode displacement field u(r) is evaluated at each spatial position r i occupied by a particle, hence it can be denoted as u i ≡ u(r i ) (which is the quantity we actually calculate to begin with).u i follows a continuum linear elastic behavior for |r i | ≫ a, yet it features significant nonlinearity and larger displacements over shorter distances.A widely used measure of spatial localization is the participation number (i.e., the participation ratio multiplied by N ) [74,75], which is O(1) in the extreme localization limit and O(N ) in the spatially extended limit.Consequently, the participation number provides an estimate for the number of particles inside the core of u i , and we estimate the dimensionless core size as Here, a 0 ≡ (V 0 /N ) 1/3 is a characteristic interparticle distance, i = 1−N runs over all particles and we used the fact that unstable plastic modes are normalized by construction, i.e., i |u i | 2 = 1.The probability distribution p(a/a 0 ) is presented in Fig. 4 for two values of r c , and under both simple shear and hydrostatic tension, for the plastic events previously analyzed in Figs.2-3.In Fig. 4a, we plot p(a/a 0 ) under simple shear for both r c = 1.5 (light blue circles) and r c = 1.2 (yellow squares).It is observed that p(a/a 0 ) is peaked at a few interparticle distances.Moreover, it shifts to larger values and becomes wider with increasing r c .Note that in the definition in Eq. (3) we do not account for the increase in the volume per particle in dilation (but rather use the fixed linear size a 0 ) because over the range of dilatational strains we consider, the implied changes are small.
In Fig. 4b, we present the corresponding results under hydrostatic tension, which are quantitatively similar to the simple shear results presented in panel (a).The inset in panel (a) shows an example of the amplitude of an unstable plastic mode ⟨|u|⟩ Ω (the solid angle Ω average of the norm of u(r)) as a function of |r|/a 0 .The estimated core size a/a 0 is marked by the horizontal green double arrow and the linear elastic far-field power-law ∼ r −2 is highlighted (red dashed-dotted line).The results of Fig. 4 show that the core size of plastic instabilities is nearly independent of the symmetry of the loading, but does depend on the interatomic interaction potential.Specifically, it becomes more compact with increasing attractive forces, i.e., decreasing r c .This trend is also consistent with observed correlations between decreasing Poisson's ratio ν and the plastic core size (sometimes termed the STZ size/volume) [57,66,76], see the variation of ν with r c in Table II of [55].
The characteristic core size a is similar to the core size ξ g of quasilocalized, nonphononic modes in glasses [77][78][79][80][81][82].Quasilocalized glassy modes, defined in the reference/undeformed glass, have been recently shown to follow a universal density of states ∼ ω 4 (distinct from Debye's density of states of low-frequency phonons, where ω is the vibrational frequency), and to feature a spatial structure similar to that of unstable plastic modes.In particular, quasilocalized modes are characterized by a nonlinear, disordered core of linear size ξ g and a linear elastic power-law decay |r| −2 in the far field, |r| ≫ ξ g , exactly like unstable plastic modes.Moreover, the trend of variation of a and ξ g with varying r c is similar (compare our results to Fig. 7i in [57]).Overall, our findings reinforce the suggestion that plastic instabilities are predominantly quasilocalized modes that are driven to a saddle-node bifurcation.

V. THE SOFTENING OF ELASTIC MODULI UNDER DILATATIONAL PLASTICITY
We provided above a basic quantitative characterization of the geometry and statistical-mechanical properties of the elementary micro-mechanical carriers of plasticity in glasses, for two end members of driving force symmetries (simple shear vs. hydrostatic tension) and variable strength of cohesive/attractive interatomic interactions.In particular, the distributions of geometric properties (quantified by the ratios R t and R p ) and a characteristic length (core size a) of plastic instabilities have been calculated.Our focus was on dilatational plasticity and its comparison to shear plasticity, both in terms of the driving forces and the material response manifested in the localization length a and the eigenstrain tensor E * .
The collective spatiotemporal accumulation of these elementary plasticity processes, including their coupling and emerging spatial organization, gives rise to largerscale plasticity in glasses.Much of these collective phenomena, in relation to dilatational plasticity in particular, remain to be explored and understood.We stress (again) that we focused on unstable plastic modes and not on the outcome of the instabilities, and also note that we do not claim to have exhaustively identified all elementary dilatational plasticity processes in glasses.For example, we have not discussed micro-cavitation, i.e., the process by which cavities on the particle scale are formed (to be distinguished from the larger-scale cavitation observed in Fig. 1b).These may be related to the micromechanical objects we identified (e.g., a subset of the outcomes of the identified plastic instabilities) or correspond to different objects.Yet, we would like to conclude this paper by demonstrating the type of new effects associated with collective dilatational plasticity.We already demonstrated a qualitative difference between collective shear and dilatational elasto-plastic dynamics in Fig. 1, manifested at the level of stress-strain curves.Here, we provide another example, focusing on the strain evolution of elastic moduli, which can serve as global proxies for the structural changes experienced by a glass during elasto-plastic deformation.In Fig. 5, we present the strain evolution of the ensemble-averaged bulk modulus K and shear modulus µ, under simple shear and hydrostatic tension.We stress that K and µ, which characterize minima of the potential energy (see precise definitions in [55]), should not be confused with the tangent moduli, which characterize derivatives of the global, fluctuations averaged stress-strain curves.The two loading protocols are applied to the very same glass ensemble, composed of 200 samples generated with N = 128K and r c = 1.2 (e.g., compared to N = 10K and r c = 1.5 in Fig. 1), resulting in strains approaching the steady-state regime in the simple shear case and going beyond the stress drop in the hydrostatic tension case (cf.Fig. 1).
It is observed that under simple shear (Fig. 5a), both elastic moduli are largely unaffected by plastic deformation.In particular, the bulk modulus K(γ) is essentially independent of γ and µ(γ) softens by ∼ 7% from the initial (as cast, undeformed) glass to the steady shear flow.The corresponding results under hydrostatic tension are presented in Fig. 5b.First, it is observed that both moduli experience a large drop associated with large-scale cavitation (see Fig. 1b), here around a cavitation strain of ϵ c ≃ 0.045.Second, both moduli significantly soften prior to ϵ c , where the softening is more pronounced for the bulk modulus.
We note that the expressions for the elastic moduli include an overall factor 1/V [55], where V is the current volume.Under simple shear, which is volume preserving, we have V (γ) = V 0 .However, under hydrostatic tension, V (ϵ) ≥ V 0 is an increasing function of the dilatational strain ϵ.It would be therefore interesting to disentangle the contribution of 1/V (ϵ) to the observed softening under hydrostatic tension from other contributions by defining the reduced moduli, μ(ϵ) ≡ µ(ϵ)V (ϵ)/V 0 and K(ϵ) ≡ K(ϵ)V (ϵ)/V 0 .The results are presented in the inset of Fig. 5b.Interestingly, by superimposing the corresponding moduli under simple shear, which are identical to those presented in Fig. 5a, it is observed that the reduced shear modulus almost coincides under both loading symmetries, prior to cavitation (in the hydrostatic tension case).On the other hand, the reduced bulk modulus under hydrostatic tension still significantly deviates from its simple shear counterpart, indicating the existence of intrinsic softening processes on top of the varying volume.This observation demonstrates qualitative differences between the shear and bulk moduli.
The softening of the elastic moduli emerges from the accumulation of dilatational plastic deformation and possibly micro-cavitation in a way that is not yet understood.Understanding these softening processes will also shed light on the qualitative differences between the shear and bulk moduli.More generally, understanding the spatiotemporal dynamics of collective dilatational plasticity upon approaching large-scale cavitation is a challenge for future work.

VI. BRIEF SUMMARY AND OUTLOOK
In this work, we studied elementary processes in glass plasticity, with a focus on dilatational plasticity and its comparison to its well-studied shear plasticity counterpart.We have extracted, using large-scale AQS computer simulations, the basic micro-mechanics, geometry and statistical-mechanical properties of unstable plastic modes (zero crossing of an eigenvalue of the glass Hes-sian) as a function of both the symmetry of the applied driving forces (simple shear vs. hydrostatic tension) and the strength of the cohesive/attractive part of the interatomic interaction.These modes feature an effective Eshelby eigenstrain tensor over a spatial scale that defines their plastic core.
In particular, we computed 3 probability distribution functions of the following quantities: (i) the plastic eigenstrain triaxility R t (a dimensionless measure of the relative magnitude of the dilatational and deviatoric parts of the eigenstrain tensor), (ii) the planarity ratio R p (of the deviatoric part of the eigenstrain tensor) and (iii) the linear core size a.We found that R t strongly depends on the symmetry of the applied driving forces and on the strength of the cohesive/attractive part of the interatomic interaction (Fig. 2).We also found that the deviatoric part of the eigenstrain tensor is generally nonplanar, and that the statistical deviation from planarity is larger for more cohesive/attractive interatomic interactions (Fig. 3).Finally, we found that the statistics of a are almost independent of the symmetry of the applied driving forces, but dependent on the strength of the cohesive/attractive part of the interatomic interaction (Fig. 4).The latter results provide additional support for the intrinsic relations between unstable plastic modes and nonphononic modes in glasses [33,82,83].
At larger scales, involving the collective, spatiallycoupled dynamics of many elementary plasticity processes, increasing hydrostatic tension leads to a cavitation instability upon which internal surfaces are formed.Cavity formation is accompanied by a large tension drop, but not an entire loss of load bearing capacity, which persists to larger dilatational strains (Fig. 1b).Interestingly, large-scale cavitation is preceded by a significant softening of the elastic moduli (Fig. 5).
These results open the way for various research directions in dilatational plasticity, which can be roughly classified into two categories.First, additional insight into elementary plasticity processes should be gained.It remains to be understood whether the initial glassy state of disorder, e.g., as controlled by the quench rate through the glass transition temperature, affects the core properties of unstable plastic modes or just their occurrence probability [34].
In this work, we studied unstable plastic modes, which are well-defined micro-mechanical objects that correspond to saddle points in the glass potential energy landscape.In this context, it is important to better understand both the origin of dilatational plastic instabilities and their outcomes, where the latter constitute the actual contribution to plastic deformation.Preliminary results (not shown here) indicate that the core properties of unstable plastic modes are correlated with the plastic strain accumulated as the glass reaches another potential energy minimum.Future work should systematically explore these correlations.In the context of the origin of plastic instabilities, the relations between quasilocalised, nonphononic modes in glasses (and other structural indi-cators [33]) to plastic instabilities under hydrostatic tension should be explored.In addition, the emergence of micro-cavitation, i.e., of particle-scale cavities, should be studied.
Second, at larger scales, future work should address the collective spatiotemporal organization of plastic deformation in glasses in the presence of hydrostatic driving forces, also going beyond the AQS limit (i.e., including finite temperatures and strain-rates).These should include the softening of the elastic moduli on the way to cavitation, demonstrated above, as well as large-scale cavity formation and subsequent dynamics.The coupling and competition between shear and dilatational plasticity should be addressed, including the relations and interplay between shear-banding and cavitation (e.g., [16]).The emerging insight about dilatational plasticity should be eventually incorporated into coarse-grained elasto-plastic models, which will enable to address the fracture toughness of glasses (i.e., irreversible deformation and damage near the edges of crack defects [84,85]) and phenomena such as ductile-to-brittle transitions [16,34,86].

ACKNOWLEDGMENTS
A.M. acknowledges support from the Minerva center on "Aging, from physical materials to human tissues".E.B. acknowledges support from the Ben May Center for Chemical Theory and Computation and the Harold Perlman Family.

Appendix A: Methodology
In this Appendix, we describe in more detail the methodology used in this work and the motivation behind it.Additional technical details are offered in [55].We employed in this work computer models to address elementary plasticity processes in glasses.Large-scale computer simulations played crucial roles in various recent developments in glass physics (see, for example, [27,87]).There are multiple reasons for this situation; first, computer glasses provide access to particle-level spatial scales that are essential for understanding glass plasticity, yet they are inaccessible through cutting-edge, real-time experimental techniques.
Second, computer simulations allow to control interatomic interactions in a way that goes well beyond current laboratory techniques.Particularly relevant for understanding dilatational plasticity is the ability to control the strength of the cohesive/attractive part of the interatomic interaction, as noted above.Specifically, we employed potential energy functions U (x), where x are the particle coordinates, composed of central force interatomic interactions of the Lennard-Jones type (see detailed formulation in [55]) in which the cohesive/attractive strength is continuously adjustable, through the dimensionless parameter r c [55-58] intro-duced above.The effect of varying r c on various mesoscopic [57,58] and macroscopic [34] quantities has been recently studied.Among these, we highlighted its effect on the fracture toughness of glasses, where reducing r c can lead to a ductile-to-brittle transition [34].
Third, computer glass simulations can be performed under athermal quasi-static (AQS) conditions [17][18][19][20], corresponding to the zero temperature and strain-rate limits.The AQS protocol is a powerful tool for studying fundamental aspects of the micro-mechanics, geometry and statistical-mechanical properties of glassy deformation.Its main merit is that it allows to exhaustively and unambiguously identify discrete plastic processes along the entire deformation path, as done in this work.Finally, computer glasses offer great advantages in implementing the deformation protocols described in Fig. 1.In the context of simple shear deformation, cf.Fig. 1a, employing periodic boundary conditions allows to eliminate surface effects and hence reach very large strains without failure.Moreover, the hydrostatic tension test, cf.Fig. 1b, allows to represent a mesoscopic portion of a macroscopic glass that experiences predominantly hydrostatic forces exerted by the surrounding material.We also note in passing that the application of isotropic dilation, which is readily accessible on the computer, is not easy to realize experimentally on the global scale (e.g., compared to the uniaxial tension test [16]), but is feasible (e.g., [88]).
Glass samples in 3D, each with a fixed number of particles N and initial volume V 0 , were generated by rapidly quenching high-temperature, equilibrium binarymixture liquids into zero temperature inherent states, as detailed in [55].While the non-equilibrium thermal history (or more generally thermo-mechanical history) of a glass has profound implications for its glassy state of disorder (e.g., [36]), and correspondingly for its material properties -plastic deformability in particular -, we did not vary it in this work.It is important to note that our glasses feature vanishingly small initial pressure (as evident in Fig. 1b), which for the fixed quenching protocol is achieved by tuning V 0 (at a given N ).This is important for revealing the intrinsic dilatational plasticity response of glasses.
Simple shear deformation (cf.Fig. 1a) in a given direction is controlled by a shear strain γ.The latter is obtained through the accumulation of AQS strain increments dγ [55], for which the deformation gradient tensor F s = I + dγ x ⊗ ŷ is applied (here I is the identity tensor in 3D, x and ŷ are Cartesian unit vectors, and ⊗ is a dyadic/outer product).The stress tensor, and in particular the simple shear stress component σ(γ), as well as the shear µ(γ) and bulk K(γ) moduli, were extracted [55].
Hydrostatic tension (pure dilation, cf.Fig. 1b) is controlled by a dilatational (volumetric) strain ϵ.The latter is obtained through the accumulation of AQS strain increments dϵ [55], for which the deformation gradient tensor F d = (1 + dϵ) I is applied.The stress tensor, and in particular the hydrostatic tension (negative pressure) −p(ϵ), as well as the shear µ(ϵ) and bulk K(ϵ) moduli, are extracted.
The potential energy U (x) of the glass is minimized during AQS deformation, controlled either by γ or ϵ.As shown above, and consistently with [43], at least part of the plastic deformation in glassy materials -independently of whether obtained under simple shear or hydrostatic tension (or more complex stress states) -occurs through the accumulation of discrete plastic rearrangements (instabilities) that correspond to a zero crossing of an eigenvalue of the Hessian M ≡ ∂ 2 U ∂x∂x (saddle-node bifurcation [17,19]).
The corresponding particle-level eigenfunctions/eigenmodes u(r) (where r is a position vector relative to the center of the eigenmode [55]), is accurately identified and extracted under AQS conditions.The plastic rearrangements (modes) u(r) feature a highly-nonlinear, disordered core of linear size a (i.e., of volume V ∝ a 3 ) and a power-law decay |r| −2 in the far field (e.g., [20]), |r| ≫ a, associated with a linear elastic continuum behavior [89].The irreversible deformation occurs inside the nonlinear core, whose averaged effect is quantified through an eigenstrain tensor E * in the framework of Eshelby's inclusions formalism [59,60].
A recently developed method [61], used in this work, allows to extract V E * .Earlier approaches, discussed in the works of [64][65][66]91] (some of which are discussed above), also invoked Eshelby's inclusions in similar contexts.We note that we analyzed the unstable modes very close to the saddle-node bifurcation, which are well defined mathematical objects, and not the outcome of instability obtained once a new energy minimum is reached.The relation between the two was briefly discussed above.Note also that we considered simple shear and hydrostatic tension as end members of a continuum of stress states, which can be studied as well.
Supplemental material for: "Elementary processes in dilatational plasticity of glasses" In this Supplemental materials file, we provide additional technical details regarding the quantities used and computations reported on in the manuscript.In particular, we provide details regarding the interaction potentials and the glass preparation protocol.We also provide details regarding the plastic instabilities detection algorithm, and contour integrals method for extracting the core properties of plastic instabilities.Finally, we present some additional results that support statements made in the manuscript.

Appendix A: Interaction potentials and preparation protocol
To study the effects of variation in interaction potential, we used binary mixtures of small and large particles, interacting via two variations of the sticky-sphere potential, introduced and studied in [56][57][58]92].The pairwise potential takes the form with ε being a microscopic energy scale, r ij ≡ |r i − r j | is the distance between particles i and j, x min = 2 1/6 is the dimensionless location of the minimum of the potential, x c is a modified cutoff, and λ ij is an interaction length, expressed in terms of the "small-small" interaction length λ small small , λ large small = 1.18λ small small , and λ large small = 1.4λ small small .The parameters a, b, {c 2ℓ } detailed in Table I are set such that the attractive and repulsive parts of ϕ ij , as well as its first two derivatives are continuous at x min and x c .We refer to the two sticky-sphere variations used in our work by defining r c ≡ x c /x min .When plotting the potentials in the inset of Fig. 2a in the manuscript, we have used r = r ij /λ ij , and the normalized r = r/x min such that the minimum occur at r = 1 [57,58].
We have initialized 200 samples of 50:50 binary mixtures with N = 128K and equilibrated them at a high parent temperature T p = 4.0 employing the Berendsen thermostat [93].Then, we quenched the samples instantaneously to zero temperature, followed by overdamped dynamics.To obtain the illustrative stress-strain curves of Fig. 1 in the manuscript, we used 108 independent systems of N = 10K particles and performed athermal quasi-static deformations to large strains, i.e, beyond the yielding transition.
Most studies in the literature focus on simple shear deformations, where little attention is given to the initial pressure resulting from the glass preparation protocol.However, as we are interested in applying both simple shear and hydrostatic tension to the same glass realizations, we aim at ensuring that the initial pressure is small, preferably vanishingly small compared with the as-cast bulk modulus K.This is important for revealing the intrinsic dilatational plasticity response of glasses, which is commonly probed in the laboratory at ambient (hence, small) pressure (obviously varying the initial pressure may affect the subsequent response to hydrostatic tension).This is achieved by tuning the initial number density for the different r c values (for our fixed quenching protocol).Using the densities stated in Table II before the quenching procedure, we attain the stated pressures at the end of the quench, without any additional processing.
Details regarding the initial densities, and resulting pressures appear in Table II.In particular, we have p(ϵ = 0)/K(ϵ = 0) ≃ O(10 −3 ), where K(ϵ = 0) is the initial bulk modulus (to be discussed in Sec.C), for both r c 's and across samples.As explained in the manuscript, we apply to our glass samples athermal quasi-static (AQS) deformation, either in simple shear or in hydrostatic tension.Under hydrostatic tension, periodic boundary conditions in all space directions are imposed.Under simple shear, we employ the Lees-Edwards boundary conditions in the shearing direction [94].Below, we specify how we apply deformations and detect instabilities under simple shear deformations, and will use γ to denote the strain.We use the same procedure for loading and instability detection under hydrostatic tension (unless stated otherwise), such that γ and ϵ are interchangeable.
The AQS deformation is performed conventionally, where at each step we apply a global affine transformation to the particle positions x, quantified through a deformation gradient tensor F , followed by energy minimization that results in non-affine displacements.As explained in the manuscript, we either use F s = I +dγ x⊗ ŷ for simple shear or F d = (1+dϵ) I for hydrostatic tension.Here I is the identity tensor in 3D, x and ŷ are Cartesian unit vectors, and ⊗ is a dyadic/outer product.As stated FIG.S1.An illustration of the first two stages of the plastic instability detection algorithm.In the first stage, the algorithm starts from the initial state (green) and imposes athermal quasi-static deformation increments (blue dashed upper parabolas).It results in a coarse-level estimate of the instability (blue-shaded region).In the second stage, we adaptively decrease the strain increment (inset, red dashed lower parabolas), resulting in a refined estimate of the instability (red-shaded region).The adaptive refinement is applied iteratively until an accurate estimate of the instability strain is obtained, as described in detail in the text.
above, we hereafter use dγ to symbolically denote strain increments either under simple shear or hydrostatic tension.
Plastic instabilities under AQS conditions, which are the main micro-mechanical objects discussed and analyzed in the manuscript, correspond to vanishing eigenvalues of the glass Hessian M ≡ ∂ 2 U ∂x ∂x evaluated at x(γ), i.e.M is a function of the imposed strain γ.Consequently, one's goal is to calculate M(γ) and its lowest eigenvalue (i.e., diagonalize the Hessian) along the entire deformation path x(γ).This poses several serious technical challenges; first, it is computationally expensive to perform these computations along the entire deformation path.Second, it is crucial to very accurately detect the strain at the saddle point γ s at which plastic instabilities take place.In particular, an a priori choice of strain increments dγ would be most likely insufficient, i.e., too coarse.
The procedure we invoked to meet this challenge is illustrated in Fig. S1.We highlight therein 3 states, an initial state (which either corresponds to γ = 0 or to a final state resulting from a previous instability), a saddle (that we aim at accurately detecting) and a final state (emerging as the outcome of the instability).The crucial point is the selection of the strain increments dγ, as mentioned above, which is obtained using 3 stages (the first two stages were already used in [95]): • In the first stage, the coarsest one, we initially set dγ = 10 −4 , marked by the blue dashed upper parabolas in Fig. S1, and identify pairs of consecutive states x(γ) and x(γ + dγ).We then compute the non-affine displacement between each pair as ∆(γ) ≡ T naf [x(γ +dγ)−x(γ)]/ f (γ).Here, the operator T naf removes the affine displacements and accounts for boundary effects (e.g., particles crossing periodic boundaries), and f (γ) is the particle-averaged force induced by the affine transformation.The largest component of ∆(γ) is denoted as ∆ max (γ).We repeat this process for each pair of states and consider the ratio ∆ max (γ + dγ)/∆ max (γ).Once the latter surpasses a threshold value Θ 0 = 2, i.e., once the non-affine displacements significantly increase with strain, we deduce that a plastic instability took place.This is illustrated by the last blue upper parabola in Fig. S1, implying that -on this coarse level -the instability is estimated to occur somewhere in between, represented by the blue-shaded region therein.• In the second stage, which aims to refine the identification of the instability strain within the blue-shaded region, we apply an adaptive and iterative refinement process, as schematically illustrated in the inset of Fig. S1 by the red dashed lower parabolas.We consider a state obtained a few coarse strain steps earlier (before the coarselevel identification of the instability).We then redefine our threshold to be Θ ≡ ∆(γ + dγ) max /∆(γ) max > Θ 0 , halve the strain increment and repeat the deformation procedure of the first stage at this finer level.We repeat this process of halving the strain increment and increas-ing the adaptive threshold ratio Θ as long as the latter is surpassed in each step (indicating that an instability is approached), until the strain increment drops below 10 −6 .The resulting state corresponds to our refined estimate of the saddle strain γs , represented by the red vertical line in the inset of Fig. S1 and the associated red-shaded region.We denote the corresponding state as x s .In addition, we push the system over the instability, which results in the final state (outcome of instability) marked in Fig. S1.We denote the corresponding state as x fin .
• In the third stage, the most refined one, we first solve for the lowest eigenvalue of the eigenvalue problem M(γ s )ψ = ω 2 ψ, where ψ is the eigenmode associated with the minimal frequency ω. ψ, being a harmonic eigenmode, may suffer from hybridization with phononic plane waves, as extensively discussed in the literature [82,96,97].We bypass this potential difficulty by employing a recently developed cubic nonlinearmodes framework that cleanly and effectively disentangles plastic instabilities from the phononic background, as detailed in [98,99].We use ψ as an input to the nonlinear mode calculation, resulting in a displacement field u(r), which is the basic quantity used in the manuscript to extract the properties of plastic instabilities.An added value of the nonlinear-modes framework, as discussed in [98,99], is that it converges more accurately to the instability, and hence effectively provides us with an even better estimate of γ s (compared to γs ), which was our overall goal.
In addition, we are interested in the sign of u(r), which allows to distinguish between isotropic expansion/dilation, ϵ * dil > 0, and isotropic core contraction, ϵ * dil < 0 (also rendering the plastic eigenstrain triaxiality R t in Eq. ( 1) and Fig. 2 in the manuscript a signed quantity).Therefore, we fix the overall sign of u(r) such that u•(x fin − x s ) > 0. The 3-stage procedure detailed above allows us to obtain a set of initial and unstable (saddle) states, where the latter are analyzed in the manuscript.Finally, the shear strain γ (and dilatational strain ϵ), as appearing on the x axis in Figs. 1 and 5 in the manuscript, correspond to the accumulated effect of the discrete strain steps involved in the procedure (the coarse and the adaptively refined ones).
Appendix C: Extracted quantities: elastic moduli and core properties through contour integrals For each initial state, we compute the shear and bulk moduli according to [57,100] Here, derivatives are taken with respect to γ, in Eq. (S1a), and to ϵ, in Eq. (S1b).Both γ and ϵ have the same geometric meaning as in the corresponding global affine transformations used in the actual deformation, though here they are used as putative motions invoked to evaluate linear response coefficients, independently of whether the actual deformation is simple shear or hydrostatic tension.d is the spatial dimension, which in our work is set to d = 3.
The elastic moduli are subsequently employed in the contour integrals method developed in [61], used to extract the the core properties V E * of each unstable plastic mode u(r).The contour integrals is presented in great detail in [61] and we apply it in the present work as is.For completeness, we briefly note that V E * are obtained via a specifically designed projection of the mode u(r)  For simple shear, we used 4039 modes for the early strain regime and 5728 for the developed strain regime, and for hydrostatic tension, we used 5007 for the early strain regime and 4992 for the developed strain regime.Overall, we do not observe any significant change in distributions, except for a mild change in the histogram in (b) for Rt under dilation, which is somewhat anticipated.
onto the spherical harmonics [61] I 0 (r) ≡ 2 √ π λ + 2µ 3λ + 2µ S u(r) • r Y 0 0 (Ω) r dΩ , (S2a) where λ = K − 2µ/3, the surface integral S is performed on a sphere of radius r, Y 0 0 (Ω) and Y m 2 (Ω) are the real orthogonal spherical harmonics of the zeroth and second degree and orders zero and m = −2, −1, 0, 1, 2 [101,102] and Ω is the solid angle.To obtain V E * we then construct the matrix Ī2 (r) as detailed in [61], and extract the core properties at a distance r where the largest eigenvalue of Ī2 (r) in absolute value attains its maximum.As explained in the manuscript, there is some variability in the number of plastic events used for the statistical analysis between different r c values and loading symmetries.The smallest data set corresponds to r c = 1.2 under hydrostatic tension, in which we have collected 3345 For simple shear, we used 4544 modes for the early strain regime and 5456 for the developed strain regime, and for hydrostatic tension, we used 544 for the early strain regime and 2801 for the developed strain regime.Overall, we do not observe any significant change in distributions, except for a mild change in the histogram in (b) for Rt under dilation, which is somewhat anticipated.
events for analysis.In this context, it is important to verify the statistical convergence of the statistical distributions by considering smaller subsets of each ensemble, as demonstrated in Fig. S2.
We also examined the plastic eigenstrain triaxiality ratio R t , planarity ratio R p , and core size a/a 0 as the strain increased throughout the deformation path.Figure S3 shows the obtained distributions for R t , R p , and a/a 0 originating from analyzing modes at the initial and largest-strain phases of loading (γ < 0.02 and γ > 0.02 under simple shear, and ϵ < 0.02 and ϵ > 0.02 under hydrostatic tension, respectively) for r c = 1.5.It is observed that most distributions are strain independent, except for a weak dependence of p(R t ) under hydrostatic tension (panel (b)).Figure S4 shows the same distributions for r c = 1.2.The results are similar to those for r c = 1.5, indicating weak strain dependence, yet with a somewhat larger effect in p(R t ) under hydrostatic tension (panel (b)).In general, we expect strain dependence to emerge under hydrostatic tension for larger deformation level, an issue to be addressed in future work.

FIG. 1 .
FIG. 1.(a) A shear stress-strain σ(γ) of a 3D computer glass (of N = 10K particles, averaged over 108 realizations) deformed under simple shear.σ is normalized by the initial shear modulus, µ(γ = 0).The left inset shows the binarymixture glass (blue and yellow particles correspond to the two species [55]) prior to deformation and the arrows illustrate the subsequent application of simple shear.The right inset shows a side view of the same glass realization in the steady-state flow regime (here γ ≃ 0.25).(b) The corresponding dilatational stress-strain −p(ϵ) for the same glass ensemble (p is the pressure) deformed under hydrostatic tension.−p is normalized by the initial bulk modulus, K(ϵ = 0).The left inset shows the very same undeformed glass realization as in panel (a), but the arrows highlight the subsequent application of hydrostatic tension (pure dilation).The right inset shows a large-scale cavity (light blue region) inside the glass right after the observed abrupt stress drop.See text for a discussion of the results.

FIG. 2 .
FIG. 2. (a) The probability distribution p(Rt) of the plastic eigenstrain triaxiliaty Rt defined in Eq. (1) for rc = 1.5 glasses, under both simple shear and hydrostatic tension (see legend).The inset presents the interatomic pair interaction ϕ(r) vs. r, where r is the scaled distance between interacting particles and r is further normalized such that the minimum of each curve occurs at r = 1 [55, 57].Curves for the two values of rc employed in this work (see legend) are shown.(b) The same as panel (a), but for rc = 1.2.See text for an extensive discussion of the results.The visual insets present iso-surfaces of the magnitude of plastic modes with two values of Rt.The left one (red), corresponds to Rt = 0 (i.e., ϵ * dil = 0) and a planar deviatoric eigenstrain tensor with J * 2 = 1.The effective dimensionality of the deviatoric eigenstrain tensor (quantifying its degree of planarity) is discussed in relation to Fig. 3.The arrows (white) indicate the direction of the displacement (with length that is consistent with its magnitude).The right visual inset (green) corresponds to ϵ * dil = 0.5 and again a planar deviatoric eigenstrain tensor with J * 2 = 1, resulting in Rt ≃ 0.29.The arrows (black) indicate the direction of the displacement (with length that is consistent with its magnitude).

FIG. 3 .
FIG.3.The probability distribution p(Rp) of the planarity ratio Rp of the deviatoric eigenstrain tensor, cf.Eq. (2), in glasses under simple shear, for two values of rc (see legend in panel (b)).The visual insets present iso-surfaces of the magnitude of the deviatoric part of plastic modes with two values of Rp.The left one (blue), corresponds to Rp = 0 (i.e., the purely planar limit), and is identical to the left visual inset in Fig.2b.The right visual inset (yellow), corresponds to Rp = 1 (i.e., the least planar).In both cases, the arrows indicate the direction of the displacement (with length that is consistent with its magnitude).(b) The same as panel (a), but under hydrostatic tension.See text for a discussion.

1 FIG. 4 .
FIG. 4. (a)The probability distribution p(a/a0) of the dimensionless linear size a/a0 of the core of unstable plastic modes, as defined in Eq. (3), for two values of rc (see legend in panel (b)) under simple shear.(inset) The amplitude of an unstable plastic mode ⟨|u|⟩ Ω (the solid angle Ω average of the norm of u(r)) vs. |r|/a0 in an rc = 1.5 glass.The estimated core size a/a0 is marked by the horizontal green double arrow and the linear elastic far-field power-law ∼ 1/r 2 is highlighted (red dashed-dotted line).(b) The same as panel (a), but under hydrostatic tension.See text for a discussion.

30 FIG. 5 .
FIG. 5. (a)The bulk K(γ) (dashed-dotted line) and shear µ(γ) (solid line) moduli under simple shear for an ensemble of rc = 1.2 glasses.The moduli, whose precise definition is given in[55], are reported in units of the interatomic interaction energy scale divided by the atomic volume a 3 0 .(b) The same as panel (a), but under hydrostatic tension.(inset) The reduced moduli, μ(ϵ) ≡ µ(ϵ)V (ϵ)/V0 and K(ϵ) ≡ K(ϵ)V (ϵ)/V0, under hydrostatic tension (green labels and curve).The corresponding moduli under simple shear (red curves and x label), whose values are identical to those presented in panel (a) since V (γ) = V0, are superimposed for comparison.See text for a discussion of the results.

8 FIG
FIG.S3.Distributions of (a-b) Rt, (c-d) Rp, and (e-f) a/a0, for rc = 1.5 at early and late stages of loading under simple shear (left column) and hydrostatic tension (right column).For simple shear, we used 4039 modes for the early strain regime and 5728 for the developed strain regime, and for hydrostatic tension, we used 5007 for the early strain regime and 4992 for the developed strain regime.Overall, we do not observe any significant change in distributions, except for a mild change in the histogram in (b) for Rt under dilation, which is somewhat anticipated.
Appendix D: Statistical convergence and weak strain dependence

2 FIG
FIG.S4.Distributions of (a-b) Rt, (c-d) Rp, and (e-f) a/a0, for rc = 1.2 at early and late stages of loading under simple shear (left column) and hydrostatic tension (right column).For simple shear, we used 4544 modes for the early strain regime and 5456 for the developed strain regime, and for hydrostatic tension, we used 544 for the early strain regime and 2801 for the developed strain regime.Overall, we do not observe any significant change in distributions, except for a mild change in the histogram in (b) for Rt under dilation, which is somewhat anticipated.

TABLE II .
Initial density and resulting pressure (normalized by the bulk modulus K) and Poisson's ratio ν in the quenched samples for the two rc values.