Pre-yield non-affine fluctuations and a hidden critical point in strained crystals
PPre-yield non-affine fluctuations and a hidden critical point in strainedcrystals † Tamoghna Das, a , b Saswati Ganguly, b Surajit Sengupta ∗ c and Madan Rao d ‡ a Collective Interactions Unit, OIST Graduate University,1919-1 Tancha, Onna-son, Okinawa, Japan -904-0495 b Centre for Advanced Materials, Indian Association for the Cultivation of Science, Jadavpur, Kolkata700032, India c TIFR Centre for Interdisciplinary Sciences, 21 Brundavan Colony, Narsingi, Hyderabad 500075, In-dia. email: [email protected] d Raman Research Institute, C.V. Raman Avenue, Bangalore 560080, India
A crystalline solid exhibits thermally induced localised non-affine droplets in the absence of external stress. Here we show thatupon an imposed shear, the size of these droplets grow until they percolate at a critical strain, well below the value at whichthe solid begins to yield. This critical point does not manifest in bulk thermodynamic or mechanical properties, but is hidden and reveals itself in the onset of inhomogeneities in elastic moduli, marked changes in the appearance and local properties ofnon-affine droplets and a sudden enhancement in defect pair concentration. Slow relaxation of stress and an-elasticity appear asobservable dynamical consequences of this hidden criticality. Our results may be directly verified in colloidal crystals with videomicroscopy techniques but are expected to have more general validity.Mechanical properties of solids [1, 2], especially mecha-nisms of yielding in response to external stress, continue to en-gage the attention of materials scientists and engineers [3]. Pre-yield phenomena such as anelasticity, occurring at non-zerotemperatures and below the yield stress, in both crystals [4]and metallic glasses [5, 6] is somewhat less understood, how-ever. Within this anelastic regime, a crystalline solid undergoesrecoverable strain but with a long relaxation time and is accom-panied by conspicuous production and reorganisation of latticedefects [7, 8] marked by the onset of significant departure frominstantaneous and linear, “Hooke’s law”, elasticity. Mechani-cal response becomes heterogeneous ; the regions with defectsbehaving differently from the rest of the solid. In amorphoussolids too, anelasticity is marked by the appearance of similarheterogeneous mechanical response confined to localised clus-ters of particles undergoing large non-affine deformation or nonaffine droplets [9, 10, 11].In this paper, we focus on defect nucleation and associatedpre-yield phenomena in two (2D) and three-dimensional (3D),initially homogeneous, crystals of interacting particles. Weprovide strong evidence associating the proliferation of crys-talline defects during deformation at non-zero temperatures,with the critical behaviour of thermally generated non-affine † Electronic Supplementary Information (ESI) available: [attached file: Sup-plementary text] ‡ Also at National Centre for Biological Sciences (TIFR), Bellary Road, Ban-galore 560065, India droplets. These non-affine displacements exist in crystallinesolids even in the absence of an applied strain and have a purelythermal origin [12]. Upon shearing the crystalline solid, theseclusters grow and eventually percolate at a hidden mechanicalcritical point; concomitantly, defects suddenly proliferate. Sig-nificantly, this happens well before the yield point of the crys-talline solid. Apart from its intrinsic novelty, this viewpointsheds insight onto a broader set of issues. First, we find thatdislocation nucleation is stochastic and seeded from regionswith large non-affine displacements. Second, we have iden-tified the localised non-affine regions with droplet excitationsfrom nearby metastable liquid/glass [13]. Thus, our finding thatthe percolation of non-affine droplets coincides with the dis-appearance of the metastable liquid-glass spinodal, brings outthe significance of metastable configurations in understandingmechanical properties of solids. Last, this perspective mightprovide a language bridging ideas concerning the mechanicalresponse of crystals and amorphous solids.The existence of a unique reference in crystals allows usto identify localized non-affine regions at non-zero tempera-tures, even at zero applied strain. In Ref.[13] we studied thestatistics of shape and size of non-affine clusters and their lo-cal thermodynamic properties using molecular dynamics (MD)simulations of a two-dimensional Lennard-Jones (2D-LJ) solid[14, 15]. We summarise the main results below for complete-ness.For every configuration, particles i undergoing large non-1 a r X i v : . [ c ond - m a t . m t r l - s c i ] F e b ig. 1: (color-online) a - c Configurations from a part of our constant
NV T , MD simulation cell of 22500 2D-LJ particles at T = . ρ = . ε = a ). Shading (light to dark) tracks local χ . Non-affine particles ( χ > χ cut ), organised into droplets, are black. Strain increasesnon-affine particles; ε = .
04 ( b ) and droplets percolate at ε = .
052 ( c ). Percolating droplet in yellow (dark grey) with a compact node circledgreen (light grey). affine displacements are identified using the parameter χ i viz.the least square error incurred in trying to fit an instantaneouslocal volume, Ω , to an affine distortion of the same volumein the, reference, undistorted lattice [11]. A cutoff criterion χ i > χ cut then eliminates trivial harmonic distortions of the lat-tice (see Methods and supplementary information SI). Theseparticles are observed to cluster together into droplets. Thedroplets are characterized by a distribution of the local den-sity ρ c and excess pressures, ∆ p c ≡ p c − p . Here p is the meanpressure of the surrounding solid of density ρ and p c is com-puted by averaging the virial over the n c particles of the clus-ter. For any temperature T the mean scaled excess pressure¯ p = ∆ p c / T vs ρ c curve is non-monotonic and resembles a vander Waals loop seen in usual liquid-gas transitions. The two sta-ble branches with ∂ ¯ p / ∂ ρ c > “unstable” branch where ∂ ¯ p / ∂ ρ c <
0. In the unstrained solid the stablebranches were associated with inflated, compact clusters or de-flated, string-like clusters depending on whether ¯ p is negativeor positive respectively [16, 17]; the unstable branch containsonly ramified “branched-polymer” clusters [18]. By studyingthe local equation of state and density-correlations, one mayassociate compact and string-like clusters with droplet fluctua-tions from nearby liquid and glass -like metastable free-energyminima respectively [13]. As T is increased, the van der Waalsloop as well as the distinction between compact and string-like droplets vanish above a metastable critical point. Thisbehaviour is quite general and exists in other model systemsadmitting non-crystalline phases, (see discussion in SI).Here, we extend this earlier analysis [13] to the mechanicalresponse of LJ solids in two and three dimensions over a rangeof ρ , and for a few T . In what follows, the bulk of our quan-titatively detailed results, unless otherwise stated, are explic-itly demonstrated for the 2D-LJ solid for which computations,as well as experimental verification using real-time video mi- croscopy [19], are relatively cheap. Key results for the 3D-LJcase [20] are also presented showing that our main conclusionscarry over to higher dimensions. Results
Percolation of clusters
As the crystal is subject to a quasi-static (pure) shear strain ε = ε xx − ε yy (Fig.12a-c) (see Methods ), localized non-affinedeformations grow and increase in number. At a critical valueof the strain ε ∗ ( e . g . ≈ .
05 for ρ = . , T = . ) , regions con-taining non-affine particles begin to percolate. Further increaseof the strain beyond ε ∗ , eventually leads to yielding of the solid(at ε ≈ . ρ , T ) .To ascertain the nature of this percolation transition, we plotthe ratio of the number of particles in the largest cluster to thetotal number of non-affine particles, f φ in the 2D-LJ solid as afunction of the fraction of non-affine particles φ ( ρ , ε ) in Fig.2a.The collapse of f φ onto a single curve for all ρ and ε , is astrong indication of a true critical percolation transition at avalue of φ ∗ (cid:39) .
4, close to the known value of site percolationin 2D (0 .
34) [21]. We plot the probability distribution P ( n c ) ofobtaining a non-affine cluster of size n c for increasing valuesof ε in Fig.2b. As one approaches the percolation transition, P ( n c ) ∼ n − τ c with an exponent τ ≈ . ig. 2: a Scaling plot of f φ vs. φ (see text) for various ρ and ε showing data collapse for the 2D-LJ solid. The percolation threshold φ ∗ (cid:46) . Inset
Plots of φ vs ε for various ρ . b Probability P ( n c ) for ρ = .
99 and ε = ,. ,.
048 and .
052 showing crossover fromexponential to power-law behaviour at percolation ε ∗ ≈ . τ = − .
67 of P ( n c ) = n − τ c shown by dotted straight line.Fig. 3: a Plot of ¯ p vs. ρ c for various φ (key on right) and n c = φ have been shiftedin ¯ p -axis for clarity. Note that the loop disappears for φ > φ ∗ . b Plotof the spinodal obtained for n c = ,
15 and 20 in the φ − ρ c planeclearly indicating a metastable critical point coinciding with φ ∗ . Local thermodynamics
Bulk thermodynamic quantities such as equation of state etc.show no abrupt change across the transition. We, therefore,focus on measuring local thermodynamic variables within thenon-affine clusters and monitoring them across the percolationtransition. We characterise droplets using the local quantities¯ p and ρ c as in [13] (see also Methods ). The local tempera-tures, on the other hand, are close to the fixed temperature ofthe solid and are irrelevant to the subsequent analysis. Insteadof temperature here strain, ε , drives the percolation transition.We anticipate that systems with various ε and ρ behave sim-ilarly for equal φ . We test this as follows. From configurationswith fraction of non-affine particles between φ and φ + ∆ φ ,where ∆ φ is a suitable bin-size, we obtain the mean ¯ p and ρ c corresponding to droplets with fixed number of particles 10 < n c <
100 ; collecting together data from all solids with0 . ≤ ρ ≤ .
99 and ε below the yield point. This data is plot-ted in Fig.3a for droplets of size n c =
15 for each value of φ .Similar plots are obtained for droplets with other values of n c too.Over a large range of φ values, the plots are non-monotonic,showing prominent van der Waals loops with two stablebranches. We show below that these corresponding to com-pact and string-like droplets. As φ increases, the loops tendto vanish. Within the accuracy of our computations, the valueof φ above which the loops vanish is the same as φ ∗ , i.e. thenon-affine fraction at which the clusters percolate. To makethis identification quantitative, we obtain the values of ρ c forall φ at the spinodal ∂ ¯ p / ∂ ρ c = n c = ,
15 and 20. When extrapolated, the spinodal lines in-tersect at a metastable, droplet critical point. This metastablecritical point and the droplet shape transition survives for evenhigher values of n c , though for very large n c , the statistics be-comes sparse. Fig. 4: a Plot showing ∆ p c as a function of the local density ρ c for“non-affine” droplets in the strain free LJ solid at T = .
4. The dif-ferent colors show the van der Waals loop for various sizes of clusters10 < n c <
30. The data becomes sparse for much higher values of n c ,leading to worse statistics. The van der Waals loops however show sta-bilisation with larger n c . We have made this quantitative by obtainingan effective barrier height ∆ F c for this transition by integrating thepressure difference ∆ p c along ρ c for each of these curves. The resultis plotted in b and shows that ∆ F c increases with n c as expected of afirst order transition. The solid line is a linear fit. How does the extent of the van der Waals loop in ¯ p vs. ρ c depend on the size of the cluster n c ? In Fig. 4a we show acomparison of the van der Waals loops for various values of n c in the ε → p vs ρ c which gives ameasure of the free energy F c expended by the solid to pro-duce these fluctuations has a convex up region, showing that aninterpretation in terms of a metastable first order transition isnot inconsistent. A metastable first order transition implies, atleast for compact droplets, that the surface free energy barriersshould scale as l d − c , where l c is a typical linear size associatedwith the clusters and d (= ) is the dimensionality [22]. The3ctual form of the scaling may be more complicated because ofthe complex shape of the clusters. Nevertheless, in Fig.4b weshow a plot of the barrier height ∆ F c , obtained by integratingFig.4a, as a function of the size of the fluctuation. A monotoni-cally increasing curve indicates a positive surface energy of thedroplets at ε =
0. At the metastable critical point we expect thissurface energy to vanish, although large fluctuations and strongfinite size effects makes an explicit evaluation of this quantitynear φ = φ ∗ computationally impractical. The shape transition
We have shown that the van der Waals loop in ¯ p vs ρ c van-ishes as φ ≥ φ ∗ where the non-affine droplets percolate. Thevan der Waals loop corresponds to a metastable transition be-tween droplets which come in two stable shapes, either com-pact or string-like. We show now that this distinction betweenthese shapes also vanishes for φ > φ ∗ . Above this transition,string-like droplets connect the adjacent compact “nodes” in asingle giant percolating network. This network has an overall,branched polymer shape so that isolated compact and string-like droplets cease to exist beyond φ ≥ φ ∗ (see Fig.12 c). Fig. 5: The radius of gyration R g as a function of the size n c of thedroplets where we have factored out the ∼ n ν c , with ν = .
64 (branchedpolymer) dependence. Below the percolation transition, a , the dropletsfor the ¯ p > p < ν ∼
1) and compact ( ν ∼ /
2) configurations respectively. Abovepercolation b , the morphology of droplets changes drastically and alldroplets become branched polymers. For the unstrained solid [13] the ¯ p < p > R g ( n c , ∆ p c , T ) = n ν c F ( ¯ pn ν c ) , where F ( x ) is a crossoverfunction. The value of ν is 0 .
64 for branched polymers andthe − ve and + ve ¯ p branches crossover to ν = . a we recover this crossover of the dropletshapes for φ ≤ φ ∗ for the two corresponding branches of thevan der Waals loop as in the unstrained solid. As φ increases, however, the giant percolating cluster dominates the dropletconfiguration space and it is impossible to obtain compact orstring-like droplets; all droplets now show branched polymerbehaviour (Fig.5 b ) with R g ∼ n . c typical of a system abovea critical point. This disappearance of the transition in shapescoincides with the disappearance of the van der Waals loop in¯ p − ρ c plane. Nucleation of defects
Non-affine droplets are associated with coordination numberchanging deformations. To show this, we track the local con-centration of defect pairs by counting the number of nearestneighbours of particles using a local Delaunay triangulation(see Fig.6a inset). In Fig.6a, we plot the density of defect pairs ρ d as a function of the external strain for a LJ solid. While asmall number of dislocation pairs (dipoles) always exist withina solid at finite temperatures[23], there is a sudden increase in ρ d as soon as ε crosses a critical value. This increase in defectpair concentration coincides with the critical percolation transi-tion of the non-affine clusters. That this is not a mere numericalcoincidence is obvious from the scaling collapse of ρ d for all ρ and ε onto a single curve when plotted against the non-affinenumber fraction φ in Fig.6b.We have obtained identical results also at T = .
35 with datacollapse occurring for both ρ and T implying that Fig. 6b rep-resents an universal relation between ρ d and φ for all ρ , T and ε for the 2D-LJ solid. A similar phenomenon occurs for the3D-LJ solid as well. We discuss these results later in this work.Dislocations are preferentially nucleated within the perco-lating cluster. Also, each defect pair is dressed with an ex-tended region of non-affine particles contributing to the system-spanning, non-affine cluster. We show this by introducing a dis-location pair deleting a row of 20 atoms in a triangular crystalof 100 ×
100 LJ atoms at ρ = .
97. After relaxation to obtaina reference configuration, χ is calculated during a MD run at T = . Elastic nonlinearity, heterogeneity, anelasticity
How does this percolation transition affect mechanical proper-ties of the solid ? Increased non-affine fluctuations at defectsites is expected to reduce the local elastic modulus [24]. Theaverage elastic modulus therefore involves ρ d , the defect con-centration and therefore decreases significantly from the linearHooke’s law behaviour as soon as ρ d becomes large [25]. Thisresults in three observable features in the mechanical responseof the solid.4 ig. 6: a Plot of the defect concentration ρ d (particles with 7 Delau-nay neighbours) vs. ε in the 2D-LJ solid at ρ = .
99 showing sharpgrowth of ρ d for ε > ε ∗ . Inset shows a particle neighbourhood attwo different values of ε = .
039 (open circles) and .
068 (filled graycircles), showing the nucleation of a defect pair (black filled circles).Dashed and the bold lines shows change in topology of the local De-launay neighbourhood with strain. b Plot of ρ d vs. φ for all ρ and ε showing data collapse onto a single curve with ρ d increasing sharply at φ ∗ . Inset shows the same data in semi-log axes to emphasise the small φ region. Symbols: same as in Fig.2a. c Particles coloured (colourkey) according to the value of the local non affine parameter χ for anisolated dislocation-anti-dislocation pair in a triangular crystal. d Thesame configuration as in c showing non-affine particles with χ > χ cut (red filled circles). Majority of non-affine particles are associated withthe dislocation dipole (black arrows). The first concerns the appearance of a slight nonlinearity inthe stress-strain relation of the bulk solid. In Fig.7a we showquasi static stress-strain curves for the LJ solid at T = . . < ρ < .
99 scaled by σ ∗ and ε ∗ the stress and strain re-spectively at percolation. Notice that for all ε ≤ ε ∗ , the curvescollapse to the trivial, linear, response while for larger strainsthe response becomes non-linear. The non-linearity ∆ definedas the deviation of the stress from a linear fit to the small ε behaviour becomes significant at precisely ε ∗ where the non-affine droplets percolate. Similar behaviour, described later, isobserved for the 3D-LJ solid beyond the percolation point.The other mechanical feature that manifests at the percola-tion transition, is the onset of elastic heterogeneity in the solid.As the number of defects, characterised by pairs of 5 − Fig. 7: a Scaled stress σ / σ ∗ vs. ε / ε ∗ in the 2D-LJ solid at T = . ρ . Inset shows non-linearity ∆ vs. ε . The vertical dashedline marks ε ∗ . b Distribution of local stress, P ( σ B ) obtained from40 ×
40 sub-blocks. At small strains the distribution is Lorentzian.As ε crosses ε ∗ , however, spatially correlated dislocation structuresdevelop making P ( σ B ) assymetric and broad. Inset shows variation ofskewness (triangles) µ = (cid:104) ( σ − ¯ σ ) (cid:105) / Σ and kurtosis (circles) µ = (cid:104) ( σ − ¯ σ ) (cid:105) / Σ − P ( σ B ) with ε . Note that µ changes sign closeto ε = ε ∗ ; µ = µ = σ and Σ are the mean andstandard deviation of P ( σ B ) .Fig. 8: a. The pair distribution function of 7 − coordinatedparticles for a ρ = .
99 2DLJ solid at T = . φ ≈ . b. Anelastic re-sponse of the solid at ρ = .
99. A few strained configurations( ε = . , . , . , . , . , . σ ( t ) over time t . The stress responds by first decreas-ing rapidly to a negative value whose magnitude increases with theinitial ε and then slowly relaxing to zero from below. For the min-imally strained solid, the relaxation is essentially exponential. If theinitial strain is around ε ∗ , the dynamical response shows a second re-laxation which is complex and non-monotonic. The inset shows thedependence of the average stress σ a during the relaxation process onthe value of the initial ε before σ ( t ) decreases to 0. While σ a growslinearly for small ε , its growth becomes rapid as ε crosses ε ∗ percolation with defects getting nucleated preferentially withinthe percolating cluster, one expects that defects also get spa-tially correlated. Correlated defects cause the local elastic re-sponse to be different from the bulk. This may be seen from a5ocal elastic analysis, Fig.7b. Here we plot the local stress dis-tribution P ( σ B ) by dividing our simulation cell into rectangularsub-blocks [26] and obtaining local stresses σ B by averagingthe virial over the block. While the stress distribution of a per-fect crystal is expected to be a Gaussian, that arising from asmall concentration of spatially uncorrelated ideal dislocationdipoles is known to be a Lorentzian, P ( σ ) = π − σ / ( σ + σ ) with σ = D ρ dip / π where D is a constant proportional to thestrength of the dipoles, ρ dip is the dipole density[27]. Our solidat T = . ρ d ∼ − evenat ε =
0. As ε increases, so does ρ d making P ( σ B ) broader. At ε > ε ∗ , however, there is a dramatic increase of ρ d and largenumbers of correlated defect pairs are produced. The transfor-mation of P ( σ B ) can be quantified by calculating the 3 rd and4 th order moments which show that at the percolation point thestress distribution changes from being a sharp Lorentzian to aflat, “sub-Gaussian” form. Position correlation between defectscan also be seen directly from the pair distribution function of7 − coordinated particles for solids. We have plotted this quan-tity for three values of φ across the transition in FIg.8. Theincrease of correlations is obvious.Lastly, crystalline solids, as well as bulk metallic glasses, of-ten show complex, time-dependent, stress relaxation behaviourat values of external strains much below the yield point with-out permanent plastic deformation. To study the an-elasticresponse of the LJ solid, we quench a few configurations ofour LJ solid, which were equilibrated at several values of ε ,by rescaling the particle coordinates and the boundary to zerostrain. Once initial transients die out, the solid enters a regimeof slow relaxation. The nature of this relaxation of the stress, σ ( t ) , which is driven by the motion and subsequent annihila-tion of defect pairs of opposite sign, depends on the numberof such defects and therefore on the initial strain. If the ini-tial strain (and consequently, the defect density) is small, thisregime is not very prominent. On the other hand, if the initial ε is large, there are a large number of defects which reorgan-ise over a much longer timescale. In this case, therefore, thesecond relaxation, as shown in Fig.8b is complex and showsnon-monotonicity. We characterise this relaxation by plottingthe time average of the stress σ a = τ − (cid:82) τ dt σ ( t ) for variousvalues of the initial strain and observe a sharp change of slopeat ε ∗ . Since defect pairs are formed at the percolation criticalpoint, the origin of non-monotonic relaxation is traced to thepercolation transition at φ ∗ . The anelastic regime, grows withlowering ρ (and increasing T ) as expected. Effect of temperature
The percolation threshold depends on φ the fraction of non-affine particles and not separately on temperature, density andstrain. The main thesis of our work is that mechanical proper-ties, such as the defect nucleation threshold, is determined by percolation of non-affine droplets. Therefore, mechanical prop-erties of the solid should depend on external parameters such astemperature, density and strain only through the dependence of φ on these parameters. We have shown this in Fig. 2a and bfor the combination of density and strain at a fixed tempera-ture. It is clear that φ decreases with increasing density andincreases with strain. The percolation threshold at fixed tem-perature therefore moves to higher strain with increasing den-sity and lower density with increasing strain. If, on the otherhand, the temperature is lowered at fixed density one should goto higher strains to achieve percolation. Similarly, the percola-tion threshold, for fixed strain, should moves to lower densitiesat lower temperatures. To put to test the validity of our results Fig. 9: a. Plot of the non-affine fraction φ ( T , ρ ) for two temperatures T = . . ρ = . .
99 for the 2DLJ solid. b. Datacollapse plot of ρ d , the defect density, plotted against the non-affinefraction φ . Note that the data for T = . at other temperatures, we have computed ρ d for a 2DLJ solidwith ρ = .
93 and 0 .
99 at both T = .
35 and T = .
40. Weshow that φ ( T , ρ ) in Fig. 9a behaves exactly according to ourexpectations. We also show in Fig. 9b that one can collapsethe data for ρ d as the solid is quasi-statically strained at differ-ent temperatures and densities using the non-affine fraction φ .The relation shown in Fig. 6b is therefore, indeed, an univer-sal curve, valid for all densities, strains and temperatures. Weshow below that the 3D-LJ solid behaves similarly. Extension to three dimensions
We have shown that in the 2D-LJ solid at any density andtemperature, defect nucleation happens at a particular value ofstrain ε ∗ ( ρ , T ) where non-affine clusters percolate. The me-chanical response of the solid σ ( ε ) is trivial below ε ∗ and thestress-strain curves collapse onto each other when scaled by ε ∗ and σ ∗ . The cause for nontrivial mechanical response is thenucleation of defects which increase rapidly for ε > ε ∗ whereclusters containing non-affine particles percolate. Are our re-sults valid only for the 2D-LJ solid or does this have more gen-eral validity? To check this, we have looked at the 3D-LJ solid6 ig. 10: a - c Non affine particles in 3D-LJ, FCC solid of 16384 particles at T = . ρ = . ε = . ε = .
04 and ε = .
06 respectively.Particles with less than 2 non-affine neighbours have been made translucent for clarity. Percolation occurs for ε (cid:38) .
04 at this density where φ (cid:39) . which has a close packed FCC structure. Since fluctuationsgenerally decrease with increase in the number of dimensions,a na¨ıve argument may suggest that non-affine fluctuations maynot have much of a role in determining the mechanical responseof the 3D-LJ solid. The situation is however more subtle, sincein higher dimensions there are more modes that are non-affine .For example, choosing a neighbourhood which includes onlynearest neighbours, the triangular lattice with 6 neighbours,features 4 affine (volume, uniaxial extension, shear and rota-tion) and 8 non-affine modes [12], while in three dimensions,the FCC lattice with 12 neighbours has 9 affine and 27 non-affine modes. It is possible that the larger number of availablenon-affine modes more than compensate for the reduction offluctuations due to the increased dimensionality.To check this we simulate the 3D-LJ solid (see Methods ) atLJ reduced temperatures of T = . . ρ = . .
2. The fraction of non-affine particles φ is obtained using aprocedure analogous to what we use in 2D. The solids are sub-jected to deviatoric strain ε = ε xx − ε yy − ε zz and the conjugatestress σ , φ and f φ obtained as in 2D. As expected, φ increaseswith ε for fixed ρ and T and decreases at fixed ε as either ρ is decreased or T increased similar to the 2D case; f φ vs. φ is an universal curve with a percolation transition at φ ≈ . ε ∗ and σ ∗ , the crit-ical strain and conjugate stress at percolation and use this toscale the stress-strain curve shown in Fig. 11b. Again, as in2D, the stress-strain response becomes non-trivial only beyondthe point that non-affine droplets percolate.Before we end this section, we point out a remarkable as-pect of the data shown in Fig.11. In 3D, identification of thedominant topological defect is non-trivial. In the FCC lattice,lattice dislocations, frequently decompose into partials produc-ing stacking faults, making a computation of ρ d from latticeconfigurations much more difficult that the simple counting ofnear neighbour coordination which we did in 2D. The nucle- ation of stacking faults governed by their energy determines thedominant deformation mechanism in FCC crystals. Thereforeconsiderable effort is needed to define and compute stackingfault energies [3]. Our results may offer an alternative. Wehave shown that the non-affine fraction φ provides an equallygood measure for significant departure of the lattice from thetopology of the ideal reference and may be used to understandpre-yielding phenomena in higher dimensional solids too. Ofcourse, a computation of φ and f φ needs data for instantaneousparticle positions which may be available only in computer sim-ulations or from colloidal solids. Fig. 11: a The percolation transition for non affine particles in the3D-LJ, FCC solid at two values of T = . . ρ = . . Inset shows the corresponding φ ( ρ , T ) curves. b Scaled stress-strainrelations for the sheared solids are shown in a ; symbols have the samemeaning as in a . Non-affine particles percolate the solid at ε / ε ∗ = Inset shows the σ vs ε curve for T = . ρ = .
5. The arrowmarks ε ∗ and the solid line marks the linear Hooke’s law behaviour. Discussion
In summary, we have unearthed a hidden mechanical criti-cal point associated with the percolation of non-affine droplets7hich is intimately tied to the onset of complex mechanical re-sponse in a crystalline solid. This finding is quite unexpected ina system as familiar as a crystalline solid. This transition doesnot manifest in mechanical or thermodynamic properties in thebulk, but subtly reveals itself in at least four ways :(i) Dislocations, nucleated upon shear, are associated witha spatially extended non-affine cloud . The spatial scale is de-termined by two coarse-graining parameters χ cut and Λ (see Methods ). It is this cloud of non-affinity that eventually perco-lates across the sample, the percolation transition being fairlyrobust to small changes in χ cut and Λ .(ii) The mechanical critical point is not apparent in bulk ther-modynamics but in the restricted (finite size) thermodynam-ics of the localized droplets. The droplet fluctuations repre-sent regions where the solid explores nearby minima in the freeenergy corresponding to metastable glassy or liquid configu-rations [28]. At an algorithmic level, since the droplets aresmall and transient, we require a special algorithm to distin-guish them from the overwhelmingly large contribution of nor-mal fluctuations in the equilibrium solid.(iii) This mechanical critical point is not reflected in the bulkelastic properties of the solid, but shows up in the distribution oflocal elastic stresses and local elastic constants and the emer-gence of shows strong spatial heterogeneity invisible in bulkelastic experiments.(iv) Finally, the mechanical critical point is linked to theonset of a strong nonlinear dynamical response of the crys-talline solid. Hysteretic response of the stress during strain cy-cling, slow creep relaxation of the strain under constant loadand relaxation of stress below yield point are all manifesta-tions of so-called “an-elastic” behaviour [29] caused by non-linearities in the elastic response, elastic heterogeneities andtime-dependent, slow, defect reorganization [4]. These featuresare precisely what characterizes the hidden mechanical criticalpoint.The fact that the emergence of plastic behaviour in both crys-talline and amorphous solids is associated with the percolationof localised non-affine deformations[30], suggests that theremight be a common language which describes the mechani-cal response of solids in general. Indeed a similar percolationtransition for quasi-statically strained amorphous solids in theathermal ( T =
0) limit has been recently described in [31].The hidden critical point betrays itself in our MD simula-tions only through local properties of non-affine droplets iden-tified using carefully chosen cut-offs and thresholds. Before weend, we speculate on whether the critical point may be revealedby changing system parameters or by introducing novel forcesso that it begins to affect bulk thermodynamics [32]. We haveshown in Ref. [12] that an external field h χ which couples to χ may be tuned to control the equilibrium value of χ and hencethe defect density. This field may be introduced in MD sim-ulations as well as realised experimentally in a colloidal solid using holographic optical tweezer techniques[33]. The resultspresented here for h χ = h χ − σ space occurring for h χ (cid:54) = Methods
Simulation details
Two dimensions:
We prepare a triangular arrangement of22 ,
500 identical particles interacting pairwise via the LJ potential.Equilibration of our 2D model system is carried out at the fixed T using the dissipative particle dynamics (DPD) protocol implementedwithin the LAMMPS package (http://lammps.sandia.gov) and pre-serving both volume and number. Solids at each ρ and T are equi-librated with 10 molecular dynamics time-steps (MDS) where eachtime-step equals 10 − LJ time units. Temperature fluctuations mea-sured at equilibrium is of the order of 1 in 10 . Pure shear is thenapplied to the system following the protocol: L (cid:48) x = ( + e ) L x and L (cid:48) y = ( − e ) L y ; L x and L y denoting the box lengths in X - and Y -directions respectively with primed quantities representing the sameafter shear, e is the shear step set to 0 . × MDS. Configurations for the analy-sis is stored in the last 1 × MDS in regular intervals of 20 MDS.This cumulative shearing process is continued till the solid fails at theyield point ε c . The strain rate associated with the process thereforeequals to 10 − per LJ unit. The sound velocity propagating throughthe solid c s = (cid:112) G / ρ where G is the elastic modulus for pure shearand ρ is the density of the medium. For ρ = .
99, ( L x = .
96 and L y = .
51 in LJ units), sound at velocity c s ∼ .
12 in LJ units takesapproximately 27 LJ time unit to travel across the simulation box in X -direction. Considering this, we refer to our shear protocol as quasi-static . The model solid at various densities is strained using the samequasi-static shear protocol. Three dimensions:
The number of particles is 16384 i.e. 16 × ×
16 cubic unit cells of the FCC lattice, each containing 4 basisatoms. The same LJ interactions and the same ensemble with fixednumber, volume and T (using the DPD thermostat) are used. Thesolid is in pure shear with L x being extended by ε , and L y and L z shortened such that the volume (and density) is preserved. The strainstep is 0 .
002 as before. At every strain step the solid is equilibratedfor 1 . × MD steps with a time step of 1 . × − in LJ units. Datais collected for the next 5 . × at every 500 steps. Identification of droplets
We identify non-affine droplets using the a prescription described in[13]. Briefly, at a fixed state point (fixed T , ρ and ε ), we calculate χ for each particle within a neighbourhood Ω , defined by r i j ≤ Λ forall configurations using the T = ρ and ε as reference. We next obtain the distribution P ( χ ) . We identifya fixed cutoff χ cut above which the fluctuations are deemed anhar-monic. All particles in a single snapshot having χ > χ cut are tagged nd tagged particles residing within the first nearest neighbor shell ofanother tagged particle belong to the same cluster; any cluster hav-ing at least 7 particles is identified as a non-affine droplet. Note that Λ and χ cut are parameters which are chosen such that Λ = . . χ cut excludes 90% of the weight of the first peak in P ( χ ) . Small deviations of these parameters do not change our resultssubstantially. Local thermodynamic quantities
We describe below our procedure for obtaining the local thermody-namic quantities for the droplets.
Local density ρ c : Once a droplet has been identified, weinclude one extra layer of particles so that the Voronoi near-est neighbors of all the particles within the cluster are in-cluded (whether non-affine or not). Then we count thenumber of Delaunay triangles using
GEOMPACK (code avail-able at http://people.sc.fsu.edu/˜jburkardt/f_src/geompack1 ). The correct density of the droplet is then given byhalf of the ratio between number of particles n c in the droplet and thetotal area of all the triangles ρ c ≡ n c / ( A c ) . Local pressure p c : This is computed using the virial (cid:104) F i j · r i j (cid:105) where F i j and r i j are the nearest neighbor forces and distances respec-tively for particles i and j belonging to the droplet. The average (cid:104) ... (cid:105) isover all n c . The droplets are characterized by a distribution of ρ c andexcess pressures, ∆ p c ≡ p c − p , where p is the mean pressure of thesurrounding solid. Block stress σ B : To compute σ B , we divide the simulation cellinto sub-blocks B of equal size. The stress is computed by averagingthe virial for all particles N B in the blocks, over all the blocks and fi-nally, over many independent configurations. The components of thestress distribution are calculated as ( σ B ) αγ = ( / N B ) ∑ k , l ∈ B ( x ( l ) α − x ( k ) α ) f ( kl ) γ where k , l are particle indices while α , γ denotes the carte-sian components of the stress. Acknowledgments
We thank F. Spaepen, P. Sollich and S. Kar-makar for discussions. TD acknowledges support from the Col-lective Interactions Unit of the Okinawa Institute of Scienceand Technology Graduate University. SG thanks CSIR (India)for support from a Senior Research Fellowship.
Author contributions
MR and SS conceived of the project andwrote the paper. TD did the necessary calculations, designedthe specific algorithms and helped in writing the paper. SGchecked the calculations and helped in designing algorithmsand obtaining the data.
Additional information
References
Mechanical Metallurgy , (McGraw-Hill, NewYork, 1961)2 Cahn R. W. & Haasen P.,
Physical Metallurgy t h Ed.(North-Holland, Amsterdam, 1996)3 Jo, M. et. al. Theory for plasticity of face-centered cubic met-als., Proc. Natl. Acad. Sc. , 6560-6565 (2014).4 Nowick, A. S. & Berry, B. S.
Anelastic relaxation in crys-talline solids , (Academic Press, New York, 1972)5 Luborsky, F. E. eds.
Amorphous Metallic Alloys (Butter-worths, London, 1983)6 Tomida, T. & Egami, T. Molecular-dynamics study of struc-tural anisotropy and anelasticity in metallic glasses, Phys.Rev. B , 3048-3057 (1993).7 Hirth, J. P. & Lothe, J. Theory of Dislocations , (Kreiger, Mal-abar, 1982)8 Nowick, A. S. & Heller, W. R. Dielectric and anelastic relax-ation of crystals containing point defects., Adv. in Phys. ,101-166 (1965)9 Argon, A. Plastic deformation in metallic glasses., Acta Met. , 47-58 (1979).10 Spaepen, F. A microscopic mechanism for steady state in-homogeneous flow in metallic glasses., Acta Met. , 407-415(1977).11 Falk, M. L. & Langer, J. S. Dynamics of viscoplastic de-formation in amorphous solids., Phys. Rev. E, , 7192-7205(1998).12 Ganguly, S. Sengupta, S. Sollich, S. P. & Rao, M. Nonaffinedisplacements in crystalline solids in the harmonic limit.,Phys. Rev. E , 042801(2013).13 Das, T. Sengupta S. & Rao, M. Nonaffine heterogeneitiesand droplet fluctuations in an equilibrium crystalline solid.,Phys. Rev. E , 041115 (2010).14 Frenkel, D. & Smit B. Understanding Molecular Simula-tions (Academic Press, San Deigo 2002).15 Barker, J. A. Henderson, D. & Abraham, F. F. Phase di-agram of the two-dimensional Lennard-Jones system; Ev-idence for first-order transitions., Physica , 226-238(1981).16 Leibler, S. Singh, R. R. P. & Fisher, M. E. Thermodynamicbehavior of two-dimensional vesicles., Phys. Rev. Lett. ,1989-1992 (1987).17 Maggs, A. C. Leibler, S. Fisher, M. E. & Camacho, C. J.Size of an inflated vesicle in 2 dimensions., Phys. Rev. A ,691-695 (1990).18 Hsu, H.-P. Nadler, W. & Grassberger, P. Statistics of latticeanimals., Comput. Phys. Commun. , 114 (2005).19 Zahn, K. Wille, A. Maret, G. Sengupta, S. & Nielaba, P.Elastic Properties of 2D Colloidal Crystals from Video Mi-croscopy., Phys. Rev. Lett. , 155506 (2003).20 Hansen, J.-P. & Verlet, L. Phase Transitions of the Lennard-9ones System., Phys. Rev. , 151-161 (1969).21 Stauffer, D. & Aharony, A. Introduction to Percolation The-ory nd Edition, (Taylor and Francis, London, 1994).22 Binder, K. Theory of First-Order Phase Transitions., Rep.Prog. Phys. , 783-859 (1987).23 Sengupta, S. Nielaba, P. & Binder, K. Elastic moduli, dis-location core energy, and melting of hard disks in two dimen-sions., Phys. Rev. E , 6294-6301 (2000).24 Zhou, C. et al. Dislocation-induced anomalous softening ofsolid helium., Phil. Mag. Lett. , 608-616 (2012).25 Hashin, Z. & Shtrikman, S. A variational approach to thetheory of the elastic behaviour of multiphase materials., J.Mech. Phys. Sol. , 127-140 (1963).26 Rovere, M. Nielaba, P. & Binder, K. Simulation studies ofgas-liquid transitions in 2 dimensions via a subsystem-block-density distribution analysis, Z. Phys. , 215-228 (1993).27 Csikor F. F. & Groma, I. Probability distribution of internalstress in relaxed dislocation systems.,Phys. Rev. B , 064106(2004).28 Yukalov, V. I. Phase-transitions and heterophase fluctua-tions., Phys. Rep. , 395-489, (1991).29 Laurson L. & Alava, M. J. Dynamic Hysteresis in CyclicDeformation of Crystalline Solids., Phys. Rev. Lett. ,155504 (2012).30 Schall, P. Weitz, D. A. & Spaepen, F. Structural rearrange-ments that govern flow in colloidal glasses, Science ,1895-1899 (2007).31 Chikkadi, V. et al., Percolating Plastic Failure as aMechanism for Shear Softening in Amorphous Solids.,arXiv:1312.4136.32 ten Wolde P. R. & Frenkel, D. Enhancement of protein crys-tal nucleation by critical density fluctuations., Science, ,1975-1978 (1997).33 Irvine, W. T. M. Hollingsworth, A. D. Grier D. G. &Chaikin, P. M. Dislocation reactions, grain boundaries, andirreversibility in two-dimensional lattices using topologicaltweezers. Proc. Nat. Acad. Sc., , 15544-15548 (2013).34 Karmakar, S. Lemaitre, A. Lerner, E. and Procaccia, I.Predicting Plastic Flow Events in Athermal Shear-StrainedAmorphous Solids., Phys. Rev. Lett., , 215502 (2010). 10 Supplementary Information for:Pre-yield non-affine fluctuations and a hidden critical point in strained crystals
Tamoghna Das, a , b Saswati Ganguly, b Surajit Sengupta ∗ c and Madan Rao d ‡ a Collective Interactions Unit, OIST Graduate University,1919-1 Tancha, Onna-son, Okinawa, Japan - 904-0495 b Centre for Advanced Materials, Indian Association for the Cultivation of Science, Jadavpur, Kolkata 700032, India c TIFR Centre for Interdisciplinary Sciences, 21 Brundavan Colony, Narsingi, Hyderabad 500075, In-dia. email: [email protected] d Raman Research Institute, C.V. Raman Avenue, Bangalore 560080, India
A.1 The non-affine parameter χ , identification of non-affine clusters and Probability P ( χ , ε ) of non-affine displacements as a function of strain In this section, we give details of how we identify non-affine droplets in the 2D-LJ and 3D-LJ solids under pure shear [1]. Ourprocedure follows closely that used in Ref.[2]. Since all the particles in the ideal solid are identical, non-affine fluctuations incrystals arise solely from thermal effects. Non-affine regions therefore are composed of those particles which suffer somewhatlarge random deviations from their ideal lattice positions due to thermal motion. We identify particles in non-affine regions bycomputing the instantaneous non-affine parameter χ .To obtain the non-affine parameter χ for a particle i we use a metric neighbourhood Ω around particle i defined using a cutoffsize Λ = . χ = ∑ j [ r i − r j − ( + ε )( R i − R j )] with respectto choices of the strain tensor ε . Here r i and R i correspond to the instantaneous and reference positions of i , respectively andthe sum is over the neighbours of i . A similar procedure was used in Ref.[3] to obtain a measure of non-affine displacements inamorphous solids. However, we modify the definition used in [3] in two important ways.Firstly, we compute χ using a fixed reference { R i } . For the case where pure shear strain is zero (as in Ref.[2]) we have usedthe T = ε for our reference configuration. Secondly, the number of neighbours in Ω betweenthe instantaneous and reference neighbourhoods is allowed to be different, accounting for inclusions within Ω which changesthe metric coordination number. If there are no inclusions, the computed χ and ε are identical in meaning to that in [3]. If,however, there are inclusions, then extra contributions to both ε and χ arise from a change in coordination. In this case, the fitted ε looses its meaning as a local strain and χ computes non-affineness in both displacements and the number of inclusions. For theparameters that we study, the number of inclusions is always at most one.In Fig.12 we plot the probability distribution P ( χ ) obtained for a solid of ρ = .
99 for strains (pure shear) 0 . ≤ ε ≤ . P ( χ ) is an unimodal distribution[4]. When the number of neighbours is allowed to vary betweenthe reference and the deformed configurations then P ( χ ) has additional peaks corresponding to each inclusion. As ε increases,the weight of P ( χ ) shifts to the peaks at large χ . From the distribution P ( χ ) we determine a cut-off χ cut above which particledisplacements are defined to be non-affine . This is defined as in Ref.[2] by obtaining P ( χ ) for an equivalent harmonic systemand equating χ cut with the value above which P ( χ ) for the harmonic lattice has negligible weight. This effectively discards 90%of the weight of the first peak in P ( χ ) (Fig.12 a ) as coming from trivial harmonic fluctuations. All particles in a single snapshothaving χ > χ cut are tagged and tagged particles residing within the first nearest neighbour shell of another tagged particle belongto the same cluster; any cluster having at least 7 particles is identified as a non-affine droplet.Our results are robust with respect to small variations of the cut-offs Λ and χ cut . Increasing Λ increases the non-affine fractionwhile increasing χ cut decreases it. Large deviations of these numbers therefore shifts the focus from the fluctuations of interestnamely all non-harmonic, non-affine fluctuations occurring within a neighbourhood of the dimensions of the interaction volumeof the system.A similar situation exists in 3D. We present the normalised distribution of χ as a 3D-LJ solid (face centred cubic) is strainedvia pure shear protocol (Fig.13 (a) & (b)) as in 2D. In this case, we use a metric neighbourhood with Λ = . χ similarto the 2D-LJ case. As the system under consideration is deep in the solid region ( ρ = . , T = . P ( χ ) forthe unstrained case ( ε = .
0) is unimodal and breaks up into a bimodal distribution once sheared. First peak of the distributiondecreases upon increasing shear at cost of higher, co-ordination number changing, non-affineness emerging in the system which11 ✏ = 0 . ✏ = 0 . ✏ = 0 . ✏ = 0 . ✏ = 0 . ✏ = 0 . ✏ = 0 . P ( ) P ( ) ⇢ = 0 . a b Fig. 12: Normalized distribution of non-affine parameter P ( χ ) is shown for the case of ρ = .
99 at T = .
4. The two different figures a & b ,shows two different ranges of the same distribution. Note the difference in the scale of the y -axis. P ( χ ) is essentially bimodal for the unstrainedcase as in Ref.[2]. Applied strain ε enhances non-affine displacements as observed by the increase in height of the second peak (2 . ≤ χ ≤ . χ ≤ . Λ = .
5, the interaction range, the subsequent peaks arise when there areinclusions within the reference volume Λ changing the number of neighbors. contribute to the second peak. The situation is qualitatively in fair accord with its 2D counterpart. These distributions have beenused to set the cut-off χ using the same criterion as in 2D and thus to define the non-affine clusters presented in the main text forthe 3D case. Fig. 13: Increase in non-affineness upon shearing a 3D LJ solid (FCC) at ρ = . T = . a & b . Both of them are part of thesame distribution, but shown with separate axes for clarity. A.2 Local thermodynamics of non-affine droplets in other model two dimensional solids
In the present paper [1] and our earlier work [2] we have related non-affine displacement fluctuations in the Lennard Jones (LJ)solid to droplet fluctuations from “nearby” metastable states such as the liquid and glass. How general is this association andwhat are the limits of its validity ? To test the generality of this result, we carry out our analysis on another model solid where theparticles interact with a different potential (for which it is known that both liquid and glassy free energy minima exist). To testthe limits of this result, we study a model solid where, by construction, there are no locally stable phases except for the crystallinesolid. 12 ig. 14: a A scatter plot showing the local pressure difference ∆ p c as a function of the local density ρ c for “non-affine” droplets in theharmonic triangular net obtained using a lower cut-off χ cut = .
2. The different colours correspond to solids at different densities obtainedby applying a hydrostatic compression (or expansion) on the T = K which is the same for all the solids. b Plot of local thermodynamicquantities ∆ p c vs. ρ c for non-affine droplets in the WCA system at a T ∗ = . n c = ,
15 and 20. The prominent van derWaals loop in this quantity is similar to the corresponding plots (see Fig. 2 in [1] or Fig. 7(b) in [2]) for the LJ system.
A.2.1 Harmonic triangular net
The harmonic triangular net where particles in a two dimensional triangular lattice areconnected to their nearest neighbors by harmonic springs (spring constant, K ) is the simplest model of a network solid. Thepotential energy for this system is given by, V = ∑ (cid:104) i j (cid:105) K ( | r i − r j | − r ) . (1)Here, r i is the position vector of the i − th particle, r is the equilibrium bond length and the sum extends over the nearest neighborshell. Since defects, which change the connectivity of the lattice are not allowed by construction, such a solid never melts and thefree energy surface in configuration space is particularly simple having just one quadratic minimum corresponding to the T = P ( χ ) in this system using the same procedure asused in the LJ system and using a neighborhood which includes the first two neighbor shells. The control parameter in our modelis the lattice parameter l p = (cid:113) / √ ρ where ρ is the density of the solid and we hold the temperature T , the equilibrium bondlength r and K fixed at unity. The reference configuration is the ideal triangular lattice with the same lattice parameter. As inthe LJ solid, the P ( χ ) is bimodal and resembles Fig.12 and Fig. 2 in [2] for the unstrained solid at values of K low enoughsuch that particles enter or leave the reference volume Ω . When non-affine droplets are defined using the same cutoff χ cut asin [2], the local thermodynamics of these droplets as given by the ∆ p c vs. ρ c scatter plots are entirely unremarkable and arenormally distributed around the average density of the solid. This is shown in Fig.14 a . Further, this behaviour is, in this system,independent of the cut-off χ cut used. We do not see any evidence of a van der Waals loop in this case.This is an important null result which shows that systems where there are no liquid like or glassy minima, do not have nontriviallocal thermodynamics of the non-affine droplets and these droplets therefore do not correspond to fluctuations of any metastableminima. A.2.2 The Weeks-Chandler-Anderson solid
In contrast to the harmonic solid, the system of particles interacting with thepurely repulsive Weeks, Chandler and Anderson (WCA) potential [5] does undergo melting [6] and glass [7] transitions. Thispotential is closely related to the LJ potential and is given by, V ( r ) = ( r − − r − ) + r ≤ / (2) = N = P ( χ ) . Local thermodynamics of these droplets then yield ∆ p c vs. ρ c curves which are qualitatively similar to the LJ case in[2].We have thus shown that there is strong correlation between the appearance of non-trivial local thermodynamics of the non-affine droplets and the presence of metastable liquid or glass like phases. While not a complete proof, taken together, our resultsstrongly suggest that non-affine fluctuations in a crystalline solid are droplet fluctuations from nearby metastable minima. References , 041115 (2010).3 M. L. Falk and J. S. Langer, Phys. Rev. E , 7192 (1998).4 S. Ganguly, S. Sengupta, P. Sollich, M. Rao, Phys. Rev. E , 042801 (2013).5 J. D. Weeks, D. Chandler, and H. C. Andersen, J. Chem. Phys. , 5237 (1971).6 A. Ahmed and R. J. Sadus, Phys. Rev. E. , 061101 (2009).7 T. Kawasaki, T. Araki and H. Tanaka, Phys. Rev. Lett.99