Fluctuation-dissipation relation between shear stress relaxation modulus and shear stress autocorrelation function revisited
FFluctuation-dissipation relation between shear stress relaxationmodulus and shear stress autocorrelation function revisited
J.P. Wittmer, H. Xu, ∗ O. Benzerara, and J. Baschnagel Institut Charles Sadron, Universit´e de Strasbourg & CNRS,23 rue du Loess, 67034 Strasbourg Cedex, France LCP-A2MC, Institut Jean Barriol, Universit´e de Lorraine & CNRS, 1 bd Arago, 57078 Metz Cedex 03, France (Dated: August 18, 2015)The shear stress relaxation modulus G ( t ) may be determined from the shear stress ˆ τ ( t ) afterswitching on a tiny step strain γ or by inverse Fourier transformation of the storage modulus G (cid:48) ( ω )or the loss modulus G (cid:48)(cid:48) ( ω ) obtained in a standard oscillatory shear experiment at angular frequency ω . It is widely assumed that G ( t ) is equivalent in general to the equilibrium stress autocorrelationfunction C ( t ) = βV (cid:104) δ ˆ τ ( t ) δ ˆ τ (0) (cid:105) which may be readily computed in computer simulations ( β beingthe inverse temperature and V the volume). Focusing on isotropic solids formed by permanent springnetworks we show theoretically by means of the fluctuation-dissipation theorem and computationallyby molecular dynamics simulation that in general G ( t ) = G eq + C ( t ) for t > G eq being thestatic equilibrium shear modulus. A similar relation holds for G (cid:48) ( ω ). G ( t ) and C ( t ) must thusbecome different for a solid body and it is impossible to obtain G eq directly from C ( t ). I. INTRODUCTION
Background.
A central rheological property charac-terizing the linear response of liquids and solid elasticbodies is the shear relaxation modulus G ( t ) [1–3]. As-suming for simplicity an isotropic system, the shear re-laxation modulus G ( t ) = δτ ( t ) /γ may be obtained fromthe stress increment δτ ( t ) = (cid:104) ˆ τ ( t ) − ˆ τ (0 − ) (cid:105) for t > | γ | (cid:28) t = 0 as sketched in panel (a) of Fig. 1. The instan-taneous shear stress ˆ τ ( t ) at time t may be determined ex-perimentally by probing the forces acting on the walls ofthe shear cell and in a numerical study, as shown in panel(b) for a sheared periodic simulation box, from the im-posed model Hamiltonian and the particle positions andmomenta [4–7]. It is well known that the componentsof the Fourier transformed relaxation modulus G ( t ), thestorage modulus G (cid:48) ( ω ) and the loss modulus G (cid:48)(cid:48) ( ω ), aredirectly measurable in an oscillatory shear strain exper-iment [1] as shown in panel (d) and panel (e). Usingeither G ( t ) or G (cid:48) ( ω ) one obtains in the static limit (boldhorizontal solid lines) the shear modulus [1, 8] G eq = lim t →∞ G ( t ) = lim ω → G (cid:48) ( ω ) . (1)The shear modulus is an important order parameter char-acterizing the transition from the liquid/sol ( G eq = 0) tothe solid/gel state ( G eq >
0) where the particle permuta-tion symmetry of the liquid state is lost on the time win-dow probed [2, 9]. Examples of current interest for thedetermination of G eq include crystalline solids [10], glassforming liquids and amorphous solids [11–26], colloidalgels [27], polymeric networks [1, 28–30], hyperbranchedpolymer chains with sticky end-groups [31] or bridgedequilibrium networks of telechelic polymers [32]. ∗ Electronic address: [email protected]
G’( ) ω σ F µ A σ F L affine response A t G eq G(t)
G(t) non − a ff i n e µ eq G t=0 =0 γ y γ τ F C(t) C(t) C(t)=0 σ >0 switched on γ (b)xL wall (c)(a)(d) (e) ω G eq C’( ) ω t t > 0 τ ^ γ FIG. 1: Sketch of notations and addressed problem: (a) experimental setup, (b) plain shear with periodic bound-ary conditions, (c) sketch of key equation, Eq. (2), with σ F = βV (cid:104) δ ˆ τ (cid:105) characterizing the stress fluctuations at γ = 0, (d) imposed sinusoidal shear strain γ ( t ) (solid line) and mea-sured shear stress ˆ τ ( t ) (dashed line), (e) comparison of C (cid:48) ( ω )(dash-dotted line) and storage modulus G (cid:48) ( ω ), Eq. (3). Key issue.
Surprisingly, it is widely assumed [11, 20,25, 28, 29] that in the linear response limit ( γ → G ( t ) must become equiv-alent to the stress autocorrelation function [4] C ( t ) ≡ βV (cid:104) δ ˆ τ ( t ) δ ˆ τ (0) (cid:105) computed in the NV γ T-ensemble at im-posed particle number N , volume V , shear strain γ and a r X i v : . [ c ond - m a t . s t a t - m ec h ] A ug temperature T ( β = 1 /k B T denoting the inverse temper-ature) [33]. Since G ( t ) = C ( t ) is assumed to hold, G eq is supposed to be measurable from some transient “finitefrozen-in amplitude” of C ( t ) [20]. We call this belief “sur-prising” since, obviously, C ( t ) → t (cid:29) θ with θ characterizingthe typical relaxation time of a shear stress fluctuation(properly defined in Sec. V B). At variance to this beliefwe shall show by inspection of the fluctuation-dissipationtheorem (FDT) [34] that more generally G ( t ) − G eq = C ( t ) for t ≥ G ( t ) = 0 for t < G ( t ) only becomes equivalent to C ( t ) for t > G eq = 0 and2. that the shear modulus G eq is only probed by G ( t )on time scales t (cid:29) θ where C ( t ) must vanish.In principle, it is thus impossible to obtain the static shearmodulus G eq of an elastic body only from C ( t ). We shallshow that a similar relation G (cid:48) ( ω ) − G eq = C (cid:48) ( ω ) ≡ (cid:90) ∞ C ( t ) sin( ωt )d( ωt ) (3) G (cid:48)(cid:48) ( ω ) = C (cid:48)(cid:48) ( ω ) ≡ (cid:90) ∞ C ( t ) cos( ωt )d( ωt ) (4)holds in the angular frequency domain as sketched inpanel (e) of Fig. 1 for G (cid:48) ( ω ). We refer below to Eqs. (2-4) as the “key relations”. Outline.
The presented work is closely related to therecent paper [26] where we focus on the difference ofstatic fluctuations and autocorrelation functions in con-jugated ensembles and where we also discuss briefly tran-sient self-assembled networks. The present paper pro-vides complementary informations focusing on perma-nent elastic networks in the NV γ T-ensemble and on theresponse to an oscillatory shear strain γ ( t ) = γ sin( ωt ).We begin by reminding in Sec. II A the “affine” and“stress fluctuation” contributions µ A and σ F to the equi-librium shear modulus G eq = µ A − σ F and demonstratethen the key relations theoretically. In Sec. III we defineour two-dimensional spring model and give some algo-rithmic details. The construction of the network andsome properties of its athermal ground state are pre-sented in Sec. IV. Our computational results for one finitetemperature are given in Sec. V. Some static propertiesare summarized in Sec. V A. We illustrate numerically inSec. V B that Eq. (2) holds. We focus in Sec. V C onthe storage and loss moduli G (cid:48) ( ω ) and G (cid:48)(cid:48) ( ω ) computeddirectly by applying a sinusoidal shear strain and com-pare this result to the Fourier transformations of G ( t )and C ( t ). We summarize our work in Sec. VI where webriefly comment on the generalization of the key relationsto linear response functions of other intensive variables. II. THEORETICAL CONSIDERATIONSA. Reminder of some static properties
Affine canonical displacements.
Let us consider firstan infinitessimal pure shear strain γ of the periodic sim-ulation box in the ( x, y )-plane as sketched in panel (b)of Fig. 1. We assume that not only the box shape ischanged but that the particle positions r (using the prin-cipal box convention) follow the macroscopic constraintin an affine manner according to r x → r x + γr y for | γ | (cid:28) . (5)Albeit not strictly necessary for the demonstration of thekey relations we focus on, we shall assume, moreover,that this shear transformation is also canonical [5, 35].This implies that the momenta must transform as p x → p x − γp y for | γ | (cid:28) . (6)All other coordinates of the positions and momenta re-main unchanged in the presented case [36]. We emphasizethe negative sign in Eq. (6) which assures that Liouville’stheorem is obeyed [35]. Instantaneous shear stress and affine shear elasticity.
Let ˆ H ( γ ) denote the system Hamiltonian of a given state s written as a function of the shear strain γ . The instan-taneous shear stress ˆ τ and the instantaneous affine shearelasticity ˆ µ A may be defined as the expansion coefficientsassociated to the energy change δ ˆ H ( γ ) /V = ˆ τ γ + ˆ µ A γ / | γ | (cid:28) γ = 0 being the reference, i.e. [37]ˆ τ ≡ ˆ H (cid:48) ( γ ) /V | γ =0 and (8)ˆ µ A ≡ ˆ H (cid:48)(cid:48) ( γ ) /V | γ =0 . (9)With ˆ H id ( γ ) and ˆ H ex ( γ ) being the standard ideal kineticand the (conservative) excess interaction contributions tothe total Hamiltonian ˆ H ( γ ) = ˆ H id ( γ ) + ˆ H ex ( γ ), this im-plies similar relations for the corresponding contributionsˆ τ id and ˆ τ ex to ˆ τ = ˆ τ id + ˆ τ ex and for the contributions ˆ µ A , id and ˆ µ A , ex to ˆ µ A = ˆ µ A , id + ˆ µ A , ex . The ideal contributionsˆ τ id and ˆ µ A , id are then due to the change of ˆ H id ( γ ) im-posed by the momentum transform Eq. (6), the excesscontributions ˆ τ ex and ˆ µ A , ex due the change of the ˆ H ex ( γ )imposed by the strained particle positions, Eq. (5). Ideal contributions.
Using Eq. (6) and p i = m i v i with m i being the mass of particle i and v i its velocity, thekinetic energy of the strained system becomes ˆ H id ( γ ) = (cid:80) i m i [( v i,x − γv i,y ) + v i,y . . . ] /
2. This implies thatˆ τ id = − V N (cid:88) i =1 m i v i,x v i,y and (10)ˆ µ A , id = 1 V N (cid:88) i =1 m i v i,y (11)for the ideal contributions to the shear stress and theaffine shear elastiticy. Note that the minus sign for theshear stress is due to the minus sign in Eq. (6) requiredfor a canonical transformation. Excess contributions.
We focus below on pairwise ad-ditive excess energies ˆ H ex = (cid:80) l u ( r l ) with u ( r ) being apair potential and where the running index l labels theinteraction between two particles i < j . (The poten-tial u ( r ) may in addition explicitly depend on the in-teraction l .) Due to Eq. (5) a squared particle distance r = x + y + . . . becomes r ( γ ) = ( x + γy ) + y + . . . Straightforward application of the chain rule [21] showsfor the excess contributions to ˆ τ and µ A thatˆ τ ex = 1 V (cid:88) l r l u (cid:48) ( r l ) n l,x n l,y and (12)ˆ µ A , ex = 1 V (cid:88) l (cid:0) r l u (cid:48)(cid:48) ( r l ) − r l u (cid:48) ( r l ) (cid:1) n l,x n l,y + 1 V (cid:88) l r l u (cid:48) ( r l ) n l,y (13)with n l = r l /r l being the normalized distance vector r = r j − r i between the particles i and j . Interest-ingly, Eq. (12) is strictly identical to the correspondingoff-diagonal term of the Kirkwood stress tensor [4]. Thermodynamic averages.
Let us assume an isotropicelastic body at imposed particle number N , constant vol-ume V , shear strain γ and mean temperature T (NV γ T-ensemble). Note that the (intensive) shear strain γ corresponds thermodynamically to an extensive variable X = γV . We write F ( γ ) for the free energy and Z ( γ ) = exp( − βF ( γ )) = (cid:88) s exp( − β ˆ H ( γ )) (14)for the corresponding partition function with (cid:80) s beingthe sum over all accessible system states. Thermal av-erages are given by (cid:104) . . . (cid:105) = (cid:80) s . . . exp( − β ˆ H ( γ )) /Z ( γ ).Interestingly, the definition of the instantaneous shearstress ˆ τ given above, Eq. (7), is consistent with the ther-modynamic mean shear stress τ ≡ (cid:104) ˆ τ (cid:105) , while the averageaffine shear elasticity µ A ≡ (cid:104) ˆ µ A (cid:105) only corresponds to an upper bound to the shear modulus G eq . To see this let usfirst note for convenience that [37] ∂ log( Z ( γ )) ∂γ = Z (cid:48) ( γ ) Z ( γ ) (15) ∂ log( Z ( γ )) ∂γ = Z (cid:48)(cid:48) ( γ ) Z ( γ ) − (cid:18) Z (cid:48) ( γ ) Z ( γ ) (cid:19) (16)and Z (cid:48) ( γ ) = − (cid:88) s β ˆ H (cid:48) ( γ ) e − β ˆ H ( γ ) (17) Z (cid:48)(cid:48) ( γ ) = (cid:88) s (cid:16) ( β ˆ H (cid:48) ( γ )) − β ˆ H (cid:48)(cid:48) ( γ ) (cid:17) e − β ˆ H ( γ ) . (18) It then follows using Eq. (15) and Eq. (17) that the meanshear stress is indeed [21] τ = ∂F ( X ) ∂X (cid:12)(cid:12)(cid:12)(cid:12) X =0 (19)= (cid:88) s ˆ H (cid:48) ( γ ) /V | γ =0 exp( − β ˆ H ) Z = (cid:104) ˆ τ (cid:105) . (20)Using Eq. (16) and Eq. (18) one verifies also that [38] G eq = V ∂ F ( X ) ∂X (cid:12)(cid:12)(cid:12)(cid:12) X =0 (21)= µ A − σ F with σ F ≡ βV (cid:10) δ ˆ τ (cid:11) (22)characterizing the shear-stress fluctuations at γ = 0.Since σ F ≥ µ A is only an upper bound to G eq .We emphasize that Eq. (22) is a special case of thegeneral stress-fluctuation relations for elastic moduli[12, 24, 39, 40]. It provides a computational convenientmethod to obtain G eq for systems in the NV γ T-ensembleused in many recent numerical studies [12–14, 21, 25, 26].We have demonstrated Eq. (22) without introducing a lo-cal displacement field as in Ref. [39]. We note en passant that the averages τ and µ A are “simple averages”, i.e.no fluctuations, and are thus identical for any ensemblegiven that the same state point is sampled [4, 21]. Lebowitz-Percus-Verlet transform.
Equation (22) canalternatively be obtained from the general transforma-tion relation for a fluctuation (cid:104) δ ˆ Aδ ˆ B (cid:105) of two observables A and B due to Lebowitz, Percus and Verlet [41] (cid:68) δ ˆ Aδ ˆ B (cid:69)(cid:12)(cid:12)(cid:12) I = (cid:68) δ ˆ Aδ ˆ B (cid:69)(cid:12)(cid:12)(cid:12) X + ∂ ( βI ) ∂X ∂ (cid:104) ˆ A (cid:105) ∂ ( βI ) ∂ (cid:104) ˆ B (cid:105) ∂ ( βI ) (23)with X = V γ denoting again the extensive variable and I = τ the conjugated intensive variable [42]. This gives σ F | τ = σ F | γ + G eq (24)i.e. the thermodynamic shear modulus G eq compares theshear stress fluctuations in the conjugated ensembles atconstant mean shear stress τ and imposed shear strain γ .The latter formula can be made more useful for compu-tational studies by rewriting the shear stress fluctuations σ F | τ at constant shear stress τ . Note that the normal-ized weight of a state s in the NV τ T-ensemble is given by p ( γ ) ∼ exp[ − β ( ˆ H ( γ ) − V γτ )]. Using the instantaneousshear stress ˆ τ defined in Eq. (8) we thus have p (cid:48) ( γ ) = − βV [ˆ τ ( γ ) − τ )] p ( γ ) . (25)Using integration by parts it is then readily seen [21] thatthis leads to σ F | τ = (cid:104) ˆ H (cid:48)(cid:48) ( γ ) /V (cid:105)| τ = (cid:104) ˆ H (cid:48)(cid:48) ( γ ) /V (cid:105)| γ = µ A in agreement with Eq. (9). This confirms Eq. (22) [33]. Simplifications.
Thermal averaging of Eq. (10) andEq. (11) implies τ id = µ A , id = P id with P id being the idealnormal pressure. We note further that σ F = σ F , id + σ F , ex may again be rewritten as the sum of an ideal con-tribution σ F , id ≡ βV (cid:104) δ ˆ τ (cid:105) = P id and an excess term σ F , ex ≡ βV (cid:104) δ ˆ τ (cid:105) . All ideal contributions to G eq thuscancel and one may rewrite Eq. (22) as G eq = µ A , ex − σ F , ex . (26)Since an ideal gas must have a vanishing shear modulus,this simplification is, of course, expected and could havebeen used from the start [21]. Using Eq. (10) and the factthat in an isotropic system all coordinates are equivalentit is seen that the average shear elasticity reduces to [21] µ A , ex = µ B − P ex with µ B = 1 V (cid:88) l (cid:10)(cid:0) r l u (cid:48)(cid:48) ( r l ) − r l u (cid:48) ( r l ) (cid:1) n l,x n l,y (cid:11) (27)being the Born-Lam´e term [40, 43] and P ex the excesscontribution to the normal pressure P = P id + P ex . B. Demonstration of Eq. (2)
Static limits.
As shown by the dash-dotted line inpanel (c) of Fig. 1, by definition C ( t ) → σ F for t → C ( t ) → t → ∞ [34]. Equation (2) thus impliesthat the relaxation modulus becomes G ( t ) → σ F + ( µ A − σ F ) = µ A for t → + , (28)which is consistent with an affine canonical shear strainimposed at t = 0, Eq. (7), and G ( t ) → G eq for t → ∞ asexpected from Eq. (1). Fluctuation Dissipation Theorem.
We show next thatEq. (2) must hold for all times. Using Boltzmann’s su-perposition principle for an arbitrary strain history γ ( t )[1] the shear stress becomes [3] τ ( t ) = (cid:90) t −∞ d s G ( t − s ) dγ ( s ) ds (29)= G ( t − s ) γ ( s ) | t −∞ − (cid:90) t −∞ d s dG ( t − s ) ds γ ( s )using integration by parts in the second step. Introduc-ing the “after-effect function” or “dynamic response func-tion” [37] χ ( t ) ≡ − G (cid:48) ( t ) = G (cid:48) ( − t ) this may be rewrittenfor a step strain imposed at t = 0 as G ( t ) = G eq + (cid:90) ∞ t χ ( s ) d s (30)being consistent with the expected G ( t ) → G eq for t → ∞ . Since according to the FDT as formulatedby Eq. (7.6.13) of Ref. [34], the after-effect function is χ ( t ) = − C (cid:48) ( t ), this yields the claimed relation Eq. (2). Alternative demonstration.
It is of some importancethat Eq. (2) may be also obtained from the general trans-formation relation Eq. (23) with A = τ ( t ) and B = τ (0).It is assumed here that the shear-barostat imposing amean shear stress τ = 0 is sufficiently slow such thatthe system trajectory is not altered over the time scales used for the determination of the correlation functions[4, 41, 44]. Generalizing the relation between the staticfluctuations Eq. (24) into the time domain, this yieldsimmediately C ( t ) | τ = C ( t ) | γ + G eq . (31)As before for the static shear stress fluctuations σ F | τ onecan show that C ( t ) | τ = (cid:28) ∂ ˆ τ ( t ; γ ) ∂γ (cid:29)(cid:12)(cid:12)(cid:12)(cid:12) τ = G ( t ) for t > τ ( t ; γ ) − τ ][ˆ τ (0; γ ) − τ ] p ( γ ) using Eq. (25). In the second stepwe have used that within linear response G ( t ) does notdepend on γ . Using Eq. (32) and C ( t ) ≡ C ( t ) | γ , Eq. (31)implies again Eq. (2) [33]. C. Oscillatory shear
Experimentally, the relaxation modulus G ( t ) is, ofcourse, commonly sampled in a linear viscoelastic mea-surement [1, 2] using an oscillatory shear imposing, e.g., asinusoidal shear strain γ ( t ) = γ sin( ωt ) of amplitude γ and angular frequency ω as shown in panel (d) of Fig. 1.This implies an average shear stress τ ( t ) = τ ref + γ ( G (cid:48) ( ω ) sin( ωt ) + G (cid:48)(cid:48) ( ω ) cos( ωt )) (33)with τ ref being the (not necessarily vanishing) referenceshear stress at γ = 0. The storage modulus G (cid:48) ( ω ) andthe loss modulus G (cid:48)(cid:48) ( ω ) may thus be determined as theFourier coefficients G (cid:48) ( ω ) = 2 pT ω (cid:90) pT ω sin( ωt ) τ ( t ) γ d t (34) G (cid:48)(cid:48) ( ω ) = 2 pT ω (cid:90) pT ω cos( ωt ) τ ( t ) γ d t (35)with T ω = 1 / πω being the period of the oscillation.Averages over instantaneous shear stresses ˆ τ ( t ) are per-formed over all periods p sampled. As stated, e.g., byEq. (7.149) and Eq. (7.150) of Ref. [1] both moduli areon the other side quite generally given by the Fourier-Sineand Fourier-Cosine transforms of G ( t ) G (cid:48) ( ω ) − G eq = ω (cid:90) ∞ ( G ( t ) − G eq ) sin( ωt )d t, (36) G (cid:48)(cid:48) ( ω ) = ω (cid:90) ∞ ( G ( t ) − G eq ) cos( ωt )d t. (37)The latter two relations together with Eq. (2) imply thekey relations Eq. (3) and Eq. (4) announced in the Intro-duction. We shall pay special attention to the low- ω andhigh- ω limits of G (cid:48) ( ω ) and G (cid:48)(cid:48) ( ω ) at the end of Sec. V C. FIG. 2: Snapshot of a small subvolume of linear length 10containing about 100 verticies of the elastic network consid-ered in this work assuming Eq. (38). The total periodic box oflinear length L ≈ . N = 10 verticies and 9956harmonic springs. The light grey circles indicate the positionsand diameters of the quenched polydisperse LJ bead systemwhich was used for the construction of the elastic network asdescribed in Sec. IV. The lines represent the quenched forcesof the athermal ( T = 0) reference configuration. The dark(black) lines indicate repulsive forces between the verticies,while the light (red) lines represent tensile forces. The linewidth is proportional to the tension/repulsion. III. MODEL AND ALGORITHMIC DETAILS
Hamiltonian.
To illustrate our key relations wepresent in Sec. V numerical data obtained using a pe-riodic two-dimensional ( d = 2) network of harmonicsprings. The model Hamiltonian is given by the sumˆ H = ˆ H id + ˆ H ex of a kinetic energy contribution ˆ H id = m (cid:80) i v i / m ) and anexcess potentialˆ H ex = (cid:88) l u l ( r l ) with u l ( r ) = 12 K l ( r − R l ) (38)where K l denotes the spring constant, R l the referencelength and r l = | r i − r j | the length of spring l . Thesum runs over all springs l connecting pairs of beads i and j at positions r i and r j . A small subvolume of thenetwork considered is represented in Fig. 2. An experi-mentally relevant example for such a permanent networkis provided by endlinked or vulcanized polymer networks[1, 28, 29]. Since the network topology is by construction permanently fixed , the shear response G ( t ) must becomefinite for t → ∞ for all temperatures at variance to sys-tems with plastic rearrangements as considered, e.g., inRef. [10]. Note that the particle mass m and Boltzmann’sconstant k B are set to unity. Lennard-Jones (LJ) units[4] are assumed throughout this paper. Computational methods, parameters and observables.
The construction and the characterization of the refer-ence network at zero temperature is presented in Sec. IV.As discussed in Sec. V, this network is then investi-gated numerically by means of standard molecular dy-namics (MD) simulation [4, 5] at constant particle num-ber N = 10 , box volume V ≈ . and a small, butfinite mean temperature T = 0 . δt MD = 10 − . The temperature T is fixedusing a Langevin thermostat with a friction constant ζ ranging from ζ = 0 .
01 up to ζ = 5 as specified below. Animportant property sampled in this study is the instan-taneous stress tensor [4, 37]ˆ σ αβ = − mV N (cid:88) i =1 v i,α v i,β + 1 V (cid:88) l r l u (cid:48) l ( r l ) n l,α n l,β (39)where the Greek indices α and β denote the d = 2 spa-tial directions x and y , v i,α the α -component of the ve-locity of particle i and n l,α the α -component of the nor-malized vector r l . The first term in Eq. (39) gives theideal gas contribution, the second term corresponds tothe Kirkwood excess stress for pair interactions. Thetrace of the stress tensor yields the instantaneous nor-mal pressure ˆ P = − (ˆ σ xx + ˆ σ yy ) /d , the off-diagonal el-ements ˆ σ xy = ˆ σ yx are consistent with the shear stressˆ τ = ˆ τ id + ˆ τ ex , Eqs. (10,12), obtained in Sec. II A. Theshear strain γ is set to zero for the computation of theautocorrelation function C ( t ) = βV (cid:0) (cid:104) ˆ τ ( t ) (cid:105) − (cid:104) ˆ τ (cid:105) (cid:1) . Atiny step strain γ (cid:28) t = 0 in order tocompute G ( t ), as discussed in Sec. V B, and a sinusoidalstrain γ ( t ) = γ sin( ωt ) to measure directly the storageand loss moduli G (cid:48) ( ω ) and G (cid:48)(cid:48) ( ω ) from the sine and co-sine transforms of the instantaneous shear stress ˆ τ ( t ),Eqs. (34,35). As already described in the first paragraphof Sec. II A we perform in both cases affine, canonical [35]and (essentially) infinitessimal strain transformations ofthe box shape and the particle positions and momenta,Eqs. (5,6). Note that the transformation of the momentais, however, not important for the simulations presentedhere where we focus on low temperatures and mainly onhigh values of the Langevin friction constant. Moreover,as shown in Ref. [26], similar results, especially concern-ing Eq. (2), are also obtained using Brownian dynamicsor Monte Carlo simulations [4, 5]. IV. REFERENCE CONFIGURATION
Construction of network.
As explained in detail inRef. [21], our network has been constructed using the dy-namical matrix M of a polydisperse LJ bead glass com-prising N = 10 particles. Prior to forming the networka LJ bead system has been quenched to T = 0 using aconstant quenching rate and imposing a normal pressure P = 2. The original LJ beads are represented in Fig. 2 bygrey polydisperse circles, the permanent spring network FIG. 3: Energy difference per volume e after an affine strain γ is applied (circles) and after the network has relaxed to thenew ground state (squares). Main panel: Double-logarithmicrepresentation of e vs. γ . The thin solid line indicatesEq. (40), the bold solid line Eq. (41). Inset: Half-logarithmicrepresentation of reduced energy increment ( e − τ γ ) / ( γ / vs. γ where the residual shear stress τ = 0 . G eq ≈ . created from the quenched bead system by lines betweenverticies. The network has a number density ρ ≈ . L ≈ . P ≈ τ ≈ . µ A ≈ . τ ofthe reference for all the correlation functions considered.We also note that the force network, Fig. 2, is stronglyinhomogeneous with zones of weak attractive links em-bedded within a strong repulsive skeleton as discussed inRefs. [13, 14]. Affine displacements.
As shown in Fig. 3, it is possi-ble to obtain numerically the shear stress τ and the affineshear elasticity µ A from the energy increment per volume e ≡ ( ˆ H ex ( γ ) − ˆ H ex (0)) /V caused by an affine shear strain γ (spheres). In agreement with Eq. (7), the energy in-creases as e = τ γ + µ A γ / γ are indicatedby, respectively, the dashed line and the dash-dotted line.We have used for the coefficients τ and µ A in Eq. (40) FIG. 4: Snapshot of non-affine displacement field after anaffine shear strain with γ = 0 .
01 is applied and the systemis allowed to relax. The arrow length is proportionel to thenon-affine particle displacements. The displacement field iscorrelated over distances much larger than the typical particledistance. All displacements are strictly linear in γ . the values obtained using Eq. (12) and Eq. (13). This ismerely a self-consistency check since both properties areactually defined assuming a virtual affine strain transfor-mation as reminded in Sec. II A [3, 21]. As shown in theinset, an accurate verification of µ A over essentially thefull γ -range is obtained by plotting in half-logarithmiccoordinates the reduced energy ( e − τ γ ) / ( γ / γ an even more precise value of the substractedresidual shear stress τ is required to collapse even thesedata points on the dash-dotted horizontal line. Non-affine displacements.
The forces f i acting on theverticies i of an affinely strained network do not vanishin general. Following Ref. [14] the network is first re-laxed by steepest descend, i.e. imposing displacementsproportional to the force, and then by means of the con-jugate gradient method [6]. These non-affine displace-ments lower the final energy of the strained system asindicated by the squares in Fig. 3. As shown by the boldsolid line in the main panel, these final energies scale as e = τ γ + G eq γ / . (41)This scaling is similar to the affine strain energy, Eq. (40),having the same linear term but with µ A being replacedby the shear modulus G eq . As shown in the inset, theseenergies can thus be used to determine G eq ≈ . γ = 0. It follows from Eq. (40) and Eq. (41) thatdue to the non-affine displacements a substantial fraction1 − G eq /µ A ≈ / FIG. 5: Determination of shear modulus G eq . Inset: Low-est eigenfrequencies ω p obtained by diagonalization of thedynamical matrix. The bold horizontal lines correspond to G eq /ρ ≈
16 according to Eq. (43) fitting the eigenvalues for p = 3 , , ,
6. Main panel: Fourier transformed transversedisplacement field for a temperature T = 0 .
001 as a func-tion of wavevector q . The horizontal line corresponds againto G eq ≈
16. The open symbols indicate data obtained us-ing the ground state positions as reference, the filled symbolsassume the average vertex position as reference. for large γ in agreement with Ref. [14]. A snapshot ofthese non-affine displacements is given in Fig. 4. Notethat the non-affine displacements are correlated over dis-tances much larger than the typical particle distance andone thus expects deviations from standard continuummechanics if similar length scales are probed. Eigenstates of the reference.
Following Refs. [13, 14]the shear modulus G eq of the network at T = 0 mayalternatively be computed from the lowest non-trivialeigenfrequencies ω p . (The running index p increaseswith frequency.) This is shown in the inset of Fig. 5.The eigenfrequencies ω p have been determined by di-agonalization of the dynamical matrix M by means ofLanczos’ method [6]. It follows from continuum elastic-ity [8] that the eigenfunctions must be planewaves withwavevectors q quantified by the boundary conditions, i.e.( q x , q y ) = (2 π/L ) ( n, m ) with n, m = 0 , ± , . . . being twoquantum numbers. An example for such a planewave isgiven in Fig. 6 for p = 4 and the pair of quantum numbers( n, m ) = (1 , λ p = 2 π/ | q | = L/ (cid:112) n + m (42)and for the eigenfrequencies of the transverse modes ω p = 2 πc ⊥ L (cid:112) n + m with c ⊥ = (cid:113) G eq /ρ (43)being the transverse wavevelocity. As shown by the hori-zontal lines, we obtain by fitting Eq. (43) to the frequen-cies for p = 3 , , , n + m = 1) that FIG. 6: Eigenvector for p = 4 ( n = 1 , m = 0) correspond-ing to a transverse planewave with wavector q pointing inthe horizontal direction, i.e. with displacements essentially inthe vertical direction, and wavelength λ = L . Interestingly,heterogeneous deviations from the planewave solution of con-tinuum mechanics are even visible for this low eigenmode. G eq ≈
16 and c ⊥ ≈
4. Interestingly, the degeneracy ofthe eigenvalues expected from continuity elasticity is al-ready lifted for p = 7 , , ,
10 ( n + m = 2), i.e. the boxsize L does not allow a precise determination of G eq usingthese eigenvalues. Deviations from the planewave solu-tion are even visible from the eigenvector for p = 4 pre-sented in Fig. 6. Continuum elasticity must break downin any case if wavelengths of order of about the particledistance are probed, λ ≈
1. This implies that Eq. (43)can at best hold up to a frequency ω = 2 πc ⊥ /λ ≈ V. FINITE-TEMPERATURECOMPUTATIONAL RESULTSA. Some static properties
Introduction.
Since the temperature T = 0 .
001 israther small, one expects all static properties such as thepressure P or the elastic modulus G eq to be similar totheir ground state values. As we have checked comparingvarious methods one confirms indeed that P ≈ P ex ≈ τ ≈ τ ex ≈ . µ A ≈ µ A , ex ≈ σ F ≈ σ F , ex ≈ G eq ≈
16. The same applies in fact to all smalltemperatures T (cid:28) Displacement correlations.
A finite- T method forcomputing G eq , which is also of experimental relevance[20], is presented in the main panel of Fig. 5. Usingan ensemble of 10 configurations sampled over a totaltime interval 10 with ζ = 5 we obtain first the dis- FIG. 7: Determination of G eq ≈
16 using Eq. (22) as a func-tion of the sampling time t at T = 0 . µ A ≈
34 converges immediately (small filled squares),the stress fluctuations σ F ( t ) are seen to increase monotonouslyto the large- t limit σ F ≈
18 (bold solid line). This is shownfor four different friction constants ζ = 5 , , . .
01. Thethin solid line compares the data for ζ = 5 with the integratedcorrelation function C ( t ) using Eq. (45). placements u i for each vertex particle i of the networkeither by taking the position of the ground state net-work as reference for defining the displacement (opensquares) or alternatively the average particle positionin the ensemble (filled squares). The displacement field u ( r ) = (cid:80) i u i δ ( r − r i ) is then Fourier transformed accord-ing to u ( q ) = (cid:80) Ni =1 exp( iq · r ) u ( r ) / √ N [20]. Note thatthe wavevector q must be commensurate to the squaresimulation box of linear length L , Eq. (42). The compo-nent of u ( q ) perpendicular to q corresponds to the trans-verse component u ⊥ ( q ) of the Fourier transformed dis-placement field. Using that according to the equiparti-tion theorem every independent elastic mode correspondsto an average kinetic or potential energy k B T /
2, contin-uum mechanics implies that [8, 20] y ≡ q (cid:10) | u ⊥ ( q ) | (cid:11) → k B T ρ/G eq for q → (cid:104) . . . (cid:105) is taken over the ensemble ofstored configurations. As can be seen from Fig. 5, y becomes rapidly constant below q ≈ G eq ≈
16 as indicated by the bold horizontal line. Incidentally,both definitions of the displacement field yield identicalresults. Since only the second definition using the averageparticle position can normally be used in experimentalstudies, this is rather reassuring. Note also that the firstpeak for large wavevectors corresponds to a wavelength λ = 2 π/q ≈ Stress fluctuations.
While it is thus possible to get G eq from the displacement correlations, it is from thecomputational point of view more convenient to deter-mine the modulus using the stress-fluctuation formula, FIG. 8: Stress relaxation modulus G ( t ) and autocorrelationfunction C ( t ) sampled by means of a Langevin thermostat offriction constant ζ = 5. An affine step strain γ = 0 .
01 isapplied at t = 0 to determine G ( t ). Eq. (22). This is shown in Fig. 7 where the affineshear elasticity µ A (small filled squares) and the stress-fluctuation term σ F ( t ) are presented as functions of thesampling time t . (The notation σ F without time argu-ment refers to the static thermodynamic large- t limit,while σ F ( t ) indicates that this property has been deter-mined using a finite time window t .) While the sim-ple average µ A is obtained immediately, σ F ( t ) is seen toincrease monotonously from zero to the large- t plateau σ F ≈
18 (horizontal bold solid line). Friction effects areimportant only in the intermediate sampling time win-dow between t ≈ . ζ = 0 . .
01 are essentially identical.The time dependence of σ F ( t ) for different friction con-stants is further analysed at the end of Sec. V B. B. Time domain of key relations
High friction limit.
We turn now to the numericalverification of the key relations stated in the Introduc-tion focusing first on results obtained in the time domainwhich have already been reported in Ref. [26]. The datapresented in Fig. 8 have been computed using a Langevinthermostat with one high friction constant ζ = 5. Thisis done merely for presentational reasons since in thislimit the kinetic degrees of freedom can be disregardedand oscillations are suppressed. The stress relaxationmodulus G ( t ) given has been computed as indicated inpanel (b) of Fig. 1 by applying an affine canonical shearstrain γ = 0 .
01 and by averaging over 1000 runs start-ing from independent configurations at t = 0 − . Due tothe strong damping G ( t ) decreases monotonously from G (0 + ) = µ A to G eq while C ( t ) decays from C (0) = σ F to zero. Confirming Eq. (2), only after vertically shifting FIG. 9: Stress relaxation modulus G ( t ) obtained using a stepstrain (open symbols) and the rescaled autocorrelation func-tion G eq + C ( t ) at γ = 0 (filled symbols) for several frictionconstants ζ . Inset: Shear stress relaxation time θ ( ζ ) measuredusing θ /e , θ and θ defined in the main text. C ( t ) → C ( t ) + G eq one obtains a collapse on the directlycomputed modulus G ( t ). This is the most important nu-merical result of the present work. Friction effects and stress relaxation time.
A simi-lar scaling collapse of G ( t ) and C ( t ) + G eq has beenalso obtained for different friction constants ζ as maybe seen from Fig. 9. As one expects, the MD data de-cays more rapidly with decreasing ζ and reveals anti-correlations for the lowest ζ probed. The effect of thefriction constant may be compactly described using thestress relaxation time θ shown in the inset of Fig. 9. Wecompare here three different definitions. The definition C ( t = θ /e ) = σ F /e , indicated by spheres, describesthe time where the correlation functions first becomessmall. This simple definition underestimates correlationsat larger times. Alternatively, one may attempt to define θ by the ratio c /c of the first to the zeroth moment of C ( t ) with c n ≡ (cid:82) ∞ t n C ( t )d t (not shown). Unfortunately, c is found to converge badly and we have not been ableto obtain reliable values over the full ζ -range. A nu-merically well-behaved alternative is based on an exactrelation between the autocorrelation function C ( t ) andthe monotoneously increasing fluctuation term σ F ( t ) de-termined as a function of the sampling time t . Assumingtime-translational invariance it can be shown [26] that y ( t ) ≡ − σ F ( t ) σ F = 2 t (cid:90) t d s (1 − s/t ) C ( s ) σ F , (45)i.e. the dimensionless ratio y ( t ) characterizes the dif-ference of the zeroth and the first moment of the au-tocorrelation function integrated up to t . That this re-lation holds can be seen from the thin line presented inFig. 7. We note en passant that y ( t ) must ultimately van-ish slowly as 1 /t since the integral over C ( s ) in Eq. (45) becomes constant [26]. (This 1 /t -decay is consistent withthe general finite-sampling time corrections for fluctua-tions [7].) Since σ F ( t ) can be determined directly, onemay define θ ( ζ ) using a fixed ratio y ( θ ). We have pre-sented in Fig. 9 the constants y ( θ ) = 5% (squares)and y ( θ ) = 1% (triangles). Note that θ correspondsto the time where σ F ( t ) begins to saturate, while θ indicates the time where the stress-correlations becomeneglible and σ F ( t ) thus allows a good estimation of G eq using the stress-fluctuation formula. C. Oscillatory shear
Direct determination of G (cid:48) ( ω ) and G (cid:48)(cid:48) ( ω ) . As alreadystressed in Sec. II C, the relaxation modulus G ( t ) iscommonly determined experimentally by inverse Fouriertransformation of the storage modulus G (cid:48) ( ω ) and/or theloss modulus G (cid:48)(cid:48) ( ω ) obtained in a linear viscoelastic mea-surement imposing an oscillatory shear [1, 2]. Motivatedby this we perform non-equilibrium MD simulations inthe linear response limit by imposing a sinusoidal shearstrain of frequency ω with an amplitude γ = 0 . γ it has been checked that all reported valuesare in the linear regime.) This is done by performingevery time step δt MD = 10 − an affine canonical strain,Sec. II A. We use production runs of total length pT ω with T ω = 2 π/ω and at least p = 100 oscillation peri-ods, i.e. the computational load increases as p/ω withdecreasing frequency. This sets a lower limit ω ≈ − for the angular frequencies we have been able to sam-ple. Two production runs are compared to rule out tran-sient behavior. The instantaneous shear stress ˆ τ ( t ) issampled each δt MD . Using Eq. (34) and Eq. (35) oneobtains G (cid:48) ( ω ) and G (cid:48)(cid:48) ( ω ) as indicated by open trianglesin Fig. 10, where we focus for simplicity on one value ζ = 5 of the friction constant, and large open symbols inFig. 11, where data for different ζ are presented. Test of key relations.
These values of G (cid:48) ( ω ) and G (cid:48)(cid:48) ( ω ) are compared in Fig. 10 with the Sine-Fourier andCosine-Fourier transforms of the shear modulus G ( t ) ob-tained from the step strain experiment (circles) and ofthe shear stress autocorrelation function C ( t ). To reducethe well-known (but considering our tiny time step δt MD surprisingly severe) numerical problems at large frequen-cies ( ω (cid:29) G (cid:48) ( ω ) and G (cid:48)(cid:48) ( ω ). The crucial point is here that it is G eq + C (cid:48) ( ω ) which collapses on G (cid:48) ( ω ), while C (cid:48) ( ω ) does not being much too small considering that G eq ≈
16, asseen in panel (a). We have verified that a similar data col-lapse for G eq + C (cid:48) ( ω ) onto G (cid:48) ( ω ) and C (cid:48)(cid:48) ( ω ) onto G (cid:48)(cid:48) ( ω ) isalso obtained for other friction constants as may be seenfrom Fig. 11. Having settled this fundamental scalingissue let us turn finally to the description of the generalshape and asymptotic behavior of both moduli.0 FIG. 10: Oscillatory shear moduli (a) G (cid:48) ( ω ) using a log-linearrepresentation and (b) G (cid:48)(cid:48) ( ω ) using a double-logarithmic rep-resentation for a large friction constant ζ = 5. We comparethe Fourier coefficients (triangles) obtained by applying a si-nusoidal shear strain with amplitude γ = 0 .
001 to the Fouriertransformed step-strain relaxation modulus G ( t ) (circles) andshear autocorrelation function C ( t ) (squares). The bold solidlines indicate the expected asymptotic limits for small ω , thedashed lines the large- ω asymptotics, the thin vertical linesthe frequencies ω ≈
25 and ω ≈
15 of the maxima of bothmoduli. The latter values coincide with the breakdown ofcontinuum elasticity at a planewave frequency ω = 2 πc ⊥ /λ with λ ≈ (b) :Phase angle δ with maximum δ ≈ at ω ≈ Low-frequency limit.
The low- ω limit of G (cid:48) ( ω ) and G (cid:48)(cid:48) ( ω ) are readily obtained by expanding sin( x ) in Eq. (3)and cos( x ) in Eq. (4). As one expects [1], one obtains toleading order that G (cid:48) ( ω ) → G eq , Eq. (1), and G (cid:48)(cid:48) ( ω ) → ωc as indicated by bold solid lines in Fig. 10 and Fig. 11.We have fitted G (cid:48)(cid:48) ( ω ) for ζ = 5 using the independentlydetermined value c ≈ Intermediate peaks.
In the intermediate frequencyregime one observes a striking peak at ω ≈
25 for G (cid:48) ( ω ),i.e. at the planewave frequency corresponding to the typ- FIG. 11: Effect of Langevin friction constant ζ with open sym-bols referring to data obtained using oscillatory shear withamplitude γ = 0 . (a) Storage modulus G (cid:48) ( ω ). The fric-tion constant ζ only matters in an intermediate frequencyrange between ω ≈ . ω ≈
10 below the peak at ω ≈ (b) Loss modulus G (cid:48)(cid:48) ( ω ). The dot-dashed line is a phe-nomenological fit for large ω and small ζ assuming Eq. (47). ical particle distance, Eq. (43), and at a slightly smallervalue ω ≈
15 for G (cid:48)(cid:48) ( ω ). The phase angle δ characteriz-ing the ratio G (cid:48)(cid:48) ( ω ) /G (cid:48) ( ω ) ≡ tan( δ ) is presented in theinset (in grad). It reveals a maximum at ω ≈ High-frequency limit.
In the high-frequency limit itbecomes increasing difficult to relax the imposed affinestrain by non-affine displacements and, hence, to dissi-pate the work done on the system. This implies that G (cid:48) ( ω ) = G eq + C (cid:48) ( ω ) → ( µ A − σ F ) + σ F = µ A (46)for ω → ∞ (bold dashed horizontal lines), while G (cid:48)(cid:48) ( ω )decays with a non-universal (see below) sharp cutoff(dashed and dash-dotted lines). This asymptotic behav-ior can be also understood on mathematical grounds byexpanding Eq. (3) and Eq. (4) by integration by parts.This leads to two series-expansions in terms of even and1odd derivatives of C ( t ) taken at t = 0 for, respectively, G (cid:48) ( ω ) and G (cid:48)(cid:48) ( ω ). We remind that C ( t ) is an even func-tion for classical systems [4, 34] and that in a MD simula-tion C ( t ) must become analytic at sufficiently short times( δt MD < t (cid:28) /ζ ) where thermostat effects are negligi-ble. (Strict non-analyticity of C ( t ) at t = 0 is formallyfound, e.g., for the Maxwell model or the polymer Rousemodel [3].) Analyticity at t = 0 implies that the (con-verging) expansion of G (cid:48) ( ω ) is dominated by its first termfor ω → ∞ . This demonstrates Eq. (46). The effectiveodd derivatives of C ( t ) must, however, become negligiblewith decreasing friction ζ , increasing angular frequency ω and decreasing time step δt MD . Hence, G (cid:48)(cid:48) ( ω ) mustrapidly vanish. This point is further investigated in thefinal paragraph of this section. Friction effects.
The friction effects presented inFig. 11 are rather small for G (cid:48) ( ω ) due to the large staticconstants G eq and µ A dominating the storage modulus.As one would expect, they are more marked for the lossmodulus G (cid:48)(cid:48) ( ω ), especially at low frequencies where thethermostat has more time to dissipate the applied me-chanical energy, but also at the large- ω cutoff. For ζ = 5one might be tempted to fit an 1 /ω -decay (dashed lines)which indicates an effective non-analytical behavior of C ( t ) at short times. With decreasing ζ the loss modulusis, however, better described by the exponential decay G (cid:48)(cid:48) ( ω ) = σ F (cid:112) π/ x exp( − x /
2) with x = ω ˜ θ (47)as indicated by the dash-dotted line in panel (b) ofFig. 11. This phenomenological fit is the Cosine-Fourier transform assuming a Gaussian C ( t ) = σ F exp( − ( t/ ˜ θ ) /
2) with the parameter ˜ θ ≈ . ω behaviorof is also obtained by assuming a Lorentzian C ( t ) = σ F / (1+( t/ ˜ θ ) ). We note finally that if smaller time stepsand larger frequencies could be computed, one expectseven for higher friction constants a similar non-algebraiccutoff, albeit with a different parameter ˜ θ ( ζ ). VI. CONCLUSION
Summary.
We have studied in the present work sim-ple isotropic solids formed by permanent spring networks.Some relevant static properties have been reviewed andcharacterized in Secs. II A, IV and V A. More impor-tantly, we have reconsidered theoretically (Sec. II) andnumerically (Sec. V) the linear-response relation betweenthe shear stress relaxation modulus G ( t ) and the shearstress autocorrelation function C ( t ) both in the time(Sec. V B, Figs. 8-9) and the frequency domain (Sec. V C,Figs. 10 and 11). According to our key relations, Eqs. (2-4), G ( t ) and C ( t ) must become different in the solid limit( G eq >
0) and it is thus impossible to determine G eq from C ( t ) or its Fourier transforms C (cid:48) ( ω ) and C (cid:48)(cid:48) ( ω ).The short-time behavior of C ( t ) and the large- ω limit of C (cid:48) ( ω ) only yield the stress-fluctuation contribution σ F to the equilibrium modulus G eq = µ A − σ F (Sec. II A) asalready stressed elsewhere [17, 26]. The peak frequen-cies observed for G (cid:48) ( ω ) and G (cid:48)(cid:48) ( ω ) (Figs. 10 and 11)are qualitatively expected [13, 14] from the breakdownof continuum mechanics for large wavevectors q (cid:29) Discussion.
It is obviously common and may oftenbe helpful to describe a plateau of C ( t ) at short or inter-mediate times (or equivalently at large or intermediatefrequencies for C (cid:48) ( ω )) in terms of a finite shear modulus G M of a dynamical relaxation model, such as the Maxwellmodel for viscoelastic fluids or the reptation model ofentangled polymer melts [1, 3]. However, such a modelallowing the theoretical interpretation of the data shouldnot be confused with the proper measurement procedure ,Eq. (1), and a finite model parameter G M not with theequilibrium modulus G eq of the system which, inciden-tally, vanishes both for a Maxwell fluid or a linear poly-mer melt. In this sense different operational “static” and“dynamical” definitions of the shear modulus are used inthe literature for describing glass-forming liquids closeto the glass transition [17–20, 22]. This may explainwhy qualitatively different temperature dependences —cusp singularity [17, 22] vs. finite jump [11, 18, 20] —have been predicted recently. Hence, while our recent at-tempts to determine G eq ( T ) for two glass-forming modelsystems [21] are consistent with a continuous cusp, this isnot necessarily in contradiction with a jump singularityfor G M ( T ) determined from an intermediate shoulder of C ( t ) [20, 25]. Outlook.
We note finally that generalizing our key re-lations one obtains readily that M ( t ) = C ( t ) + M eq for t > M ( t ) of any continuous inten-sive variable I with M eq = ∂I/∂X being the equilibriummodulus and C ( t ) = βV (cid:104) δ ˆ I ( t ) δ ˆ I (0) (cid:105) the correspondingautocorrelation function computed at a constraint ther-modynamically conjugated extensive variable X . In thefrequency domain this leads to M (cid:48) ( ω ) = M eq + C (cid:48) ( ω )and M (cid:48)(cid:48) ( ω ) = C (cid:48)(cid:48) ( ω ). A natural example is provided bythe relaxation modulus M ( t ) = K ( t ) associated to thenormal pressure I = P . While C ( t ) and C (cid:48) ( ω ) must ob-viously vanish in the static limit for, respectively, largetimes and small frequencies, K ( t ) and K (cid:48) ( ω ) approach afinite compression modulus M eq = K eq for stable (non-critical) thermodynamic systems. Acknowledgments
H.X. thanks the IRTG Soft Matter for financial sup-port. We are indebted to H. Meyer (ICS, Strasbourg)and M. Fuchs (Konstanz) for helpful discussions.2 [1] M. Rubinstein and R. Colby,
Polymer Physics (OxfordUniversity Press, Oxford, 2003).[2] T. Witten and P. A. Pincus,
Structured Fluids: Polymers,Colloids, Surfactants (Oxford University Press, Oxford,2004).[3] M. Doi and S. F. Edwards,
The Theory of Polymer Dy-namics (Clarendon Press, Oxford, 1986).[4] M. Allen and D. Tildesley,
Computer Simulation of Liq-uids (Oxford University Press, Oxford, 1994).[5] D. Frenkel and B. Smit,
Understanding Molecular Sim-ulation – From Algorithms to Applications (AcademicPress, San Diego, 2002), 2nd edition.[6] J. Thijssen,
Computational Physics (Cambridge Univer-sity Press, Cambridge, 1999).[7] D. P. Landau and K. Binder,
A Guide to Monte CarloSimulations in Statistical Physics (Cambridge UniversityPress, Cambridge, 2000).[8] L. D. Landau and E. M. Lifshitz,
Theory of Elasticity (Pergamon Press, 1959).[9] S. Alexander, Physics Reports , 65 (1998).[10] F. Sausset, G. Biroli, and J. Kurchan, J. Stat. Phys. ,718 (2010).[11] W. G¨otze,
Complex Dynamics of Glass-Forming Liquids:A Mode-Coupling Theory (Oxford University Press, Ox-ford, 2009).[12] J.-L. Barrat, J.-N. Roux, J.-P. Hansen, and M. L. Klein,Europhys. Lett. , 707 (1988).[13] J. P. Wittmer, A. Tanguy, J.-L. Barrat, and L. Lewis,Europhys. Lett. , 423 (2002).[14] A. Tanguy, J. P. Wittmer, F. Leonforte, and J.-L. Barrat,Phys. Rev. B , 174205 (2002).[15] L. Berthier, G. Biroli, J.-P. Bouchaud, L. Cipelletti, D. E.Masri, D. L’Hˆote, F. Ladieu, and M. Pierno, Phys. Rev.Lett. , 1797 (2005).[16] L. Berthier, G. Biroli, J.-P. Bouchaud, W. Kob,K. Miyazaki, and D. Reichman, J. Chem. Phys. ,184503 (2007).[17] H. Yoshino and M. M´ezard, Phys. Rev. Lett. , 015504(2010).[18] G. Szamel and E. Flenner, Phys. Rev. Lett. , 105505(2011).[19] H. Yoshino, J. Chem. Phys. , 214108 (2012).[20] C. Klix, F. Ebert, F. Weysser, M. Fuchs, G. Maret, andP. Keim, Phys. Rev. Lett. , 178301 (2012).[21] J. P. Wittmer, H. Xu, P. Poli´nska, F. Weysser, andJ. Baschnagel, J. Chem. Phys. , 12A533 (2013).[22] A. Zaccone and E. Terentjev, Phys. Rev. Lett. ,178002 (2013).[23] J. P. Wittmer, H. Xu, P. Poli´nska, C. Gillig, J. Hellferich,F. Weysser, and J. Baschnagel, Eur. Phys. J. E , 131(2013).[24] H. Mizuno, S. Mossa, and J.-L. Barrat, Phys. Rev. E ,042306 (2013).[25] E. Flenner and G. Szamel, Phys. Rev. Lett. , 105505(2015).[26] J. P. Wittmer, H. Xu, and J. Baschnagel, Phys. Rev. E (2015), in print.[27] E. del Gado and W. Kob, J. Non-Newtonian Fluid Mech. , 28 (2008).[28] E. Duering, K. Kremer, and G. S. Grest, Phys. Rev. Lett. , 67 (1991).[29] E. Duering, K. Kremer, and G. S. Grest, J. Chem. Phys. , 101 (1994).[30] S. Ulrich, X. Mao, P. Goldbart, and A. Zippelius, Euro-physics Lett. , 677 (2006).[31] C. Tonhauser, D. Wilms, Y. Korth, H. Frey, andC. Friedrich, Macromolecular Rapid Comm. , 2127(2010).[32] A. Zilman, J. Kieffer, F. Molino, G. Porte, and S. A.Safran, Phys. Rev. Lett. , 2003 (2003).[33] For simplicity of the notation we write σ F for σ F | γ and C ( t ) for C ( t ) | γ . We only note σ F | γ and σ F | τ or C ( t ) | γ and C ( t ) | τ if fluctuations at imposed γ (NV γ T-ensemble)and τ (NV τ T-ensemble) are compared.[34] J. Hansen and I. McDonald,
Theory of simple liquids (Academic Press, New York, 2006), 3nd edition.[35] H. Goldstein, J. Safko, and C. Poole,
Classical Mechanics (Addison-Wesley, 2001), 3th edition.[36] A more general affine canonical transformation consistsin changing both the particle coordinates and the conju-gated momenta using a linear matrix [4, 40].[37] With the exception of the storage modulus G (cid:48) ( ω ) and theloss modulus G (cid:48)(cid:48) ( ω ) a prime denotes generally a deriva-tive of a function with respect to its argument.[38] H. B. Callen, Thermodynamics and an Introduction toThermostatistics (Wiley, New York, 1985).[39] D. R. Squire, A. C. Holt, and W. G. Hoover, Physica ,388 (1969).[40] J. F. Lutsko, J. Appl. Phys , 2991 (1989).[41] J. L. Lebowitz, J. K. Percus, and L. Verlet, Phys. Rev. , 250 (1967).[42] For the simplicity of the notation we have assumed inEq. (23) that X is not the internal energy U . For a moregeneral theoretical description it is necessary to define the“entropic intensive variable” J ≡ ∂S ( X ) /∂X with S ( X )being the entropy [38]. If X (cid:54) = U , one has J = − I/T [38].These entropic intensive variables are used in Ref. [41].Note that expressing Eq. (23) in terms of I , rather thanin terms of J , changes the signs.[43] M. Born and K. Huang, Dynamical Theory of CrystalLattices (Clarendon Press, Oxford, 1954).[44] We emphasize that the dynamical Lebowitz-Percus-Verlet transform [41] relies on the condition that only thedistribution of start points of the trajectories depends onthe ensemble, but not the relaxation pathways themselves[4]. While this does not hold for extensive variables, thisis generally the case for fluctuations of instantaneous in-tensive variables. Interestingly, a similar approach basedon Ref. [41] has been used for the four-point dynamicsusceptibility χ ( tt