Shear-strain and shear-stress fluctuations in generalized Gaussian ensemble simulations of isotropic elastic networks
SShear-strain and shear-stress fluctuations in generalizedGaussian ensemble simulations of isotropic elastic networks
J.P. Wittmer, ∗ I. Kriuchevskyi, J. Baschnagel, and H. Xu 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: October 7, 2018)Shear-strain and shear-stress correlations in isotropic elastic bodies are investigated both theoreti-cally and numerically at either imposed mean shear-stress τ ( λ = 0) or shear-strain γ ( λ = 1) and formore general values of a dimensionless parameter λ characterizing the generalized Gaussian ensem-ble. It allows to tune the strain fluctuations µ γγ ≡ βV (cid:10) δ ˆ γ (cid:11) = (1 − λ ) /G eq with β being the inversetemperature, V the volume, ˆ γ the instantaneous strain and G eq the equilibrium shear modulus. Fo-cusing on spring networks in two dimensions we show, e.g., for the stress fluctuations µ ττ ≡ βV (cid:10) δ ˆ τ (cid:11) (ˆ τ being the instantaneous stress) that µ ττ | λ = µ A − λG eq with µ A = µ ττ | λ =0 being the affine shear-elasticity. For the stress autocorrelation function C ττ ( t ) ≡ βV (cid:104) δ ˆ τ ( t ) δ ˆ τ (0) (cid:105) this result is then seen(assuming a sufficiently slow shear-stress barostat) to generalize to C ττ ( t ) | λ = G ( t ) − λG eq with G ( t ) = C ττ ( t ) | λ =0 being the shear-stress relaxation modulus. PACS numbers: 05.70.-a,05.20.Gg,05.10.Ln,65.20.-w
I. INTRODUCTION
Background.
A “simple average” A = (cid:104) ˆ A (cid:105) of an ob-servable A [1] does not depend on which thermodynamicensemble it is measured in, at least not if the systemis sufficiently large and if each ensemble samples indeedthe same thermodynamic state point [1–6]. An exam-ple for such a simple average is the affine shear-elasticity µ A characterizing the energy change of an affinely shear-strained elastic body [7–10] as properly defined below insect. III. As one may verify numerically [8], it is irrele-vant whether one computes µ A in the NV γ T-ensembleat constant particle number N , volume V , shear strain γ and temperature T ≡ /β or in the conjugated NV τ T-ensemble where the strain is allowed to fluctuate sub-jected to the constraint that the internal mean shearstress τ is imposed by the external applied shear stress τ ext [3]. In contrast to this, the fluctuation (cid:104) δ ˆ Aδ ˆ B (cid:105) oftwo observables A and B may depend on whether an ex-tensive variable X or its conjugated intensive variable I is imposed [1–3]. Focusing in the present work on shear-strained isotropic elastic networks, as sketched in panel(a) of fig. 1, the relevant extensive variable is the rescaledshear strain X = V γ and the conjugated intensive vari-able the shear stress I = τ . Using the Lebowitz-Percus-Verlet (LPV) transform [2] it is seen for the (rescaled)mean-squared fluctuation µ ττ ≡ βV (cid:10) δ ˆ τ (cid:11) of the instan-taneous shear stress ˆ τ that [8, 9] µ ττ | λ =0 = µ ττ | λ =1 + G eq (1)with G eq being the equilibrium shear modulus. For laterconvenience the NV τ T-ensemble is indicated by “ λ = 0” ∗ Electronic address: [email protected] -20 -10 0 10 20 30 40 50 G ext -4-3-2-101 λ NV τ T NV γ T G eq =16.3 -G eq i n s t a b ili t y (c) (a) (b) γ γ U /V ext γ ext xy L τ ext ex t − τ γ springexternal FIG. 1: Sketch of addressed problem: (a)
Shear strain ex-periment in the ( x, y )-plane with shear strain γ and appliedexternal shear stress τ ext . (b) External potential U ext ( γ ) /V according to eq. (3) with the bold line indicating the externalspring of spring constant G ext . (c) Dimensionless parameter λ = G ext / ( G eq + G ext ) comparing the external spring con-stant G ext and the equilibrium shear modulus G eq . We set G eq = 16 . and the NV γ T-ensemble by “ λ = 1”. Generalizing eq. (1)in the time-domain it has been shown for the stress-stresscorrelation function C ττ ( t ) ≡ βV (cid:104) δ ˆ τ ( t ) δ ˆ τ (0) (cid:105) that [9, 10] C ττ ( t ) | λ =0 = C ττ ( t ) | λ =1 + G eq . (2)This assumes that the shear-barostat needed to samplethe NV τ T-ensemble must act very slowly on the time-scales used to probe C ττ ( t ) | λ =0 . This may be realizedequivalently by averaging over an ensemble of configu-rations with quenched strains ˆ γ distributed according tothe NV τ T-ensemble [9]. It is due to the ergodicity break-ing generated by the averaging over the quenched ensem- a r X i v : . [ c ond - m a t . s t a t - m ec h ] A ug ble that C ττ ( t ) | λ =0 → G eq stays finite for t → ∞ , while C ττ ( t ) | λ =1 must vanish for large times (as one commonlyexpects for all correlation functions in ergodic systems[11]). By integration by parts it is seen that C ττ ( t ) | λ =0 is equivalent to the experimentally important shear re-laxation modulus G ( t ) [9]. The transform eq. (2) thusimplies G ( t ) = G eq + C ττ ( t ) | λ =1 , i.e. the relaxation mod-ulus may be measured using the equilibrium modulus G eq and the correlation function C ττ ( t ) | λ =1 both determinedin the NV γ T-ensemble [9].
Generalized Gaussian ensemble.
As sketched in panel(b) of fig. 1, it is straightforward to interpolate betweenthe NV γ T-ensemble and the NV τ T-ensemble by impos-ing an external field U ext (ˆ γ ) /V = − τ ext (ˆ γ − γ ext ) + 12 G ext (ˆ γ − γ ext ) (3)with G ext being the spring constant of the external har-monic spring. The standard NV τ T-ensemble at τ = τ ext is recovered by setting G ext = 0. Note that our approachis conceptually similar to the so-called “Gaussian ensem-ble” proposed some years ago by Hetherington and others[12, 13] generalizing the Boltzmann weight of the canon-ical ensemble by an exponential factor U ext ( ˆ E ) ∝ ˆ E ofthe instantaneous energy ˆ E . A similar external spring po-tential has been also used in the recent “elastic bath” ap-proach by Workum and de Pablo [14]. Choosing the refer-ence strain γ ext equal to the mean strain γ of the NV τ T-system at vanishing shear stress τ = (cid:104) ˆ τ (cid:105) = τ ext ≡ same ther-modynamic state, i.e. all first derivatives of the energy orthe free energy [3] and all simple averages are identical.As sketched in panel (c) of fig. 1, G ext is not necessarilypositive, reducing the strain fluctuations, but may evenbecome negative, which thus amplifies the fluctuations.It has been argued [14] that this may allow a more con-venient determination of the elastic modulus. Since theexternal spring is parallel to the system, the combinedsystem and external device have an effective spring con-stant G eff = G eq + G ext . Defining the dimensionless pa-rameter λ ≡ G ext / ( G eq + G ext ) the strain fluctuations ofthe combined system are thus given by µ γγ ≡ βV (cid:10) δ ˆ γ (cid:11) = 1 /G eff = (1 − λ ) /G eq , (4)i.e. must vanish linearly with λ . NV τ T-ensemble statis-tics is expected for λ →
0, while NV γ T-statistics shouldbecome relevant in the opposite limit for G ext → ∞ , i.e. λ →
1. The system must become unstable in the limit G ext → − G eq , i.e. λ → −∞ . Key results.
The aim of the present study is to gen-eralize the relations for static fluctuations, eq. (1), anddynamical correlation functions, eq. (2), to the more gen-eral transformations between Gaussian ensembles char-acterized by the continuous parameter λ ≤
1. In thisway we want to make manifest that these transformation relations are generated by the continuous change of theconstraint imposed on the fluctuations of the extensivevariable. Focusing on spring networks in two dimensionswe show, e.g., for the stress-stress fluctuations µ ττ in dif-ferent λ ensembles that µ ττ | λ = µ ττ | λ =0 − λG eq (5)with µ ττ | λ =0 being given by the affine shear-elasticity µ A mentioned above. Assuming a very slowly acting shear-barostat, which is irrelevant for the system evolution forshort times t (cid:28) τ (cid:63) ( λ ), the above result is then seen togeneralize in the time-domain for the stress-stress corre-lation function C ττ ( t ) | λ = C ττ ( t ) | λ =0 − λG eq for t (cid:28) τ (cid:63) ( λ ) . (6)Since G ( t ) = C ττ ( t ) | λ =0 , eq. (6) allows the determina-tion of the relaxation modulus G ( t ) from the equilibriummodulus G eq and the stress-stress correlation function C ττ ( t ) | λ for any λ and t (cid:28) τ (cid:63) ( λ ). The upper time limit τ (cid:63) ( λ ) is seen to be set by the diffusion time T γγ of theinstantaneous strain ˆ γ ( t ) over the typical strain fluctua-tions for λ . As a consequence from eq. (5) and eq. (6)it appears that it is the equilibrium modulus G eq of thesystem, which generates both transformsdd λ µ ττ | λ = dd λ C ττ ( t ) | λ = − G eq (7)with µ ττ | λ = µ A and C ττ ( t ) | λ = G ( t ) for λ = 0. Similarrelations, allowing also to determine the creep compliance J ( t ) [16, 17], are obtained for strain-strain and strain-stress correlations functions. Outline.
The stated key relations eq. (5) and eq. (6)are justified theoretically in sect. II. The static fluctua-tions are discussed in sect. II A before we address the dy-namical correlation functions in sects. II B, II C and II Dand the macroscopic linear response in sect. II E. Algo-rithmic details are given in sect. III where the specificmodel system considered is introduced in sect. III A. Thecanonical affine plane shear transformations used in thiswork are specified in sect. III B, the instantaneous shearstress ˆ τ and the instantaneous affine shear-elasticity ˆ µ A in sect. III C. The zero-temperature ground state prop-erties of the two-dimensional elastic network are summa-rized in sect. III D. Some technical details related to thefinite-temperature simulation of the Gaussian λ -ensembleusing a Metropolis Monte Carlo (MC) scheme as a func-tion of the maximum strain displacement rate κ , the sec-ond key operational parameter of this study, are given insect. III E. Section IV presents the numerical results ob-tained by molecular dynamics (MD) simulations at onefinite (albeit small) temperature T and one relativelyhigh value of the friction constant ζ of the Langevinthermostat used [1]. The relevant static properties aredescribed in sect. IV A before we turn to the dynamicalstrain-strain (sect. IV B), strain-stress (sect. IV C) andstress-stress (sect. IV D) correlation functions. Our workis summarized in sect. V. II. THEORETICAL CONSIDERATIONSA. Static fluctuations
Lebowitz-Percus-Verlet transform.
We begin by re-stating the LPV transform in a convenient form assumingthat the relevant extensive variable is the (rescaled) shearstrain X = V γ and the conjugated intensive variable theshear stress I = τ . Following Lebowitz et al. [1, 2], oneverifies (see also sect. II.A of ref. [8]) that to leadingorder (cid:68) δ ˆ Aδ ˆ B (cid:69)(cid:12)(cid:12)(cid:12) λ =0 = (cid:68) δ ˆ Aδ ˆ B (cid:69)(cid:12)(cid:12)(cid:12) λ =1 + G eq βV ∂A∂τ ∂B∂τ (8)for the transformation of (cid:104) δ ˆ Aδ ˆ B (cid:105) with δ ˆ A ≡ ˆ A − (cid:104) ˆ A (cid:105) and δ ˆ B ≡ ˆ B − (cid:104) ˆ B (cid:105) from the NV γ T-ensemble ( λ = 1) tothe NV τ T-ensemble ( λ = 0). The more general transfor-mation between arbitrary λ -ensembles can be found byreworking the saddle-point approximation of ref. [2] tak-ing into account that the fluctuations around the peak ofthe distribution of the extensive variable is now not setby the modulus G eq of the system but by the effectivemodulus G eff = G eq + G ext . How this may be done canbe seen in sect. 2.5 of ref. [18] for the volume X = V asextensive variable and the (negative) pressure I = − P as intensive variable. Rewriting the latter result to thepresent case we get the generalized LPV transform (cid:68) δ ˆ Aδ ˆ B (cid:69)(cid:12)(cid:12)(cid:12) λ = (cid:68) δ ˆ Aδ ˆ B (cid:69)(cid:12)(cid:12)(cid:12) λ =1 + (1 − λ ) G eq βV ∂A∂τ ∂B∂τ (9)which reduces for λ = 0 to eq. (8). Strain-strain fluctuations.
Let us check the ensem-ble dependence of the rescaled strain-strain fluctuations µ γγ ≡ βV (cid:10) δ ˆ γ (cid:11) [19]. The generalized Gaussian ensem-ble corresponds to replacing the shear modulus G eq ofthe system by the effective modulus G eff = G eq + G ext .As already stated above, eq. (4), this leads to [20] µ γγ = 1 /G eff = (1 − λ ) /G eq . (10)One verifies readily that the postulated LPV transformfor general λ , eq. (9), is consistent with this result. Tosee this one sets A = B = V γ and uses the fact that thestrain fluctuations must vanish in the NV γ T-ensemble,i.e. that the first term on the right hand-side of eq. (9)must vanish.
Strain-stress fluctuations.
Setting A = V γ and B = τ and using again that the NV γ T-term in eq. (9) mustvanish, it is seen from the LPV transform that the strain-stress fluctuations µ γτ ≡ βV (cid:104) δ ˆ γδ ˆ τ (cid:105) [19] should scale as µ γτ = 1 − λ. (11)This result can be also obtained directly by replacing inthe definition of µ γτ the fluctuation δ ˆ τ by G eq δ ˆ γ andusing then the strain-strain relation, eq. (10). Stress-stress fluctuations.
We turn to the most im-portant stress-stress fluctuation µ ττ ≡ βV (cid:10) δ ˆ τ (cid:11) [19]. Weset now A = B = τ in the LPV transform. Since thestress fluctuations do not vanish in the NV γ T-ensemblethe corresponding term must now be included. Thisyields [20] µ ττ | λ = µ ττ | λ =1 + (1 − λ ) G eq . (12)Interestingly, the contribution (1 − λ ) G eq may be rewrit-ten using the notation (cid:104) f (ˆ γ ) (cid:105) γ ≡ (cid:82) dˆ γ p (ˆ γ ; λ ) f (ˆ γ ) forthe strain-average of a property f (ˆ γ ) with p (ˆ γ ; λ ) be-ing the normalized strain-distribution for the consid-ered λ -ensemble. The total mean stress τ of the en-semble is thus given by τ = (cid:104) τ (ˆ γ ) (cid:105) γ with τ (ˆ γ ) beingthe average shear-stress of all configurations of shear-strain ˆ γ . Using that p (ˆ γ ; λ ) is a Gaussian and that δτ (ˆ γ ) ≡ τ (ˆ γ ) − τ = G eq (ˆ γ − γ ) with γ being the maxi-mum of the distribution, it is then seen that(1 − λ ) G eq = G G eff = βV (cid:10) δτ (ˆ γ ) (cid:11) γ . (13)Using eq. (12) this implies in turn µ ττ | λ = µ ττ | λ =1 + βV (cid:10) δτ (ˆ γ ) (cid:11) γ (14)as one expects to lowest order from a standard saddle-point approximation. Stress-fluctuation formula for shear modulus.
Sub-stracting the transform for λ = 0 from eq. (12) confirmsthe key result eq. (5). As shown in ref. [10], µ ττ | λ =0 canbe reduced by integration by parts to the affine shear-elasticity µ A . Since the latter expression is a simple av-erage, i.e. the same value µ A is obtained in any ensemble,eq. (5) may be further simplified as [20] µ ττ | λ = µ A − λG eq . (15)It is thus sufficient to compute the material constants µ A and G eq in any ensemble to obtain the stress-stressfluctuations µ ττ as a function of λ . Note that the spe-cial case for λ = 1 corresponds to the well-known stress-fluctuation formula [4, 8, 21, 22] G eq = µ A − µ ττ | λ =1 (16)used in numerous numerical studies to compute themodulus G eq conveniently in the NV γ T-ensemble [8–10, 18, 22–26].
B. Dynamics: Definitions and general relations
Introduction.
We define now several dynamical ob-servables of interest, which will be investigated numer-ically in sect. IV, and remind some well-known gen-eral relations [1, 11, 27]. Time translational symmetry( t ↔ t + δt ) and time reversal symmetry ( t ↔ − t ) are as-sumed. We use ˆ a ( t ) and ˆ b ( t ) for either the instantaneousshear strain ˆ γ or the shear stress ˆ τ , a ≡ (cid:104) ˆ a (cid:105) and b ≡ (cid:104) ˆ b (cid:105) for their thermodynamic averages and δ ˆ a ( t ) ≡ ˆ a ( t ) − a and δ ˆ b ( t ) ≡ ˆ b ( t ) − b for their time-dependent fluctuations. Mean-square displacements.
We define (generalized)mean-square displacements (MSD) by g ab ( t ) ≡ βV (cid:68) (ˆ a ( t ) − ˆ a (0))(ˆ b ( t ) − ˆ b (0)) (cid:69) (17)where the prefactor βV / g ab ( t ) = g ba ( t ), g ab ( t ) = g ab ( − t )as may be seen using the above-mentioned symmetries, g ab ( t ) → t → g ab ( t ) → µ ab ≡ βV (cid:68) δ ˆ a δ ˆ b (cid:69) for t (cid:29) T ab (18)with T ab being the characteristic time needed to reachthis thermodynamic limit [28]. Correlation functions.
Similarly, we define dynamiccorrelation functions by C ab ( t ) ≡ βV (cid:104) δ ˆ a ( t ) δ ˆ b (0) (cid:105) [19].Obviously, C ab (0) = µ ab , C ab ( t ) = C ab ( − t ) and C ab ( t ) = C ba ( t ). The latter identity is seen from [3, 11] (cid:68) ˆ a ( t )ˆ b (0) (cid:69) = (cid:68) ˆ a (0)ˆ b ( − t ) (cid:69) = (cid:68) ˆ b ( t )ˆ a (0) (cid:69) (19)where the time translational invariance is used in the firststep and the time reversal symmetry in the second step.Using again both symmetries one verifies that [11, 27] g ab ( t ) = C ab (0) − C ab ( t ) = µ ab − C ab ( t ) . (20)This implies that C ab ( t ) → t (cid:29) T ab [28].Albeit g ab ( t ) and C ab ( t ) thus contain the same informa-tion it will be sometimes better for theoretical or numer-ical reasons to focus on either g ab ( t ) or C ab ( t ). Since g ab ( t ) → t → g ab ( t ) to clarifythe power-law scaling of the correlations at short times. Finite-time dependence of fluctuation estimation.
Asdefined in sect. II A, µ ab is a thermodynamic ensembleaverage, hence, a time-independent property. In practice, µ ab may, however, often be estimated using a time series( a k , b k ) with a finite number n of more or less correlatedentries [1, 5]. It is supposed here that these series aresampled with a constant time interval of length δt , i.e. n corresponds to a time t = ( n − δt . With A k ≡ n (cid:80) nk =1 A k denoting such a finite average, one samplesthe n -dependent observable [19] µ ab ( n ) ≡ βV (cid:28) (ˆ a k − ˆ a k )(ˆ b k − ˆ b k ) (cid:29) (21)with (cid:104) . . . (cid:105) standing for an additional ensemble averageover different time series of n subsequent data points. Wehave µ ab ( n ) = 0 for n = 1 and µ ab ( n ) → µ ab for n → ∞ in an ergodic system. Interestingly, the detailed time-dependence of µ ab ( n ) can be obtained from a weightedintegral of the correlation function C ab ( t ) [5, 9]. To seethis we note first the identity(ˆ a k − ˆ a k )(ˆ b k − ˆ b k ) = 12 n n (cid:88) k,l =1 (ˆ a k − ˆ a l )(ˆ b k − ˆ b l ) (22) which can be verified by straightforward algebra. (Seesect. 2.4 of ref. [27].) Using time translational invariancethis implies µ ab ( n ) = 2 n n (cid:88) k =1 ( n − k ) g ab ( δt k ) (23)where the weight ( n − k ) stems from the finite length ofthe trajectory used. We change now to continuous timevariables, k → s = ( k − δt , and replace the discretesum by the time integral µ ab ( t ) = 2 t (cid:90) t d s (1 − s/t ) g ab ( s ) . (24)Using eq. (20) this yields the general relation1 − µ ab ( t ) µ ab = 2 t (cid:90) t d s (1 − s/t ) C ab ( s ) C ab (0) (25)which can be used to obtain µ ab ( t ) if the correlation func-tion C ab ( t ) is known. Relaxation time θ ab . Interestingly, the time integralin eq. (25) must become constant in all cases consideredbelow where C ab ( t ) vanishes ultimately for t → ∞ . (Thisapplies if a finite shear-barostat is switched on.) Definingthe characteristic relaxation time θ ab ≡ lim t →∞ (cid:90) t d s (1 − s/t ) C ab ( s ) µ ab , (26)this leads to1 − µ ab ( t ) /µ ab → θ ab /t for t → ∞ (27)as one expects for the Poisson statistics of uncorrelatedevents [5]. One may thus determine the relaxation time θ ab from the large-time asymptotics of µ ab ( t ) [29]. C. Dynamics: Additional model assumptions
Let us suppose that the MSD g ab ( t ) is diffusive forshort times, i.e. g ab ( t ) = D ab t/ t (cid:28) T ab with D ab being a diffusion constant. Matching this short timeregime with g ab ( t ) = µ ab for t (cid:29) T ab gives the possibilityto operationally define T ab as the crossover time T ab ≡ µ ab /D ab . (28)Different short-time dynamics may suggest of course adifferent operational definition. Let us further assumean exponentially decaying correlation function C ab ( t ) = µ ab exp( − x ) with x = t/T ab . Using eq. (20) this is seento be consistent with a diffusive short-time regime andeq. (28). It follows by integration from eq. (25) that1 − µ ab ( t ) µ ab = f Debye ( x ) ≡ x [exp( − x ) − x ] (29)using the Debye function well-known in polymer science[27]. For large reduced times x (cid:29) − µ ab ( t ) /µ ab → T ab /t , i.e. by comparison with eq. (27) wehave T ab = θ ab for an exponentially decaying correlationfunction.As we shall see in sect. IV D, it may also occur thatthe correlation function is more or less constant betweena local (barostat independent time) τ A up to a (barostatdependent) time τ (cid:63) , i.e. essentially C ab ( t ) ≈ cH ( τ (cid:63) − t )with c being a constant and H ( x ) the Heaviside function.It follows from eq. (26) that1 − µ ab ( t ) µ ab ≈ (cid:26) c for τ A (cid:28) t (cid:28) τ (cid:63) c τ (cid:63) /t for τ (cid:63) (cid:28) t, (30)i.e. θ ab ≈ τ (cid:63) c/ D. Dynamics: Ensemble effects
Introduction.
We have just stated various general re-lations applying to all ensembles. These relations do,however, not allow to relate the dynamical correlationsbetween different λ . Since some ensembles are more read-ily computed than others, it would be useful to have atransformation relation such as the LPV transform forstatic fluctuations considered in sect. II A. We emphasizethat the correlation functions depend in general on theensemble, i.e. the parameter λ , and on the dynamics ofthe shear-barostat. Since in the end we want to describethe intrinsic relaxation dynamics of the system, it is suf-ficient to focus on the limit where the barostat becomesvery slow such that it becomes essentially irrelevant forthe system evolution below an upper time τ (cid:63) . We con-sider the limit where this time τ (cid:63) is be much larger thanany intrinsic relaxation time of the system. Since we wantstill to consider meaningful thermal averages with respectto the chosen ensemble, we need either trajectories of atime interval t traj much larger than τ (cid:63) to allow the fullsampling of the phase space or, equivalently, we need toaverage over independent start configurations represent-ing the ensemble. Under these two assumptions usefultransformation relations can be formulated by reworkingthe generalization of the LPV transform for dynamicalcorrelations replacing the system modulus G eq by the ef-fective modulus G eff . Being beyond the scope of thispaper, we proceed by postulating the central scaling re-lation for the MSD g ab ( t ) and argue briefly why this rela-tion is natural. We discuss then in turn the strain-straincorrelations C γγ ( t ), the strain-stress correlations C γτ ( t )and the stress-stress correlations C ττ ( t ) for different λ .Alternative, more direct ways to derive the relations arementioned. Scaling relations.
We postulate that under the twoassumptions made above the MSD g ab ( t ) does not dependon the ensemble g ab ( t ) ∼ λ for t (cid:28) τ (cid:63) ( λ ) , (31) i.e. the MSD behaves as a simple mean. This scalingpostulate is justified by two facts. Firstly, the barostat is(by construction) too weak to change for t (cid:28) τ (cid:63) the evo-lution of the system. Secondly, that the starting points ofthe trajectory at t = 0 are distributed according to theconsidered ensemble must become an irrelevant higherorder effect (vanishing rapidly with the system volume V ), since the MSD probes displacements with respect tothe starting points, not their absolute values. The fun-damental scaling, eq. (31), implies using eq. (20) that C ab ( t ) = µ ab ( λ ) − g ab ( t ) for t (cid:28) τ (cid:63) ( λ ) (32)where only the first term on the right hand-side dependson λ . This leads us finally to the general transformationrelation between correlation functions C ab ( t ) | λ = C ab ( t ) | λ =1 + µ ab | λ − µ ab | λ =1 (33)where λ = 1 stands for the NV γ T-ensemble with γ = γ ext The two static contributions on the right hand-side canbe further simplified using results from sect. II A. Thisdemonstrates the linear relation µ ab | λ − µ ab | λ =1 ∼ − λ . Strain-strain correlations.
Since ˆ γ ( t ) ≈ ˆ γ (0) for t (cid:28) τ (cid:63) , the MSD g γγ ( t ) must vanish to leading order. Usingeq. (32) this implies [20] C γγ ( t ) = µ γγ = (1 − λ ) /G eq for t (cid:28) τ (cid:63) ( λ ) . (34)This result is directly obtained by setting ˆ γ ( t ) = ˆ γ (0)in the definition of C γγ ( t ). We emphasize that it is alsoconsistent with the LPV transform, eq. (9), as one verifiesby setting I = τ , X = V γ , A = V γ ( t ) and B = V γ (0)and using that the strain fluctuation term for λ = 1 mustvanish. Strain-stress correlations.
To determine the strain-stress correlation function C γτ ( t ) one may use again thatthe MSD g γτ ( t ) must vanish since ˆ γ ( t ) ≈ ˆ γ (0) for t (cid:28) τ (cid:63) .Using eq. (32) this implies [20] C γτ ( t ) = µ γτ = 1 − λ for t (cid:28) τ (cid:63) ( λ ) . (35)This result is also obtained from the LPV transform set-ting A = V γ ( t ) and B = τ (0) and using again that (cid:104) δ ˆ Aδ ˆ B (cid:105)| λ =1 = 0. The postulated scaling eq. (31) hasthus led again to a reasonable result. Stress-stress correlations.
Interestingly, as one maysee from the NV γ T-ensemble limit considered in [9, 10],the stress-stress MSD g ττ ( t ) is not expected to simplyvanish for t (cid:28) τ (cid:63) as in the two previous cases. (Theinstantaneous stress ˆ τ fluctuates even at a fixed strainˆ γ .) However, eq. (33) still holds leading to [20] C ττ ( t ) | λ = C ττ ( t ) | λ =1 + (1 − λ ) G eq for t (cid:28) τ (cid:63) ( λ ) (36)where the static term on the right hand-side has beensimplified using eq. (5). This confirms the key relationeq. (6) announced in the Introduction. Note that thecorrelation function C ττ ( t ) | λ =1 in the NV γ T-ensemblemust vanish beyond some local time scale τ A which doesdepend on the network considered but, of course, not onthe shear-barostat. For a sufficiently slow barostat thecorrelation function thus becomes constant C ττ ( t ) = (1 − λ ) G eq for τ A (cid:28) t (cid:28) τ (cid:63) ( λ ) (37)where we have dropped | λ . Using eq. (32) this leads tothe remarkable relation g ττ ( t ) = µ A − G eq for τ A (cid:28) t (cid:28) τ (cid:63) ( λ ) (38)which must hold for all λ ≤
1. We note finally thatwhile for λ < g ττ ( t ) → µ ττ for t (cid:29) τ (cid:63) , g ττ ( t ) = µ A − G eq holds for all times t (cid:29) τ A for λ = 1according to stress-fluctuation formula, eq. (16). E. Dynamics: Macroscopic linear response
Introduction.
The experimentally important macro-scopic linear response measured by the creep compliance J ( t ) and the shear relaxation modulus G ( t ) [16, 17, 27]may be obtained conveniently in an equilibrium simula-tion at a given λ using some of the correlation functionsdiscussed above. Please note that being material prop-erties of the given state point (experimentally obtainedusing a simple average, not a fluctuation) both responsefunctions J ( t ) and G ( t ) do, of course, not depend on λ . Creep compliance.
Let us first consider the creep com-pliance J ( t ) ≡ (cid:104) δ ˆ γ ( t ) (cid:105) /δτ ext for t ≥
0. It is assumed thatfor t < τ equals the applied external stress τ ext of the NV τ T-ensemble. After imposing at t = 0 asmall increment δτ ext , the creep compliance J ( t ) mea-sures the ensuing average strain increment (cid:104) δ ˆ γ ( t ) (cid:105) . Wenote en passant that the average internal shear stress (cid:104) ˆ τ ( t ) (cid:105) does neither immediately reach the new equilib-rium value τ ext + δτ ext but shows a similar time depen-dence as the strain. Reworking the arguments put for-ward by Doi and Edwards, see eq. (3.67) of ref. [27], it isseen that J ( t ) = g γγ ( t ) | λ =0 for | δτ ext | (cid:28) , (39)i.e. the creep compliance is most readily computed usingthe strain-strain MSD g γγ ( t ) in the NV τ T-ensemble. Asdescribed in sect. III E, we shall change the shear-strainˆ γ ( t ) using a MC shear-barostat which corresponds to aperfectly viscous dynamics. One thus expects C γγ ( t ) = µ γγ exp( − t/T γγ ) for the strain-strain correlations. Using µ γγ = 1 /G eq for λ = 0 and eq. (20) this suggest J ( t ) = 1 G eq [1 − exp( − t/T γγ )] (40)in agreement with the Kelvin-Voigt model [17] represent-ing a purely viscous damper and a purely elastic springconnected in parallel. Shear relaxation modulus.
The shear relaxation mod-ulus G ( t ) ≡ (cid:104) δ ˆ τ ( t ) (cid:105) /δγ may be obtained from the stressincrement (cid:104) δ ˆ τ ( t ) (cid:105) for t > | δγ | (cid:28) t = 0. It is wellknown that the components of the Fourier transformedrelaxation modulus G ( t ), the storage modulus G (cid:48) ( ω ) andthe loss modulus G (cid:48)(cid:48) ( ω ), are directly measurable in an os-cillatory shear strain experiment [10, 17]. As seen, e.g.,by eq. (32) in ref. [10], it can be demonstrated by inte-gration by parts that G ( t ) = C ττ ( t ) | λ =0 for t (cid:28) τ (cid:63) ( λ = 0) , (41)i.e. the barostat should be irrelevant on the time scalesconsidered. Using the transformation relation betweendifferent ensembles, eq. (36), the relaxation modulus mayequivalently be obtained for other λ according to G ( t ) = C ττ ( t ) | λ + λG eq for t (cid:28) τ (cid:63) ( λ ) . (42)Since for large times G ( t ) → G eq , this implies C ττ ( t ) | λ → (1 − λ ) G eq consistently with eq. (37). For the specific case λ = 1 eq. (42) yields G ( t ) = C ττ ( t ) | λ =1 + G eq (43)which holds for all times t since the barostat becomesirrelevant for λ →
1. Note that eq. (43) may also be de-rived directly (without using the transformation relationbetween different ensembles) using Boltzmann’s super-position principle for an arbitrary strain history and thestandard fluctuation-dissipation theorem for the after-effect function [11] as shown elsewhere [9, 10]. Two im-mediate consequences of eq. (43) are (i) that G ( t ) onlybecomes equivalent to C ττ ( t ) | λ =1 for t > G eq = 0 and (ii) that theshear modulus G eq is only probed by G ( t ) on time scales t (cid:29) τ A where C ττ ( t ) | λ =1 must vanish. In principle, it isthus impossible to obtain the static shear modulus G eq ofan elastic body only from C ττ ( t ) | λ =1 as often incorrectlyassumed [26, 30]. III. ALGORITHMIC DETAILSA. Model Hamiltonian
To illustrate our key relations we present in sect. IV nu-merical data obtained using a periodic two-dimensional( d = 2) network of N l harmonic springs connecting N vertices. The model Hamiltonian is given by the sumˆ H = ˆ H id + ˆ H ex of a kinetic energy contributionˆ H id = m N (cid:88) i =1 v i , (44)with v i being the velocity of vertex i and m its (assumed)monodisperse mass, and an excess potentialˆ H ex = N l (cid:88) l =1 u l ( r l ) with u l ( r ) = 12 K l ( r − R l ) (45)where K l denotes the spring constant, R l the referencelength and r l = | r i − r j | the length of spring l . The sumruns over all springs l connecting pairs of beads i and j with i < j at positions r i and r j . The vertex mass m andBoltzmann’s constant k B are set to unity and Lennard-Jones (LJ) units [1] are assumed. B. Canonical affine shear transformations
While the box volume V is kept constant throughoutthis work, we shall frequently change the shape of thesimulation box. As sketched in panel (a) of fig. 1, weperform plane shear transformations of the instantaneousshear strain ˆ γ → ˆ γ + δγ with an essentially infinitesimalstrain increment δγ . We assume that not only the boxshape is changed but that the particle positions r (us-ing the principal box convention [1]) follow the imposedmacroscopic constraint in an affine manner according to r x → r x + δγ r y for | δγ | (cid:28) canonical [4, 31]. This implies that the x -component of the velocity must transform as [10] v x → v x − δγ v y with | δγ | (cid:28) . (47)We emphasize the negative sign in eq. (47) which assuresthat Liouville’s theorem is obeyed [10, 31]. C. Shear stress and affine shear-elasticity
Let ˆ H ( δγ ) = ˆ H id ( δγ ) + ˆ H ex ( δγ ) denote the systemHamiltonian of a configuration originally at ˆ γ strainedusing a canonical affine transformation to ˆ γ + δγ com-pactly written as a function of the strain increment δγ .The instantaneous shear stress ˆ τ and the instantaneousaffine shear-elasticity ˆ µ A may be defined as the expansioncoefficients associated to the energy change δ ˆ H ( δγ ) /V = ˆ τ δγ + ˆ µ A δγ / | δγ | (cid:28) γ being the reference, i.e. [32]ˆ τ ≡ ˆ H (cid:48) ( δγ ) /V | ˆ γ and (49)ˆ µ A ≡ ˆ H (cid:48)(cid:48) ( δγ ) /V | ˆ γ . (50)The derivatives ˆ H (cid:48) ( δγ ) and ˆ H (cid:48)(cid:48) ( δγ ) with respect to δγ may be computed as shown in sect. 2.1 of [10]. Similarrelations apply 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 . Using eq. (44) this implies [10] ˆ τ id = − V N (cid:88) i =1 m i v i,x v i,y and (51)ˆ µ A , id = 1 V N (cid:88) i =1 m i v i,y (52)for the ideal contributions to the shear stress and theaffine shear-elasticity. Note that the minus sign for theshear stress is due to the minus sign in eq. (47) requiredfor a canonical transformation. For the excess contribu-tions one obtains [10]ˆ τ ex = 1 V (cid:88) l r l u (cid:48) ( r l ) n l,x n l,y and (53)ˆ µ 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 (54)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. (53) is strictly identical to the correspondingoff-diagonal term of the Kirkwood stress tensor [1]. Thelast term in eq. (54) automatically takes into account thefinite normal pressure of the system. D. Groundstate characterization
Specific network.
As explained elsewhere [8, 10], thespecific elastic network used in this work has been con-structed using the dynamical matrix of a quenched poly-disperse LJ bead glass, i.e. at low temperatures our net-work has exactly the same mechanical and vibrationalproperties as the original discrete particle system. Priorto forming the network the bead system was cooled downto T = 0 using a constant quenching rate and imposinga normal pressure P = 2. The original LJ beads arerepresented in the snapshot shown in fig. 2 by grey poly-disperse circles, the permanent spring network createdfrom the quenched bead system by lines between vertices.The dark (black) lines indicate repulsive forces betweenthe vertices, while the light (red) lines represent tensileforces. The line width is proportional to the tension (re-pulsion). Note that the force network is strongly inho-mogeneous with zones of weak attractive links embeddedwithin a strong repulsive skeleton as already discussedin refs. [24, 25]. Only a small subvolume of the networkis represented. The total periodic box of linear length L ≈ . N = 10 vertices and N l = 9956springs. The monomer density ρ is close to unity. Finite shear stress τ . Due to the construction of thenetwork the total force acting on each vertex of the ref-erence network must vanish at T = 0. As seen in thesnapshot, this does not imply that the repulsive and/or -0.005 0.000 0.005 γ -0.20-0.100.000.100.20 τ ( γ ) affinerelaxed τ + µ A γτ +G eq γ µ A =34.2 τ ≅ G eq =16.3 γ ≅ - . r e l ax r e l ax (a)(b) FIG. 2: Groundstate of elastic network model considered inthis work assuming eq. (45): (a)
Snapshot of a small subvol-ume of linear length 10 containing about 100 vertices. Thelines represent the quenched forces of the athermal ( T = 0)reference configuration. (b) Shear stress τ ( γ ) assuming anaffine strain (squares) and after energy minimization (cir-cles). In agreement with ref. [10] we find in the first case τ ( γ ) = τ + µ A γ with µ A = 34 . τ ( γ ) = τ + G eq γ with G eq = 16 . γ = 0 where τ = τ = 0 . γ = γ = − . tensile forces transmitted along the springs must alsovanish. Due to the periodic boundary conditions andthe constant-strain constraint ( γ = 0) the shear stress τ does not necessarily vanish. For a pair potential such aseq. (45) the relevant excess contribution ˆ τ ex of the shearstress is readily computed using the Kirkwook expression,eq. (53). As shown in the main panel of fig. 2 (horizon-tal arrow), it turns out that for the specific network weuse throughout this work we have a finite shear stress τ ≡ τ ( γ = 0) = 0 . Affine shear-elasticity µ A . Let us consider a smallaffine shear strain according to eq. (46). As show inpanel (b) of fig. 2, one may now compute using eq. (53)the shear stress τ ( γ ) for different γ (squares). As oneexpects from the definition of the affine shear-elasticitycoefficient µ A , eq. (48), this yields the linear relation τ ( γ ) = τ + µ A γ with µ A ≈ . µ A is also obtained di-rectly from the unstrained configuration using eq. (54). Equilibrium shear modulus G eq . The forces f i actingon the vertices i of an affinely strained network do notvanish in general and the system is normally not at me-chanical equilibrium. As described elsewhere [10, 25], werelax these forces by first applying a steepest descend al-gorithm, i.e. by imposing displacements proportional tothe force, and then by means of the conjugate gradientmethod [6]. The ensuing non-affine displacements of thevertices decrease the energy [10] and the magnitude of the shear stress τ ( γ ) as may be seen from the large cir-cles indicated in fig. 2. As indicated by the bold solidline these final stresses τ ( γ ) scale again linearly as τ ( γ ) = τ + G eq γ = G eq ( γ − γ ) (55)with G eq ≈ . γ = − τ /G eq ≈ − . µ A has thus been replaced by the muchsmaller equilibrium shear modulus G eq of the ground-state network. Note that shear stress vanishes at a strain γ = γ as indicated by the vertical arrow. If the strain γ is allowed to change freely (e.g., using a steepest de-scend scheme) without any external force τ ext applied,the system relaxes to γ = γ . E. Finite temperature simulations
Molecular dynamics scheme.
As discussed in sect. IV,this network is investigated numerically basically bymeans of a molecular dynamics (MD) simulation [1, 4]at constant particle number N = 10 , box volume V ≈ . and a small, but finite mean temperature T = 0 . δt MD =10 − . The temperature T is fixed using a Langevin ther-mostat with a relatively large friction constant ζ = 1.This was done to suppress artificial long-range correla-tions in the two-dimensional periodic simulation box. Thermal averages.
Using the measured instantaneousshear stress ˆ τ and affine shear-elasticity ˆ µ A we compute,e.g., the thermal averages τ ≡ (cid:104) ˆ τ (cid:105) , µ ττ ≡ βV (cid:10) δ ˆ τ (cid:11) and µ A ≡ (cid:104) ˆ µ A (cid:105) . It can be shown for the respective idealcontributions that [10] τ id = µ ττ, id = µ A , id = P id = T ρ (56)with P id being the ideal normal pressure contribution.This holds irrespective of the considered λ -ensemble. Theideal contributions are thus negligible. Due to the equiv-alence of the different axes it is also seen for the excesscontributions that [10] µ 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) (57)being the Born-Lam´e coefficient [7, 22–25, 33] and P ex the excess part of the normal pressure P = P id + P ex . Monte Carlo shear-barostat.
The MD algorithm forthe particles is coupled for λ < γ → ˆ γ + δγ (sect. III B).First, a strain increment δγ is randomly chosen froma uniformly distributed interval [ − δγ max , δγ max ]. In or-der to determine the Metropolis weight [5] we computenext the energy change δ ˆ H = ˆ H (ˆ γ + δγ ) − ˆ H (ˆ γ ) ofthe network which comprises both an excess contribu-tion due to eq. (46) and an ideal contribution due to G ext λ µ γγ µ ττ A η T γγ T γτ θ ττ -10 -1.59 0.159 60.1 0.997 0.0081 9.2 4.3 13-6 -0.58 0.097 43.7 0.996 0.0103 5.6 2.7 6-3 -0.23 0.075 37.9 0.996 0.0117 4.3 2.1 40 0 0.061 34.2 0.996 0.0130 3.5 1.7 216.3 0.50 0.031 26.1 0.994 0.0184 1.8 0.85 1100 0.86 0.009 20.1 0.993 0.0348 0.5 0.24 0.31000 0.98 0.001 18.1 0.979 0.1030 0.06 0.03 0.310000 0.998 - 17.9 0.935 0.3233 0.007 0.003 0.4TABLE I: Summary of some properties as a function of theexternal spring constant G ext : λ = G ext / ( G eq + G ext ) with G eq = 16 .
3, strain-strain fluctuation µ γγ (fig. 5), stress-stressfluctuation µ ττ (fig. 7), acceptance rate A of MC shear-barostat (fig. 3), η ( λ, κ )-parameter defined by eq. (61), strain-strain crossover time T γγ ≈ θ γγ (fig. 8), strain-stress crossovertime T γτ ≈ θ γτ (fig. 10), stress-stress relaxation time θ ττ (fig. 15). The dynamical properties (columns 5-9) are onlygiven for κ = 10 − . The time scales T γγ , T γτ and θ ττ shouldbe compared to the intrinsic relaxation time τ A ≈ .
13 of thestress-stress correlations for κ → eq. (47). Since our system is subjected to an externalfield U ext (ˆ γ ), eq. (3), this gives rise to an additional con-tribution δU ext = U ext (ˆ γ + δγ ) − U ext (ˆ γ ). The suggestedstrain move δγ is accepted if ξ ≤ exp[ − β ( δ ˆ H + δU ext )] (58)with ξ being a uniformly distributed random variablewith 0 ≤ ξ < Operational parameters λ and κ . We assume for theexternal field U ext (ˆ γ ) that τ ext ≡ γ ext ≡ γ ,i.e. the average shear stress τ is imposed to vanishfor all ensembles studied and all ensembles compare thesame thermodynamic state point. (This was explicitlychecked.) The only remaining operational parameterfrom the static point of view is the external spring con-stant G ext or, equivalently, λ ≡ G ext / ( G eq + G ext ) assketched in panel (c) of fig. 1. As seen in table I orfig. 3, we vary G ext from −
10, i.e. λ ≈ − .
59, over G ext = λ = 0 (NV τ T-ensemble) up to G ext = 10000,i.e. λ ≈ . γ T-ensemble, i.e. the strain fluctua-tions become irrelevant for most properties. The secondoperational parameter of this study is the maximum at-tempted strain displacement δγ max which determines theimpact of the shear-barostat on the relaxation dynamics.Since the Metropolis MC move for the strain is performedevery MD time step of length δt MD , it is convenient touse instead of δγ max the maximum strain increment rate κ ≡ δγ max /δt MD . (Since all simulations are performedwith the same δt MD , δγ max and κ are strictly equivalent.)We compare below the dynamical strain and stress cor-relations for the five rates κ = 1 , − , − , − and10 − for a broad range of λ . For λ = 0 . κ = 3 , . , . , .
003 and
FIG. 3: Acceptance rate A of the Metropolis MC shear-barostat, eq. (58), for a large range of λ and κ as a functionof the scaling variable η ( λ, κ ) as defined by eq. (61). . κ = 10 − are summarizedin the columns 5-9 of table I. The simulations becomeincreasingly time consuming with decreasing κ and dataobtained for κ = 10 − have to be taken with care. Free strain diffusion limit.
The impact of κ can bebest judged from the simple limit where neither the sys-tem nor the external field restricts the barostat. This isthe case for sufficiently small κ depending on λ . In thislimit one expects the free diffusion of the strain ˆ γ ( t ), i.e. g γγ ( t ) = D γγ t/ t (cid:28) T γγ . Since the (attempted andaccepted) strain increment δγ is uniformly distributedin [ − δγ max , δγ max ], we have a mean-square strain step (cid:104) δγ (cid:105) = δγ / δt MD . This implies a diffusionconstant D γγ = βV δγ δt MD ∼ κ (59)in the free-diffusion limit. Using eq. (28) and eq. (4) thisyields in turn the corresponding crossover time T γγ = 6 µ γγ βV δt MD δγ = 6 δt MD /η (60)where we have defined the dimensionless variable η ( λ, κ ) ≡ (cid:113) βV δγ /µ γγ = κ δt MD (cid:112) (cid:104) δ ˆ γ (cid:105) . (61)Since G eq and δt MD are kept constant in all our simula-tions, the parameter η determines the dynamical regimefor a system of operational parameters λ and κ . As seenin fig. 3 using η as a scaling variable the acceptance rate A of the MC shear-barostat of a broad range of λ and κ can be brought to collapse. For 1 (cid:28) η (cid:28)
10 theMC barostat is most efficient for static properties. Fordynamical properties we shall focus below on κ -valueswhere η (cid:28) . A ≈ FIG. 4: Normalized histogram p ( x ) of x = ˆ γ − γ ext for different λ . The line indicates for λ = − .
59 ( G ext = −
10) a Gaussianwith (cid:10) δ ˆ γ (cid:11) = 1 /βV G eff ≈ . .FIG. 5: Second moment of strain fluctuations: (a) Dimen-sionless y = µ γγ G eq with G eq ≡ . λ showing the linear decay (bold line) expected from eq. (4). (b) Shear modulus G γγ , eq. (62), confirming G eq ≈ . (c) MSD g γγ ( t ) vs. time t for several λ and κ = 10 − . The dash-dotted line indicates the free-diffusionlimit, eq. (59), for short times, the horizontal lines the long-time limit µ γγ for each λ . Also indicated is for λ = 0 and κ = 1 (filled circles) an example with η (cid:29) . IV. COMPUTATIONAL RESULTSA. Static properties
Strain-strain fluctuations.
The equilibrium fluctua-tions of the strain ˆ γ at a finite temperature T = 0 .
001 fordifferent values of G ext or, equivalently, λ = G ext / ( G eq + G ext ) with G eq ≡ . FIG. 6: Strain-stress correlations: (a)
Pre-averaged shearstress τ ( x ) as a function of x = ˆ γ − γ ext . All data collapse onthe linear slope with G eq = 16 . (b) Scat-ter plot of 100 data points (ˆ γ, ˆ τ ) for λ = − .
23 ( G ext = − δt = 10 − (dashed line) and squares to a large δt = 10 +2 (bold solidline) between subsequent configurations. The given data refer to a maximum strain rate κ = 10 − for the Metropolis MC step ˆ γ → ˆ γ + δγ attempted aftereach MD time step of δt MD = 10 − (sect. III E). That thisleads indeed to the desired Gaussian distributions can beseen from fig. 4 where the normalized histogram p ( x ) isplotted as a function of x = ˆ γ − γ ext with γ ext = γ .As one expects, p ( x ) → δ ( x ) for λ →
1. The rescaledsecond moment µ γγ ≡ βV (cid:104) δ ˆ γ (cid:105) of the distribution is fur-ther analyzed in fig. 5. As one may see from panel (a), y ( λ ) ≡ µ γγ G eq decreases linearly with λ (bold solid line)and vanishes for λ = 1 as it should according to eq. (10).In turn this implies that one may determine G eq from G γγ ≡ /µ γγ − G ext . (62)As can be seen from panel (b) of fig. 5, using this relationone confirms that the groundstate value G eq ≈ . T (cid:28)
1. The ac-curate determination of G γγ becomes of course more del-icate with increasing λ since the MC acceptance rate be-comes eventually too small. The data quality in this limitcan be readily improved (not shown) by either increasingthe sampling time t traj or by decreasing κ ∼ η . (Thesampling becomes again inefficient if η gets too small.) Strain-stress correlations.
Since the shear modulus G eq is finite, this implies that shear strain and shearstress fluctuations must be correlated [8]. These corre-lations are addressed in fig. 6. The main panel presentsthe pre-averaged stress τ ( x ) for an instantaneous strain x = ˆ γ − γ ext for several λ [34]. Irrespective of x or λ all data collapse on the linear slope indicated by thebold line with G eq = 16 .
3. Naturally, it follows fromfig. 4 that the statistics must strongly decrease for x (cid:29) FIG. 7: Determination of shear modulus G eq ≈ . G γγ for λ < G γτ for λ < G ττ for λ (cid:54) = 0. Also indicated are the affine shear-elasticity µ A ( λ ) ≈ . µ ττ whichis seen to decay linearly (bold dash-dotted line) according toeq. (5) down to µ ττ | λ =1 ≈ µ γτ ≡ βV (cid:104) δ ˆ γδ ˆ τ (cid:105) , i.e. a larger linear-slope window is vis-ible in fig. 6 for smaller λ . Panel (b) of fig. 6 presentstwo scatter plots of (ˆ γ, ˆ τ ) for λ = − .
23 ( G ext = − δt = δt MD . On these time scales the configuration hasno time to relax the affine displacements imposed by theshear-barostat. As emphasized by the dashed slope thisleads to a linear stress-strain relation with a coefficient µ A = 34 .
2. Other short-time sequences yield similar,but horizontally shifted linear slopes (not shown). A δt -independent representation of the static strain-stresscorrelations for large times is obtained if (ˆ γ, ˆ τ ) is indi-cated for 100 configurations with a time interval δt = 10 (squares). The scatter plot of the latter data is already innice agreement with a coefficient G eq = 16 . all data tuples (ˆ γ, ˆ τ ) of a given λ -ensembleone verifies that the shear modulus G eq is accurately de-termined using the linear regression coefficient G γτ ≡ (cid:104) δ ˆ γδ ˆ τ (cid:105)(cid:104) δ ˆ γ (cid:105) = µ γτ µ γγ for λ < . (63)Since, as shown in sect. II A, µ γγ = (1 − λ ) /G eq and µ γτ = 1 − λ , this ratio must yield G eq for all λ < G eq = 16 . G γτ for different λ are identical to G γγ as shown in fig. 7. Stress-stress fluctuations.
The stress-stress fluctua-tions µ ττ ≡ βV (cid:104) δ ˆ τ (cid:105) presented in fig. 7 (circles) decreaselinearly with λ as stated in the Introduction, eq. (5).Being a simple mean µ A (small squares) is found tobe strictly λ -independent. As expected, µ ττ → µ A for λ →
0. Using eq. (5) the generalized stress-fluctuation
FIG. 8: Diffusion of instantaneous strain ˆ γ ( t ): (a) MSD g γγ ( t ) for λ = 0 (open symbols) compared to the explicitlycomputed creep compliance J ( t ) = (cid:104) δ ˆ γ ( t ) (cid:105) /δτ ext for a stepstress δτ ext = 0 .
01 (filled symbols). The dash-dotted lines in-dicate the regime used to determine the diffusion coefficient D γγ , the bold horizontal line the expected asymptotic limit µ γγ and the thin solid line the Kelvin-Voigt model, eq. (40),for κ = 10 − . (b) Double-logarithmic representation of therescaled crossover time T γγ / δt MD vs. η . The dash-dottedline indicates the expected power law, eq. (60). formula for G eq reads G ττ ≡ µ A − µ ττ − ( µ A − µ ττ ) /G ext for λ (cid:54) = 0 . (64)For λ →
1, i.e. G ext → ∞ , this reduces to eq. (16). Equa-tion (64) thus generalizes the stress-fluctuation formulafor the NV γ T-ensemble [4, 8–10, 21–23, 33] to general λ . As seen from fig. 7 (diamonds), it is thus possible todetermine G eq from µ A and µ ττ for all λ (cid:54) = 0. B. Dynamics: Strain-strain correlations
Strain-strain MSD.
The strain histograms p ( x ) fordifferent λ and their second moments shown, respectively,in fig. 4 and fig. 5 have been rapidly sampled by aver-aging over short time series from trajectories of length t traj = 10 . This is demonstrated in panel (c) of fig. 5presenting the MSD g γγ ( t ) of the instantaneous strainˆ γ ( t ) for κ = 10 − . The statistics is improved by perform-ing a gliding average over the 10 data tuples sampled [1].It is seen that g γγ ( t ) converges rapidly after a time T γγ of order unity to its asymptotic limit µ γγ ( λ ) indicatedby horizontal lines. Free diffusion is observed for shorttimes, i.e. g γγ ( t ) ≈ D γγ t/ t (cid:28) T γγ (dash-dottedline), and this essentially irrespective of the ensemble inagreement with eq. (59). Note that free strain diffusionis observed for all λ and sufficiently small κ if the param-eter η ( λ, κ ) (cid:28)
1. As shown in the main panel of fig. 8for different κ and λ = 0 (open symbols), this behavior2 FIG. 9: Determination of equilibrium shear modulus G eq as afunction of the sampling time t for G γγ ( t ) (thin lines), G γτ ( t )(filled symbols) and G ττ ( t ) (large open symbols) for κ = 10 − .All fluctuation formulae are equivalent for t > G eq for large times. Note that G γτ ( t ) for all λ < G ττ ( t ) for λ = 1. may be used to determine first the short-time diffusioncoefficient D γγ for systems with sufficiently small η andthen using eq. (28) the crossover time T γγ of the strainfluctuations. The rescaled times T γγ / δt MD are repre-sented in panel (b) of fig. 8 as a function of the scalingvariable η for several λ . The dash-dotted line indicatesthe power-law slope, eq. (60). A perfect data collapse onthis line is observed for all η (cid:28)
1. A successful data col-lapse can be also obtained for g γγ ( t ) over a broad rangeof λ and κ by plotting g γγ ( t ) /µ γγ as a function of thereduced time x = t/T γγ (not shown). As seen by thethin solid line in the main panel of fig. 8, one verifiesthat the scaling function f ( x ) = g γγ ( t ) /µ γγ is given by f ( x ) = 1 − exp( − x ). This is, of course, consistent with anexponentially decaying strain-strain relaxation function C γγ ( t ) = µ γγ exp( − x ) as we have also checked directly. Creep compliance J ( t ) . As discussed in sect. II E,a MSD g γγ ( t ) computed in the NV τ T-ensemble corre-sponds to a creep compliance J ( t ) = (cid:104) δ ˆ γ ( t ) (cid:105) /δτ ext mea-sured by imposing at t = 0 a step stress increment | δτ ext | (cid:28) τ ext applied to thesystem. Using δτ ext = 0 .
01 and averaging over 100 config-urations this yields the data shown by the filled symbolsin panel (a) of fig. 8 for two values of κ . An excellent datacollapse g γγ ( t ) ≈ J ( t ) is observed. While it is natural touse the NV τ T-ensemble for obtaining the creep compli-ance from the equilibrium fluctuations, it is worthwhileto note that due to the scaling g γγ ( t ) /µ γγ = f ( x ), J ( t )may be also determined from other λ -values. Time-dependent strain-strain fluctuations.
As ex-pected from sect. II B, the strain-strain fluctuations µ γγ ( t ) computed according to eq. (22) by gliding aver-age from a finite time series increase monotonously with FIG. 10: Scaling of strain-stress MSD g γτ ( t ) for different λ : (a) scaling collapse of y = g γτ ( t ) /µ γτ vs. x = t/T γτ for κ =10 − with the dash-dotted line indicating the diffusive regimeand the bold solid line the long-time limit, (b) crossover time T γτ , eq. (28), rescaled as T γτ / δt MD vs. the scaling variable η using the same symbols as in the main panel. The dash-dotted line indicates the prediction eq. (65). time t from zero (if only one configuration is sampled) tothe asymptotic thermodynamic limit µ γγ (not shown).Note that one determines first the strain fluctuations inan interval [ t , t = t + t ] and averages then over allpossible t . If the estimate G γγ , eq. (62), of the modulus G eq is determined using µ γγ ( t ) instead of µ γγ , it becomesa time dependent quantity called G γγ ( t ). As can be seenfrom fig. 9 for three different λ values, G γγ ( t ) thus be-comes a monotonously decreasing function of time (thinlines). Interestingly, G γγ ( t ) appears not to depend on λ .We note finally that using either the large-time asymp-totics eq. (27) or, more conveniently, the Debye relationeq. (29), one may determine the relaxation time θ γγ ofthe strain-strain fluctuations µ γγ ( t ). One verifies that θ γγ ≈ T γγ holds as expected for an exponentially decay-ing correlation function (sect. II C). C. Dynamics: Strain-stress correlations
Strain-stress MSD.
The scaling of the strain-stressMSD g γτ ( t ) is investigated in fig. 10. The relaxation time T γτ presented in panel (b) has been determined, as be-fore for T γγ , by matching the short-time diffusive regimewith the long-time limit g γτ ( t ) → µ γτ . As indicated bythe dash-dotted line, the data are well described by T γτ = G eq µ A T γγ ∼ /κ for η (cid:28) . (65)That this scaling must hold can be seen by replacingin the definition of g γτ ( t ) the stress fluctuation δ ˆ τ ( t ) by µ A δ ˆ γ ( t ) which implies D γτ = µ A D γγ for the diffusion3coefficients and in turn eq. (65) using again eq. (28). Westress that in agreement with the inset of fig. 6 the systemhas not enough time to relax the affine strain imposed bythe shear-barostat for the short times t (cid:28) τ A ≈ . µ A and not G eq which has to beused as the linear slope coefficient. The main panel (a)presents the collapse of y = g γτ ( t ) /µ γτ as a function of x = t/T γτ for κ = 10 − . Minor deviations are visible for x ≈
1. The quality of the collapse improves by decreas-ing η ( λ, κ ), whereas the crossover becomes more sudden(even a hump may occur) if η is too large (not shown). Time-dependent correlation coefficient G γτ ( t ) . Sincethe strain-shear correlations are dominated for t (cid:28) τ A by the affine shearing, one expects G γτ ( t ) computed us-ing the ratio of µ γτ ( t ) and µ γγ ( t ) to be similar to µ A for t →
0. As one sees from G γτ ( t ) presented in fig. 9(filled symbols), this is indeed the case. Interestingly, G γτ ( t ) is seen to be independent of λ for all times t , i.e.the λ -dependences of µ γτ ( t ) and µ γγ ( t ) do cancel for alltimes just as they cancel for the ratio of the equilibriummoments µ γτ /µ γγ . D. Dynamics: Stress-stress correlations
Introduction.
We turn finally to the more intricatecharacterization of the stress-stress correlations as a func-tion of our two operational parameters λ and κ . We focusagain on the regime with η ( λ, κ ) (cid:28)
1. As above we beginby discussing the MSD g ττ ( t ). Since according to eq. (20) g ττ ( t ) = C ττ (0) − C ττ ( t ) with C ττ (0) = µ ττ , (66)the MSD g ττ ( t ) and the correlation function C ττ ( t ) con-tain in principal the same information. From the presen-tational point of view g ττ ( t ) has the advantage that theshort-time power-law behavior of the correlations can bemade manifest using a double-logarithmic plot as shownin fig. 11. We describe then the large time behavior ofthe correlation function C ττ ( t ) and demonstrate that τ (cid:63) is given by T γγ . Finally, we turn to the time-dependentstress fluctuation µ ττ ( t ) which is used to determine thestress-stress relaxation time θ ττ . Stress-stress MSD.
Let us concentrate first on thedata for λ = − .
59 presented by open symbols in fig. 11.If κ is large, all internal dynamics is destroyed and wefind immediately g ττ ( t ) ≈ µ ττ as seen for κ = 1 (smallcircles). For smaller, but not too small κ one observesa free-diffusion regime with g ττ ( t ) = D ττ t/ g γτ ( t ) in sect. IV C, δ ˆ τ ( t ) in the definition of g ττ ( t ) canbe replaced by µ A δ ˆ γ ( t ). This argument yields the dif-fusion coefficient D ττ = µ D γγ successfully used in theplot for κ = 10 − and κ = 10 − . If κ is further reduced, g ττ ( t ) becomes κ -independent for short times. As seenfor κ = 10 − and 10 − the dynamics becomes “determin-istic” in the sense that g ττ ( t ) = µ ττ t/ ˜ θ ) ∼ κ λ for t (cid:28) τ A (67) FIG. 11: Stress-stress MSD g ττ ( t ) focusing on λ = − . κ (open symbols). Additionally, wepresent λ = 0 for κ = 1 (small filled circles) and κ = 10 − , λ = 0 . κ = 10 − and the NV γ T-ensemble (stars). Thedot-dashed lines indicate the free diffusive behavior for inter-mediate η , the thin solid line the analytic short-time limit forsmall η , the horizontal dashed line the thermodynamic limit µ ττ for λ = − .
59 and the bold solid line the expected inter-mediate plateau, eq. (38). Note that for sufficiently small κ all data approach the κ - and λ -independent lower envelope in-dicated by the solid lines in agreement with the fundamentalscaling postulate eq. (31). with ˜ θ ≈ .
16 as shown by the thin solid line. Albeit µ ττ depends on λ , the ratio µ ττ / ˜ θ does not, in agreementwith the fundamental postulate eq. (31). This impliesthat ˜ θ depends somewhat on λ . The t -scaling in eq. (67)is expected [10] since the associated correlation function C ττ ( t ) must be a continuous and symmetric function at t = 0 which, moreover, should be an analytic function ifbarostat and thermostat effects can be ignored [35]. Thisimplies that g ττ ( t ) must be an even expansion in termsof t which to leading order leads to eq. (67). Since thereare several dynamical regimes it is not meaningful to de-termine a crossover time T ττ using eq. (28) as before for T γγ and T γτ over the full range of λ and κ . We shall seeat the end of this section how a longest relaxation time θ ττ for the stress-stress correlations might be obtained.For times t (cid:29) τ A an intermediate plateau appears (boldsolid line) as predicted by eq. (37). As can be seen forfour different λ values this intermediate plateau does not depend on the ensemble provided that κ is sufficientlyweak. (As one expects from fig. 3 deviations from thisasymptotic limit arise again for large κ if η ( λ, κ ) (cid:29) . λ and κ wehave simulated. The MSD g ττ ( t ) behaves thus as anensemble-independent simple average as we have arguedin sect. II D. One may operationally define the crossover4 FIG. 12: Relaxation modulus G ( t ) obtained from the stressresponse to an applied strain δγ = 0 .
01 (stars) compared tothe equilibrium correlation function C ττ ( t ) | λ for several λ ob-tained by averaging over 1000 quenched-strain configurations. time τ A by matching eq. (67) and eq. (38). This yields τ A = ˜ θ (cid:115) − G eq /µ A − λG eq /µ A ≈ .
13 (68)as indicated by the vertical line. Note that τ A = ˜ θ √ λ = 1. We stress that while ˜ θ is a function of λ , τ A is an ensemble-independent constant for the givenelastic network. (Moreover, it can be shown that itbarely depends on the friction constant ζ of the Langevinthermostat [10].) Using eq. (68) one may reformulateeq. (67) in a manifest λ -independent form: g ττ ( t ) =( µ A − G eq ) ( t/τ A ) . The MSD ultimately approaches µ ττ ( λ ) for even larger times t (cid:29) τ (cid:63) as indicated by thedashed horizontal line for λ = − .
59. Remembering that µ ττ = µ A − G eq for λ = 1 it is clear that the NV γ T-datastays constant for all times.
Relaxation modulus.
The large- t scaling is further ad-dressed in fig. 12 and fig. 13 presenting the stress-stresscorrelation function C ττ ( t ). The most central result ofthis work stated by eq. (6) is demonstrated in fig. 12.We compare the directly measured relaxation modulus G ( t ) with the equilibrium correlation function C ττ ( t ) as-suming an asymptotically slow barostat ( τ (cid:63) = ∞ ). Av-eraging over 1000 independent configurations from anNV γ T-ensemble with γ = 0 for t <
0, the relaxationmodulus G ( t ) = (cid:104) ˆ τ ( t ) (cid:105) /δγ measures the stress responseafter a (canonical and affine) shear strain δγ = 0 .
01 wasbeen applied at t = 0 (stars). The correlation functions C ττ ( t ) for λ < λ -ensemble as in-dicated. Switching off the shear-barostat the strain ˆ γ ofeach configuration is quenched. Note that C ττ (0) = µ ττ holds due to the ensemble averaging. As expected fromeq. (37), C ττ ( t ) decreases monotonously from µ ττ down FIG. 13: Effect of the barostat parameter κ for C ττ ( t ) with λ = 0 .
5. The large stars correspond to the quenched γ -ensemble ( τ (cid:63) = ∞ ) and all other symbols to a switched-on barostat of finite κ . The vertical arrows indicate T γγ for κ = 10 − , − , − and 10 − confirming that roughly T γγ ≈ τ (cid:63) . A much longer run is warranted for κ = 10 − . to (1 − λ ) G eq for t (cid:29) τ A . This corresponds to the λ -independent intermediate plateau g ττ ( t ) = µ A − G eq in fig. 11. It is seen that G ( t ) = C ττ ( t ) for λ = 0(large circles) and that for all λ one may obtain G ( t ) byvertically shifting the correlation functions C ττ ( t ) | λ → C ττ ( t ) | λ + λG eq = G ( t ) using the already determinedshear modulus G eq . Note also that C ττ ( t ) | λ =1 → C ττ ( t ) remain finite. As al-ready emphasized at the end of sect. II E, it is thus impos-sible to obtain the modulus G eq solely from C ττ ( t ) | λ =1 .Please note that the observed short-time value µ ττ ≈ λ = 1 is quite different from theequilibrium modulus G eq ≈
16 (bold solid line).
Finite- κ effects. Focusing on one ensemble with λ =0 . κ = 0, τ (cid:63) = ∞ ). As one expects, the decay of C ττ ( t ) be-comes systematically slower with decreasing κ . For large κ > − we have computed C ττ ( t ) by gliding averagingover a single trajectory of t traj = 10 . For smaller κ itbecomes increasingly harder to have a sufficiently longtrajectory of { ˆ τ k } probing the phase space, i.e. for which C ττ (0) = µ ττ , and to store at the same time the shearstresses with a sufficiently fine time resolution. The datafor κ ≤ − have thus been obtained by using an ensem-ble of 100 configurations of different ˆ γ and averaging overthe trajectories obtained for each configuration using afinite κ . As expected from eq. (37), for sufficiently small κ all correlation functions have an intermediate plateaufor τ A (cid:28) t (cid:28) τ (cid:63) (bold solid line). That the barostat isnot completely switched off can be seen from the final de-cay. The vertical arrows indicate the crossover times T γγ for the four smallest κ . As one expects, this time scaleappears to coincide roughly (up to a constant prefactor)with the upper limit τ (cid:63) of the respective plateau. In other5 FIG. 14: Shear-stress fluctuation µ ττ ( t ) | λ as a function of thesampling time. Also indicated is the affine shear-elasticity µ A obtained for the NV τ T-ensemble (small filled squares)which is seen to immediately converge to the long-time limit.The thin lines correspond to the integrated shear-stress auto-correlation function C ττ ( t ) | λ according to eq. (25). words, the plateau disappears when the ergodicity break-ing associated with the averaging over a quenched-strainensemble is lifted by the shear-barostat. Time-dependent stress-fluctuation formula.
Let usreturn to fig. 9. According to the generalized stress-fluctuation formula G ττ , eq. (64), the shear modulus G eq may be obtained from the affine shear-elasticity µ A andthe equilibrium stress-stress fluctuation µ ττ for λ ≤ λ (cid:54) = 0 (fig. 7). If instead of µ ττ a time-averagedfluctuation µ ττ ( t ) is used, this leads again to a time-dependent monotonously decreasing estimate G ττ ( t ) asseen from the open symbols in fig. 9. No shear-barostatis of course applied for the NV γ T-ensemble indicated bythe large squares. We emphasize the remarkable findingthat within numerical accuracy G ττ ( t ) | λ =1 ≈ G γτ ( t ) | λ< (69)for all times and a broad range of κ for G γτ ( t ). As seen for λ = − .
58, for small λ and not too short times G ττ ( t ) and G γγ ( t ) are generally similar. Since µ ττ ( t ) must vanish forsmall times, this implies that G ττ ( t ) must become finite G ττ ( t ) → µ A − µ A /G ext for t → G ττ ( t ) → µ A for G ext (cid:29) µ A . This limit isindicated by the bold dashed line for λ = 1. Time-dependence of stress-stress fluctuations.
Thetime-dependence of µ ττ ( t ) for several λ is explicitly rep-resented in fig. 14 assuming κ = 10 − for λ <
1. Whilethe simple mean µ A converges immediately with time, µ ττ ( t ) increases monotonously from zero to the large-time limit µ ττ (horizontal lines). As already pointed outin sect. II B, assuming time-translational invariance thetime-dependence of µ ττ ( t ) is a consequence of the time-dependence of the C ττ ( t ). To emphasize this point we FIG. 15: Stress-stress relaxation time θ ττ : (a) Determina-tion of θ ττ from the decay of y = 1 − µ ττ ( t ) /µ ττ accord-ing to eq. (27) for λ = 0 . κ . The thin dashedlines indicate eq. (27), the bold horizontal line the expectedintermediate plateau (1 − λ ) G eq /µ ττ for 10 τ A (cid:28) t (cid:28) τ (cid:63) . (b) Relaxation time θ ττ for three not too large λ plotted as θ ττ / δt MD vs. η ( λ, κ ). The dash-dotted line indicates thepower law expected for η (cid:28) have integrated C ττ ( t ) as suggested by eq. (25). As seenby the thin lines indicated in fig. 14, this leads to con-sistent results. If one has thus characterized C ττ ( t ), thisimplies µ ττ ( t ) and then in turn G ττ ( t ). Note that for all λ the asymptotic limit is reached at about 100 τ A whichis of the same order as T γγ ≈ . λ = 0 and κ = 10 − . Relaxation time θ ττ . The convergence of µ ττ ( t ) is sys-tematically analyzed in fig. 15. The main panel presentsthe dimensionless deviation y = 1 − µ ττ ( t ) /µ ττ for several κ and λ = 0 .
5. As emphasized by the thin dashed lines,a correlation time θ ττ may be measured from the large-time 1 /t -decay following eq. (27). (Unfortunately, thesampling time is yet insufficient for our smallest κ .) Withdecreasing κ the decay becomes systematically slowerand a horizontal plateau with y = (1 − λ ) G eq /µ ττ (boldsolid line) appears in an intermediate time window. Thisplateau is expected from eq. (30) and eq. (37). Followingthe discussion in sect. II B, this implies that for suffi-ciently small κ and not too large λ the relaxation time θ ττ must scale (up to a prefactor of order unity) as θ ττ ≈ (1 − λ ) G eq µ ττ τ (cid:63) ∼ (1 − λ ) /κ for η (cid:28) τ (cid:63) ≈ T γγ ∼ /κ . That the κ -scaling holds can be seen inpanel (b) of fig. 15. (For λ = 0 . θ ττ also for the additional values κ = 3 , . , . , . . θ ττ this is for tworeasons not a rigorous scaling plot over the full rangeof operational parameters considered. Firstly, due to the6(1 − λ ) G eq /µ ττ -factor in eq. (71) there is an additional λ -dependence even for η (cid:28) λ . (This canbe accounted for using a different representation.) Moreimportantly, for λ ≈ τ A essentiallysets θ ττ and not τ (cid:63) ≈ T γγ . We thus observe θ ττ ≈ τ A λ κ for all λ > . κ = 10 − . Basically,since we have two time scales τ A and T γγ , it is not possibleto put all θ ττ ( λ, κ ) on one master curve. V. CONCLUSION
Static properties.
Focusing on two-dimensional elas-tic networks (sect. III) we have characterized theoreti-cally (sect. II) and numerically (sect. IV) the static anddynamical fluctuations of shear-strain and shear-stressin generalized Gaussian ensembles as a function of thedimensionless parameter λ = G ext / ( G eq + G ext ), fig. 1.Monitoring various properties we interpolate between theNV γ T-ensemble ( λ = 1) and the NV τ T-ensemble ( λ = 0)and consider even negative λ . The static strain-strainfluctuations µ γγ (fig. 5), the strain-stress fluctuations µ γτ and the stress-stress fluctuations µ ττ (fig. 7) have beenshown to decrease linearly with increasing λ [20]. As aconsequence, the static shear modulus G eq may be ob-tained (fig. 7) from either the strain-strain fluctuationformula G γγ , eq. (62), or from the strain-stress formula G γτ , eq. (63), for λ < G ττ , eq. (64), for λ (cid:54) = 0. The latter formulais a generalization of the well-known stress-fluctuationformula eq. (16) for λ = 1. When sampled from a fi-nite time series (ˆ γ k , ˆ τ k ), G γγ ( t ), G γτ ( t ) and G ττ ( t ) be-have similarly with time t for all λ (fig. 9). Interestingly, G ττ ( t ) | λ =1 and G γτ ( t ) | λ< are identical for all times. Dynamic properties.
The influence of the parameter κ of the MC shear-barostat (sect. III E) has been investi-gated for various dynamical properties, especially for thestress-stress correlation function C ττ ( t ) (fig. 13). The up-per time limit τ (cid:63) ( λ ), below which barostat effects can beignored, is set by the strain-strain crossover time T γγ ( λ )(fig. 8). While a rapid barostat may be of advantagefor static properties, a sufficiently slow barostat witha large τ (cid:63) corresponds to the fundamentally importantlinear-response limit where the barostat is only relevantfor the distribution of the start points of the trajectories (and, hence, C ab (0) = µ ab ) but not for their dynamicalpathways for t (cid:28) τ (cid:63) . We have argued that this impliesthe fundamental scaling g ab ( t ) ∼ λ κ for t (cid:28) τ (cid:63) ( λ ) asdemonstrated explicitly for g ττ ( t ) in fig. 11. Assumingsuch a slow barostat one may obtain G eq using eq. (6)from C ττ ( t ) | λ ≈ (1 − λ ) G eq for λ < G eq alone from the autocorrelation function C ττ ( t ) | λ =1 . The experimentally important shear-stressrelaxation modulus G ( t ) may be most readily determinedat λ = 1 calculating first G eq = µ A − µ ττ | λ =1 and then G ( t ) = C ττ ( t ) | λ =1 + G eq [36]. We have commented brieflyon the creep compliance J ( t ). Being well-described bythe Kelvin-Voigt model, eq. (40), J ( t ) is best obtainedfrom the strain-strain MSD g γγ ( t ) at λ = 0 (fig. 8). More general thermodynamic variables.
Most of thetheoretical relations and numerical techniques discussedabove, especially the key formulae eq. (5) and eq. (6),generalize readily for any pair of continuous extensive andintensive variables X and I with M eq = ∂I/∂X being theequilibrium modulus. With M ext being the spring con-stant of the external potential controlling the fluctuationsof the extensive variable and λ = M ext / ( M eq + M ext ) oneobtains, e.g., the relaxation modulus M ( t ) = C ( t ) | λ =0 = C ( t ) | λ + λM eq for t (cid:28) τ (cid:63) ( λ ) (72)with C ( t ) ≡ βV (cid:104) δ ˆ I ( t ) δ ˆ I (0) (cid:105) being the relevant au-tocorrelation function. The associated MSD g ( t ) ≡ ( βV / (cid:104) ( ˆ I ( t ) − ˆ I (0)) (cid:105) is expected to be strictly λ -independent in the same time window. These propertiesmust be computed again either using a slow “barostat”with a large, albeit finite τ (cid:63) (cid:28) t traj or, equivalently, byaveraging over an equilibrium ensemble of configurationswith frozen ˆ X . For times t (cid:29) τ (cid:63) a switched-on “baro-stat” ultimately restores the ergodicity and C ( t ) mustvanish while g ( t ) approaches its finite thermodynamicvalue βV (cid:104) δ ˆ I (cid:105) . Acknowledgments
H. Xu thanks the IRTG Soft Matter for financial sup-port. We are indebted to O. Benzerara (Strasbourg) andJ. Helfferich (Freiburg) for helpful discussions. [1] M. Allen and D. Tildesley,
Computer Simulation of Liq-uids (Oxford University Press, Oxford, 1994).[2] J. L. Lebowitz, J. K. Percus, and L. Verlet, Phys. Rev. , 250 (1967).[3] H. B. Callen,
Thermodynamics and an Introduction toThermostatistics (Wiley, New York, 1985).[4] D. Frenkel and B. Smit,
Understanding Molecular Sim-ulation – From Algorithms to Applications (AcademicPress, San Diego, 2002), 2nd edition. [5] D. P. Landau and K. Binder,
A Guide to Monte CarloSimulations in Statistical Physics (Cambridge UniversityPress, Cambridge, 2000).[6] J. Thijssen,
Computational Physics (Cambridge Univer-sity Press, Cambridge, 1999).[7] M. Born and K. Huang,
Dynamical Theory of CrystalLattices (Clarendon Press, Oxford, 1954).[8] J. P. Wittmer, H. Xu, P. Poli´nska, F. Weysser, andJ. Baschnagel, J. Chem. Phys. , 12A533 (2013). [9] J. P. Wittmer, H. Xu, and J. Baschnagel, Phys. Rev. E , 022107 (2015).[10] J. P. Wittmer, H. Xu, O. Benzerara, and J. Baschnagel,Mol. Phys. , 10.1080/00268976.2015.1023225 (2015).[11] J. Hansen and I. McDonald, Theory of simple liquids (Academic Press, New York, 2006), 3nd edition.[12] J. Hetherington, J. Low Temp. Phys. , 145 (1987).[13] M. Costeniuc, R. Ellis, H. Touchette, and B. Turkington,Phys. Rev. E , 026105 (2006).[14] K. van Workum and J. de Pablo, Phys. Rev. E , 011505(2003).[15] The focus of Hetheringtons’s Gaussian ensemble [12], asof related generalizations [13], is on the transformationbetween the microcanonical ensemble , characterized bythe (possibly non-concave) entropy as a function of theenergy, and the (generalized) canonical ensemble , char-acterized by the free energy as function of the inversetemperature β , i.e. different pairs of conjugated variablesare considered compared to the present work.[16] T. Witten and P. A. Pincus, Structured Fluids: Polymers,Colloids, Surfactants (Oxford University Press, Oxford,2004).[17] M. Rubinstein and R. Colby,
Polymer Physics (OxfordUniversity Press, Oxford, 2003).[18] J. P. Wittmer, H. Xu, P. Poli´nska, C. Gillig, J. Helfferich,F. Weysser, and J. Baschnagel, Eur. Phys. J. E , 131(2013).[19] The observable becomes T -independent at low temper-atures due to the prefactor β = 1 /T and system sizeindependent due to the prefactor V .[20] Using the slightly modified definitions ˜ µ γγ = µ γγ G eq ,˜ µ γτ = µ γτ and ˜ µ ττ = µ γτ /G eq , one obtains the generaltransform ˜ µ ab | λ = ˜ µ ab | λ =0 − λ for all fluctuations. Simi-larly, the rescaling ˜ C γγ ( t ) = C γγ ( t ) G eq , ˜ C γτ ( t ) = C γτ ( t )and ˜ C ττ ( t ) = C ττ ( t ) /G eq leads to the compact transfor-mation relation ˜ C ab ( t ) | λ = ˜ C ab ( t ) | λ =0 − λ for the corre-lation functions. Equation (7) may be generalized asdd λ ˜ µ ab | λ = dd λ ˜ C ab ( t ) | λ = − λ = 0.[21] D. R. Squire, A. C. Holt, and W. G. Hoover, Physica ,388 (1969).[22] H. Mizuno, S. Mossa, and J.-L. Barrat, Phys. Rev. E ,042306 (2013).[23] J.-L. Barrat, J.-N. Roux, J.-P. Hansen, and M. L. Klein,Europhys. Lett. , 707 (1988).[24] J. P. Wittmer, A. Tanguy, J.-L. Barrat, and L. Lewis,Europhys. Lett. , 423 (2002).[25] A. Tanguy, J. P. Wittmer, F. Leonforte, and J.-L. Barrat,Phys. Rev. B , 174205 (2002).[26] E. Flenner and G. Szamel, Phys. Rev. Lett. , 105505(2015).[27] M. Doi and S. F. Edwards, The Theory of Polymer Dy-namics (Clarendon Press, Oxford, 1986).[28] This assumes that the system is ergodic and does notapply if one samples, e.g., configuration ensembles withquenched shear-strain.[29] As seen in fig. 15, in some cases rather long trajecto-ries are needed to determine clearly a relaxation time θ ab from the 1 /t -decay of y ( t ) ≡ − µ ab ( t ) /µ ab . As shown inref. [10] one may instead characterize a relaxation timeusing y ( t = θ f ) = f with f (cid:28) f = 1 /
100 [10]. If f is sufficiently small, eq. (27)implies θ ab = θ f f/ , 178301 (2012).[31] H. Goldstein, J. Safko, and C. Poole, Classical Mechanics (Addison-Wesley, 2001), 3nd edition.[32] A prime denotes generally a derivative of a function withrespect to its argument.[33] J. F. Lutsko, J. Appl. Phys , 2991 (1989).[34] For each x one only takes into account (ˆ γ, ˆ τ ) with ˆ γ cor-responding to a finite δx -bin around x .[35] A good fit for the short time behavior is given by theGaussian C ττ ( t ) = µ ττ exp( − ( t/ ˜ θ ) /
2) [10]. To leadingorder this yields eq. (67).[36] See ref. [10] for the determination of G (cid:48) ( ω ) and G (cid:48)(cid:48) ( ω )using an oscillatory shear strain of frequency ωω