Shear stress relaxation and ensemble transformation of shear stress autocorrelation functions revisited
SShear stress relaxation and ensemble transformationof shear stress autocorrelation functions revisited
J.P. Wittmer, ∗ H. Xu, 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 7, 2018)We revisit the relation between the shear stress relaxation modulus G ( t ), computed at finiteshear strain 0 < γ (cid:28)
1, and the shear stress autocorrelation functions C ( t ) | γ and C ( t ) | τ computed,respectively, at imposed strain γ and mean stress τ . Focusing on permanent isotropic spring networksit is shown theoretically and computationally that in general G ( t ) = C ( t ) | τ = C ( t ) | γ + G eq for t > G eq being the static equilibrium shear modulus. G ( t ) and C ( t ) | γ thus must become differentfor solids and it is impossible to obtain G eq alone from C ( t ) | γ as often assumed. We comment brieflyon self-assembled transient networks where G eq ( f ) must vanish for a finite scission-recombinationfrequency f . We argue that G ( t ) = C ( t ) | τ = C ( t ) | γ should reveal an intermediate plateau set bythe shear modulus G eq ( f = 0) of the quenched network. I. INTRODUCTION
Shear stress relaxation.
The static equilibrium shearmodulus G eq [1–7] is an important order parameter [8–10]characterizing the transition from the liquid/sol ( G eq =0) to the solid/gel state ( G eq >
0) where the particlepermutation symmetry of the liquid state is lost for thetime window probed [6, 7]. Examples of current inter-est for the determination of G eq include crystalline solids[11], glass-forming liquids and amorphous solids [5, 12–27], colloidal gels [28], permanent polymeric networks[2, 29–31], hyperbranched polymer chains with stickyend-groups [32] or networks of telechelic polymers [33].As emphasized by the thin horizontal line in Fig. 1, theshear modulus of an isotropic solid may be determinedexperimentally from the long-time limit [2, 34] G eq ≡ lim t →∞ G ( t ) (1)of the shear stress relaxation modulus (bold solid line) de-fined as G ( t ) ≡ δτ ( t ) /γ . It measures the stress increment δτ ( t ) = (cid:104) ˆ τ ( t ) − ˆ τ (0 − ) (cid:105) due to a step strain 0 < γ (cid:28) t = 0. Here ˆ τ ( t ) denotes the instantaneousshear stress which may be measured experimentally fromthe forces acting on the walls of the shear cell [2]. Correlation functions.
A quantity related to G ( t ) isthe shear stress autocorrelation function [4, 5] C ( t ) ≡ βV (cid:104) δ ˆ τ ( t ) δ ˆ τ (0) (cid:105) ≡ ˜ C ( t ) − βV (cid:104) ˆ τ (cid:105) (2)with β = 1 /k B T being the inverse temperature and V thevolume. We write C ( t ) | γ or C ( t ) | τ if C ( t ) is computed,respectively, in the NV γ T-ensemble at imposed particlenumber N , volume V , shear strain γ and temperature T or in the conjugated NV τ T-ensemble where instead of γ the mean shear stress τ is imposed. The effect of the ∗ Electronic address: [email protected] γ C(t)| =0C(t)| γ C(t)| γ σ | F γ σ | F γ affine response A t G eq G(t)
G(t) non − a ff i n e µ eq G t=0 γ >0 switched on=0 γ FIG. 1: Schematic comparison of the shear relaxation mod-ulus G ( t ) (bold solid line) and the shear stress autocorrela-tion function C ( t ) | γ computed in the NV γ T-ensemble (dash-dotted line). Note that G (0 + ) = µ A = σ F | τ and C (0) | γ = σ F | γ with µ A being the affine Born-Lam´e contribution to theshear modulus G eq = µ A − σ F | γ with σ F ≡ βV (cid:104) δ ˆ τ (cid:105) character-izing the static shear stress fluctuations [12, 14, 17, 23, 35, 36]. latter constraint is assumed to be arbitrarily slow, suchthat γ ( t ) barely changes over the time window probed.This separation of time scales implies˜ C ( t ) | τ = (cid:90) d γ p ( γ ) ˜ C ( t ) | γ (3)with p ( γ ) being the normalized distribution of strains γ in the NV τ T-ensemble. The conceptionally importantuniversal limit, Eq. (3), may be realized experimentallyusing an overdamped external force or computationallyby either using a strong frictional Langevin force addedto a standard molecular dynamics (MD) “shear-barostat”[37–41] imposing an average shear stress τ or (as usedbelow) a Monte Carlo (MC) scheme with a low attempt-frequency for an affine canonical δγ -change [23]. Key issue.
Interestingly, it is often assumed [5, 21, 29,30, 37] that G ( t ) and C ( t ) | γ become generally equivalentin the linear response limit ( γ → G ( t ) = C ( t ) | γ , a r X i v : . [ c ond - m a t . s t a t - m ec h ] A ug the equilibrium shear modulus G eq , Eq. (1), may thenbe identified with some transient plateau G P or “finitefrozen-in amplitude” of C ( t ) | γ [21] and, hence, with the“nonergodicity parameter” of the mode-coupling theoryfor glass-forming liquids [5]. Here we raise concerns withsuch an identification. It will be shown that in fact G ( t ) = C ( t ) | τ = C ( t ) | γ + G eq for t > G ( t ) = 0 for t < intensive variables [37, 42, 43], this is thekey relation we want to stress in this paper. Two im-portant consequences of Eq. (4) are that (i) G ( t ) onlybecomes equivalent to C ( t ) | γ for t > G eq = 0 and that (ii) a finite shear modulus G eq isonly probed by G ( t ) on time scales where C ( t ) | γ actuallyvanishes. While the static shear modulus G eq can be ob-tained from C ( t ) | τ , this is not possible using only C ( t ) | γ without making additional model-specific assumptions. Outline.
We recall in Sec. II A the “affine” contribu-tion µ A and the “stress fluctuation” contribution σ F | γ tothe equilibrium shear modulus G eq = µ A − σ F | γ . Thekey relation Eq. (4) is then demonstrated theoreticallyin Sec. II B using several (albeit not completely inde-pendent) lines of thought. If the stress fluctuation con-tribution σ F ( t ) is determined numerically over a finitetime window t , it must systematically underestimate thevalue σ F for asymptotically long sampling times [44]. Itis seen in Sec. II C how σ F ( t ) is quite generally relatedto the correlation function C ( t ). We briefly commenton self-assembled transient elastic networks in Sec. II D.The specific model system considered numerically is in-troduced in Sec. III. A well-defined solid with finite equi-librium shear modulus G eq for t → ∞ is assumed. Forthis reason we replace the Lennard-Jones (LJ) interac-tions of a quenched bead system by a permanent elasticspring network corresponding to its dynamical matrix atzero temperature [13, 14, 23]. Some static properties andmeasurement procedures are summarized in Sec. IV A.Using our simple model Hamiltonian the key relation isconfirmed numerically in Sec. IV B by means of molecu-lar dynamics (MD), Brownian dynamics (BD) and MonteCarlo (MC) simulations [37–40]. This work is summa-rized in Sec. V. We finally state the generalization ofEq. (4) for autocorrelation functions of other intensivevariables and comment briefly on ongoing simulations ofself-assembled transient networks. II. THEORETICAL CONSIDERATIONSA. Static properties
Static stress fluctuations.
We begin by reminding [23]that the shear modulus G eq of a solid body may be ob- tained in principle from σ F | τ = σ F | γ + G eq (5)by comparing the (reduced) shear stress fluctuations σ F ≡ C ( t = 0) ≡ βV (cid:10) δ ˆ τ (cid:11) (6)at constant mean shear stress τ (NV τ T-ensemble) withthe fluctuations at imposed strain γ (NV γ T-ensemble).This relation is obtained directly from the Lebowitz-Percus-Verlet transformation for a fluctuation (cid:104) δ ˆ Aδ ˆ B (cid:105) oftwo observables A and B [23, 37, 43] (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 ) (7)with X = V γ being in our case the extensive variable, I = τ the conjugated intensive variable and ˆ A = ˆ B = ˆ τ [8]. For the simplicity of the notation we have assumed inEq. (7) 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 [8]. If X (cid:54) = U , one has J = − I/T [8].These entropic intensive variables are used in Ref. [43].Note that expressing Eq. (7) in terms of I , rather thanin terms of J , changes the signs. From fluctuations to simple means.
From the compu-tational point of view it is important that Eq. (5) can befurther simplified. With ˆ H ( γ ) = ˆ H id ( γ ) + ˆ H ex ( γ ) beingthe Hamiltonian of a given state of the system param-eterized in terms of an affine strain γ [23, 35, 36, 41],its normalized weight in the NV τ T-ensemble is given by p ( γ ) ∼ exp[ − β ( ˆ H ( γ ) − V γτ )]. We thus have p (cid:48) ( γ ) = − βV [ˆ τ ( γ ) − τ )] p ( γ ) with ˆ τ ( γ ) ≡ ˆ H (cid:48) ( γ ) /V (8) defining the instantaneous shear stress [23]. (A primedenotes a derivative of a function with respect to its ar-gument.) For small γ it follows that ˆ τ id ≡ ˆ H (cid:48) id ( γ ) /V reduces to the standard instantaneous ideal shear stressand ˆ τ ex ≡ ˆ H (cid:48) ex ( γ ) /V for pair potential interactions to theKirkwood virial expression of the shear stress [23, 37, 45].By integration by parts the stress fluctuation σ F | τ can beexpressed as the “simple average” [23, 24] σ F | τ = 1 V (cid:68) ˆ H (cid:48)(cid:48) ( γ ) (cid:69) = (cid:104) ˆ τ (cid:48) ( γ ) (cid:105) | τ ≡ µ A (9)which can be directly computed in any ensemble assum-ing that the same state point is sampled. The “affineshear elasticity” µ A characterizes the mean total (kineticand excess) energy change µ A V γ / homo-geneous affine shear transformation of the system as itmay be done in a computer experiment by changing themetric of system [35, 37, 38, 41]. For pair potentials µ A can be further reduced to µ A = µ B − P ex + P id (10)with µ B being the well-known Born-Lam´e coefficient, P ex the excess pressure and P id the ideal pressure contribu-tion. We have thus rewritten σ F | τ as a simple average ofmoments of first and second derivatives of the potentialplus P id . (Since second derivatives are considered, impul-sive corrections must be taken into account for truncatedand shifted potentials as stressed in Ref. [22].) The shearmodulus can hence be conveniently computed by meansof the stress-fluctuation formula G eq = G F ≡ µ A − σ F | γ (11)in the NV γ T-ensemble [12, 14, 17, 23, 27, 36]. Since for aplain shear strain at constant volume the ideal free energycontribution does not change, the explicit kinetic energycontributions must be irrelevant for G eq . (An ideal gascannot elastically support a finite shear stress.) As onethus expects, the kinetic contributions µ A , id = σ F , id | γ = P id to µ A and σ F | γ cancel and can be dropped when G eq is determined using Eq. (11). B. Demonstration of key relation
Asymptotic limits.
As shown by the dash-dotted linein Fig. 1, by definition C ( t ) | γ → σ F | γ for t → C ( t ) | γ → t → ∞ [4]. Equation (4) thus impliesthat G ( t ) → σ F | γ + ( µ A − σ F | γ ) = µ A for t → + —which is consistent with the affine shear strain imposedat t = 0 — and G ( t ) → G eq for t → ∞ as it should.We note also that by definition C ( t ) | τ → σ F | τ = µ A for t →
0. Interestingly, the autocorrelation function C ( t ) | τ does not vanish in general in the large- t limit.This is a direct consequence of the time scale separationmentioned above, Eq. (3), from which it is seen that C ( t ) | τ = (cid:90) d γ p ( γ ) C ( t ) | γ + C ∞ with (12) C ∞ ≡ βV (cid:90) d γ p ( γ ) (cid:16) (cid:104) ˆ τ (cid:105)| γ − τ (cid:17) . (13)The first contribution to C ( t ) | τ in Eq. (12) vanishes for t → ∞ . Note that C ∞ differs from σ F | τ = µ A due to theunderlined term in Eq. (13). Using that p ( γ ) is Gaussianand (cid:104) δγ (cid:105) = k B T / ( V G eq ), it is seen that C ( t ) | τ → C ∞ = βV G (cid:10) δγ (cid:11) = G eq for t → ∞ . (14)We show now that Eq. (4) must hold for all times. First equality of the key relation.
Generalizing Eq. (9)one shows for the shear stress fluctuations at constantstress (assuming a slow shear-barostat) that C ( t ) | τ = (cid:28) ∂ ˆ τ ( t ; γ ) ∂γ (cid:29)(cid:12)(cid:12)(cid:12)(cid:12) τ = G ( t ) for t > . (15)To show this, we have reexpressed in the first step[ˆ τ ( t ; γ ) − τ ][ˆ τ (0; γ ) − τ ] p ( γ ) using Eq. (8) and integra-tion by parts. In the second step we have used thatwithin linear response G ( t ) does not depend on γ . Thisdemonstrates the first equality stated in Eq. (4). Second equality of the key relation.
Using Boltz-mann’s superposition principle the shear stress τ ( t ) foran arbitrary strain history γ ( t ) may be written [2, 3] τ ( t ) = (cid:90) t −∞ d s G ( t − s ) dγ ( s ) ds (16)= G ( t − s ) γ ( s ) | t −∞ − (cid:90) t −∞ d s dG ( t − s ) ds γ ( s )using integration by parts. Since γ ( t ) is a step func-tion and introducing the “after-effect function” χ ( t ) ≡− G (cid:48) ( t ) = G (cid:48) ( − t ) [4] this gives G ( t ) = G (0) − (cid:90) t χ ( s ) d s = G eq + (cid:90) ∞ t χ ( s ) d s (17)where G eq appears as an integration constant. Sinceaccording to the FDT as formulated by Eq. (7.6.13)of Ref. [4], the after-effect function is given by χ ( t ) = − C (cid:48) ( t ) | γ , this demonstrates G ( t ) = G eq + C ( t ) | γ as statedby the second equality in Eq. (4) [46]. Alternatively, fromEq. (12) and Eq. (14) one obtains directly C ( t ) | τ → C ( t ) | γ + G eq for V → ∞ (18)using steepest-descent, (cid:82) d γ p ( γ ) C ( t ) | γ → C ( t ) | γ , with γ corresponding to the maximum of p ( γ ). Together withEq. (15) this confirms again our key relation. Dynamical Lebowitz-Percus-Verlet transform.
Inter-estingly, Eq. (18) may be also obtained by generalizingthe Lebowitz-Percus-Verlet transformation, Eq. (7), intothe time domain with ˆ A = ˆ τ ( t ) and ˆ B = ˆ τ (0). We re-mind that this transform relies on the condition that onlythe distribution of start points of the trajectories dependson the ensemble, but not the relaxation pathways them-selves [37]. While this does not hold for extensive vari-ables if the same extensive variable if imposed, this isgenerally the case for fluctuations of instantaneous in-tensive variables which we focus on here. Interestingly, asimilar approach based on Ref. [43] has been used for thefour-point dynamic susceptibility χ ( t ) comparing its de-cay at constant temperature and constant energy [15, 16]. C. Time dependence of stress fluctuations
Introduction.
We have seen in Sec. II A that the shearmodulus G eq may be obtained by measuring the staticstress fluctuations σ F = βV (cid:104) δ ˆ τ (cid:105) . As for any fluctua-tion measured along a trajectory [37, 40] one expects thestress fluctuations σ F ( t ) computed over a too short timewindow t to yield only a time-dependent lower bound tothe true asymptotic long-time limit σ F [44]. (This re-mains even true if as a second step one averages over in-dependent trajectories.) This may seriously restrict theuse of the stress-fluctuation formula, Eq. (11), as will beseen at the end of Sec. IV A below. It is thus importantthat for systems with time translational symmetry σ F ( t )can be rewritten as an integral over the stress autocorre-lation function C ( t ). Correlated trajectories.
Let us consider N (cid:29) x n ≡ √ βV ˆ τ ( n ) with n = 1 , . . . , N stored at equidistant time steps δt over a total time in-terval t = N δt . Using similar steps as for the calcu-lation of the radius of gyration of polymer chains (Sec.2.4 of Ref. [3]) one may rewrite the expectation value (cid:104) ( x n − (cid:104) x n (cid:105) ) (cid:105) of the shear stress fluctuations as σ F ( N ) = (cid:42) N N (cid:88) n,m =1 ( x n − x m ) (cid:43) (19)where the average is performed over different trajecto-ries. Defining a “mean-square displacement” g ( s ) ≡(cid:104) ( x m = n + s − x n ) (cid:105) this allows to rewrite Eq. (19) as σ F ( N ) = 1 N N (cid:88) s =1 ( N − s ) g ( s ) . (20)where the weight ( N − s ) stems from the finite tra-jectory length. Using the correlation function C ( s ) = (cid:104) x m = n + s x n (cid:105) one verifies that g ( s ) = 2( C (0) − C ( s )) [3].Since (cid:80) Ns =1 ( N − s ) ≈ N / σ F ( N ) = C (0) − N N (cid:88) s =1 ( N − s ) C ( s ) . (21)Using that C (0) = σ F and rewriting the discrete sum asa continuous time integral this yields σ F ( t ) = σ F − t (cid:90) t d s (1 − s/t ) C ( s ) (22)independent of whether γ or τ are imposed. It followsthat the stress-fluctuation formula G F ≡ µ A − σ F | γ maybe rewritten quite generally as [44] G F ( t ) = G eq + 2 t (cid:90) t d s (1 − s/t ) C ( s ) | γ (23)= 2 t (cid:90) t d s (1 − s/t ) G ( s ) (24)where we have used the key relation, Eq. (4), in thesecond step. For large times C ( t ) | γ → C ( t ) | γ becomes constant. As expected forgeneral finite-sampling time corrections for fluctuations[40], the second term in Eq. (23) vanishes thus extremelyslowly as 1 /t [47]. As seen from Eq. (24), G F ( t ) and G ( t ) are different in general albeit closely related. Theyhave the same asymptotic limits G F (0) = G (0) = µ A and G F ( ∞ ) = G ( ∞ ) = G eq . Intermediate plateau.
It is of some interest to con-sider briefly the case of model systems where the stressautocorrelation function C ( t ) | γ reveals a broad interme-diate plateau, C ( t ) | γ = G P , extending over several ordersof magnitude up to a time τ P . It is readily seen usingEq. (23) or Eq. (24) that G ( t ) = C ( t ) | τ ≈ G eq + G P ≈ G F ( t ) for t (cid:28) τ P , (25)i.e. G ( t ) and G F ( t ) may become identical and constantfor a finite time window. D. Digression: Self-assembled transient networks
While the present work focuses on permanent elasticnetworks let us mention that this study can be extendednaturally on self-assembled transient networks as hy-perbranched polymer chains with sticky end-groups [32]or microemulsions bridged by telechelic polymers [33].Such networks may be modeled using purely repulsive LJbeads representing the oil droplets of the microemulsionwhich are connected reversibly by ideal springs similarlyas in MC simulations of equilibrium polymers [48]. Thetopological rearrangement of the network may be done byrandomly choosing a spring with a scission-recombinationfrequency f and making a hopping attempt to reconnectit with other neighboring beads subject to a standardMetropolis criterion [40]. If one freezes an equilibratednetwork, i.e. if one sets f = 0, the network must behaveexactly as the permanent solids we focus on in this work,i.e. Eq. (4) should hold with a finite G eq ( f = 0). If oneconsiders a very small, but finite frequency f and verylong sampling times, one expects G F ( t ) to show an inter-mediate plateau G P up to the relaxation time of the net-work τ P ( f ) ∼ /f and to decay for larger times. ( τ P ( f )characterizes the time needed to restore the particle per-mutation symmetry of the liquid state.) The plateau G P of G ( t ) should be set by the modulus G eq ( f = 0) of thequenched network, since for small times t (cid:28) τ P ( f ) thescission-recombination events become irrelevant. Since G eq ( f ) = 0 for finite f , this implies according to Eq. (4)that G ( t ) = C ( t ) | τ = C ( t ) | γ for all times. Using nowEq. (23) and Eq. (25) these arguments suggest confirm-ing Ref. [11] that G ( t ) = C ( t ) | γ = G F ( t ) = G P = G eq ( f = 0) (26)for intermediate times t (cid:28) τ P ( f ). In this sense C ( t ) | γ may indeed measure a shear modulus. It is however not the shear modulus G eq ( f ) of the system computed at fi-nite f (which must vanish) but of the quenched referencenetwork at f = 0. We return after this digression tosolids formed by permanently connected springs. III. COMPUTATIONAL MODEL
To illustrate our key relation we present in Sec. IV nu-merical data obtained using a periodic two-dimensionalnetwork of ideal harmonic springs of interaction energyˆ H ex = 12 (cid:88) l K l ( r l − R l ) (27)with K l being the spring constant, R l the reference lengthand r l = | r i − r j | the length of spring l . The sum runsover all springs l between topologically connected vertices i and j of the network at positions r i and r j . Note thatthe mass of the particles is set to unity and LJ units areassumed throughout. As explained in detail in Ref. [23],our network has been constructed using the dynamicalmatrix M of a strongly polydisperse LJ bead glass com-prising N = 10 particles. (An experimentally more rel-evant example for such permanent networks is providedby endlinked or vulcanized polymer networks [2, 29, 30].)Prior to forming the network the latter bead system hadbeen quenched down to T = 0 using a constant quench-ing rate and imposing a relatively large normal pressure P = 2. This yields systems of number density ρ ≈ . L ≈ .
3. Since thenetwork topology is by construction permanently fixed ,the shear response G ( t ) must become finite for t → ∞ for all temperatures at variance to systems with plasticrearrangements as considered, e.g., in Ref. [11], or thetransient networks mentioned in Sec. II D. IV. COMPUTATIONAL RESULTSA. Static properties
Ground state values.
Following Refs. [13, 14] one maycompute the shear modulus G eq of the ground state ofthe model at T = 0 from the lowest non-trivial four-fold degenerated eigenfrequencies ω p associated to trans-verse eigenmodes. (The running index p increases withfrequency.) Such eigenmodes can be determined by nu-merical diagonalization of the dynamical matrix M bymeans of Lanczos’ method [39]. For planar transversemodes one expects from continuum theory [1] that ω p = 2 πc T L (cid:112) n + m with c T = (cid:113) G eq /ρ (28)being the transverse wave velocity and n, m = 0 , , . . . two quantum numbers. One thus obtains for p = 3 , , , G eq ≈
16. By applying an affine shear strain to thesystem or by using Eq. (9) one determines an affine shearelasticity µ A ≈
34. In turn this implies a shear stressfluctuation σ F | γ = µ A − G eq ≈
18, i.e. about half of theenergy implied by an affine strain is relaxed by non-affinedisplacements as discussed in more detail in Ref. [14].
Static properties at small finite temperatures.
We fo-cus below on systems with a finite temperature T = 10 − .Since this temperature is rather small, one expects allstatic properties such as the pressure P or the elasticmodulus G eq to be barely changed. As we have checkedcomparing various methods one confirms indeed that P ≈ P ex ≈ µ A ≈ σ F | γ ≈
18 and G eq ≈
16 andthe same applies to all small temperatures T (cid:28) Convergence of stress-fluctuation formula.
How theshear modulus is obtained using the stress-fluctuatingformula, Eq. (11), can be seen from Fig. 2 where wepresent data obtained by standard velocity-Verlet MDsimulation [37, 38] using a (rather cautious) time step δt MD = 10 − . The temperature T = 10 − is imposedby means of a Langevin thermostat with a large frictionconstant ζ = 5. We have used here one long productionrun over a time t run = 10 for one equilibrated start con-figuration. Various properties, such as the instantaneous FIG. 2: Affine shear elasticity µ A ( t ), shear stress fluctua-tions σ F ( t ) | γ and their difference G F ( t ) = µ A − σ F ( t ) | γ as afunction of the measurement time t for one network of idealsprings in two dimensions (NV γ T-ensemble). The data havebeen sampled by MD simulation with time step δt MD = 10 − ,temperature T = 10 − and friction constant ζ = 5. The solidlines present a consistency check of σ F ( t ) | γ and G F ( t ) with thecorrelation function C ( t ) | γ confirming Eq. (22) and Eq. (23). values of the shear stress ˆ τ ( t ) or the affine shear elasticityˆ µ A ( t ), Eq. (10), have been written down at equidistanttime steps δt = 10 − . The data correspond to averagestaken first over a given time interval [ t , t = t + t ],i.e. using 1 + t/δt entries, and taking then in a secondstep gliding averages over all times t possible for t [37].(Naturally, the error bars thus increase somewhat with t .)The horizontal axis in Fig. 2 indicates the interval length t . As one expects, the simple average µ A ( t ) becomes im-mediately constant (filled squares), i.e. µ A ( t ) = µ A ≈ σ F ( t ) | γ are seen to in-crease monotonously from zero to the asymptotic plateau σ F | γ ≈
18. Interestingly, this plateau is only reached forsurprisingly large times t (cid:29) . The stress-fluctuationestimate G F ( t ) ≡ µ A − σ F ( t ) | γ (diamonds) of the shearmodulus G eq decreases thus monotonously from µ A toits large- t limit G eq ≈
16 indicated by the bold horizon-tal line. A too short production run thus leads to anoverestimation of G eq . The two thin solid lines present aconsistency check for σ F ( t ) | γ and G F ( t ) integrating theshear stress autocorrelation function C ( t ) | γ as suggestedby Eq. (22) and Eq. (23). The slow convergence of thestress-fluctuation relation G F ( t ) noted in Refs. [19, 23]can thus be traced back to the sluggish 1 /t -decay of thesecond term in Eq. (23). We turn now to the descriptionof C ( t ) | γ and other correlation functions. FIG. 3: Stress relaxation modulus G ( t ) and stress autocorre-lation functions C ( t ) | τ (triangles) and C ( t ) | γ (squares) sam-pled by MD simulation. An affine strain γ = 10 − is appliedto determine G ( t ). The filled triangles correspond to C ( t ) | τ computed using one single trajectory with a slow MC shear-barostat with δγ max = 10 − , the crosses to a too large value δγ max = 10 − for which C ( t ) | τ is seen to decay rapidly dueto additional relaxation pathways. B. Computational test of key relation
Introduction.
Having shown in Fig. 2 how a finiteshear modulus G eq ≈
16 may be determined in theNV γ T-ensemble using the stress fluctuation formula,Eq. (11), we now demonstrate numerically our key rela-tion, Eq. (4), by comparing the explicitly computed out-of-equilibrium stress relaxation modulus G ( t ) with theequilibrium autocorrelation functions C ( t ) | γ and C ( t ) | τ .As before we show first in Fig. 3 data obtained by MDsimulations using a high friction constant ζ = 5, whichsimplifies the data by enforcing a monotonic decay of thecorrelations. We discuss then results obtained using dif-ferent friction constants and computational schemes. Stress relaxation and autocorrelation functions.
Thestress relaxation modulus G ( t ) presented in Fig. 3 hasbeen computed from the shear stress increment δ ˆ τ ( t )measured after an affine shear strain γ = 10 − was im-posed at t = 0 [41]. We average over 10 runs start-ing from independent reference configurations at t = 0 − .The shear stress relaxation modulus G ( t ) decreases (dueto the strong damping) monotonously from G (0 + ) = µ A to a finite G eq . In contrast to this C ( t ) | γ (open squares)decays from σ F | γ to zero. Confirming Eq. (4), the verti-cally shifted autocorrelation function C ( t ) | γ + G eq (filledsquares) is seen to collapse onto G ( t ). The autocorrela-tion function C ( t ) | τ (open triangles) has been obtainedby preparing first an NV τ T-ensemble of mean stress τ = 0 containing 10 independent start configurations.We sample ˜ C ( t ) | γ and (cid:104) ˆ τ (cid:105)| γ for each configuration keep-ing γ constant and average then over all configurations,Eq. (12). Confirming Eq. (15) we observe G ( t ) ≈ C ( t ) | τ . FIG. 4: Stress relaxation modulus G ( t ) (open symbols) andrescaled autocorrelation function C ( t ) | γ + G eq (filled symbols)for MD simulations for several friction constants ζ , BD simu-lations with ζ = 1 and local MC moves. Shifting horizontallythe MC data by a factor 1 / Keeping the shear-barostat switched on.
As shown bythe small filled triangles in Fig. 3, the same result is alsoobtained by sampling ˆ τ ( t ) for every δt = 10 − using anextremely slow shear-barostat for one single trajectoryup to t = 10 . This large time is needed for a sufficientensemble sampling. Otherwise, σ F ( t ) | τ would remain be-low its asymptotic large- t limit µ A = σ F | τ [44]. We haveused here a hyprid MD-MC scheme where after everyMD step a Metropolis MC attempt was made to changethe metric [23, 41] by a small amount | δγ | ≤ δγ max =10 − . Additional (non-universal) relaxation pathwaysbecome important if γ ( t ) changes too strongly. If in-stead δγ max = 10 − is used (all other parameters keptconstant) this naturally leads to a rapid decay of C ( t ) | τ (crosses). The conceptionally important point is herethat in the limit of sufficiently small δγ max Eq. (3) holds.The quenched NV τ T-ensemble average C ( t ) | τ (open tri-angles) is then obtained without completely switchingoff the shear-barostat. The detailed description of thepresumably non-universal scaling of the additional relax-ation pathways at larger δγ max is of course also of inter-est. This should be addressed in future work. Different numerical schemes.
The scaling collapse of G ( t ) and C ( t ) | γ + G eq has been also obtained for differ-ent temperatures T (not shown) and friction constants ζ as may be seen from Fig. 4. As one expects, the MDdata decay more rapidly with decreasing ζ and revealanti-correlations and oscillations for the lowest ζ probed.Also included in Fig. 4 is data obtained by (overdamped)BD simulations with a friction constant ζ = 1 and MCsimulations with local monomer jump attempts unifor-mally distributed in a disk of radius 0 .
01 [37]. Both datasets for each simulation type are again found to collapse.Note that it is possible to collapse additionally the BDand MC data by shifting the MC data horizontally.
Gauge freedom for the instantaneous stress.
A techni-cal point should finally be mentioned. While for our MDsimulations the instantaneous shear stress ˆ τ = ˆ τ id + ˆ τ ex comprises both an ideal contribution ˆ τ id and an excesscontribution ˆ τ ex and correspondingly µ A = µ A , id + µ A , ex and σ F | γ = σ F , id | γ + σ F , ex | γ take ideal contributions µ A , id = σ F , id | γ = P id , this is not possible for BD andMC simulations. (Note that for the low temperatureconsidered contributions of order P id ≈ − are in anycase negligible.) Within the so-called “MC-gauge” [24],we thus set ˆ τ id ≡ µ A , σ F | γ , C ( t ) | τ , C ( t ) | γ and G ( t ) are measured consistently . V. CONCLUSION
Main results.
Focusing on permanent isotropic net-works in thermal equilibrium (Sec. III) we have revis-ited theoretically (Sec. II B) and numerically (Sec. IV B)the linear-response relation between the shear stress re-laxation modulus G ( t ) and the shear stress autocorre-lation functions C ( t ) | γ and C ( t ) | τ computed, respec-tively, at imposed strain γ and mean stress τ . It hasbeen demonstrated that according to Eq. (4) or Eq. (18) G ( t ) = C ( t ) | τ and C ( t ) | γ must become different in thesolid limit for G eq >
0. While G ( t ) may be determinednumerically directly from C ( t ) | τ using either a quenchedNV τ T-ensemble or an asymptotically slow shear-barostatfor which Eq. (3) holds (Fig. 3), this is not possible alonefrom C ( t ) | γ . Digression.
More briefly, we have commented on self-assembled transient elastic networks characterized bya scission-recombination frequency f for the springs(Sec. II D). For a finite, but small frequency f the shearmodulus G eq ( f ) must vanish for long sampling times.Following Ref. [11] we have argued that G ( t ) = C ( t ) | τ = C ( t ) | γ = G F ( t ) should reveal an intermediate plateau G P and that this plateau is set by the finite shear modulusof the quenched network, G P = G eq ( f = 0). Discussion.
More generally, it is obviously often help-ful to describe an observed intermediate plateau or shoul- der of G ( t ) in terms of a phenomenological shear modu-lus G P of a dynamical model, such as the Maxwell modelfor viscoelastic fluids or the reptation model of entangledpolymer melts [2, 3]. However, such a model allowing thetheoretical interpretation of the data should not be con-fused with the proper measurement procedure and themodel parameter G P should not be identified with thethermodynamic equilibrium modulus G eq of the system.Note that the shear modulus G eq of a Maxwell fluid or alinear polymer melt must vanish while the phenomeno-logical parameter G P describing the short or intermediatetime stress response is finite. Since in this sense differ-ent operational “static” and “dynamical” definitions ofthe shear modulus are used for describing glass-formingliquids close to the glass transition [17, 18, 20, 21, 25],this may explain why qualitatively different temperaturedependences (cusp singularity [17, 25] vs. finite jump[5, 18, 21, 26]) have been predicted. Hence, while our re-cent attempts to determine G eq ( T ) for two glass-formingmodel systems [23] are consistent with a continuous cusp,this is not necessarily in contradiction with a jump sin-gularity for G P ( T ) determined from C ( t ) | γ [21, 26]. Outlook.
It should be noted that generalizing Eq. (4)one obtains readily that M ( t ) = C ( t ) | I = C ( t ) | X + M eq for t > M ( t ) being the relaxation modulus of an intensivevariable I , M eq = ∂I/∂X the associated equilibriummodulus and C ( t ) = βV (cid:104) δ ˆ I ( t ) δ ˆ I (0) (cid:105) the correspondingautocorrelation function for any (continuous) intensive variable ˆ I ( t ). We note finally that we are currently simu-lating transient elastic networks formed by dense, purelyrepulsive beads which are reversibly connected by har-monic springs. The preliminary results support Eq. (26)suggested in Sec. II D, i.e. it is seen that G ( t ), C ( t ) | γ and G F ( t ) approach with decreasing, but finite scission-recombination frequency f , i.e. G eq ( f ) = 0, an interme-diate plateau given by the shear modulus G eq ( f = 0) ofthe quenched reference network [49]. Acknowledgments
H.X. thanks the IRTG Soft Matter for financial sup-port. We are indebted to O. Benzerara and J. Farago(both ICS, Strasbourg) and M. Fuchs (Konstanz) forhelpful discussions. [1] L. D. Landau and E. M. Lifshitz,
Theory of Elasticity (Pergamon Press, 1959).[2] M. Rubinstein and R. Colby,
Polymer Physics (OxfordUniversity Press, Oxford, 2003).[3] M. Doi and S. F. Edwards,
The Theory of Polymer Dy- namics (Clarendon Press, Oxford, 1986).[4] J. Hansen and I. McDonald,
Theory of simple liquids (Academic Press, New York, 2006), 3nd edition.[5] W. G¨otze,
Complex Dynamics of Glass-Forming Liquids:A Mode-Coupling Theory (Oxford University Press, Ox- ford, 2009).[6] S. Alexander, Physics Reports , 65 (1998).[7] T. Witten and P. A. Pincus,
Structured Fluids: Polymers,Colloids, Surfactants (Oxford University Press, Oxford,2004).[8] H. B. Callen,
Thermodynamics and an Introduction toThermostatistics (Wiley, New York, 1985).[9] D. Chandler,
Introduction to Modern Statistical Mechan-ics (Oxford University Press, New York, 1987).[10] P. M. Chaikin and T. C. Lubensky,
Principles of con-densed matter physics (Cambridge University Press,1995).[11] F. Sausset, G. Biroli, and J. Kurchan, J. Stat. Phys. ,718 (2010).[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] B. Schnell, H. Meyer, C. Fond, J. P. Wittmer, andJ. Baschnagel, Eur. Phys. J. E , 97 (2011).[20] H. Yoshino, J. Chem. Phys. , 214108 (2012).[21] C. Klix, F. Ebert, F. Weysser, M. Fuchs, G. Maret, andP. Keim, Phys. Rev. Lett. , 178301 (2012).[22] H. Xu, J. Wittmer, P. Poli´nska, and J. Baschnagel, Phys.Rev. E , 046705 (2012).[23] J. P. Wittmer, H. Xu, P. Poli´nska, F. Weysser, andJ. Baschnagel, J. Chem. Phys. , 12A533 (2013).[24] J. P. Wittmer, H. Xu, P. Poli´nska, C. Gillig, J. Hellferich,F. Weysser, and J. Baschnagel, Eur. Phys. J. E , 131(2013).[25] A. Zaccone and E. Terentjev, Phys. Rev. Lett. ,178002 (2013).[26] M. Ozawa, T. Kuroiwa, and A. Ikeda, Phys. Rev. Lett. , 205701 (2012).[27] H. Mizuno, S. Mossa, and J.-L. Barrat, Phys. Rev. E ,042306 (2013).[28] E. del Gado and W. Kob, J. Non-Newtonian Fluid Mech. , 28 (2008).[29] E. Duering, K. Kremer, and G. S. Grest, Phys. Rev. Lett. , 67 (1991).[30] E. Duering, K. Kremer, and G. S. Grest, J. Chem. Phys. , 101 (1994).[31] S. Ulrich, X. Mao, P. Goldbart, and A. Zippelius, Euro-physics Lett. , 677 (2006).[32] C. Tonhauser, D. Wilms, Y. Korth, H. Frey, andC. Friedrich, Macromolecular Rapid Comm. , 2127(2010).[33] A. Zilman, J. Kieffer, F. Molino, G. Porte, and S. A.Safran, Phys. Rev. Lett. , 2003 (2003).[34] Numerically, it might be better to measure the average shear stress τ ( γ ), to fit a spline to this equation of stateand to determine G eq ( γ ) from the fit by taking the deriva-tive with respect to the shear strain γ .[35] M. Parrinello and A. Rahman, J. Chem. Phys. , 2662(1982).[36] J. F. Lutsko, J. Appl. Phys , 2991 (1989).[37] M. Allen and D. Tildesley, Computer Simulation of Liq-uids (Oxford University Press, Oxford, 1994).[38] D. Frenkel and B. Smit,
Understanding Molecular Sim-ulation – From Algorithms to Applications (AcademicPress, San Diego, 2002), 2nd edition.[39] J. Thijssen,
Computational Physics (Cambridge Univer-sity Press, Cambridge, 1999).[40] D. P. Landau and K. Binder,
A Guide to Monte CarloSimulations in Statistical Physics (Cambridge UniversityPress, Cambridge, 2000).[41] A canonical affine transformation consists in changingboth the particle coordinates r ⇒ h r and the conjugatedvelocities v → h − v using a linear matrix h [35, 36].For a uniform shear in two dimensions this amounts to r x → r x + γr y and v x → v x − γv y for the transformation ofthe x -components of the particle positions and velocities.[42] M. S. Green, Phys. Rev. , 829 (1960).[43] J. L. Lebowitz, J. K. Percus, and L. Verlet, Phys. Rev. , 250 (1967).[44] The notation σ F , Eq. (6), or G F , Eq. (11), without timeargument refer to static thermodynamic properties deter-mined using asymptotically long sampling times. In orderto avoid additional notations we write σ F ( t ) or G F ( t ) toindicate that these properties have been determined us-ing a finite time window t . Please note also that G F ( t )should not be confused with the response function G ( t )albeit both properties are related as discussed in Sec. II C.[45] In general there is a considerable freedom for definingan instantaneous intensive variable ˆ I as long as its aver-age I = (cid:104) ˆ I (cid:105) does not change. The experimentally naturalchoice for ˆ I = ˆ τ is given by the instantaneous forcesacting on the shear cell boundaries. For theory and sim-ulation, ˆ τ ≡ ˆ H (cid:48) ( γV ), i.e. the energy change with respectto an assumed affine strain, is the most common choicedue to Eq. (8).[46] It is instructive to obtain G ( t ) = G eq + C ( t ) | γ directlyby rewriting the derivation of the FDT given in Ref. [3]for a strain γ switched on at t = 0. This shows that the G eq -term stems naturally from the residual finite shearstress of the strained system at equilibrium.[47] For an exponentially decaying stress-autocorrelationfunction C ( t ) = C (0) exp( − x ) with x = t/τ one ob-tains by integration σ F ( t ) /σ F = 1 − f Debye ( x ) with f Debye ( x ) ≡ − x ) − x ) /x being the Debye func-tion well-known in polymer science [3]. For large times x (cid:29) σ F ( t ) /σ F = 1 − /x .[48] J. P. Wittmer, A. Milchev, and M. E. Cates, J. Chem.Phys. , 834 (1998).[49] All simple means are independent of f , especially µ A .Interestingly, this is different for the shear stress fluctua-tions σ F | γ , i.e. a well-defined static four-point correlationfunction, which shows a singular behavior for asymptoti-cally long sampling times. While σ F | γ = µ A − G eq ( f = 0)for the quenched network, one obtains σ F | γ = σ F | τ = µ A for finite f (consistent with G eqeq