Elasticity and Plasticity in Stiff and Flexible Oligomeric Glasses
Oleg Gendelman, H. George E. Hentschel, Pankaj K. Mishra, Itamar Procaccia, Jacques Zylberg
EElasticity and Plasticity in Stiff and Flexible Oligomeric Glasses
Oleg Gendelman, H. George E. Hentschel ∗ , Pankaj K. Mishra, Itamar Procaccia and Jacques Zylberg Department of Chemical Physics, The Weizmann Institute of Science, Rehovot 76100, Israel ∗ Department of Physics, Emory University, Atlanta GA, USA (Dated: October 13, 2018)In this paper we focus on the mechanical properties of oligomeric glasses (waxes), employing amicroscopic model that provides, via numerical simulations, information about the shear modulusof such materials, the failure mechanism via plastic instabilities and about the geometric responsesof the oligomers themselves to a mechanical load. We present a microscopic theory that explainsthe numerically observed phenomena, including an exact theory of the shear modulus and of theplastic instabilities, both local and system spanning. In addition we present a model to explain thegeometric changes in the oligomeric chains under increasing strains.
I. INTRODUCTION
A polymer is a macromolecule that consists of a largenumber of monomer subunits [1]. Polymeric glasses aresolids composed of a large number of such polymericunits. Subjected to homogeneous strain such solids canexhibit a variety of interesting phenomena including craz-ing instabilities, shear banding, strain hardening etc [2].Considerable effort was expended to describe these phe-nomena on the microscopic level using theory and sim-ulations [3–5]. Under tensile strains cavities may nucle-ate in a hitherto homogeneous polymeric glass. It wasargued that the formation of cavities takes place in re-gions of local low elastic modulus [6]. Polymeric glassessubjected to large strains exhibit strain hardening; thismay suppress strain localization and consequent crazing,necking, shear banding etc. Strain hardening is presum-ably caused by ordering the polymer beyond a certainstrain threshold. The microscopic origin of strain hard-ening was studied using molecular dynamic simulationsin Ref. [7, 8], finding that the origin of this phenomenon isrelated to plastic rearrangements of the monomers. Thisalso leads to short-range ordering. In spite of the abovementioned efforts a first-principles theory of these inter-esting phenomena is still incomplete. In particular inthis paper we propose a microscopic theory that relatesmacroscopic observables with the conformational defor-mation of the oligomers under pure shear.In recent years there has been great progress in under-standing the mechanical properties of amorphous solidsfrom first principles [9–11]. This progress was based onidentifying elementary plastic events as the loss of me-chanical stability when a Hessian eigenvalue hits zero [9–11]. This event is connected to a saddle node bifurcationin the generalized energy landscape. It was demonstratedalso that these elementary events can aggregate and con-catenate to yield shear localization and eventually shearbands [12, 13]. The aim of this paper is to extend this an-alytic approach to plasticity from simple Lennard-Jonesglasses (and recently some glasses with magnetic prop-erties) [14–16] to the realm of short oligomeric (or wax)glasses. These are amorphous solids whose constituentsare short chains of the order of 10-30 monomers, where the full impact of polymeric entanglement is still not cru-cial [17]. Nevertheless the existence of fairly long chainsof connected monomers introduces a hierarchy of newlength scales and energy scales related to valence bonds,valence angles and inter-oligomer interactions. In partic-ular the persistence length (cid:96) p of the oligomer turns out tobe crucial. Thus a variety of new phenomena and ques-tions arise, calling for a careful numerical simulation andanalytic assessment. Among the issues arising we willprovide a microscopic theory for the shear modulus ofthese materials, for the failure mechanism through plas-ticity (both local and system spanning) and shed lighton the geometric characteristic of the oligomers undermechanical yield.The outline of the paper is as follows: In Sect. II wedescribe the atomistic model used in further simulations.The model employs Lennard-Jones, angular and FENEinteractions (and see below for details). Sect. III presentsfirstly the results of numerical simulations for the stressvs. strain curves, the energy budget, characteristic of theoligomeric chains like end-to-end distance etc. For ana-lytic transparency we perform the simulation in quasi-static athermal conditions to highlight the plastic eventswithout any thermal fluctuations or strain rate effectsthat mask the fundamental physics. The same sectionprovides some theory of these characteristics. In Sect. IVwe present a theory for elementary plastic events. Nextin Sect. V we discuss the failure mechanism involvingshear localization and eventually shear bands. The fol-lowing section VI presents the analytic calculation of theshear modulus and a comparison with the numerics. II. DESCRIPTION OF THE MODEL
We consider a system composed of N p chains each com-prising n monomers (oligomers). Thus the total numberof particles in our system is N = N p × n . The interactionbetween monomers belonging to the same or to differentoligomers is different. Inter-oligomer interaction are sim-ply given by a truncated and smoothed Lennard-Jonespotential φ LJ , see below in Eq. (3). Within a givenoligomer the interactions have three contributors. First,all monomers within the Lennard-Jones cutoff range r co a r X i v : . [ c ond - m a t . s o f t ] J u l exert a force on each other which is derived from thepotential φ LJ . Secondly, a contribution χ is added tothe energy of any two successive monomers within thepolymer (to mimic the valence bond interaction). Thethird contribution to the energy is an angular potentialto constrain the value of the valence angle θ determinedby three successive monomers within a oligomer. Thisinteraction is denoted below ψ ( θ ). Thus the total energycan be written as [18] U = U LJ + U FENE + U Angle (1) U LJ = N (cid:88) (cid:104) ij (cid:105) φ ij LJ , U FENE = N p (cid:88) k =1 n − (cid:88) i =1 χ i,i +1 k U Angle = N p (cid:88) k =1 n − (cid:88) i =2 ψ i − ,i,i +1 k (2)The notation is such that successive particles are i and i +1 within a oligomer chain and ψ ik stands for the angularcontribution formed by any three successive particles ( i − , i, i + 1) within the k (cid:48) th oligomer where i is the vertex.The truncated and smoothed potential Lennard-Jonespotential is defined as: φ ij LJ = 4 ε (cid:34)(cid:18) λr ij (cid:19) − (cid:18) λr ij (cid:19) (cid:35) , r ij ≤ r min (3) φ ij LJ = ε (cid:34) a (cid:18) λr ij (cid:19) − b (cid:18) λr ij (cid:19) + (cid:88) (cid:96) =0 c (cid:96) (cid:16) r ij λ (cid:17) (cid:96) (cid:35) r min < r ij < r co , (4) φ ij LJ = 0 , r ij ≥ r co . (5)Here r min /λ is the length where the potential attains itsminimum, and r co /λ is the cut-off length for which thepotential vanishes. The coefficients a, b and c (cid:96) are cho-sen such that the repulsive and attractive parts of thepotential are continuous with two derivatives at the po-tential minimum and the potential goes to zero contin-uously at r co /λ with two continuous derivatives as well.The unit of length λ = 1 . ε is the unit of energy andthe Boltzmann constant k B = 1.For any two successive particles within the k (cid:48) th chainthere is the Finite Non-Elastic Elongation (FENE) po-tential with finite length r which is defined as: χ i,i +1 k ( r ) = (cid:40) − ηr ln[1 − ( r/r ) ] ; r < r ∞ ; r ≥ r (6)where r ≡ r i,i +1 /λ and η is a parameter with units offorce per unit length.Finally, for any three successive monomers within the k (cid:48) th oligomer with vertex i there is an angle constraintaround a chosen equilibrium angle ϕ eq and is defined as: ψ i − ,i,i +1 k ( ϕ i ) = κ [cos ϕ ik − cos ϕ eq ] α = κ [1 + cos ϕ ik ] α . (7) Below we will also employ the angle θ were θ ≡ π − ϕ .Thus for a stiff polymer ϕ ≈ π while θ is close to zero.We distinguish between two cases, that of a stiffoligomer with α = 1 and a semi-flexible oligomer with α = 2. The meaning of the words “stiff” and “semi-flexible” will be made clear in the sequel. The valuesof all the parameters used in the simulation are given inTable I. III. NUMERICAL SIMULATIONS
We prepare a 2-dimensional system consisting of 256polymers having 20 monomers in a chain. The initialdensity ρ = 0 . τ α LJ time units at pressure P=1.0 (LJ units), where τ α is the alpha relaxation time. After equilibration thepolymer melt is coupled to a heat-bath at temperatureT=0.01 (LJ units) and constant pressure (P=1). Thesystem is then further equilibrated for another 100 LJtime units. Finally the glass sample is taken to the near-est inherent minimum state using a conjugate gradientscheme. This protocol is referred to as “infinitely fast”quench. Below we also consider samples prepared by fi-nite quench rates.Having prepared the oligomeric glass sample it is sub-jected to an Athermal Quasi-Static strain (AQS) as de-scribed in detail in ref. [11]. In brief, each monomer isfirst displaced by the affine transformation x i → x i + δγy i ; y i → y i , (8)where r i ≡ ( x i , y i ) is the initial position of the i th monomer and δγ is the strain step applied during eachaffine transformation. The above transformation leadsto non-zero resultant forces on the monomers. Theseforces are annulled by a non-affine transformation r i → r i + u i , where u i is displacement of the monomer neces-sary to return to mechanical equilibrium. The non-affinedisplacement is computed using the conjugate gradientscheme. The strain step is chosen for the present studyis δγ = 10 − . The simulation is performed under peri-odic boundary condition along each direction of the boxusing the Lees-Edward formalism.We consider both stiff and semi-flexible polymers inour studies (see Eq. (7)). For stiff polymers the resultswill be presented for κ = 2 , , , and 15; for the semi-flexible case we consider κ = 2 , ,
8, 10. Unless statedspecifically the results reported below will refer to thestiff oligomer case with α = 1. TABLE I: The parameters used in the simulation.a b c0 c2 c4 c6 η r ϕ eq × − -0.0207 0.10691 -0.143794 30 1.5 π γ σ xy FIG. 1: A typical stress vs strain curve obtained by AQSstraining of 256 polymers of chain length 20 with stiffnessparameter κ = 2. Smooth (linear) increases in the strainare punctuated with sharp drops. The trajectory of stress vsstrain is reversible only until the first drop. The sharp dropsare plastic events as explained in the text. Inset: a blow-upof the first few plastic drops. A. Mechanical response of the polymer
A typical stress vs. strain curve that results in theAQS protocol is shown in Fig. 1 for a single realizationof a stiff oligomeric glass with κ = 2. The stress growslinearly at first with the strain, and the protocol can bereversed to return to initial state. Upon increasing thestrain the stress vs strain trajectory gets punctuated withsharp drops, these are irreversible, and after the occur-rence of the first one we cannot return to the initial stateby reversing the protocol. After each plastic drop thestress rises again linearly with the applied strain (butnot necessarily with the same slope) until the next plas-tic drop takes place. Generally speaking both the stressvs. strain and the energy vs. strain curves reach even-tually a kind of steady state in which the average stressand energy do no longer change even though they stillexperience elastic increases and plastic drops.At first, when the external strain is still small, the en-ergy drops associated with the plastic events are small,and do not increase with the system size. These plasticenergy drops are associated with localized events as is log(N) l o g ( ∆ σ x y ) data ( κ =2)−0.489496X+0.429642 log(N) l o g ( ∆ U ) data ( κ =2)0.487055X+0.567501 FIG. 2: log-log plots of the average magnitude of the stressdrops ∆ σ and energy drops ∆ U in the steady state (after theyield strain had been passed) as a function of the number ofparticles in the system explained in the next section. On the other hand, whenthe external strain is increased, at a threshold value ofthe external strain (also known as the yield strain σ Y )much bigger energy drops become possible. Once theyield stress has been achieved, there is a quantitativechange in the nature of the plastic drops since they be-come system-size dependent. We can examine the statis-tics of the magnitude of the energy and stress drops inthe steady state. In Fig. 2 we show the average magni-tude of energy (cid:104) ∆ U (cid:105) and stress drops (cid:104) ∆ σ (cid:105) for systemsof increasing number of particles N . It appears that thedata support the scaling laws (cid:104) ∆ U (cid:105) ∼ N α , α ≈ . , (9) (cid:104) ∆ σ (cid:105) ∼ N β , β ≈ − . . (10)In Ref. [11] it was shown that these exponents satisfy ascaling relation α − β = 1 as these exponents do. Sincethe system spanning events are confined to linear struc-tures one is not surprised with the exponent α = 1 / β = − / α = 0 . α = 1 /
2. The data −4 −2 ∆ U P ( ∆ U ) N=2000N=2560N=3000N=3600N=5120 −4 −2 ∆ U/N P ( ∆ U / N . ) FIG. 3: Raw and re-scaled pdf’s for the energy drops in thesteady state. The data collapse in the tails of the distributionssupports the scaling laws presented in Eqs. (9) and (10). collapse of the pdfs in the tails shows that the exponentis adequate. Note that the re-scaling does not collapsethe data for small drops, these continue to be system sizeindependent.A theoretical discussion of the localized and the subex-tensive plastic events is provided in Sect. IV. Neverthe-less the reader should note that a continuum descriptionof the stress vs strain curves in our open system is stillunder debate, even in the case of simpler examples like bi-nary Lennard-Jones glasses. Here the quantitative theoryof energy input by mechanical strain, including the sharetaken by stress vs. oligomeric conformation changes onthe one hand, and energy dissipated to the heat bath onthe other hand is still unavailable. Such an understand-ing is prerequisite to any continuum theory.
B. Stress and energy averaged over realizations asa function of the stiffness parameter
In addition to the measuring the plastic drops in singlerealization of the glass, it is interesting to examine theenergy and the stress averaged over many realizations.Such graphs should be closer to what is expected in thethermodynamic limit when N p → ∞ . In particular wecan examine the dependence on the stiffness parameter κ . In Fig. 4 we see the stress vs. strain and the energyvs. strain averaged over 40 independent realizations as afunction of κ . It is interesting to see that both the en-ergy and the stress appear to reach the same steady statefor all κ (cid:54) = 0, but not for κ = 0. To underline the factthat the attainment of the same steady state is not at alltrivial, we show in Fig. 5 the dependence on the strainof the various contributions to the energy coming from γ σ x y γ U κ = 0 κ = 2 κ = 4 κ = 5 κ = 8 κ = 10 κ = 15 FIG. 4: Upper panel: Stress Vs strain for different κ for thestiff polymer ( α = 1). The stress reaches to the same steadystate for the finite κ . Bottom Panel: the variation of totalinternal potential energy with strain. The data averaged over40 independent realizations are shown. the different terms in the Hamiltonian. It is quite evi-dent that the various contributions to the total energy doNOT reach the same steady state, and the result shownin Fig. 4 is the consequence of an interesting and sub-tle cancellation that needs to be explained. Currently wehave no explanation to this observation. To be more con-fident in the correctness of the observation we changedthe parameter η in Eq. 6 and repeated the measurements;the observation remains invariant. C. Changes in geometry of polymers with appliedstrain
In addition to the energy and the stress in the system,the oligomeric glass presents also interesting responses toexternal strains in the resulting geometry of the chains.Of course, the configuration of the oligomers in the glassdepends on the stiffness of the chains. As the chainsbecome stiffer the oligomer chain is easier to bend (sincein our convention the straight chain is the minimal energystate). In order to characterize the configuration of theoligomer chains we compute the end-to-end length R ee of the chain and follow how it changes with the appliedstrain. Fig. 6 shows the variation of R ee for the stiffcase as the applied strain is increased. We see that thetendency is different for small and large value of κ . Forsmall κ the chains start from a coiled state, with R ee being of the order of √ n . Then the action of the straintends to straighten the chains to increase R ee until a κ -dependent steady state. On the other hand for large κ one starts with almost straight chains, such the R ee is of U L J U F E N E γ U A n g l e κ = 0 κ = 2 κ = 4 κ = 5 κ = 8 κ = 10 κ = 15 FIG. 5: Variation of contribution in the total potential energydue to different interaction with strain. Upper panel: U LJ Vsstrain for different κ . The potential energy due to the LJinteraction increases on increase of κ . Middle panel: U FENE
Vs strain for different κ . Bottom panel: U angle Vs strain. Thepotential energy U angle decreases on increase of the stiffnessof the chain. the order of n ; straining now leads to bending, increasingthe energy of the system, reaching again a κ -dependentsteady state. The process described can be seen directlyin snapshots of the system under strain. This is shownfor κ = 0 and κ = 10 in Fig. 7 D. Theoretical remarks on the end-to-end distanceand the average angular distribution
The first observation that needs to be explained is theend-to-end distance in equilibrium (at γ = 0). It turnsout that just using the angular potential is sufficient togive a good estimate of this distance. The reason is thatbecause the oligomers are fairly stiff, the Lennard Jonesterm does not lead to strong short-range particle-particlerepulsion, while the main effect of the FENE term is sim-ply to re-normalize the individual inter-bond distances.Thus using the angular potential we can simply calculatethe average angle of the oligomer chain which given by: (cid:104) cos θ (cid:105) = (cid:82) π dθ cos θ exp (cid:16) − U Angle T (cid:17)(cid:82) π dθ exp (cid:16) − U Angle T (cid:17) . (11) γ R ee κ = 0 κ = 2 κ = 4 κ = 10 κ = 15 FIG. 6: Variation of end to end length R ee of the stiff polymerswith the applied strain γ as a function of κ .FIG. 7: Snapshot of the polymer in the box for κ = 0 (upperpanel) and κ = 10 (lower panel). From left to right (a) γ = 0,(b) γ = 1, and (c) γ = 5. For κ = 0 the polymers are coiledfor zero-strain and they stretch on average upon increasing ofthe strain. The opposite occurs for κ = 10. Here the temperature T should be taken to be of the or-der of the fluid melt from which the glass was quenched.Below we take T = 1. This integral can be performedexactly, and its value is I ( κ/T ) /I ( κ/T ) where I and I are the modified Bessel function of order 1 and 0 re-spectively. In the upper panel of Fig. 8 we compare thetheoretical evaluation of (cid:104) cos θ (cid:105) to its numerically com-puted counterpart, and conclude the comparison is good.Using the average angle we can write the average value (cid:104) R ee (cid:105) [19] as: (cid:104) R ee (cid:105) = n (cid:18) (cid:104) cos θ (cid:105) − (cid:104) cos θ (cid:105) − n (cid:104) cos θ (cid:105) (1 − (cid:104) cos θ (cid:105) n )(1 − (cid:104) cos θ (cid:105) ) (cid:19) (12)Taking the square-root of this expression we plot it inthe middle panel of Fig. 8 and compare it with the nu- Κ X cos Θ \ Κ R ee H L Κ FIG. 8: The average angular measure (cid:104) cos θ (cid:105) (upper panel),the average equilibrium end-to-end length R ee (0) (middlepanel) and the average persistence length ( (cid:96) p ) (lower panel)as a function of κ at zero strain ( γ = 0). The red dots witherror bars are the data obtained from numerical simulationsand black curves are the theoretical estimates obtained usingthe theory discussed in Sec.III D. merically calculated value of R ee at γ = 0, averaged over40 different initial conditions. The agreement is quiteacceptable.Finally, it is advantageous to define ‘persistence length’ (cid:96) p using the relationship (cid:104) cos θ (cid:105) ≡ exp( − /(cid:96) p ) . (13)The resulting (cid:96) p for γ = 0 is shown in the lower panelof Fig.8. We note that the persistence length becomes ofthe order of R ee when the latter is about 10.Returning to Fig. 6 one notes three interesting features:1. For small values of κ the end-to-end distance riseswith increasing strain.2. For large values of κ the end-to-end distance de-creases with increasing strain. 3. In either case the end-to-end distance attains a κ -dependent asymptotic value for large γ which isnevertheless not the fully stretched state.For κ small there is a simple estimate of the asymptoticvalue of R ee which involves a balance between the forcedue to straining which tends to stretch the oligomer andthe entropic force which tends to keep the oligomer coiled.Estimating the force due to stress as σR and the entropicforce as TR ee ( R − R ee ) [20]. Balancing the two expressionswe predict that R ee ( σ ) = R ee ( σ = 0)1 − σR ee ( σ = 0) /T . (14)Indeed, the observed increase in the end-to-end distanceat small values of κ is in accordance with this prediction.Of course for larger values of σ the FENE terms need tobe invoked to cure the apparent divergence in Eq. (14).Once the persistence length is of the order of the initialvalue of R ee we can assume that the oligomers are entirelystretched. Then the effect of the shear strain is opposite,in reducing the end-to-end distance. This stems simplyfrom the fact that any inclined stretched polymer willbend under the action of shear, since its two ends moveat different speeds. The reader can see this phenomenonoccurring in the lower panel of Fig. 7. Thus the effect ofincreasing γ will initially decrease R ee as is observed inFig. 6.In both cases the estimate of the asymptotic value of R ee is not easy, and we leave it for future research. IV. THEORY OF PLASTIC EVENTS
The stability of amorphous solids is determined by theHessian matrix which is made of second derivatives of theHamiltonian with respect to all the degrees of freedom.This matrix is always symmetric and real and thereforediagonalizable. As long as all the eigenvalues are positive,the system is mechanically stable. Plastic instabilities arecharacterized by an eigenvalue going to zero signaling theloss of mechanical stability.
A. Calculation of the Hessian matrix
To calculate the Hessian matrix for the oligomeric glasswe recognize the three contributions to the potential en-ergy φ LJ , χ and ψ . These contributions result in threesub matrices that need to be summed up to yield the fullHessian. We denote the sub matrices as H LJ , H FENE and H Angle : H = H LJ + H FENE + H Angle (15)We begin with H LJ : H LJ ( i, j ; α, β ) = ∂ φ ij LJ ∂x jβ ∂x iα = ∂ φ ij LJ ∂ ( r ij ) ∂r ij ∂x jβ ∂r ij ∂x iα + ∂φ ij LJ ∂r ij ∂ r ij ∂x jβ ∂x iα , (16)and for i = j it is: H LJ ( i, i ; α, β ) = (cid:88) j (cid:54) =i −H LJ ( i, j ; α, β ) . (17)Note that unless otherwise stated latin letters (e.g., i,j, etc.) will be used for the particle’s coordinate andgreek letters (e.g., α, β , etc.) will be used to the denotethe displacement coordinate of the particles. In order tocompute the terms used in the above equation (Eq. 16)explicitly, we take the advantage of the identities: ∂r m(cid:96) ∂x iα = r m(cid:96)α r m(cid:96) ( δ (cid:96)i − δ mi ) , (18)and ∂ r m(cid:96) ∂x jβ ∂x iα = (cid:32) δ αβ r m(cid:96) − r m(cid:96)α r m(cid:96)β ( r m(cid:96) ) (cid:33) (cid:0) δ (cid:96)j − δ mj (cid:1) (cid:0) δ (cid:96)i − δ mi (cid:1) . (19)Now we consider the terms in which the bond of apolymer connecting particles k and (cid:96) contributes to the H FENE . These terms are written as: H FENE ( k, l ; α, β ) = ∂ χ k(cid:96) ∂x kα ∂x kβ ∂ χ k(cid:96) ∂x kα ∂x (cid:96)β ∂ χ k(cid:96) ∂x (cid:96)α ∂x kβ ∂ χ k(cid:96) ∂x (cid:96)α ∂x (cid:96)β , (20)where α and β stand for all coordinates and thus thedimensions of A are 2 d × d . The four entries of matrix A are not necessarily adjacent in H FENE ; depending on the values of k and (cid:96) they are positioned at ( dk, dk ),( dk, d(cid:96) ), ( d(cid:96), dk ) and ( d(cid:96), d(cid:96) ) respectively where d is thedimensionality of the system.Further we consider the terms related to the i th va-lence angle within a polymer chain. This angle is definedby three successive particles with indices k , (cid:96) and m . Thecontribution of these terms to the Hessian H Angle is ex-pressed as: H Angle = ∂ ψ i ∂x kα ∂x kβ ∂ ψ i ∂x kα ∂x (cid:96)β ∂ ψ i ∂x kα ∂x mβ ∂ ψ i ∂x (cid:96)α ∂x kβ ∂ ψ i ∂x (cid:96)α ∂x (cid:96)β ∂ ψ i ∂x (cid:96)α ∂x mβ ∂ ψ i ∂x mα ∂x kβ ∂ ψ i ∂x mα ∂x (cid:96)β ∂ ψ i ∂x mα ∂x mβ , (21)where α and β stand for all coordinates and thus thedimensions of H Angle are 3 d × d . The four entries of H Angle are not necessarily adjacent, depending on thevalues of k , (cid:96) and m they are positioned at ( dk, dk ),( dk, d(cid:96) ), ( dk, dm ), ( d(cid:96), dk ), ( d(cid:96), d(cid:96) ), ( d(cid:96), dm ), ( dm, dk ),( dm, d(cid:96) ) and ( dm, dm ), respectively, where d is dimen-sionality of the system.Formally the angular contribution to the Hessian canbe expressed as: H Angle ( i, j ; α, β ) = ∂ ψ (cid:96) ∂x jβ ∂x iα = ∂ ψ (cid:96) ( ∂ cos ϕ (cid:96) ) ∂ cos ϕ (cid:96) ∂x iα ∂ cos ϕ (cid:96) ∂x jβ + ∂ψ (cid:96) ∂ cos ϕ (cid:96) ∂ cos ϕ (cid:96) ∂x jβ ∂x iα , (22)where the cosine of the valence angle is defined as:cos ϕ l = − r l − ,lγ r l,l +1 γ r l − ,l r l,l +1 , with r klγ = r lγ − r kγ . Now inserting the formula for cosinein the Eq. 22, one obtains: ∂ cos ϕ l ∂x mα = − (cid:20) r l,l +1 γ r l,l +1 ∂∂x mα (cid:18) r l − ,lγ r l − ,l (cid:19) + r l − ,lγ r l − ,l ∂∂x mα (cid:18) r l,l +1 γ r l,l +1 (cid:19)(cid:21) ∂ cos ϕ l ∂x mα ∂x qβ = − r l,l +1 γ r l,l +1 ∂ ∂x mα ∂x qβ (cid:18) r l − ,lγ r l − ,l (cid:19) + r l − ,lγ r l − ,l ∂ ∂x mα ∂x qβ (cid:18) r l,l +1 γ r l,l +1 (cid:19) ++ ∂∂x mα (cid:18) r l − ,lγ r l − ,l (cid:19) ∂∂x qβ (cid:18) r l,l +1 γ r l,l +1 (cid:19) + ∂∂x qβ (cid:18) r l − ,lγ r l − ,l (cid:19) ∂∂x mα (cid:18) r l,l +1 γ r l,l +1 (cid:19) The other auxiliary expressions are given by: ∂∂x mα (cid:18) r klγ r kl (cid:19) = (cid:18) δ αγ r kl − r klα r klγ ( r kl ) (cid:19) (cid:0) δ lm − δ km (cid:1) ∂ ∂x mα ∂x qβ (cid:18) r klγ r kl (cid:19) = (cid:18) r klα r klβ r klγ ( r kl ) − δ αβ r klγ + δ αγ r klβ + δ βγ r klα ( r kl ) (cid:19) (cid:0) δ lm − δ km (cid:1) (cid:0) δ lq − δ kq (cid:1) . Combining the above expressions together we have: ∂ cos ϕ l ∂x mα = − (cid:20) r l,l +1 γ r l,l +1 ∂∂x mα (cid:18) r l − ,lγ r l − ,l (cid:19) + r l − ,lγ r l − ,l ∂∂x mα (cid:18) r l,l +1 γ r l,l +1 (cid:19)(cid:21) == − (cid:20)(cid:0) δ l,m − δ l − ,m (cid:1) (cid:18) r l,l +1 α r l,l +1 r l − ,l − r l − ,lα r l − ,lγ r l,l +1 γ r l,l +1 ( r l − ,l ) (cid:19) + (cid:0) δ l +1 ,m − δ l,m (cid:1) (cid:18) r l − ,lα r l,l +1 r l − ,l − r l,l +1 α r l − ,lγ r l,l +1 γ r l − ,l ( r l,l +1 ) (cid:19)(cid:21) , ∂ cos ϕ l ∂x mα ∂x qβ = − r l,l +1 γ r l,l +1 ∂ ∂x mα ∂x qβ (cid:18) r l − ,lγ r l − ,l (cid:19) + r l − ,lγ r l − ,l ∂ ∂x mα ∂x qβ (cid:18) r l,l +1 γ r l,l +1 (cid:19) ++ ∂∂x mα (cid:18) r l − ,lγ r l − ,l (cid:19) ∂∂x qβ (cid:18) r l,l +1 γ r l,l +1 (cid:19) + ∂∂x qβ (cid:18) r l − ,lγ r l − ,l (cid:19) ∂∂x mα (cid:18) r l,l +1 γ r l,l +1 (cid:19) == − (cid:18) r l − ,lα r l − , β r l − ,lγ r l,l +1 γ r l,l +1 ( r l − ,l ) − [ δ αβ r l − ,lγ r l,l +1 γ + r l,l +1 α r l − , β + r l − , α r l,l +1 β ] r l,l +1 ( r l − ,l ) (cid:19) (cid:0) δ lm − δ l − ,m (cid:1) (cid:0) δ lq − δ l − ,q (cid:1) ++ (cid:18) r l,l +1 α r l, β r l − ,lγ r l,l +1 γ r l − ,l ( r l,l +1 ) − [ δ αβ r l − ,lγ r l,l +1 γ + r l,l +1 α r l − , β + r l − , α r l,l +1 β ] r l − ,l ( r l,l +1 ) (cid:19) (cid:0) δ l +1 ,m − δ l,m (cid:1) (cid:0) δ l +1 ,q − δ l,q (cid:1) ++ (cid:18) δ αβ r l − ,l r l,l +1 − r l − ,lα r l − , β r l,l +1 ( r l − ,l ) − r l,l +1 α r l, β r l − ,l ( r l,l +1 ) + r l − ,lα r l, β r l − ,lγ r l,l +1 γ ( r l − ,l ) ( r l,l +1 ) (cid:19) (cid:0) δ l,m − δ l − ,m (cid:1) (cid:0) δ l +1 ,q − δ l,q (cid:1) ++ (cid:18) δ αβ r l − ,l r l,l +1 − r l − ,lα r l − , β r l,l +1 ( r l − ,l ) − r l,l +1 α r l, β r l − ,l ( r l,l +1 ) + r l,l +1 α r l − , β r l − ,lγ r l,l +1 γ ( r l − ,l ) ( r l,l +1 ) (cid:19) (cid:0) δ l +1 ,m − δ l,m (cid:1) (cid:0) δ l,q − δ l − ,q (cid:1) . B. Elementary plastic events
Having calculated the Hessian matrix we can now ex-amine the elementary plastic events that occur at smallvalues of γ . As said above, the mechanical stability is lostwhen an eigenvalue of the Hessian goes to zero. This isoccurring via a saddle-node bifurcation in which the min-imum in which the system resides collides with a saddleof the global energy surface. During a saddle node bifur-cation the approach of the eigenvalue to zero is generic,following a square-root singularity [11] λ p ∼ (cid:112) γ p − γ , (23)where λ P is the eigenvalue that reaches zero at γ = γ P .An example of this square-root singularity for a stiffoligomeric glass with κ = 2 is shown in Fig. 9. As theinstability is approached the non-affine response becomescloser to the eigenvector of the Hessian matrix that is as-sociated with λ P , denoted as Ψ P . This phenomenon isdemonstrated in Fig. 10. It is important to stress thatthe square-root singularity is generic and characteristicto saddle-node bifurcations. It should be therefore inde-pendent of the system parameters and even the nature ofthe system. In our case we demonstrate this universalityby changing form stiff to semi-flexible and measuring theeigenvalue λ P for two values of κ as shown in Fig. 11. V. FORMATION OF SHEAR BANDS
Oligomeric glasses, like simple binary glasses and themuch more complex metallic glasses, exhibit, in addi-tion to localized plastic events also a second class ofsystem spanning, shear localizing events. These eventsare precursors to shear banding, and they need a finiteamount of stress or strain to accumulate before they be-come possible. In previous analysis it was shown that γ λ P −18 −17 −16 −15 −14−7−6−5 log ( γ P − γ ) l o g ( λ P ) FIG. 9: The variation of the smallest eigenvalue λ P as γ isincreased, In the upper panel we see the eigenvalue dips tozero, then recovers after the instability is over, and again dipsto zero at the next instability. In the lower panel we chooseto blow up the region of the first instability to demonstratethe approach of the eigenvalue to zero with a square-root sin-gularity Eq. (23). shear localizing events occur when the strain exceeds avalue γ Y which depends on the Poisson ratio of the ma-terial but is usually around 5-7% [12]. It appears thatthe present oligomeric glasses are not much different inthis respect. We begin to see shear localizing instabilitieswhen γ is of the order of 10% or less. The shear local-ization event is rather dramatic; even though we shearhomogeneously with our affine transformation the sys-tem chooses to respond by localizing all the shear overa small band of the size of the core of the Eshelby so-lution, and see Refs. [12, 13] for details. It was shownin Refs. [12, 13] that this solution minimizes the energy Eigen vector
Non−affine field
FIG. 10: The eigenvector Ψ P and the non-affine displacementfield associated with the first plastic instability as λ → −18 −17.5 −17 −16.5 −16 −15.5−7−6−5 log( γ P − γ ) l o g ( λ P ) data ( κ =2)0.480980X+2.517667 −18 −17 −16 −15 −14−8−6−4 log( γ P − γ ) l o g ( λ P ) data ( κ =8)0.491688X+2.515884 FIG. 11: Log-log plot of the eigen value of the plastic mode λ P Vs γ P − γ for κ = 2 (top panel) and for κ =8 (bottompanel) near the first elementary plastic event for semi-flexible( α = 2) polymer. The exponent is approximately 0 . compared to a random array of elementary plastic events.The nature of the shear localizing events is similar towhat had been seen previously: an eigenvalue of the Hes-sian matrix dips to zero, but now instead of a singlequadrupolar structure a whole string of those, concate-nated along a line in 2-dimensions [12] or on a planein 3-dimensions [13], appear simultaneously. They havea global connection now, with the outgoing direction ofone quadrupolar structure connecting immediately to theincoming direction of the next quadrupole, thus arrang-ing the displacement field to go in two different directionabove and below the line (or plane). For pure shear theline (or plane) is in 45 o to the principal stress axis. Otherangles are possible for uniaxial loading [21].The best way to demonstrate the phenomenon is todisplay the eigenfunction or the displacement field asso-ciated with the event. In Fig. 12 we show both, theeigenfunction in the upper panel and the directly sim-ulated non-affine displacement field at the instability inthe lower panel. Both images show how the shear is nowconcentrated over a narrow band, with the displacementfield pointing to the “right” above the band and to the“left” below the band. In a stress control rather than astrain controlled experiment such an event would lead tomacroscopic failure.In the next section we will present a theoretical for-malism to compute the shear modulus for the oligomericglasses. VI. SHEAR MODULUS µ The shear modulus that is a measure of linear responseof the material under the applied strain characterizes themechanical behaviour of the system. Here we provide thetheory that relates the shear modulus to the microscopicvariables like Hessian, non-affine displacements, etc.We recall that for homogeneous shear strain the shearmodulus is defined as the second derivatives of the po-tential energy with respect to the applied strain γ , i.e., µ = 1 V d U ( r , r , ...., r N ; γ ) dγ . (24)In this expression the second derivative contains two con-tributions: one coming from the affine part and anotherfrom the non-affine motion of the monomers. Thus wehave [22] ddγ = ∂∂γ + ∂∂ u i · ∂ u i ∂γ ≡ ∂∂ r i · ∂ u i ∂γ , (25)where the second equality follows from the relation: d r i = d u i . Now the expression for shear modulus has the form: µ = ∂ U∂γ + ∂u i ∂γ ∂U∂r i ∂γ . (26)Further we note that the affine step is followed by thenon-affine step that returns the system to the equilibriumstate. For the equilibrium state0 Eigen vector
Non−affine field
FIG. 12: The eigenfunction Ψ P and the non-affine displace-ment field associated with a shear localizing plastic instabilityas λ →
0. Note the global connection between the series ofquadrupoles arranged along the line, such that the displace-ment field is pointing right above and left below the line. ThisIS the phenomenon of shear localization. df i dγ ≡ − ddγ ∂U∂r i = 0 , (27)where f i is the force on the i th particle. As we use theEq. 25 in the above equation (Eq. 27) we obtain d u i dγ = − H − ij · Ξ j , (28)where H ij is the Hessian and Ξ j = ∂ U∂γ∂r j is the non-affine −5 0 5 10 15 2010152025 κ µ numericaltheoretical FIG. 13: Comparison of the theoretically calculated and thenumerically estimated shear modulus for various values of κ .The relatively large error bars stem from the relative small-ness of the system, in which different realizations give a spreadof values of the shear modulus. Nevertheless the agreementbetween theory and simulations is quite satisfactory. force. Now putting back Eq.(28) into Eq.(26) we obtainthe expression for shear modulus as µ = 1 V ∂ U ( r , r , ....r n ; γ ) ∂γ − V (cid:88) i,j Ξ i · H − ij · Ξ j . (29)The first term in the above expression represents con-tribution in the shear modulus as a result of the affinedisplacement (also called as Born term), while the secondone is the contribution due to the non-affine responses.The Born term is computed analytically in Appendix A.The so called “non-affine force” Ξ is calculated directlyfrom the knowledge of the potential, see the Appendix,and then we solve the inverted equation (28) H · d u dγ = Ξ using conjugate gradient minimization. Having at handthe non-affine velocity d u /dγ we can get the non-affinecontribution to the shear modulus using Eq. (26).A comparison between the theoretically calculatedshear modulus and the one estimated directly from thestress vs. strain curves at very small γ is provided in Fig.13. VII. SUMMARY AND CONCLUDINGREMARKS
In this paper we discussed the mechanics of oligomericglasses, also known as waxes, with a special attentionto the stress and energy vs strain, the characteristics ofthe oligomeric chains and their changes under strain, theshear modulus, and the plastic failure modes. We pro-posed a microscopic outlook which extends the availabletheory for simple binary glasses to this much more com-plex oligomeric example. This resulted in an exact the-ory for the shear modulus, and a full understanding of1the plastic failure, both in the localized and the extendedmodes.There are a few open problems that call for further the-oretical and numerical considerations. The most relevantare:1. A continuum theory of the stress vs strain and en-ergy vs strain is lacking. To be realistic, this isa hard task, and even for the simpler case of bi-nary glasses such a theory is still under hard debate[23, 24]. Understanding the energy budget will becrucial in achieving progress along these lines.2. A theory of the conformational changes of theoligomeric chain under strain is missing. We haveprovided above a theory of the end-to-end distancefor the case γ = 0 but not for finite γ .3. The extension of the approach to three dimensionsis highly desirable. There one can expect interest-ing effects of oligomer interpenetration, trappingand reptation, especially with longer oligomers andunder higher strains.At least the last of these open issues is under activestudy in our laboratory, and we hope to present it in thenear future. Appendix A: Analytic computation of Born term
For small strain field the potential energy can be ex-pressed as: U = U + ∂U∂(cid:15) αβ (cid:15) αβ + 12 ∂ U∂(cid:15) αβ ∂(cid:15) ην (cid:15) αβ (cid:15) ην + O ( (cid:15) ) . (A1)Also for the simple shear with affine transformation h we have: h T h = (cid:32) γ (cid:33) · (cid:32) γ (cid:33) = (cid:32) γγ γ (cid:33) = 2 (cid:15) + I . (A2)Thus the strain field (cid:15) can be written as: (cid:15) = (cid:32) γ/ γ/ γ / (cid:33) (A3)The Born contribution to the shear modulus µ B canbe expressed as: µ B ≡ d Udγ = ∂ U∂(cid:15) xy + ∂U∂(cid:15) yy (A4)At this stage we need to compute the first and secondderivative of the strained potential: ∂U∂(cid:15) αβ ; ∂ U∂(cid:15) αβ ∂(cid:15) ην . (A5)Let us recall that for the polymer case the potentialenergy is given by: U = N (cid:88) (cid:104) ij (cid:105) φ ij LJ + N p (cid:88) k =1 n − (cid:88) i =1 χ ik + N p (cid:88) k =1 n − (cid:88) i =2 ψ ik (A6)Regarding the pair-wise interactions we have: ∂U LJ or FENE ∂(cid:15) αβ = ∂φ ij LJ ∂r ij ∂r ij ∂(cid:15) αβ or, ∂χ ij ∂r ij ∂r ij ∂(cid:15) αβ (A7)In order to compute ∂r ij /∂(cid:15) αβ we define the change in r ij using ˆ r ijα = h αβ r ijβ . Therefore,ˆ r ij = (cid:113) (ˆ r ijλ ) (A8)= (cid:113) r ijα h T hr ijβ (A9) ≈ (cid:113) ( r ij ) + 2 (cid:15) αβ r ijα r ijβ (A10)= r ij (cid:115) (cid:15) αβ r ijα r ijβ ( r ij ) (A11) ≈ r ij (cid:32) (cid:15) αβ r ijα r ijβ ( r ij ) −
12 ( (cid:15) αβ r ijα r ijβ ) ( r ij ) + O ( (cid:15) ) (cid:33) (A12)= r ij + (cid:15) αβ r ijα r ijβ r ij −
12 ( (cid:15) αβ r ijα r ijβ ) ( r ij ) + O ( (cid:15) ) . (A13)Considering the coefficient of first order term we get: ∂r ij ∂(cid:15) αβ = r ijα r ijβ r ij , (A14)and from the the second order term we have: ∂ U LJ ∂(cid:15) ην ∂(cid:15) αβ = ∂∂(cid:15) ην (cid:18) ∂φ ij ∂r ij ∂r ij ∂(cid:15) αβ (cid:19) (A15)= ∂ φ ij ∂ ( r ij ) ∂r ij ∂(cid:15) ην ∂r ij ∂(cid:15) αβ + ∂φ ij ∂r ij ∂ r ij ∂(cid:15) ην ∂(cid:15) αβ . (A16)Using the aforementioned definition of the change in r ij we obtain: ∂ r ij ∂(cid:15) ην ∂(cid:15) αβ = 2 ( r ijα ) ( r ijβ ) ( r ij ) . (A17)We now turn for the computation of the contributionin the shear modulus coming from the angular part ofthe potential ψ , which is: ∂ψ (cid:96) ∂(cid:15) αβ = ∂ψ (cid:96) ∂ cos ϕ (cid:96) ∂ cos ϕ (cid:96) ∂(cid:15) αβ , (A18)2 ∂ ψ (cid:96) ∂(cid:15) ην ∂(cid:15) αβ = ∂∂(cid:15) ην (cid:18) ∂ψ (cid:96) ∂ cos ϕ (cid:96) ∂ cos ϕ (cid:96) ∂(cid:15) αβ (cid:19) = ∂ ψ (cid:96) ∂ (cos ϕ (cid:96) ) ∂ cos ϕ (cid:96) ∂(cid:15) ην ∂ cos ϕ (cid:96) ∂(cid:15) αβ + ∂ψ (cid:96) ∂ cos ϕ (cid:96) ∂ cos ϕ (cid:96) ∂(cid:15) ην ∂(cid:15) αβ . (A19)The terms that remained to computed are : ∂ cos ϕ (cid:96) ∂(cid:15) αβ ; ∂ cos ϕ (cid:96) ∂(cid:15) αβ ∂(cid:15) ην (A20) Using the definition of the cosine as:cos ϕ l = − r l − ,lγ r l,l +1 γ r l − ,l r l,l +1 , r klγ = x lγ − x kγ , we have ∂ cos ϕ l ∂(cid:15) αβ = − (cid:20) r l,l +1 γ r l,l +1 ∂∂(cid:15) αβ (cid:18) r l − ,lγ r l − ,l (cid:19) + r l − ,lγ r l − ,l ∂∂(cid:15) αβ (cid:18) r l,l +1 γ r l,l +1 (cid:19)(cid:21) ∂ cos ϕ l ∂(cid:15) ην ∂(cid:15) αβ = − r l,l +1 γ r l,l +1 ∂ ∂(cid:15) ην ∂(cid:15) αβ (cid:18) r l − ,lγ r l − ,l (cid:19) + r l − ,lγ r l − ,l ∂ ∂(cid:15) ην ∂(cid:15) αβ (cid:18) r l,l +1 γ r l,l +1 (cid:19) + ∂∂(cid:15) αβ (cid:18) r l − ,lγ r l − ,l (cid:19) ∂∂(cid:15) ην (cid:18) r l,l +1 γ r l,l +1 (cid:19) + ∂∂(cid:15) ην (cid:18) r l − ,lγ r l − ,l (cid:19) ∂∂(cid:15) αβ (cid:18) r l,l +1 γ r l,l +1 (cid:19) where: ∂∂(cid:15) αβ (cid:18) r k(cid:96)γ r k(cid:96) (cid:19) = r k(cid:96) ∂r k(cid:96)γ ∂(cid:15) αβ − r k(cid:96)γ ( r k(cid:96) ) ∂r k(cid:96) ∂(cid:15) αβ ∂ ∂(cid:15) ην ∂(cid:15) αβ (cid:18) r k(cid:96)γ r k(cid:96) (cid:19) = − r k(cid:96) ) ∂r k(cid:96) ∂(cid:15) ην ∂r k(cid:96)γ ∂(cid:15) αβ + r k(cid:96) ∂ r k(cid:96)γ ∂(cid:15) ην ∂(cid:15) αβ − r k(cid:96) ) ∂r k(cid:96) ∂(cid:15) αβ ∂r k(cid:96)γ ∂(cid:15) ην + r k(cid:96)γ ( r k(cid:96) ) ∂r k(cid:96) ∂(cid:15) ην ∂r k(cid:96) ∂(cid:15) αβ − r k(cid:96)γ ( r k(cid:96) ) ∂ r k(cid:96) ∂(cid:15) ην ∂(cid:15) αβ New derivatives that need to be defined are: ∂r k(cid:96)x ∂(cid:15) αβ = α = x, β = xr k(cid:96)y ; α = x, β = yr k(cid:96)y ; α = y, β = x ≈ α = y, β = y (A21) ∂r k(cid:96)y ∂(cid:15) αβ = 0 (A22) ∂ r k(cid:96)γ ∂(cid:15) ην ∂(cid:15) αβ = 0 (A23)Plugging the latter into cosine derivations we have: ∂ cos ϕ (cid:96) ∂(cid:15) yy = − (cid:34) r (cid:96),(cid:96) +1 γ r (cid:96),(cid:96) +1 ∂∂(cid:15) yy (cid:32) r (cid:96) − ,(cid:96)γ r (cid:96) − ,(cid:96) (cid:33) + r (cid:96) − ,(cid:96)γ r (cid:96) − ,(cid:96) ∂∂(cid:15) yy (cid:32) r (cid:96),(cid:96) +1 γ r (cid:96),(cid:96) +1 (cid:33)(cid:35) = − (cid:34) r (cid:96),(cid:96) +1 x r (cid:96),(cid:96) +1 ∂∂(cid:15) yy (cid:18) r (cid:96) − ,(cid:96)x r (cid:96) − ,(cid:96) (cid:19) + r (cid:96),(cid:96) +1 y r (cid:96),(cid:96) +1 ∂∂(cid:15) yy (cid:32) r (cid:96) − ,(cid:96)y r (cid:96) − ,(cid:96) (cid:33) + r (cid:96) − ,(cid:96)x r (cid:96) − ,(cid:96) ∂∂(cid:15) yy (cid:18) r (cid:96),(cid:96) +1 x r (cid:96),(cid:96) +1 (cid:19) + r (cid:96) − ,(cid:96)y r (cid:96) − ,(cid:96) ∂∂(cid:15) yy (cid:32) r (cid:96),(cid:96) +1 y r (cid:96),(cid:96) +1 (cid:33)(cid:35) − (cid:34) r (cid:96),(cid:96) +1 x r (cid:96),(cid:96) +1 (cid:32) − r (cid:96) − ,(cid:96)x ( r (cid:96) − ,(cid:96)y ) ( r (cid:96) − ,(cid:96) ) (cid:33) + r (cid:96),(cid:96) +1 y r (cid:96),(cid:96) +1 (cid:32) − ( r (cid:96) − ,(cid:96)y ) ( r (cid:96) − ,(cid:96) ) (cid:33) + r (cid:96) − ,(cid:96)x r (cid:96) − ,(cid:96) (cid:32) − r (cid:96),(cid:96) +1 x ( r (cid:96),(cid:96) +1 y ) ( r (cid:96),(cid:96) +1 ) (cid:33) + r (cid:96) − ,(cid:96)y r (cid:96) − ,(cid:96) (cid:32) − ( r (cid:96),(cid:96) +1 y ) ( r (cid:96),(cid:96) +1 ) (cid:33)(cid:35) (A24)3 ∂ cos ϕ (cid:96) ∂(cid:15) xy = − (cid:34) r (cid:96),(cid:96) +1 γ r (cid:96),(cid:96) +1 ∂∂(cid:15) xy (cid:32) r (cid:96) − ,(cid:96)γ r (cid:96) − ,(cid:96) (cid:33) + r (cid:96) − ,(cid:96)γ r (cid:96) − ,(cid:96) ∂∂(cid:15) xy (cid:32) r (cid:96),(cid:96) +1 γ r (cid:96),(cid:96) +1 (cid:33)(cid:35) (A25)= − (cid:34) r (cid:96),(cid:96) +1 x r (cid:96),(cid:96) +1 ∂∂(cid:15) xy (cid:18) r (cid:96) − ,(cid:96)x r (cid:96) − ,(cid:96) (cid:19) + r (cid:96),(cid:96) +1 y r (cid:96),(cid:96) +1 ∂∂(cid:15) xy (cid:32) r (cid:96) − ,(cid:96)y r (cid:96) − ,(cid:96) (cid:33) + r (cid:96) − ,(cid:96)x r (cid:96) − ,(cid:96) ∂∂(cid:15) xy (cid:18) r (cid:96),(cid:96) +1 x r (cid:96),(cid:96) +1 (cid:19) + r (cid:96) − ,(cid:96)y r (cid:96) − ,(cid:96) ∂∂(cid:15) xy (cid:32) r (cid:96),(cid:96) +1 y r (cid:96),(cid:96) +1 (cid:33)(cid:35) (A26)= − (cid:34) r (cid:96),(cid:96) +1 x r (cid:96),(cid:96) +1 (cid:32) r (cid:96) − ,(cid:96)y r (cid:96) − ,(cid:96) − ( r (cid:96) − ,(cid:96)x ) r (cid:96) − ,(cid:96)y ( r (cid:96) − ,(cid:96) ) (cid:33) + r (cid:96),(cid:96) +1 y r (cid:96),(cid:96) +1 (cid:32) − r (cid:96) − ,(cid:96)x ( r (cid:96) − ,(cid:96)y ) ( r (cid:96) − ,(cid:96) ) (cid:33) + r (cid:96) − ,(cid:96)x r (cid:96) − ,(cid:96) (cid:32) r (cid:96),(cid:96) +1 y r (cid:96),(cid:96) +1 − ( r (cid:96),(cid:96) +1 x ) r (cid:96),(cid:96) +1 y ( r (cid:96),(cid:96) +1 ) (cid:33) + r (cid:96) − ,(cid:96)y r (cid:96) − ,(cid:96) (cid:32) − r (cid:96),(cid:96) +1 x ( r (cid:96),(cid:96) +1 y ) ( r (cid:96),(cid:96) +1 ) (cid:33)(cid:35) (A27) ∂ cos ϕ (cid:96) ∂(cid:15) xy ∂(cid:15) xy = − (cid:34) r (cid:96),(cid:96) +1 γ r (cid:96),(cid:96) +1 ∂ ∂(cid:15) xy ∂(cid:15) xy (cid:32) r (cid:96) − ,(cid:96)γ r (cid:96) − ,(cid:96) (cid:33) + r (cid:96) − ,(cid:96)γ r (cid:96) − ,(cid:96) ∂ ∂(cid:15) xy ∂(cid:15) xy (cid:32) r (cid:96),(cid:96) +1 γ r (cid:96),(cid:96) +1 (cid:33) +2 ∂∂(cid:15) xy (cid:32) r (cid:96) − ,(cid:96)γ r (cid:96) − ,(cid:96) (cid:33) ∂∂(cid:15) xy (cid:32) r (cid:96),(cid:96) +1 γ r (cid:96),(cid:96) +1 (cid:33)(cid:35) (A28)= − (cid:34) r (cid:96),(cid:96) +1 x r (cid:96),(cid:96) +1 ∂ ∂(cid:15) xy ∂(cid:15) xy (cid:18) r (cid:96) − ,(cid:96)x r (cid:96) − ,(cid:96) (cid:19) + r (cid:96),(cid:96) +1 y r (cid:96),(cid:96) +1 ∂ ∂(cid:15) xy ∂(cid:15) xy (cid:32) r (cid:96) − ,(cid:96)y r (cid:96) − ,(cid:96) (cid:33) + r (cid:96) − ,(cid:96)x r (cid:96) − ,(cid:96) ∂ ∂(cid:15) xy ∂(cid:15) xy (cid:18) r (cid:96),(cid:96) +1 x r (cid:96),(cid:96) +1 (cid:19) + r (cid:96) − ,(cid:96)y r (cid:96) − ,(cid:96) ∂ ∂(cid:15) xy ∂(cid:15) xy (cid:32) r (cid:96),(cid:96) +1 y r (cid:96),(cid:96) +1 (cid:33) +2 ∂∂(cid:15) xy (cid:18) r (cid:96) − ,(cid:96)x r (cid:96) − ,(cid:96) (cid:19) ∂∂(cid:15) xy (cid:18) r (cid:96),(cid:96) +1 x r (cid:96),(cid:96) +1 (cid:19) + 2 ∂∂(cid:15) xy (cid:32) r (cid:96) − ,(cid:96)y r (cid:96) − ,(cid:96) (cid:33) ∂∂(cid:15) xy (cid:32) r (cid:96),(cid:96) +1 y r (cid:96),(cid:96) +1 (cid:33)(cid:35) (A29)= − (cid:34) r (cid:96),(cid:96) +1 x r (cid:96),(cid:96) +1 (cid:18) − r (cid:96) − ,(cid:96) ) ∂r (cid:96) − ,(cid:96) ∂(cid:15) xy ∂r (cid:96) − ,(cid:96)x ∂(cid:15) xy + 2 r (cid:96) − ,(cid:96)x ( r (cid:96) − ,(cid:96) ) ∂r (cid:96) − ,(cid:96) ∂(cid:15) xy ∂r (cid:96) − ,(cid:96) ∂(cid:15) xy (cid:19) + r (cid:96),(cid:96) +1 y r (cid:96),(cid:96) +1 (cid:32) r (cid:96) − ,(cid:96)y ( r (cid:96) − ,(cid:96) ) ∂r (cid:96) − ,(cid:96) ∂(cid:15) xy ∂r (cid:96) − ,(cid:96) ∂(cid:15) xy (cid:33) + r (cid:96) − ,(cid:96)x r (cid:96) − ,(cid:96) (cid:18) − r (cid:96),(cid:96) +1 ) ∂r (cid:96),(cid:96) +1 ∂(cid:15) xy ∂r (cid:96),(cid:96) +1 x ∂(cid:15) xy + 2 r (cid:96),(cid:96) +1 x ( r (cid:96),(cid:96) +1 ) ∂r (cid:96),(cid:96) +1 ∂(cid:15) xy ∂r (cid:96),(cid:96) +1 ∂(cid:15) xy (cid:19) + r (cid:96) − ,(cid:96)y r (cid:96) − ,(cid:96) (cid:32) r (cid:96),(cid:96) +1 y ( r (cid:96),(cid:96) +1 ) ∂r (cid:96),(cid:96) +1 ∂(cid:15) xy ∂r (cid:96),(cid:96) +1 ∂(cid:15) xy (cid:33) +2 (cid:18) r (cid:96) − ,(cid:96) ∂r (cid:96) − ,(cid:96)x ∂(cid:15) xy − r (cid:96) − ,(cid:96)x ( r (cid:96) − ,(cid:96) ) ∂r (cid:96) − ,(cid:96) ∂(cid:15) xy (cid:19) (cid:18) r (cid:96),(cid:96) +1 ∂r (cid:96),(cid:96) +1 x ∂(cid:15) xy − r (cid:96),(cid:96) +1 x ( r (cid:96),(cid:96) +1 ) ∂r (cid:96),(cid:96) +1 ∂(cid:15) xy (cid:19) +2 (cid:32) r (cid:96) − ,(cid:96) ∂r (cid:96) − ,(cid:96)y ∂(cid:15) xy − r (cid:96) − ,(cid:96)y ( r (cid:96) − ,(cid:96) ) ∂r (cid:96) − ,(cid:96) ∂(cid:15) xy (cid:33) (cid:32) r (cid:96),(cid:96) +1 ∂r (cid:96),(cid:96) +1 y ∂(cid:15) xy − r (cid:96),(cid:96) +1 y ( r (cid:96),(cid:96) +1 ) ∂r (cid:96),(cid:96) +1 ∂(cid:15) xy (cid:33)(cid:35) (A30)= − (cid:34) r (cid:96),(cid:96) +1 x r (cid:96),(cid:96) +1 (cid:32) − r (cid:96) − ,(cid:96)x ( r (cid:96) − ,(cid:96)y ) ( r (cid:96) − ,(cid:96) ) + 2( r (cid:96) − ,(cid:96)x ) ( r (cid:96) − ,(cid:96)y ) ( r (cid:96) − ,(cid:96) ) (cid:33) + r (cid:96),(cid:96) +1 y r (cid:96),(cid:96) +1 (cid:32) r (cid:96) − ,(cid:96)x ) ( r (cid:96) − ,(cid:96)y ) ( r (cid:96) − ,(cid:96) ) (cid:33) + r (cid:96) − ,(cid:96)x r (cid:96) − ,(cid:96) (cid:32) − r (cid:96),(cid:96) +1 x ( r (cid:96),(cid:96) +1 y ) ( r (cid:96),(cid:96) +1 ) + 2( r (cid:96),(cid:96) +1 x ) ( r (cid:96),(cid:96) +1 y ) ( r (cid:96),(cid:96) +1 ) (cid:33) + r (cid:96) − ,(cid:96)y r (cid:96) − ,(cid:96) (cid:32) r (cid:96),(cid:96) +1 x ) ( r (cid:96),(cid:96) +1 y ) ( r (cid:96),(cid:96) +1 ) (cid:33) +2 (cid:32) r (cid:96) − ,(cid:96)y r (cid:96) − ,(cid:96) − ( r (cid:96) − ,(cid:96)x ) r (cid:96) − ,(cid:96)y ( r (cid:96) − ,(cid:96) ) (cid:33) (cid:32) r (cid:96),(cid:96) +1 y r (cid:96),(cid:96) +1 − ( r (cid:96),(cid:96) +1 x ) r (cid:96),(cid:96) +1 y ( r (cid:96),(cid:96) +1 ) (cid:33) +2 (cid:32) − r (cid:96) − ,(cid:96)x ( r (cid:96) − ,(cid:96)y ) ( r (cid:96) − ,(cid:96) ) (cid:33) (cid:32) − r (cid:96),(cid:96) +1 x ( r (cid:96),(cid:96) +1 y ) ( r (cid:96),(cid:96) +1 ) (cid:33)(cid:35) (A31)4For non-affine forces= 0: ∂ cos ϕ (cid:96) ∂x mν ∂(cid:15) xy = − ∂∂(cid:15) xy (cid:34)(cid:0) δ (cid:96),m − δ (cid:96) − ,m (cid:1) (cid:32) r (cid:96),l +1 ν r (cid:96),(cid:96) +1 r (cid:96) − ,(cid:96) − r (cid:96) − ,(cid:96)ν r (cid:96) − ,(cid:96)γ r (cid:96),(cid:96) +1 γ r (cid:96),(cid:96) +1 ( r (cid:96) − ,(cid:96) ) (cid:33) + (cid:0) δ (cid:96) +1 ,m − δ (cid:96),m (cid:1) (cid:32) r (cid:96) − ,(cid:96)ν r (cid:96),(cid:96) +1 r (cid:96) − ,(cid:96) − r (cid:96),(cid:96) +1 ν r (cid:96) − ,(cid:96)γ r (cid:96),(cid:96) +1 γ r (cid:96) − ,(cid:96) ( r (cid:96),(cid:96) +1 ) (cid:33)(cid:35) (A32)= (cid:0) δ (cid:96) − ,m − δ (cid:96),m (cid:1) (cid:32) ∂∂(cid:15) xy r (cid:96),l +1 ν r (cid:96),(cid:96) +1 r (cid:96) − ,(cid:96) − ∂∂(cid:15) xy r (cid:96) − ,(cid:96)ν r (cid:96) − ,(cid:96)γ r (cid:96),(cid:96) +1 γ r (cid:96),(cid:96) +1 ( r (cid:96) − ,(cid:96) ) (cid:33) + (cid:0) δ (cid:96),m − δ (cid:96) +1 ,m (cid:1) (cid:32) ∂∂(cid:15) xy r (cid:96) − ,(cid:96)ν r (cid:96),(cid:96) +1 r (cid:96) − ,(cid:96) − ∂∂(cid:15) xy r (cid:96),(cid:96) +1 ν r (cid:96) − ,(cid:96)γ r (cid:96),(cid:96) +1 γ r (cid:96) − ,(cid:96) ( r (cid:96),(cid:96) +1 ) (cid:33)(cid:35) (A33)= (cid:0) δ (cid:96) − ,m − δ (cid:96),m (cid:1) (cid:18) r (cid:96) − ,(cid:96) ∂∂(cid:15) xy r (cid:96),l +1 ν r (cid:96),(cid:96) +1 + r (cid:96),l +1 ν r (cid:96),(cid:96) +1 ∂∂(cid:15) xy r (cid:96) − ,(cid:96) − r (cid:96) − ,(cid:96)ν r (cid:96) − ,(cid:96)γ ( r (cid:96) − ,(cid:96) ) ∂∂(cid:15) xy r (cid:96),(cid:96) +1 γ r (cid:96),(cid:96) +1 − r (cid:96),(cid:96) +1 γ r (cid:96),(cid:96) +1 ∂∂(cid:15) xy r (cid:96) − ,(cid:96)ν r (cid:96) − ,(cid:96)γ ( r (cid:96) − ,(cid:96) ) (cid:33) + (cid:0) δ (cid:96),m − δ (cid:96) +1 ,m (cid:1) (cid:18) r (cid:96),(cid:96) +1 ∂∂(cid:15) xy r (cid:96) − ,(cid:96)ν r (cid:96) − ,(cid:96) + r (cid:96) − ,(cid:96)ν r (cid:96) − ,(cid:96) ∂∂(cid:15) xy r (cid:96),(cid:96) +1 − r (cid:96),(cid:96) +1 ν r (cid:96),(cid:96) +1 γ ( r (cid:96),(cid:96) +1 ) ∂∂(cid:15) xy r (cid:96) − ,(cid:96)γ r (cid:96) − ,(cid:96) − r (cid:96) − ,(cid:96)γ r (cid:96) − ,(cid:96) ∂∂(cid:15) xy r (cid:96),(cid:96) +1 ν r (cid:96),(cid:96) +1 γ ( r (cid:96),(cid:96) +1 ) (cid:33)(cid:35) (A34)Using previously defined expressions and: ∂∂(cid:15) αβ r k(cid:96)ν r k(cid:96)γ ( r k(cid:96) ) = r k(cid:96)γ ( r k(cid:96) ) ∂r k(cid:96)ν ∂(cid:15) αβ + r k(cid:96)ν ( r k(cid:96) ) ∂r k(cid:96)γ ∂(cid:15) αβ − r k(cid:96)ν r k(cid:96)γ ( r k(cid:96) ) ∂r k(cid:96) ∂(cid:15) αβ (A35)we have the final expression of the second derivative of cosine as: ∂ cos ϕ (cid:96) ∂x mν ∂(cid:15) xy = (cid:0) δ (cid:96) − ,m − δ (cid:96),m (cid:1) (cid:34) r (cid:96) − ,(cid:96) (cid:32) δ νx r (cid:96),(cid:96) +1 y r (cid:96),(cid:96) +1 − r (cid:96),(cid:96) +1 ν r (cid:96),(cid:96) +1 x r (cid:96),(cid:96) +1 y ( r (cid:96),(cid:96) +1 ) (cid:33) − r (cid:96),l +1 ν r (cid:96),(cid:96) +1 (cid:32) r (cid:96) − ,(cid:96)x r (cid:96) − ,(cid:96)y ( r (cid:96) − ,(cid:96) ) (cid:33) − r (cid:96) − ,(cid:96)ν r (cid:96) − ,(cid:96)γ ( r (cid:96) − ,(cid:96) ) (cid:32) r (cid:96),l +1 y r (cid:96),(cid:96) +1 − r (cid:96),l +1 γ r (cid:96),l +1 x r (cid:96),l +1 y ( r (cid:96),(cid:96) +1 ) (cid:33) − r (cid:96),(cid:96) +1 γ r (cid:96),(cid:96) +1 (cid:32) r (cid:96) − ,(cid:96)y ( δ νx r (cid:96) − ,(cid:96)γ + r (cid:96) − ,(cid:96)ν )( r (cid:96) − ,(cid:96) ) − r (cid:96) − ,(cid:96)γ r (cid:96) − ,(cid:96)ν r (cid:96) − ,(cid:96)x r (cid:96) − ,(cid:96)y ( r (cid:96) − ,(cid:96) ) (cid:33)(cid:35) + (cid:0) δ (cid:96),m − δ (cid:96) +1 ,m (cid:1) (cid:34) r (cid:96),(cid:96) +1 (cid:32) δ νx r (cid:96) − ,(cid:96)y r (cid:96) − ,(cid:96) − r (cid:96) − ,(cid:96)ν r (cid:96) − ,(cid:96)x r (cid:96) − ,(cid:96)y ( r (cid:96) − ,(cid:96) ) (cid:33) + r (cid:96) − ,(cid:96)ν r (cid:96) − ,(cid:96) (cid:32) r (cid:96),(cid:96) +1 x r (cid:96),(cid:96) +1 y ( r (cid:96),(cid:96) +1 ) (cid:33) − r (cid:96),(cid:96) +1 ν r (cid:96),(cid:96) +1 γ ( r (cid:96),(cid:96) +1 ) (cid:32) r (cid:96) − ,(cid:96)y r (cid:96) − ,(cid:96) − r (cid:96) − ,(cid:96)ν r (cid:96) − ,(cid:96)x r (cid:96) − ,(cid:96)y ( r (cid:96) − ,(cid:96) ) (cid:33) − r (cid:96) − ,(cid:96)γ r (cid:96) − ,(cid:96) (cid:32) r (cid:96),(cid:96) +1 y ( δ νx r (cid:96),(cid:96) +1 γ + r (cid:96),(cid:96) +1 ν )( r (cid:96),(cid:96) +1 ) − r (cid:96),(cid:96) +1 γ r (cid:96),(cid:96) +1 ν r (cid:96),(cid:96) +1 x r (cid:96),(cid:96) +1 y ( r (cid:96),(cid:96) +1 ) (cid:33)(cid:35) (A36) [1] M. Rubinstein and R. H. Colby, Polymer Physics, OxfordUniversity Press (2003).[2] A. M. Donald and E. J. Kramer, J. Mat. Science , 1765(1982). [3] P. H. Mott, A. S. Argon, and U. W. Suter, Phil. MagazineA, , 931 (1993).[4] M. Warren and J. Rottler, Phys. Rev. E , 031802(2007). [5] F. Varnik, L. Bocquet, and J.-L. Barrat, J. Chem. Phys. , 6 (2004).[6] A. Makke, M. Perez, J. Rottler, O. Lame, and J. L. Bar-rat, Macromol. Theory Simul. , 826 (2011).[7] R. S. Hoy and M. O. Robbins, Phys. Rev. Lett. ,117801 (2007);R. S. Hoy and M. O. Robbins, Phys. Rev.E , 031801 (2008).[8] J. L. Barrat, J. Baschnagel and A. Lyulin, Soft Matter , 3430 (2010).[9] D. L. Malandro and D. J. Lacks, J. Chem. Phys. , 353 (1999).[10] C. E. Maloney and A. Lemaˆıtre, Phys. Rev. E , 016118(2006).[11] E. Lerner, and I. Procaccia, Phys. Rev. E , 066109(2009).[12] R. Dasgupta, H. G. Hentschel and I. Procaccia, Phys.Rev. Lett. , 255502 (2012); R. Dasgupta, H. G. E.Hentschel and I. Procaccia, Phys. Rev. E , 115 (2003).[18] K. Kramer and G. S. Grest, J. Chem. Phys. (8), 5057(1990).[19] R. Auhl, R. Everaers, G. S. Grest, K. Kremer, and S. J.Plimpton, J. Chem. Phys. , 12718 (2003).[20] P. G. De Gennes, J. Chem. Phys. , 5030 (1974).[21] Ashwin J., O. Gendelman, I. Procaccia, Carmel Shor,Phys. Rev. E. , 022310 (2013).[22] S. Karmakar, E. Lerner, and I. Procaccia, Phys. Rev. E , 026105 (2010). For a fuller detailed exposition seearXiv:1004.2198.[23] M.L. Falk and J.S. Langer, Annu. Rev. Condens. Matt.Phys. , 353 (2011).[24] S. M. Fielding, P. Sollich and M. E. Cates, J. Rheol.44