Exploring the link between crystal defects and non-affine displacement fluctuations
Pankaj Popli, Sayantani Kayal, Peter Sollich, Surajit Sengupta
EExploring the link between crystal defects and non-affine displacement fluctuations
Pankaj Popli, Sayantani Kayal, Peter Sollich,
2, 3 and Surajit Sengupta Tata Institute for Fundamental Research, Centre for Interdisciplinary Sciences, 36/P Gopanapally, Hyderabad 500107, India. King’s College London, Department of Mathematics, Strand, London WC2R 2LS, UK Institute for Theoretical Physics, University of G¨ottingen, 37077 G¨ottingen, Germany (Dated: July 18, 2019)We generalize, and then use, a recently introduced formalism to study thermal fluctuations ofatomic displacements in several two and three dimensional crystals. We study both close packed aswell as open crystals with multi atom bases. Atomic displacement fluctuations in a solid, once coarse-grained over some neighborhood may be decomposed into two mutually orthogonal components. Inany dimension d there are always d affine displacements representing local strains and rotationsof the ideal reference configuration. In addition, there exists a number of non-affine localizeddisplacement modes that cannot be represented as strains or rotations. The number of these modesdepends on d and the size of the coarse graining region. All thermodynamic averages and correlationfunctions concerning the affine and non-affine displacements may be computed within a harmonictheory. We show that for compact crystals, such as the square and triangular in d = 2 and thesimple, body-centered and face-centered cubic crystals in d = 3, a single set of d − fold degeneratemodes always dominate the non-affine sub-space and are separated from the rest by a large gap.These modes may be identified with specific precursor configurations that lead to lattice defects. Inopen crystals, such as the honeycomb and kagome lattices, there is no prominent gap although softnon-affine modes continue to be associated with known floppy modes representing localized defects.Higher order coupling between affine and non-affine components of the displacements quantify thetendency of the lattice to be destroyed by large homogeneous strains. We show that this coupling islarger by almost an order of magnitude for open lattices as compared to compact ones. Deformationmechanisms such as lattice slips and stacking faults in close packed crystals can also be understoodwithin this framework. The qualitative features of these conclusions are expected to be independentof the details of the atomic interactions. PACS numbers: 62.20.D-, 63.50.Lm, 63.10.+a
I. INTRODUCTION
A crystalline solid is formed from the liquid whentranslation and orientation symmetries are broken asa result of a thermodynamic phase transformation [1].At any non-zero temperature, atomic fluctuations in thecrystalline solid tend to restore those symmetries. Thecurrent paradigm classifies atomic fluctuations as beingeither “smooth” or “singular”. The former compriselong wavelength, hydrodynamic, smooth variations of theelastic displacements and density fields [2, 3], while thelatter are defects, where the displacement becomes dis-continuous [4–6]. Such a classification has proved tobe immensely useful in understanding many commonlyobserved properties of solids [7] as well as melting [1]and failure of solids in response to external mechani-cal loads [4, 8]. This viewpoint has also been extendedwith some success towards understanding the mechanicalproperties of amorphous solids [9]. However, in contrastto crystals, the lack of long ranged structural order inamorphous solids precludes a clear distinction betweensmooth and singular displacements [10, 11]. It has beensuggested (and experimentally observed) that localized, non-singular displacement configurations play the rˆole oflattice defects in such solids [9, 12, 13].For a while now, there has been an effort by some ofus to formulate a different way of classifying atomic dis-placements in solids, which we hope, may have more gen- eral applicability [14–23]. Borrowing an idea first usedto study mechanical deformation in glasses [24], it wasshown that any set of atomic displacements of an atomand its neighbors, within a specified “coarse graining” re-gion, may be decomposed into two mutually orthogonalsub-spaces using a well defined projection formalism [14–16]. The affine component of these displacements repre-sents homogeneous linear transformations of a referenceconfiguration within the coarse graining volume. Ignor-ing trivial uniform displacements, these are isotropic ex-pansion, shear strains and local rotations. Since, in gen-eral, not all displacements can be described completelyby these linear transformations, inevitably, a non-affine component remains. By construction, the affine and thenon-affine parts of the displacements are linearly inde-pendent [14]. Thermodynamically conjugate fields maybe defined, which enhance or suppress each part indepen-dently of the other. The affine displacements couple tolocal stresses (and torques) while the non-affine compo-nent of the displacement couples to a new “non-affine”field [15, 21, 22]. Enhancing non-affine fluctuations by in-creasing temperature, applying large strains or the non-affine field leads to the creation of defects [15, 18, 22].Indeed, atomic fluctuations that act as precursors to theformation of defects such as dislocation dipoles have beenshown to be the most prevalent, though not the sole,non-affine displacement even within a small oscillation,harmonic, approximation [15, 18]. a r X i v : . [ c ond - m a t . s t a t - m ec h ] J u l This decomposition of displacements into affine andnon-affine have lead to a deeper understanding in severaldisparate contexts. For example, it has been possibleto show that rigid solids are always thermodynamicallymetastable for infinitesimal stresses. Loading a rigid solidis tantamount to quenching it across an equilibrium firstorder phase transition [22]. Yielding of a crystal underload is simply the decay of this metastable phase to thestable phase, which eliminates stress by non-affine atomicrearrangements. This radically different viewpoint nev-ertheless allows one to calculate, for the first time, strainrate dependent yielding thresholds using classical nucle-ation theory.In networked solids, where atoms are bound by strongchemical bonds, dislocations do not form. Neverthelesssuch solids may deform by special singular formationscalled pleats or “ripplocations”, which have been de-scribed within the same language of non-affine displace-ments [17, 20].It has further been shown that colloidal crystals withany given interaction may be arranged in any structurewhatsoever if it is only possible to suppress non-affine dis-placements away from this reference configuration. Spe-cial, dynamic, feedback controlled laser traps have beenproposed, though not yet experimentally realized, whichmay be able to perform this feat. Unlike static traps, thestructures stabilized by such a process are translationallyinvariant and possess all allowed zero modes [21].Finally, in a protein - a large molecule consisting ofinteracting atoms with no spatial long ranged order - ithas been shown that important conformational changes,which precede binding to ligands are always non-affine.These may be discovered by simply projecting out thelocal atomic displacements using the same projection for-malism. Regions with large susceptibility for non-affinedisplacements correlate with binding hotspots and spa-tial correlations of the magnitude of non-affine ness marksites of allosteric control [23].In this paper, we carry this program forward by re-turning to a study of small amplitude non-affine displace-ments in periodic crystals. We attempt to answer somequestions that naturally arise from our previous work. Sofar, except for the work on proteins [23], the non-affineanalysis has been mainly restricted to two dimensionalsolids for simplicity. Most experiments, however, are inthree dimensions and an extension of our work to threedimensions is therefore necessary. We have also felt theneed for a comparative study of harmonic non-affine fluc-tuations among different kinds of crystals, especially be-cause of their importance as precursors to lattice defects.Finally, we would like to test the robustness of our re-sults to variation of some parameters such as the extentof coarse graining and nature of interactions. One of theaims of the present work is to discover features that arecommon to all crystalline solids and differentiate themfrom those that depend on details of the crystal struc-ture, dimensionality and interaction parameters.Our main results are as follows. The eigenvalue spec- trum of the Hessian of the coarse-grained Hamiltoniantaken with respect to local atomic displacements and pro-jected onto the non-affine subspace [15], always shows aprominent gap between the largest eigenvalues and theothers, for all Bravais lattices with a monatomic basis.The gap increases with the size of the coarse graining vol-ume. The dominant non-affine eigenmode correspondingto the largest eigenvalue features displacements that maybe identified with defect precursors [15, 16]. For open lat-tices featuring a multi atom basis, the gap is much lessprominent, although large eigenvalue floppy modes con-tinue to resemble precursors for known defects. The rel-ative prominence of modes in the non-affine eigenvaluespectrum for open lattices is more sensitive to the na-ture of the interactions compared to those in close packedcrystals. Spatial correlation functions of the affine andnon-affine modes are similar in nature among the variouscrystal structures studied. Affine and non-affine modescouple at, higher than linear, order [14]. This couplingmeasures the susceptibility of the crystal to producingnon-affine displacements in response to small externalstress. We find that open crystals are more susceptiblethan close packed ones by almost an order of magnitude.The rest of the paper is arranged as follows. In thenext section (Section II) we introduce the projection for-malism for atomic displacements. Some of this materialhas been discussed in our earlier publications but we gen-eralize the derivation here to make it directly applicableto crystals with a multi-atom basis and in any dimen-sion. This is followed by Section III where we present abrief description of the crystal lattices studied and theinteraction models. In two dimensions we study the tri-angular, square, honeycomb and kagome lattices, whilein three dimensions we restrict ourselves to cubic crys-tals and study the simple cubic, body centered cubic andface centered cubic structures. In Section IV we presentour results for several quantities related to the statisticsof displacement fluctuations, critically comparing themamong the different lattices. In Section V we summarizeour main findings in detail by discussing our main resultsand analyzing their implications, especially with respectto indications to directions of future research.
II. NON AFFINE FLUCTUATIONS: THEPROJECTION METHOD
In this section we present the method we use to projectout local atomic displacements into affine and non-affinecomponents, extending it from earlier formulations [14–16] by taking care that the derivation is directly applica-ble to crystals with a multi-atom basis in any dimension.In the most general setting, non-affine displacements arethose atomic displacements that cannot be captured bya homogeneous deformation. For example, imagine anideal (defect free) crystal at zero temperature consistingof atoms placed at their reference position. The Cauchy-Born rule [25] (CBR) states that any external deforma-tion caused by changing the shape of the boundary ofthe solid is distributed homogeneously among the atomsof the crystal, which are shifted appropriately from theirreference position. This, of course, amounts to statingthat for an ideal crystal there are no non-affine defor-mations at zero temperature. At finite temperature, oneexpects this rule to hold on the average, with the (finitesize scaled) elastic moduli [26, 27] setting the scale forthe displacement fluctuations averaged within some localcoarse graining volume Ω. While this is more or less true,there is a subtle point here that needs attention. Indeed,atomic displacement fluctuations within Ω decomposenaturally into two mutually orthogonal subspaces. One,designated as affine, consists of fluctuations for which theCBR holds locally and instantaneously, while the otheris the set of non-affine displacement modes for which theCBR is violated. While the former may be directly con-nected to the elastic moduli, the latter represents fluctu-ations that act as precursors for defects. We show belowhow this decomposition may be carried out for a genericcrystalline solid. We consider a d dimensional lattice with N lattice sites and N b basis particles per site. The to-tal number of particles in the system is N × N b . Wetake { R iα } as the equilibrium position vector for anysite i ∈ { , , . . . N − } , where α ∈ { , , . . . N b − } represents the index of a basis particle. To distinguishaffine and non-affine displacements we consider relativedisplacements of pairs of atoms whose reference positionsare within some fixed coarse-graining distance r of eachother. Specifically, the coarse-graining region around thebasis on lattice site i is defined asΩ( i ) = { ( jγ, iα ) | < | R jγ − R iα | ≤ r ∧ ( j (cid:54) = i ∨ γ > α ) } . (1)In words, Ω( i ) contains all pairs of indices ( jγ, iα ) ofparticles within the specified dsitance r ; at least one ofthese particles has to belong to the basis around latticesite i . The last constraint in Eq.(1) merely avoids double-counting of pairs within the basis, by insisting that indexpairs of the form ( iγ, iα ) have got ordered basis indices, γ < α . We denote the number of particle pairs in Ω( i ) by N Ω and number the elements of Ω( i ) in some arbitraryfashion as Ω( i ) = { ( j n γ n , iα n ) , n = 1 . . . N Ω } We note that the lattice symmetries mean that all neigh-bourhoods Ω( i ) for different i are just translated copiesof each other.When the lattice is deformed particles undergo dis-placements and take new positions { r iα } ; we write thedisplacement from their equilibrium positions as u iα = r iα − R iα . It has been shown[14–16] that the displace-ments in a given deformed coarse-graining volume can beexpressed as a linear combination of independent affineand non-affine deformations. For a fully affine deforma-tion of the coarse-graining volume around lattice site i there is by definition a local d × d dimensional deforma- tion matrix D such that u j n γ n ,iα n = D R j n γ n ,iα n , n = 1 . . . N Ω (2)using the abbreviations u j n γ n ,iα n = u j n γ n − u iα n and R j n γ n ,iα n = R j n γ n − R iα n . In general, { u iα } will havecontributions from non-affine transformations as well. Insuch cases D is defined as the matrix that minimizes χ i = min D (cid:32) N Ω (cid:88) n =1 ( u j n γ n ,iα n − D R j n γ n ,iα n ) (cid:33) (3)Therefore χ i is a measure of the non-affinity at the givenlattice site i .We now introduce some simplified notation by rear-ranging components of u j n γ n − u iα n for all n into a col-umn vector ∆ of length N Ω d with elements∆ nµ = u µj n γ n ,iα n where µ = 1 . . . d denotes the spatial components of thedisplacement vectors. Similary a column vector e oflength d is obtaind by arranging the elements of D inorder ( D , D ... D d ... D d , D d ... D dd ) . With these defi-nitions, Eq. (3) may be written concisely as χ = min e (cid:16) [ ∆ − R e ] (cid:17) (4)Here we have introduced a block matrix R whose elementsare given by R nµ,νν (cid:48) = δ µν R ν (cid:48) j n γ n ,iα n and the entries of the vector R e are given by (cid:80) νν (cid:48) R nµ,νν (cid:48) e νν (cid:48) in the obvious way. Once Eq. (4) isminimized, we obtain the residual contribution from non-affine deformation and the “best-fit” affine strain. Thesecan be expressed as χ = ∆ T P ∆ (5) e = Q ∆ . (6)with the matrices P = I − RQ and Q = (cid:0) R T R (cid:1) − R T . Notethat P is a projection matrix ( P = P ), which, when act-ing on the space of ∆ , extracts only the non-affine com-ponent of the displacements. It can be seen that P has d zero and N Ω d − d non-zero eigenvalues correspondingto the affine and non-affine subspaces respectively. Asusual I − P will then project out the affine componentof ∆ . The elements of the best-fit affine transformationmatrix D from Eq. (6) can be written explicitly as D µν = (cid:88) ν (cid:48) ( M − ) νν (cid:48) (cid:88) n R ν (cid:48) j n γ n ,iα n ∆ nµ (7)in terms of the matrix M with elements M νν (cid:48) = (cid:88) n R νj n γ n ,iα n R ν (cid:48) j n γ n ,iα n For the lattices considered in this paper, M is diago-nal due to lattice symmetries (see also the discussion inRef. [14]) so that Eq.(7) simplifies to D µν = 1 (cid:80) n ( R νj n γ n ,iα n ) (cid:88) n R νj n γ n ,iα n ∆ nµ (8)We now obtain the statistics of ( χ, e ) in the classicalcanonical ensemble for any lattice and dimension. Forany given Hamiltonian H the canonical probability dis-tribution is P ( p , u ) = 1 Z e − β H ( p , u ) (9)Here we restrict the Hamiltonian to the harmonic ap-proximation , H = (cid:88) iα p iα m iα + 12 (cid:88) iα,jγ (cid:88) µν u µiα φ µνiα,jγ u νjγ (10)where p iα are the momenta, m iα the masses and φ µνiα,jγ are the elements of the Hessian. The Hamiltonian can beeasily diagonalized if one takes a plane wave ansatz forthe displacements. We therefore write u iα = (cid:90) d q V BZ u α ( q ) e i q · R iα and similarly u iα,jγ = (cid:90) d q V BZ u α ( q ) (cid:0) e i q · R iα − e i q · R jγ (cid:1) where the integration runs over the first Brillouine zonewith volume V BZ and q is the wave vector. The LatticeGreen’s Function (LGF) may be obtained as the inverse G ( q ) = D − ( q ) of the dynamical matrix D ( q ) with ele-ments D µναγ ( q ) = (cid:88) i φ µνiα, γ e i q · ( R iα − R γ ) (11)With the knowledge of the LGF, thermal averages of dif-ferent quantities are easy to calculate. For example thedisplacement correlator reads in Fourier space (cid:10) u µα ( q ) u νγ ( q (cid:48) ) (cid:11) = G µναγ ( q ) δ ( q + q (cid:48) ) V BZ , and in real space (cid:10) u µiα u νjγ (cid:11) = 1 β (cid:90) d q V BZ G µναγ ( q ) e i q · ( R iα − R jγ ) . (12)Along similar lines, for our coarse-graining volume onecan obtain the thermal average of any observable A ( ∆ )as (cid:104) A ( ∆ ) (cid:105) = 1 Z (cid:90) d ∆ A ( ∆ ) e − ∆ T C − ∆ (13) with the normalisation constant Z = (2 π ) N Ω d/ (cid:112) | C | .The covariance matrix C in the above equation can beobtained from the LGF via C nµ,mν = (cid:104) ∆ nµ ∆ mν (cid:105) = (cid:90) d q βV BZ (cid:20) G µνγ n γ m ( q ) e i q · R jnγn,jmγm (14) − G µνα n γ m ( q ) e i q · R iαn,jmγm − G µνγ n α m ( q ) e i q · R jnγn,iαm + G µνα n α m ( q ) e i q · R iαn,iαm (cid:21) To obtain the statistics of ( χ, e ) we make use of Eq. (13)and obtain the characteristic function Φ (k , κκκ ) for thejoint probability distribution P ( χ, e ).Φ (k , κκκ ) = exp (cid:18) − κκκ T QC ( I − PC ) − Q T κκκ (cid:19) (15) × (cid:112) | I − ik PCP | . Using the identity( I − ik PC ) − = I + ( I − ik PC ) − (2 ik PC )the last result can be written in terms of the character-istic function for the marginals as follows:Φ (k , κκκ ) = Φ χ (k) Φ e ( κκκ )e − ik κκκ T QC ( I − PC ) − PCQ T κκκ (16)where Φ χ ( k ) = 1 (cid:81) l √ − ikσ l (17) Φ e ( κκκ ) = e − κκκ T QCQ T κκκ (18)and the σ l are the eigenvalues of PCP . For κκκ = 0 and k = 0, Φ (k , κκκ ) reduces to Φ χ ( k ) and Φ e ( κκκ ) respectively.The term in the exponential governs the (non-linear) cou-pling between the non-affine and affine components ofthe displacements. Previous work has shown that thiscoupling is significant only for large uniaxial and shearstrains [14]. For smaller strains, it can largely be ignored.With the knowledge of the characteristic function,thermal averages and other higher order moments maybe computed such as (cid:104) χ (cid:105) = Tr ( PCP ) (19) (cid:10) ee T (cid:11) = QCQ T (20)From Eq. (19) it is clear that (cid:104) χ (cid:105) is a sum over the eigen-values of ( PCP ). Each eigenvalue represents the contri-bution of a specific non-affine mode to χ . It has beenshown [15] that these eigenvalues are elements of the in-verse Hessian of the free energy in the direction of thenon-affine mode in configuration space. A large eigen-value implies a small value of the local curvature of thefree energy minimum. We shall see later in this paperthat the corresponding eigenvectors are precisely thoseatomic displacement fluctuations that lead to lattice de-fects or other imperfections tending to destroy crystallineorder. The non-affine mode corresponding to the largesteigenvalue therefore has the highest contribution to thisprocess.The non-affine parameter χ depends linearly on tem-perature T and is inversely proportional to the strengthof the interaction. Due to the fact that the underly-ing distribution of ∆ is Gaussian, (cid:104) e (cid:105) vanishes unless anexternal stress is present. Finally, the leading order non-linear coupling between non-affine and affine modes isgiven by (cid:10) χ ee T (cid:11) − (cid:104) χ (cid:105) (cid:10) ee T (cid:11) = 2 QC [ P , C ] Q T . (21)Two-point distributions and spatial correlation func-tions for χ and e may also be calculated following theprocedure explained in Refs. [14–16]. Below we includea brief description for completeness.The spatial correlations of the non-affine parameter aregiven by (cid:104) χ ¯ χ (cid:105) − (cid:104) χ (cid:105)(cid:104) ¯ χ (cid:105) = 2 Tr ( P¯CP )( P¯CP ) T (22)The two point covariance ¯C between relative displace-ments within two coarse-graining neighbourhoods Ω ≡ Ω( i ) and ¯Ω ≡ Ω( k ) around lattice positions R i and R k ,respectively, is defined as¯ C nµ,mν = (cid:104) ∆ nµ ¯∆ mν (cid:105) (23)It is obtained from an expression similar to Eq. (14) (seeRef. [14] for details):¯ C nµ,mν = (cid:90) d q βV BZ (cid:20) G µνγ n γ m ( q ) e i q · R jnγn,lmγm (24) − G µνα n γ m ( q ) e i q · R iαn,lmγm − G µνγ n α m ( q ) e i q · R jnγn,kαm + G µνα n α m ( q ) e i q · R iαn,kαm (cid:21) where we have assumed that the elements of Ω( k ) arenumbered ( l m γ m , kα m ). For all simple lattices in 2 d and3 d with one particle basis, the correlations (cid:104) χ ¯ χ (cid:105) − (cid:104) χ (cid:105)(cid:104) ¯ χ (cid:105) are short ranged.Strain-strain correlation may be obtained from (cid:10) e ¯ e T (cid:11) = (cid:10) Q ∆ ¯ ∆ T Q T (cid:11) = Q¯CQ T . (25)It is often useful to consider these correlations in Fourierspace, where they can be expressed as [14] (cid:10) e ¯ e T (cid:11) ( q ) = Q¯C ( q ) Q T (26)with the Fourier transform ¯C ( q ) defined via¯ C nµ,mν = (cid:90) d q V BZ ¯ C nµ,mν ( q ) e i q · R i ,k . Comparison with Eq.(24) then shows that β ¯ C nµ,mν ( q ) isgiven directly by the term in square brackets in Eq.(14).This follows from the fact that the particle pairs in Ω( k )are just those in Ω( i ) translated by R k ,i ; e.g. in thefirst term of Eq.(24) one has after extracting the Fourierfactor e i q · R i ,k the phase term e i q · ( R jnγn,lmγm − R i ,k ) = e i q · [ R jnγn − ( R lmγm − R k ,i )] = e i q · ( R jnγn − R jmγm ) = e i q · R jnγn,jmγm Correlations of the affine displacements viz. local vol-ume change ( e v ), uniaxial or deviatoric strain ( e u ), shearstrain ( e s ) and rotation ( w ) respectively, may be obtainedusing the component forms as follow, (cid:10) e v (cid:11) ( q ) = E + E + 2 E (27) (cid:10) e u (cid:11) ( q ) = E + E − E (cid:10) e s (cid:11) ( q ) = E + E + 2 E (cid:10) w (cid:11) ( q ) = E + E − E , Here we have used the same notation for the fourthrank tensor, E = Q ¯ C ( q ) Q T as in Ref. [14], and (cid:10) e v (cid:11) ( q )etc. are strain correlators in Fourier space. Expressionsfor 3d can be obtained in similar fashion, for instance,strain correlation in 3d for volume, uniaxial and shearin the x − y plane are as follows, (cid:10) e v (cid:11) ( q ) = E + E + E (28)+ 2 ( E + E + E ) (cid:10) e u (cid:11) ( q ) = E + E + E − E + E − E ) (cid:10) e s (cid:11) ( q ) = E + E + 2 E . Other components of the shear strain (and rotation) canbe written down by analogy.
III. MODELS
In this paper, we use the methods of the last sectionto obtain the statistics of affine and non-affine displace-ments for a number of lattices in two (2d) and three (3d)dimensions. In 2d we consider lattices both with a sin-gle atom basis such as the triangular and square latticesas well as those with a multi-atom basis like the planarhoneycomb and the kagome lattices. In 3d we confineourselves to a study of cubic systems, namely, the simplecubic, body centered and face cantered cubic lattices. Inorder to keep the discussion general we model the inter-actions by harmonic springs. Our results are thereforevalid for any crystalline solid at sufficiently low temper-atures where anharmonic effects may be neglected. Atypical Hamiltonian for such interactions is H = (cid:88) iα p iα m + (cid:88) (cid:104) iα,jγ (cid:105) k iα,jγ (cid:16) u iα,jγ · ˆ R iα,jγ (cid:17) (29)Here k iα,jγ determines the spring constant acting be-tween particle pairs iα, jγ . The k iα,jγ are chosen suchthat the lattice is stable and satisfies Maxwell’s stabilitycriteria [29]. In particular, we take k iα,jγ to be equal to k for nearest neighbours and k for next nearest neigh-bours; interactions beyond the second neighbor shell areneglected. Additionally, throughout the paper, the near-est neighbour bond strength k and the lattice constant a are chosen to be unity without any loss of generality.This sets the scales for energy and length respectively.We have also, in some cases, studied the effect of in-cluding simple three body bond-angle dependent poten-tials in order to introduce an energy cost for bond bend-ing. These interactions are very well documented in theliterature mostly on the system like graphene [30, 31]. Tomodel bond bending we take the Kirkwood [36] model, H bend = k b (cid:88) (cid:104) iα,jγ,kδ (cid:105) (∆ θ iα,jγ,kδ ) . which in the small oscillation approximation can be writ-ten as H bend (cid:39) k b (cid:80) (cid:104) iα,jγ,kδ (cid:105) [cot θ ( ˆ R jγ,kδ · u jγ,kδ + ˆ R iα,kδ · u iα,kδ ) − θ ( ˆ R iα,kδ · u jγ,kδ + ˆ R jγ,kδ · u iα,kδ )] (30)where θ is the equilibrium angle and angular bracketsdenote triples of particles where iα and jγ are both near-est neighbours of kδ .In the following section we discuss our results for spe-cific lattices in 2d and 3d. In Fig. 1 we have shown theselattices schematically and indicated the bonding interac-tions that we have considered. IV. RESULTS
We are now in a position to use the methods describedin Sec. II to obtain the statistics of coarse grained non-affine and affine displacements of particles interactingthrough the Hamiltonians presented in Sec. III for acollection of 2d and 3d lattices (see Fig. 1). As discussedabove, the statistics of χ can be obtained once one hasknowledge of the matrix PCP . The projection matrix P only depends upon the reference position of particles inthe lattice and can be constructed easily. The covari-ance matrix C can be calculated using Eq. (14) once one a. b.c.d. e.f.g. FIG. 1. Schematic diagram showing the triangular (a), square(b), planar honeycomb (c), kagome (d), simple cubic (e), bodycentered cubic (f) and face centered cubic (g) lattices. Thenearest neighbor bonds are shown in bold while the next near-est neighbor bonds, whenever present, are drawn using dashedlines. The parameter a is the lattice constant, chosen to beunity. The equilibrium bond angle θ has also been markedfor the triangular and honeycomb lattices. In the 3d casesshaded color regions have been added to make the cubic ge-ometry clearer. knows the dynamical matrix D ( q ). For the harmonic in-teractions with nearest (and next nearest) neighbours wecan compute D ( q ) in a straightforward manner for alllattices; the results are listed in Appendix A. The prob-ability distribution for χ can then be obtained using theeigenvalues of PCP . We have also checked our resultsby directly simulating the model systems using standardmolecular dynamics in the canonical ensemble [33] as im-plemented in the LAMMPS simulation package [34]. Allour results scale linearly with temperature k B T = β − ,where k B is the Boltzmann constant. We have used dif-ferent temperatures for different lattices only for ease ofpresentation.Our results for P ( χ ) obtained by numerically inverting Φ χ ( k ) are shown in Fig. 2 and Fig. 3 together with theresults from direct simulations. All the averaging is doneover at least 1000 well equilibrated and uncorrelated con-figurations. Needless to say, the agreement is perfect asexpected.Once our formalism is thus established for all the sys-tems considered in 2d and 3d, we turn to each lattice indetail below. We show that using our method one canfind the most prominent non-affine displacement modes(eigenvectors of PCP ) for any lattice. Often these modesturn out to be precursors for the most commonly ob-served defect structures for a given lattice system. Wenote that the relative probabilities of different non-affinemodes also depend on the lattice and the interactionsand can be easily captured using our approach.
FIG. 2. Scaled distibution P ( χ (cid:63) ) for all 2d lattices where χ (cid:63) = χ/ (cid:104) χ (cid:105) . The solid colored lines are from our analyticcalculations and the points are simulation results (using N =1024 except for Honeycomb, where N = 512 and Kagome,where N = 300). Triangular: light pink with k b = 0; Square:sky blue, k = 0 .
5; Honeycomb: brown, k = 0 . k = 0 .
5. The distribution for square and triangleplotted here is for the smallest coarse-graining volume usedin Sections IV A and IV B. Whereas, for other lattices coarse-graining volume is same as shown in Fig 6 and 8.FIG. 3. Scaled distibution P ( χ (cid:63) ) for all 3d lattices, where χ (cid:63) = χ/ (cid:104) χ (cid:105) . The meaning of the symbols is the same as inFig. 2. The distribution is plotted for SC: sky blue, k = 0 . N = 1000. BCC : brown, k = 0 . N = 2000. FCC: purple, k = 0 . N = 4000. The coarse-graining volume is describedin the text. A. Triangular
The triangular lattice is the only close packed struc-ture observed in 2d [1]. It has just one basis particle persite. We have previously established that the most promi-nent non-affine mode, i.e. the eigenvector correspondingto the largest eigenvalue of
PCP , corresponds to the in-
FIG. 4. Non-affine modes for triangular lattice with differ-ent sizes of the coarse graining volume Ω at inverse temper-ature β = 1000 with no next nearest neighbor bonds. (a)The spectrum of the eigenvalues of PCP is shown for threedifferent choices of Ω (inset), which consists of all particleswithin the first (Ω - light blue), second (Ω - magenta) andthird (Ω - dark green) nearest neighbor shells. The refer-ence positions of particles are shown by small yellow circles.The horizontal lines show the eigenvalues. Note the largegap between the largest eigenvalue and the rest of the spec-trum. (b) The two degenerate eigenvectors corresponding tothe largest eigenvalue of PCP . Note that a nearest neighborbond is being stretched and a next near neighbor bond nearlyperpendicular to it has been shortened. This mode is sameas the one discussed in Ref. [15]. This displacement tendsto replace the six-fold neighborhoods by two five- and twoseven-fold neighbors producing a tightly bound dislocation–antidislocation pair. (c) and (d) show that increasing Ω doesnot affect the nature of this mode. cipient dissociation of a tightly packed dislocation-anti-dislocation pair [15]. To reach this conclusion we useda coarse-graining volume that included only six nearestneighbor particles (see Fig. 4). This makes
PCP a 12 × < Ω < Ω . At the sametime the gap between the first eigenvalue and the othersincreases significantly.It is interesting that as the size of Ω increases, ad-ditional vibrational modes go on to populate the lowereigenvalues, keeping the gap intact. We show later thatthis phenomenon is quite general and observed for many(but not all) lattices. We will comment further on thisobservation in the discussion (Section V). a) b) c) d) a. b.c.d. i
PCP . These modestend to shift a row of atoms relative to adjacent rows. (c) and(d) show that, as in the triangular case, increasing Ω does notaffect the nature of these modes.
B. Square
For the square lattice, we need to include both near-est and next nearest neighbor bonds in order to satisfythe Maxwell criteria for stability [29]. We have also cho-sen the smallest Ω such that all particles to which thecentral particle is bonded by nonzero interactions are in-cluded. This yields an Ω containing four neighbor andfour next neighbor particles. Hence, one obtains
PCP as a 16 ×
16 matrix. This has 16 eigenvalues with eightnon-zero values. As the size of Ω and with it the de-gree of coarse-graining, is increased one then observesthe same effect on the square lattice as in the triangularlattice (see Fig. 5). Again there is a gap in the eigen-value spectra between the largest eigenvalue and the rest;this gap increases with the coarse-graining scale chosen.Fig 5 shows, in addition to the coarse-graining volumesand the eigenvalue spectra, the softest degenerate eigen-modes. The nature of the mode corresponding to thelargest eigenvalue is somewhat different to the triangularlattice case. Instead of introducing defects, it tends toshift the middle row of atoms with respect to its neigh-boring rows in a direction parallel to the rows. One caneasily see that this corresponds to a precursor that cantake a square lattice to a triangular one by a shuffle ofalternate layers [21].
C. Planar honeycomb
The planar honeycomb (or simply honeycomb!) lat-tice occurs in many condensed matter systems, the most i0 i1 i j j0 j1 a. b.c. i
34 dimensional matrix. The dynamical matrixis a 4 × C as shown in Eq. (14). Our experience with the triangu-lar and square lattices shows that increasing the size ofΩ does not influence the nature of the most prominentnon-affine modes, although it does considerably increasecomputational complexity.The planar honeycomb structure has been studied indetail in an earlier publication [16]. We include here someof those results for completeness. The probability distri-bution of χ is shown in Fig. 2 together with the resultsof other lattices. Fig. 6 shows the non-affine eigenvaluespectra. We observe that, in contrast to the triangu-lar and square lattice, there is no clear gap between thelargest eigenvalue and the others. We suspect that thepresence of optical modes produces an eigenvalue spec-trum that does not have pronounced gaps between dif-ferent modes. It has also been shown in [16] that thenature of this spectrum remains unaffected if one soft-ens the lattices by reducing the value of spring constant k . In fact, as one softens the lattice these eigenvaluesgrow without bound, producing more non-affinity in thesystem. This is obvious because C is proportional to thelattice Green’s function which itself diverges when thespring constant vanishes.Eigenvectors of PCP corresponding to the first threelargest eigenvalues are plotted in Fig. 6. Intriguingly,the second mode represents the precursor to a Stone-Wales (SW) defect [35]. In SW a central bond flips by90 ◦ creating pentagonal and heptagonal voids. One ofthe questions that was not addressed in [16] was what, ifany, is the effect of introducing a bending rigidity [32, 36]to the bonds? We take up this issue below. FIG. 7. Plot of the eigenvalues σ i of PCP (colored lines, toppanel) as a function of the bond angle rigidity parameter k b ,showing the effect of including bending rigidity of bonds inthe triangular a. and planar honeycomb b. structures. Whilethe relative prominence of the modes is unaffected in the tri-angular lattice except for the breaking of degeneracy of somelow probability modes, in the honeycomb lattice the modecorresponding to the SW defect precursor (yellow) is stronglysuppressed with increasing k b (see text). Bond bending rigidity may be modeled as a three bodypotential given by Eq.(30). The dynamical matrix corre-sponding to bond bending can be obtained in closed form(see Appendix A). The total dynamical matrix, which isthe sum of the dynamical matrices for bond stretchingand bending, is used to calculate C .Fig. 7 shows how the spectrum of non-affine (non-zero)eigenvalues of PCP changes upon increasing the bendingconstant k b . We have included a similar calculation forthe triangular lattice for comparison. Because the trian-gular lattice is isotropic, addition of a bond bending costonly stiffens the lattice and decreases the eigenvalues andconsequently their sum, χ . We see that 5-7 defect pre-cursor modes as discussed above continue to be the mostprominent modes in the system. Less importantly, theaddition of bond bending also breaks the degeneracy ofsome of the modes corresponding to small eigenvalues.In contrast to the triangular lattice, the relative promi-nence of non-affine modes in the planar honeycomb isstrongly dependent on the value of the bending constant. Fig. 7 shows the first six eigenvalues against k b . Severalcrossovers among the different modes are visible in thesespectra. We notice that the SW mode, which earlier wasthe second most prominent mode in the system, becomesstrongly suppressed as one increases k b . This was to beexpected because the SW defect requires that nearestneighbor bonds become flexible. Our projection formal-ism is hence very general and can pick out the dominantdefect precursor modes for arbitary lattice symmetry andinteractions. FIG. 8. a. Schematic of the Kagome lattice and coarse-graining volume Ω (pink shaded region), b. the spectrumof non-affine modes and c. k = 0 . β = 100. D. Kagome
The Kagome lattice structure is found in many naturalminerals and has interested physicists and chemists be-cause of its unusual magnetic properties [1, 37]. Similarto the planar honeycomb, a kagome lattice has a trian-gular symmetry, but with three basis particles in eachcell. Each particle in the cell has four nearest neighbourand four next nearest neighbours. The dynamical matrix(Appendix A) D ( q ) becomes a 6 × P becomes a 42 ×
42 matrix. Accordingly,
PCP has 38non-zero eigenvalues corresponding to non-affine eigen-modes. The probability distribution P ( χ ) is shown inFig 2. Fig 8 shows the eigenvalue spectra. Similarly tothe planar honeycomb we notice the absence of any largegap among the eigenvalues. The non-affine modes forthe largest eigenvalues are shown in Fig. 8. These modesturn out to be the well known floppy modes [38]. If thenext nearest neighbor bonds are stiffened or bond an-gle dependent potentials are introduced, the amplitudesof these floppy modes decrease, exactly as in the honey-comb lattice.0 FIG. 9. a. Plot of the eigenvalues of the SC lattice for threedifferent values of k : 0 .
25, 0 . . b. The non-affinemode corresponding to the largest eigenvalue for k = 0 . β = 1000. Note that this is similar to what is obtained forthe square lattice. E. Simple cubic
Our discussion of lattices in three dimensions beginswith the simple cubic (SC) lattice, having a single basisatom in a cubic cell with six nearest neighbor and twelvenext nearest neighbor particles. We assume that theseparticles are connected by springs with stiffness constant k for nearest neighbors and k for next nearest neigh-bors. The dynamical matrix (Appendix A) can be cal-culated and has three acoustic branches comprising onelongitudinal and two transverse phonon modes. We pro-ceed in a similar fashion as in 2d to calculate C . Theprojection matrix has 54 eigenvalues out of which 9 arezero corresponding to nine affine modes in 3 d . Similarto the triangular and square lattices in 2d, we find thata large gap exists between the largest eigenvalue of PCP and the rest, see Fig. 9. For the SC, we find that three de-generate modes correspond to this largest eigenvalue, oneof which is shown In Fig. 9. Note that the displacementpattern in the blue shaded plane in SC is similar to thatin the square lattice. Indeed the SC lattice may be re-garded as a stacking of 2d square lattices. The other twodegenerate modes show the same movement in the othertwo orthogonal planes of the SC. This leads to the in-terpretation that the most prominent non-affine mode ofthe SC lattice simply tend to convert the stacked planesfrom square to triangular symmetry, hence generating 3dclose packed structures [1].The eigenvalue spectra in Fig. 9 also show that χ de-creases as one stiffens the lattice by increasing the stiff-ness constant k . However, this increase in stiffness doesnot affect the qualitative features of the spectrum includ-ing the continuing presence of a gap. FIG. 10. a. Plot of the non-affine eigenvalues of the BCClattice for three different values of k : 0 .
5, 0 .
75 and 1 . β = 1000. b. One of the non-affine modes with the largesteigenvalue for k = 0 . F. BCC
The body centered cubic (BCC) lattice may be thoughtof either as a system with a single-atom basis, or as a SClattice with a two atom basis [1]. For reasons of computa-tional simplicity we choose the former view to constructour coarse-graining volume: this consists of 14 particleswith 8 nearest neighbors having bond stiffness k = 1 and6 next nearest neighbors with bond stiffness k . Since Ωcomprises 14 particles, in three dimensions P becomes a42 ×
42 matrix, and has 33 non-zero eigenvalues corre-sponding to the non-affine part.After performing the projection analysis, we find thata gap below the largest eigenvalue is present regardlessof the choice of k . We notice that as for all other lat-tices discussed above, (cid:104) χ (cid:105) decreases with an increase inthe stiffness of the lattice, see Fig. 10. BCC has three de-generate modes related to the largest eigenvalue. One ofthese is shown in Fig. 10, where we notice that the centreparticle has moved along the [001] direction. The othertwo modes show a displacement of the centre particle inthe two orthogonal directions. These dominant modestogether represent the motion of the body centered par-ticle to one of the six faces of the cubic unit cell, whichcan be viewed as generating locally a single atomic planeof the FCC lattice by an atomic shuffle. G. FCC
The coarse-graining volume for the face centred cubic(FCC) lattice, we construct around a single atom ba-sis, similar to the BCC case. It consists of 12 nearestneighbor particles and 6 next neighbors. We thus have18 × PCP , of which 9 eigenval-ues representing affine deformations are zero. The eigen-value spectrum again shows a prominent gap betweenthe three largest degenerate (and mutually orthogonal)eigenmodes and the rest. It is also obvious from the spec-1 [111]
FIG. 11. a. Plot of eigenvalue spectra of FCC for differentchoice of k : 0 .
25, 0 . . β = 1000. b. The eigenmodecorresponding to the largest eigenvalue for k = 0 . c. Sameas b. but viewed from the [111] direction. Notice that thecentral particle has displaced out of plane and sits below aparticle from a different stacking layer, resulting in a stackingfault in the FCC system. tra that (cid:104) χ (cid:105) decreases as one increases the stiffness of thelattice by increasing k .One of the three non-affine modes corresponding tothe largest eigenvalue is shown in Fig. 11. We show laterthat this mode is a precursor to either a slip or a stackingfault [1, 4]. The other two degenerate modes show theanalogous deformation in orthogonal directions. H. Affine-non-affine coupling
While the affine and non-affine components of the dis-placements are orthogonal to each other by construction,they couple at higher order [14]. This has the physi-cal meaning of suggesting that at higher strains, fluctu-ations which tend to create lattice defects become moreprobable. While this has been noted in the triangularlattice [14, 22], here we undertake a systematic study in-volving many lattices.In Section II we showed that this coupling is deter-mined by the commutator [ P , C ]. We have computed thiscommutator for all the lattices considered in this paperand the results are shown in Table. I. It is interestingto see that open lattices like the planar honeycomb andkagome have larger values of this coupling than the moreclose packed ones. I. Spatial correlations
We now look at two-point spatial correlations of χ andthe affine strains e . These have been extensively stud-ied for two dimensional lattices both numerically [27, 28] Lattice Type || C || || [ P,C ] |||| C || Parameters ( k , k b )2d triangle 3.874 0.031 0.5, 02d triangle 1.290 0.013 0.5, 0.52d square 5.055 0.037 0.5, 02d honeycomb 9.362 0.160 0.5, 02d honeycomb 4.897 0.135 0.5, 0.52d kagome 9.910 0.113 0.5, 03d SC 9.361 0.024 0.5, 03d BCC 6.994 0.020 0.75, 03d FCC 7.279 0.025 0.5, 0TABLE I. The Frobenius norm (the square root of the sumof the absolute squares of the elements) of the commutator[ P , C ], made dimensionless by dividing it by the correspondingnorm of C , for a number of lattices in 2d and 3d at β = 1.Corresponding parameter values for stiffness are quoted in thelast column. Also note that the norm of P is essentially thesquare-root of total number of non-affine modes in each case.FIG. 12. Normalized χ correlations, C χχ ( ρ ) = ( (cid:104) χ (0) χ ( ρ ) (cid:105) −(cid:104) χ (0) (cid:105) ) / (cid:0) (cid:104) χ (0) (cid:105) − (cid:104) χ (0) (cid:105) (cid:1) , for several different lattices asa function of distance ρ = R · ˆ x /d nn measured in units of the nearest neighbor distance , d nn in the reference lattice along onecoordinate axis. Brown and blue are for square and trianglelattices ( d nn = a ); orange and red are for FCC ( d nn = a/ √ d nn = a √ /
2) respectively. and analytically [14–16]. The spatial correlation of theaffine strain is important because it offers a way to ob-tain elastic properties of colloidal crystals from opticalmicroscopy images [26].The spatial correlations of χ for some of the latticesconsidered in this paper are shown in Fig. 12 in a singleplot. These correlations are nearly isotropic and are plot-ted as a function of distance ρ measured in the units ofnearest neighbour distance d nn along one coordinate axis.The values of d nn for different lattices are mentioned in2Fig 12. We observe that the nature of the correlationfunction is similar for all lattices. It is a sharply de-caying function that essentially vanishes after the secondneighbor shell. More quantitatively, we observe that thecorrelations decay somewhat faster in higher dimensions.The spatial correlations for the affine strains may beobtained using the procedure outlined in Section II.These have a more non-trivial structure. They areanisotropic and can be long-ranged along particular di-rections [28]. In the q → β (cid:104) e v (cid:105) ( q ) = Q v /Q , β (cid:104) e u (cid:105) ( q ) = Q u /Q , and β (cid:104) e s (cid:105) ( q ) = Q s /Q with the ab-breviations Q = q x q y + k (cid:0) q x + q y (cid:1) + k (cid:0) q x − q y (cid:1) Q v = 2 q x q y + k (cid:0) q x − q y (cid:1) Q u = 2 q x q y + k (cid:0) q x + 6 q x q y + q y (cid:1) Q s = (cid:0) q x + q y (cid:1) + k (cid:0) q x − q y (cid:1) (31)Similar expression for the triangular lattice have al-ready been discussed in [14]. We observe that for k = 1 /
2, where the square lattice becomes elasticallyisotropic [39], the expressions in Eq. (31) differ from thosefor the triangular lattice only by an unimportant overallfactor. The general shape of these correlation functions,viz. the “butterfly pattern”, is also similar to results ob-tained in colloidal glasses using video microscopy tech-niques [40].In three dimensions, the correlation functions are con-siderably more complicated, although analytic expres-sions in Fourier space in the small wave-number limitcan still be worked out with some effort. The algebraicexpressions are given in Appendix B and they are plottedin Fig. 13.
V. DISCUSSION AND CONCLUSIONS
In this paper we have studied the nature of thermallyexcited non-affine atomic displacements for a number ofcrystalline solids in 2d and 3d. We have discovered sev-eral features that are common to many lattice systemsalthough close packed and open lattices show somewhatdifferent properties. While in close packed lattices, thecontribution to χ is dominated by a single non-affinemode (or degenerate, symmetry-related class of modes),in open lattices there is no such predominance. Further,in open lattices, the contribution of the different modesis much more sensitive to details of the interactions andthey are more strongly coupled to affine fluctuations.One of the important findings of our earlier work [15,16] was that non-affine displacement fluctuations serveas precursors to the formation of defects. In the trian-gular lattice, in the presence of strain, a dislocation-anti-dislocation pair separates and produces a slip plane [22] that has high values of χ . We end this paper by carryingout a simple exercise in order to better understand therelation between non-affine modes and defects.Accordingly, we first consider a triangular lattice wherea slip is introduced such that a part of the lattice movesa lattice spacing in a close packed direction, comparedto the rest. In the bulk, there is no contribution to χ as all atoms undergo either no motion or just a uniformtranslation. Thermal vibrations are neglected in this cal-culation and those that follow. We choose a coarse grain-ing volume corresponding to the smallest Ω as shown inSec. IV centered on an atom lying in the slip plane. TheΩ at the interface of the slipped and un-slipped regions is,of course, deformed (see Fig. 14). This deformation can-not be described by a homogeneous affine transformationof Ω and therefore contributes to χ . Including thermalcontributions would produce a P ( χ ) that is identical tothe ones calculated in the bulk (Section IV), while in thevicinity of the slip, P ( χ ) would be displaced to higher χ values [22].One can now project this deformation onto the non-affine and affine modes computed from thermal averagesof displacements to find the contribution of individualmodes to this deformation. In Fig. 14, we plot the bar-graph of the components ( c i ) obtained by projectingonto the affine and non-affine modes for this deforma-tion, where c i is the coefficient corresponding to the i th mode in the expansion of the displacement as a super-position of non-affine and affine modes. We see that thelargest contribution comes from the first two non-affinemodes as expected. There is also a non-zero contribu-tion from the affine modes, with the affine and non-affinemodes contributing equally overall. This may be eas-ily understood from the insets shown in Fig. 14 a and b . In Fig. 14 a (inset i), we show the configuration ofparticles where the two atoms belonging to the bottom-most row are displaced to the left by a lattice spacingrelative to the upper two rows. The total non-affine con-tribution is shown in inset ii of the same figure. Thisshows a relative displacements to the right of the mid-dle row consisting of three atoms. On the other hand,the affine contribution to the slip shown in Fig. 14 b (in-set) consists of homogeneous shear and local rotation asverified from the bar-graph. It is clear that the sum ofthe affine and non-affine displacements gives rise to theslipped configuration shown in Fig. 14 a (inset i). Theaffine deformation produces an internal shear stress atthe slip plane. When a crystal slips in response to an ex-ternal homogeneous shear, the internal stress cancels theexternal stress locally. By introducing a finite density ofsuch slip planes, any homogenous stress may be expelled.In Ref. [22] such an expulsion process was shown to leadto yielding of crystalline solids at any shear stress, how-ever small. Since the non-affine strains corresponding tothe largest eigenvalues do not depend on the choice of Ω(see Section IV), the mechanism described is quite gen-eral.We now turn to the square lattice. It is known [21]3 a.b.c. SC BCC FCC
FIG. 13. Iso-surfaces for the strain-strain correlation functions in Fourier space and in the q → a. deviatoric ( e xx − e yy − e zz ), b. shear ( e xy + e yx ) and c. volume ( e xx + e yy + e zz ) strains. Thevalues of the correlations at the iso-surfaces are different in each case and have been chosen for ease of presentation. They arelisted in the Appendix B along with the full algebraic expressions used to plot the iso-surfaces. The other parameters used are k = 1 / β = 1 throughout. i. ii. FIG. 14. a. Relative contribution from non-affine modeswhen a slip is introduced by translating the bottom half of atriangular lattice along a close-packed direction by a latticespacing. b. The corresponding contribution of the four affinemodes. Note that the largest contributions come from shearand rotation. The insets show the configuration of atoms inΩ after the slip ( a -i) and the separate non-affine ( a -ii) andaffine ( b ) contributions. See text for details. that a square lattice can transform to a triangular lat-tice by either a homogeneous, affine, shear or by a non- FIG. 15. Non-affine contribution in the square lattice whenthe middle row is displaced by half a lattice spacing. Notethat this transformation tends to produce a triangular latticesymmetry starting from the original square lattice. FIG. 16. Contribution from different a. non-affine and b. affine modes when one introduces a slip in the FCC lattice.The inset shows a schematic diagram of the lattice as seenfrom the [111] direction. The original position of the closepacked lattice plane is shown as a triangle with a blue dottedboundary and the final position as a blue shaded triangle.Contribution from a stacking fault from c. non-affine and d. affine modes. The inset shows, as before, the original andfinal positions of a close packed plane of atoms. affine deformation where alternate rows of atoms shift byhalf a lattice spacing together with a homogeneous relax-ation of the lattice parameters. We show this deforma-tion in Fig. V(inset) and compute the projection onto thenon-affine modes. As expected, there is an overwhelm-ingly large contribution from the non-affine mode withthe largest eigenvalue. There is no affine component forthis deformation.For the FCC lattice, we first create a slip along oneof the close-packed planes, similar to the triangular case.Fig 16 shows the original position (triangle with a dottedboundary) and the new position (blue shaded triangle)of the closed packed plane. The bar-graph of c i corre-sponding to non-affine and affine modes makes it clearthat slip in the FCC lattice behaves similarly to a slipin the triangular lattice: we also observe here that thetotal contributions of the affine and non-affine modes areequal. In the 3d FCC lattice apart from a slip, one canalso consider a stacking fault. Fig 16 also shows this de-formation where the closed packed plane is displaced byhalf a lattice parameter. The shaded blue region in Fig16 represents the new position of the closed packed plane.Again during a stacking fault it is observed that the non-affine and affine parts contribute equally and the maxi-mum contribution comes from the first three non-affine modes.This simple exercise therefore strengthens our claimthat the non-affine modes that belong to the largesteigenvalues of PCP are related to fluctuations that tendto nucleate defects. We have seen in [22] that these fluc-tuations condense under external strain to cause plasticdeformation in a 2d triangular crystal. The computationspresented here show that non-affine modes deduced fromsmall harmonic fluctuations are able to describe impor-tant processes that occur during large deformations. Wehope that this knowledge will enable us to study in detailmechanical properties of 3d cubic lattices. Preliminarycalculations are under way and will be published else-where.We believe that our work brings out an interesting as-pect concerning defects and the dynamics of deformationin crystalline solids [1, 4, 5]. Atomic fluctuations thatgenerate defects in close packed crystals are shown tobe determined by the non-affine modes with the largesteigenvalue. Representing fluctuations in crystals as con-sisting of smooth phonons and singular defects thereforeamounts to making a “largest eigenvalue approximation”and neglecting other non-affine modes that make smallercontributions to the total χ . This approximation lies atthe heart of all dislocation based theories of crystal plas-ticity [4, 5, 41, 42]. Such an approximation is excellentwhen the defect-like mode is separated from the others bya large gap in the spectrum of non-affine eigenvalues asin the triangular (Section IV A) and FCC (Section IV G)lattices. However, the approach may not work for crystallattices where such a gap does not exist or is too small,e.g. for the planar honeycomb (Section IV C) or kagome(Section IV D) structures. In amorphous matter also thisapproximation may not be so useful, even if dislocation-like structures are identifiable [11]. In such cases, newcontinuum theories of deformation that include all non-affine modes (or at least a large class of them) may beneeded. Such a theory does not exist at present and wehope that our work provides sufficient motivation to thecommunity for thinking along these lines. ACKNOWLEDGMENTS
We acknowledge useful discussions with Saswati Gan-guly, Parswa Nath, Madan Rao, Srikanth Sastry, ItamarProcaccia and J¨urgen Horbach.5
Appendix A: Dynamical Matrices
We give below the expressions for the dynamicalmatrices for all the lattices studied in this paper. Foreach case we have included nearest neighbor and nextnearest neighbor bonds, with spring constants k and k , respectively. For the triangular lattice only nearestneighbor bonds were retained ( k = 0). Separateexpressions are given for the additional terms thatresult in the triangular and honeycomb lattices from theintroduction of additional bond-bending terms.
1. Square A = 2( k + k − cos( aq x )( k + k cos( aq y ))) A = 2 k sin( aq x ) sin( aq y ) = A A = 2( k + k − cos( aq y )( k + k cos( aq x ))) D ( q ) = (cid:18) A A A A (cid:19) (A1)
2. Triangular B = k (cid:18) − aq x ) − cos (cid:16) aq x (cid:17) cos (cid:32) √ aq y (cid:33) (cid:19) B = k (cid:18) √ (cid:16) aq x (cid:17) sin (cid:32) √ aq y (cid:33) (cid:19) = B B = k (cid:18) − (cid:16) aq x (cid:17) cos (cid:32) √ aq y (cid:33) (cid:19) D ( q ) = (cid:18) B B B B (cid:19) (A2)
3. Triangular with bond bending: B (cid:48) = 3 k b (cid:20) − cos( aq x ) − (cid:16) aq x (cid:17) cos (cid:32) √ aq y (cid:33) (cid:21) = B (cid:48) B (cid:48) = 0 = B (cid:48) D ( q ) = (cid:18) B B B B (cid:19) (A3)
4. Planar honeycomb W = 3 k k − k cos (cid:18) aq x (cid:19) cos (cid:32) √ aq y (cid:33) W = √ k sin (cid:18) aq x (cid:19) sin (cid:32) √ aq y (cid:33) W = k − exp ( iaq x ) − exp (cid:0) − iaq x (cid:1) cos (cid:16) √ aq y (cid:17) W = k (cid:32) √ i exp (cid:18) − iaq x (cid:19) sin (cid:32) √ aq y (cid:33)(cid:33) W = 3 k k − k cos (cid:16) √ aq y (cid:17) − k cos (cid:18) aq x (cid:19) cos (cid:32) √ aq y (cid:33) W = k (cid:32) − (cid:18) − iaq x (cid:19) cos (cid:32) √ aq y (cid:33)(cid:33) D ( q ) = W W W W W ∗ W W W W ∗ W ∗ W W W ∗ W ∗ W ∗ W (A4)6
5. Honeycomb with bond bending: W (cid:48) = 3 k b W (cid:48) = i √ k b exp (cid:18) − iq x a (cid:19) sin (cid:32) √ q y a (cid:33) W (cid:48) = − k b exp (cid:18) − iq x a (cid:19) cos (cid:32) √ q y a (cid:33) W (cid:48) = − i √ k b exp (cid:18) − iq x a (cid:19) sin (cid:32) √ q y a (cid:33) W (cid:48) = k b (cid:32) (cid:18) q x a (cid:19) cos (cid:32) √ q y a (cid:33)(cid:33) W (cid:48) = − k b exp (cid:18) − iq x a (cid:19) (cid:32) e iqxa + cos (cid:32) √ q y a (cid:33)(cid:33) W (cid:48) = − k b (cid:32) q x a ) + exp (cid:18) iq x a (cid:19) cos (cid:32) √ q y a (cid:33) − i sin( q x a ) (cid:33) D ( q ) = W (cid:48) W (cid:48) W (cid:48) W (cid:48) W (cid:48)∗ W (cid:48) W (cid:48) W (cid:48) W (cid:48)∗ W (cid:48)∗ W (cid:48) W (cid:48)∗ W (cid:48)∗ W (cid:48) W (cid:48) W (cid:48) (A5)
6. Kagome M = k + 3 k M = 0 M = − k exp (cid:18) − i (cid:0) aq x √ aq y (cid:1)(cid:19) − k exp (cid:18) − i (cid:0) − aq x − √ aq y (cid:1)(cid:19) − k exp (cid:18) − i (cid:0) − aq x √ aq y (cid:1)(cid:19) − k exp (cid:18) − i (cid:0) aq x − √ aq y (cid:1)(cid:19) M = − √ k exp (cid:18) − i (cid:0) aq x √ aq y (cid:1)(cid:19) − √ k exp (cid:18) − i (cid:0) − aq x − √ aq y (cid:1)(cid:19) + 14 √ k exp (cid:18) − i (cid:0) − aq x √ aq y (cid:1)(cid:19) + 14 √ k exp (cid:18) − i (cid:0) aq x − √ aq y (cid:1)(cid:19) M = − k exp (cid:18) − i (cid:0) − aq x √ aq y (cid:1)(cid:19) − k exp (cid:18) − i (cid:0) aq x − √ aq y (cid:1)(cid:19) − k exp (cid:18) − i (cid:0) aq x √ aq y (cid:1)(cid:19) − k exp (cid:18) − i (cid:0) − aq x − √ aq y (cid:1)(cid:19) M = 14 √ k exp (cid:18) − i (cid:0) − aq x √ aq y (cid:1)(cid:19) + 14 √ k exp (cid:18) − i (cid:0) aq x − √ aq y (cid:1)(cid:19) − √ k exp (cid:18) − i (cid:0) aq x √ aq y (cid:1)(cid:19) − √ k exp (cid:18) − i (cid:0) − aq x − √ aq y (cid:1)(cid:19) M = 3 k + k M = −
14 3 k exp (cid:18) − i (cid:0) aq x √ aq y (cid:1)(cid:19) − k exp (cid:18) − i (cid:0) − aq x − √ aq y (cid:1)(cid:19) − k exp (cid:18) − i (cid:0) − aq x √ aq y (cid:1)(cid:19) − k exp (cid:18) − i (cid:0) aq x − √ aq y (cid:1)(cid:19) M = −
14 3 k exp (cid:18) − i (cid:0) − aq x √ aq y (cid:1)(cid:19) − k exp (cid:18) − i (cid:0) aq x − √ aq y (cid:1)(cid:19) − k exp (cid:18) − i (cid:0) aq x √ aq y (cid:1)(cid:19) − k exp (cid:18) − i (cid:0) − aq x − √ aq y (cid:1)(cid:19) M = 5 k k M = √ k − √ k M = − k exp (cid:0) i aq x (cid:1) − k exp (cid:0) − i aq x (cid:1) M = 3 k k M = − k exp (cid:0) i √ aq y (cid:1) − k exp (cid:0) − i √ aq y (cid:1) M = √ k − √ k D ( q ) = M M M M M M M ∗ M M M M M M ∗ M ∗ M M M M M ∗ M ∗ M ∗ M M M M ∗ M ∗ M ∗ M ∗ M M M ∗ M ∗ M ∗ M ∗ M ∗ M (A6)
7. SC S = − aq x ) ( k + k (cos ( aq y ) + cos ( aq z )))+ 2 k + 4 k S = 2 k sin ( aq x ) sin ( aq y ) = S S = 2 k sin ( aq x ) sin ( aq z ) = S S = − aq y ) ( k + k (cos ( aq x ) + cos ( aq z )))+ 2 k + 4 k S = 2 k sin ( aq y ) sin ( aq z ) = S S = − aq z ) ( k + k (cos ( aq x ) + cos ( aq y )))+ 2 k + 4 k D ( q ) = S S S S S S S S S (A7)8
8. BCC: G = − k cos (cid:16) aq x (cid:17) cos (cid:16) aq y (cid:17) cos (cid:16) aq z (cid:17) − k cos ( aq x ) + 83 k + 2 k G = 83 k cos (cid:16) aq z (cid:17) sin (cid:16) aq x (cid:17) sin (cid:16) aq y (cid:17) = G G = 83 k cos (cid:16) aq y (cid:17) sin (cid:16) aq x (cid:17) sin (cid:16) aq z (cid:17) = G G = − k cos (cid:16) aq x (cid:17) cos (cid:16) aq y (cid:17) cos (cid:16) aq z (cid:17) − k cos ( aq y ) + 83 k + 2 k G = 83 k cos (cid:16) aq x (cid:17) sin (cid:16) aq y (cid:17) sin (cid:16) aq z (cid:17) = G G = − k cos (cid:16) aq x (cid:17) cos (cid:16) aq y (cid:17) cos (cid:16) aq z (cid:17) − k cos ( aq z ) + 83 k + 2 k D ( q ) = G G G G G G G G G (A8)
9. FCC F = 4 k + 2 k − k cos ( aq x ) − k cos (cid:16) aq x (cid:17) (cid:16) cos (cid:16) aq y (cid:17) + cos (cid:16) aq z (cid:17)(cid:17) F = 2 k sin (cid:16) aq x (cid:17) sin (cid:16) aq y (cid:17) = F F = 2 k sin (cid:16) aq x (cid:17) sin (cid:16) aq z (cid:17) = F F = 4 k + 2 k − k cos ( aq y ) − k cos (cid:16) aq y (cid:17) (cid:16) cos (cid:16) aq x (cid:17) + cos (cid:16) aq z (cid:17)(cid:17) F = 2 k sin (cid:16) aq y (cid:17) sin (cid:16) aq z (cid:17) = F F = 4 k + 2 k − k cos ( aq z ) − k cos (cid:16) aq z (cid:17) (cid:16) cos (cid:16) aq x (cid:17) + cos (cid:16) aq y (cid:17)(cid:17) D ( q ) = F F F F F F F F F (A9) Appendix B: Strain correlation iso-surfaces for 3d lattices
We give below the equations for the iso-surfaces of the strain-strain correlation functions for the 3d lattices as shownin Fig. 13. We have used the same notation as in Section IV I.
1. SC
The iso-surfaces are given by the equations β (cid:104) e v (cid:105) ( q ) = Q SCv /Q SC = 1, β (cid:104) e u (cid:105) ( q ) = Q SCu /Q SC = 1 . β (cid:104) e s (cid:105) ( q ) = Q SCs /Q SC = 3 .
6, where, Q SC = a (cid:20) k q x q y q z + k k (cid:18) q x (cid:0) q y + q z (cid:1) + q x (cid:0) q y + 6 q y q z + q z (cid:1) + q y q z (cid:0) q y + q z (cid:1)(cid:19) + k k (cid:18) q x + 5 q x (cid:0) q y + q z (cid:1) + q x (cid:0) q y + 3 q y q z + 5 q z (cid:1) + q y + 5 q y q z + 5 q y q z + q z (cid:19) + k (cid:18) q x + 3 q x (cid:0) q y + q z (cid:1) + q x (cid:0) q y + 8 q y q z + 3 q z (cid:1) + 2 q y + 3 q y q z + 3 q y q z + 2 q z (cid:19)(cid:21) (B1a)9 Q SCv = 3 k q x q y q z + 2 k k (cid:18) q x (cid:0) q y + q z (cid:1) + q x (cid:0) q y + q z (cid:1) + q y q z (cid:0) q y + q z (cid:1)(cid:19) + k (cid:18) q x + q x (cid:0) q y + q z (cid:1) + q x (cid:0) q y + 3 q y q z + q z (cid:1) + q y + q y q z + q y q z + q z (cid:19) (B1b) Q SCu = 3 k q x q y q z + 2 k k (cid:18) q x (cid:0) q y + q z (cid:1) + q x (cid:0) q y + 8 q y q z + q z (cid:1) + q y q z (cid:0) q y + q z (cid:1)(cid:19) + k (cid:18) q x + 9 q x (cid:0) q y + q z (cid:1) + 3 q x (cid:0) q y + q y q z + 3 q z (cid:1) + q y + q y q z + q y q z + q z (cid:19) (B1c) Q SCs = k q z (cid:0) q x + q y (cid:1) + k k (cid:18) q x + q x (cid:0) q y + 4 q z (cid:1) + q x (cid:0) q y − q z (cid:1) + q y (cid:0) q y + 4 q y q z + q z (cid:1)(cid:19) + k (cid:18) q x + q x q z + 2 q x (cid:0) q y q z + q z (cid:1) +2 q y + q y q z + 2 q y q z (cid:19) (B1d)
2. BCC
The iso-surfaces are given by the equations β (cid:104) e v (cid:105) ( q ) = Q BCCv /Q BCC = 1 . β (cid:104) e u (cid:105) ( q ) = Q BCCu /Q BCC = 5 . β (cid:104) e s (cid:105) ( q ) = Q BCCs /Q BCC = 2 .
15, where, Q BCC = a (cid:20) k (cid:18) − q z (cid:0) q x + q y (cid:1) + (cid:0) q x − q y (cid:1) (cid:0) q x + q y (cid:1) − q z (cid:0) q x − q x q y + q y (cid:1) + q z (cid:19) + 3 k k (cid:18) q x + 3 q x (cid:0) q y + q z (cid:1) + 3 q x (cid:0) q y − q z (cid:1) + (cid:0) q y + q z (cid:1) (cid:19) + 9 k k (cid:0) q x + q y + q z (cid:1)(cid:18) q x (cid:0) q y + q z (cid:1) + q y q z (cid:19) + 27 k q x q y q z (cid:21) (B2a) Q BCCv = 3 k (cid:0) q x + q y (cid:1)(cid:18) k q x q y + k (cid:0) q x − q y (cid:1) (cid:19) − (cid:18) k ( k − k ) q x − (cid:0) k − k k + 9 k (cid:1) q x q y + k ( k − k ) q y (cid:19) q z − k ( k − k ) (cid:0) q x + q y (cid:1) q z + 3 k q z (B2b) Q BCCu = 3 (cid:20) k q x q y q z + 6 k k (cid:18) q x (cid:0) q y + q z (cid:1) + q y q z (cid:0) q y + q z (cid:1) + q x (cid:0) q y + 5 q y q z + q z (cid:1)(cid:19) + k (cid:18) q x + 7 q x (cid:0) q y + q z (cid:1) + (cid:0) q y − q z (cid:1) (cid:0) q y + q z (cid:1) + q x (cid:0) q y − q y q z + 7 q z (cid:1)(cid:19)(cid:21) (B2c) Q BCCs = 3 (cid:20) k (cid:0) q x + q y (cid:1) q z + 3 k k (cid:18) q x + q x (cid:0) q y − q z (cid:1) + q y (cid:0) q y + q z (cid:1) + q x (cid:0) q y + 2 q z (cid:1)(cid:19) + k (cid:18) q x − q x (cid:0) q y + 2 q z (cid:1) + (cid:0) q y − q y q z (cid:1) + q x (cid:0) − q y + 8 q y q z + q z (cid:1)(cid:19)(cid:21) (B2d)
3. FCC
The iso-surfaces are given by the equations β (cid:104) e v (cid:105) ( q ) = Q F CCv /Q F CC = 1 . β (cid:104) e u (cid:105) ( q ) = Q F CCu /Q F CC = 3 . β (cid:104) e s (cid:105) ( q ) = Q F CCs /Q F CC = 2 .
75, where, Q F CC = a (cid:20) k (cid:18) q x + 3 q x (cid:0) q y + q z (cid:1) + q x (cid:0) q y + 8 q y q z + 3 q z (cid:1) +2 q y + 3 q y q z + 3 q y q z + 2 q z (cid:19) + 4 k k (cid:18) q x + 5 q x (cid:0) q y + q z (cid:1) + q x (cid:0) q y + 3 q y q z + 5 q z (cid:1) + q y + 5 q y q z + 5 q y q z + q z (cid:19) + 16 k k (cid:18) q x (cid:0) q y + q z (cid:1) + q x (cid:0) q y + 6 q y q z + q z (cid:1) + q y q z (cid:0) q y + q z (cid:1)(cid:19) + 64 k q x q y q z (cid:21) (B3a)0 Q F CCv = 4 (cid:20) k (cid:18) q x + q x (cid:0) q y + q z (cid:1) + q x (cid:0) q y + 3 q y q z + q z (cid:1) + q y + q y q z + q y q z + q z (cid:19) + 8 k k (cid:18) q x (cid:0) q y + q z (cid:1) + q x (cid:0) q y + q z (cid:1) + q y q z (cid:0) q y + q z (cid:1)(cid:19) + 48 k q x q y q z (cid:21) (B3b) Q F CCu = 4 (cid:20) k (cid:18) q x + 9 q x (cid:0) q y + q z (cid:1) +3 q x (cid:0) q y + q y q z + 3 q z (cid:1) + q y + q y q z + q y q z + q z (cid:19) + 8 k k (cid:18) q x (cid:0) q y + q z (cid:1) + q x (cid:0) q y + 8 q y q z + q z (cid:1) + q y q z (cid:0) q y + q z (cid:1)(cid:19) + 48 k q x q y q z (cid:21) (B3c) Q F CCs = 4 (cid:20) k (cid:18) q x + q x q z + 2 q x (cid:0) q y q z + q z (cid:1) +2 q y + q y q z + 2 q y q z (cid:19) + 4 k k (cid:18) q x + q x (cid:0) q y + 4 q z (cid:1) + q x (cid:0) q y − q z (cid:1) + q y (cid:0) q y + 4 q y q z + q z (cid:1)(cid:19) + 16 k q z (cid:0) q x + q y (cid:1)(cid:21) (B3d) [1] P. Chaikin and T. Lubensky, Principles of CondensedMatter Physics (Cambridge Press, Cambridge, 1995).[2] M. Born and K. Huang,
The dynamical theory of crystallattices (Clarendon Press, Gloucestershire, 1998).[3] P. Martin, O. Parodi, P. Pershan, (1972)
Unified Hydro-dynamic Theory for Crystals, Liquid Crystals, and Nor-mal Fluids
Phys. Rev. A, 6:2401.[4] Phillips R (2004)
Crystals, defects and microstructures:Modeling across scales (Cambridge Press, Cambridge).[5] J. P. Hirth and J. Lothe,
Theory of Dislocations (McGraw-Hill, 1967).[6] F. R. N. Nabarro and M. S. Duesbery, Eds.,
Dislocationsin Solids, Vol. 10 (Elsevier, Amsterdam, 1996), p.505-594.[7] N. W. Ashcroft and N. D. Mermin,
Solid State Physics ,(Holt, Rinehart, and Winston, New York, 1976)[8] J. P. Sethna et al. (2017) Deformation of Crystals: Con-nections with Statistical Physics
Annu. Rev. Mater. Res.
Annu. Rev. Con-dens. Matter Phys.
Soft Matter ,13:5649.[11] A. Acharya and M. Widom (2017) A microscopic con-tinuum model for defect dynamics in metallic glasses,
J.Mech. Phys. Solids , 104:1[12] A. Ghosh, Z. Budrikis, V. Chikkadi, A. L. Sellerio, S.Zapperi and P. Schall (2017) Direct Observation of Per-colation in the Yielding Transition of Colloidal GlassesPhys. Rev. Lett.118:148001.[13] Zaccone A, Schall P and Terentjev EM (2014) Micro-scopic origin of nonlinear nonaffine deformation in bulkmetallic glasses. Phys. Rev. B, 90:140203(R).[14] S. Ganguly, S. Sengupta, P. Sollich, and M. Rao, Phys.Rev. E , 042801 (2013).[15] S. Ganguly, S. Sengupta, and P. Sollich, Soft Matter ,4517 (2015). [16] A. Mitra, S. Ganguly, S. Sengupta, and P. Sollich, JS-TAT, P06025 (2015).[17] S. Ganguly, P. Nath, J. Horbach, P. Sollich, S. Karmakarand S. Sengupta, J. Chem. Phys. , 124501 (2017).[18] S. Ganguly and S. Sengupta, Excess vibrational modes ofa crystal in an external non-affine field , J. Chem. Sci. , 891 (2017).[19] S. Ganguly, P. S. Mohanty, P. Schurtenberger, S. Sen-gupta, and A. Yethiraj,
Contrasting the dynamics of elas-tic and non-elastic deformations across an experimentalcolloidal Martensitic transition , Soft Matter, , 4689,(2017).[20] S. Ganguly, D. Das, J. Horbach, P. Sollich, S. Karmakar,and S. Sengupta, Plastic deformation of a permanentlybonded network: Stress relaxation by pleats , J. Chem.Phys. , 184503 (2018).[21] P. Popli, S. Ganguly and S. Sengupta,
Translationallyinvariant colloidal crystal templates , Soft Matter, , 104(2018).[22] P. Nath, S. Ganguly, J. Horbach, P. Sollich, S. Karmakarand S. Sengupta, (2018) On the existence of thermody-namically stable rigid solids, Proc. Natl. Acad. Sci. USA , , E4322 .[23] Dube Dheeraj Prakashchand, Navjeet Ahalawat, Himan-shu Khandelia, Jagannath Mondal and Surajit Sengupta, On identifying collective displacements in apo- proteinsthat reveal eventual binding pathways , PLoS Comput Biol15(1): e1006665 (2019).[24] M. L. Falk and J. S. Langer, Phys. Rev. E , 7192 (1998).[25] J. Eriksen (2008) On the Cauchy–Born Rule, Math.Mech. Solids
Elastic properties of 2D colloidal crystals from video mi-croscopy
Phys. Rev. Lett. , 155506 (2003).[27] K. Franzrahe, P. Keim, G. Maret, P. Nielaba, S. Sen-gupta, Nonlocal elastic compliance for soft solids: The-ory, simulations and experiments , Phys. Rev. E, ,026106 (2008).[28] K. Franzrahe, P. Nielaba, and S. Sengupta, Coarse-graining microscopic strains in a harmonic, two- dimensional solid: Elasticity, nonlocal susceptibilities,and nonaffine noise , Phys. Rev. E , 016112 (2010)[29] J. C. Maxwell, Phil. Mag. , 294 (1864).[30] A. K. Geim, Science, 324, 1530 (2009).[31] S. K. Jain, G. T. Barkema, N. Mousseau, C-M. Fang, andM. A. van Huis, J. Phys. Chem. C, , 9646 (2015).[32] P. N. Keating, Phys. Rev. , 637 (1966).[33] D. Frenkel and B. Smit, Understanding Molecular Simu-lations (Academic Press, San Diego, 2002).[34] S. Plimpton, J. Comp. Phys, , 1 (1995);(http://lammps.sandia.gov). [35] A. Stone and D. Wales, Chem. Phys. Lett. , 501(1986).[36] J. G. Kirkwood, J. Chem. Phys. , 506 (1939).[37] M. Mekata, Physics Today, , 073901 (2015).[39] Kim, Y., Kim, Y. & Ryu, S. J Mech Sci Technol (2018) , 2693; L Monette and M P Anderson Modelling Simul.Mater. Sci. Eng. 2 53 (1994).[40] V. Chikkadi, G. Wegdam, D. Bonn, B. Nienhuis and P.Schall, Phys. Rev. Lett. , 198303 (2011)[41] J. Lubliner Plasticity Theory , (Dover, New Yourk, 2008)[42] A. Acharya
J. Mech. Phys. Solids ,49