Compressibility and pressure correlations in isotropic solids and fluids
J.P. Wittmer, H. Xu, P. Polinsak, C. Gillig, J. Helfferich, F. Weysser, J. Baschnagel
aa r X i v : . [ c ond - m a t . s t a t - m ec h ] A ug Compressibility and pressure correlations in isotropic solids and fluids
J.P. Wittmer, ∗ H. Xu, P. Poli´nska, C. Gillig, J. Helfferich, F. Weysser, 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 FMF, University of Freiburg, Stefan-Meier-Str. 21, D-79104 Freiburg, Germany Theoretical Polymer Physics, University of Freiburg,Hermann-Herder-Str. 3, D-79104 Freiburg, Germany (Dated: September 17, 2018)Presenting simple coarse-grained models of isotropic solids and fluids in d = 1, 2 and 3 dimensionswe investigate the correlations of the instantaneous pressure and its ideal and excess contributionsat either imposed pressure (NPT-ensemble, λ = 0) or volume (NVT-ensemble, λ = 1) and for moregeneral values of the dimensionless parameter λ characterizing the constant-volume constraint. Thestress fluctuation representation F Row | λ =1 of the compression modulus K in the NVT-ensemble isderived directly (without a microscopic displacement field) using the well-known thermodynamictransformation rules between conjugated ensembles. The transform is made manifest by computingthe Rowlinson functional F Row also in the NPT-ensemble where F Row | λ =0 = Kf ( x ) with x = P id /K being a scaling variable, P id the ideal pressure and f ( x ) = x (2 − x ) a universal function. By graduallyincreasing λ by means of an external spring potential, the crossover between both classical ensemblelimits is monitored. This demonstrates, e.g., the lever rule F Row | λ = K [ λ + (1 − λ ) f ( x )]. PACS numbers: 05.20.Gg,05.70.-a,61.20.Ja,65.20.-w
I. INTRODUCTION
Strain and stress fluctuations.
Among the fundamen-tal properties of any equilibrium system are its (general-ized) elastic constants characterizing the fluctuations ofits extensive and/or conjugated intensive variables [1–9].For instance, for an isotropic solid or fluid the volumeand density fluctuations are set by the isothermal com-pression modulus K defined as [2] K ≡ V ∂ F∂V (cid:12)(cid:12)(cid:12)(cid:12) T ≡ − V ∂P∂V (cid:12)(cid:12)(cid:12)(cid:12) T = ρ ∂P∂ρ (cid:12)(cid:12)(cid:12)(cid:12) T (1)with V being the volume, N the particle number, ρ = N/V the particle density, F ( T, V ) the free energy, P the(mean) pressure and T the (mean) temperature. Albeit K may in principle be measured by fitting the pressureisotherm P ( ρ, T ) [10], it is from the computational pointof view important [11, 12] that this modulus may be ob-tained from the volume fluctuations at constant pres-sure (NPT-ensemble) and the pressure fluctuations atconstant volume (NVT-ensemble) evaluated at the same state point, i.e. at the same mean temperature, densityand pressure [2]. In the NPT-ensemble K is obtainedfrom the fluctuations δ ˆ V = ˆ V − V of the instantaneousvolume ˆ V around its mean value V = h ˆ V i using [2] K = F vol | ≡ k B T V / h δ ˆ V i (cid:12)(cid:12)(cid:12) (2) ∗ Electronic address: [email protected] with k B being Boltzmann’s constant and where | indi-cates that the average is obtained in the NPT-ensemble( λ = 0 with λ defined below). Equivalently, K may beobtained in a canonical NVT-ensemble using the stressfluctuation formula [3–5, 10, 11, 13] K = F Row | ≡ P + η Born | − η F , ex | (3)with | indicating the NVT-ensemble ( λ = 1) and otherdefinitions given immediately below. This formula wasfirst stated for liquids in the 1950s by Rowlinson [3]and later implicitly rediscovered by Squire, Hold andHoover [4] formulating the stress fluctuation formalismfor anisotropic solids as summarized in appendix C. Affine contribution.
The second term η Born in eq. (3)represents the Born approximation [14] for the interac-tion energy implied by an imposed infinitesimal strain as-suming affine microscopic particle displacements [5, 6, 9].As reminded in appendix A, for pairwise additive poten-tials this “Born-Lam´e coefficient” becomes [11] η Born ≡ d V *X l h ( r l ) + with h ( r ) ≡ r ( r u ′ ( r )) ′ (4)and d being the spatial dimension, l an index labelingthe interactions, r l the distance between two interactingparticles, u ( r ) the pair potential and a prime denoting aderivative with respect to the indicated variable. As sug-gested by eq. (4), η Born is sometimes also called “hyper-virial” [11]. We note en passant that η Born is a momentof the second derivative of u ( r ) and some care is requiredif η Born is computed using a truncated potential [15, 16].
Non-affine contribution.
In general the Born approx-imation overpredicts the free-energy change. The over-prediction is “corrected” by the stress fluctuation term η F , ex | ≡ βV D δ ˆ P E(cid:12)(cid:12)(cid:12) ≥ β = 1 /k B T being the inverse temperature and ˆ P ex the instantaneous excess pressure which for pairwise ad-ditive potentials is given by Kirkwood’s virial [11]ˆ P ex ≡ d ˆ V X l r l f l with f l = − u ′ ( r l ) (6)being the central force between two interacting parti-cles. (Although ˆ P ex is used here in Rowlinson’s formula,eq. (3), at constant volume V = ˆ V , we have written itin a slightly more general form which is necessary if vol-ume fluctuations are allowed.) In particular from ref. [5]it has become clear that stress fluctuation corrections,such as η F , ex | , do not necessarily vanish for T →
0. Aswe shall also illustrate in the present paper, this is due tothe fact that the particle displacements need not followan imposed macroscopic strain affinely [6, 17–21]. Howimportant the non-affine motions are, depends on thesystem under consideration [20]. While the elastic prop-erties of crystals with one atom per unit cell are givenby the Born term only, stress fluctuations are significantfor crystals with more complex unit cells [5]. They be-come pronounced for polymer-like soft materials [10] andamorphous solids [6, 8, 17–25].
Fluctuations in conjugated ensembles.
Focusing onthe compression modulus K we emphasize in this reportthat the numerically more convenient stress fluctuationformalism may be obtained directly using the well-knownthermodynamic transformation rules between conjugatedensembles [26, 27]. This point is crucial if the formalismis used in situations where no meaningful microscopicdisplacement field can be defined [13, 28]. ComputingRowlinson’s F Row for NPT-ensembles the general trans-form behind the formalism can be made manifest. Elab-orating a short comment [29], we show that F Row | ≡ P + η Born | − η F , ex | = K f ( x ) with f ( x ) = x (2 − x ) (7)being the universal scaling function, x ≡ P id /K the scal-ing variable and P id the ideal pressure contribution. Generalized λ -ensembles. It is straightforward to in-terpolate between the NPT- and the NVT-ensemble byimposing an external spring potential U ext ( ˆ V ) = k ext (cid:16) ˆ V − V ext (cid:17) with K ext ≡ V k ext (8)being the associated compression modulus introduced forconvenience [30]. (Our approach is conceptually simi-lar to the so-called “Gaussian ensemble” proposed someyears ago by Hetherington and others [31, 32] generaliz-ing the Boltzmann weight of the canonical ensemble byan exponential factor U ext ( ˆ E ) ∝ ˆ E of the instantaneousenergy ˆ E .) Throughout this work it is assumed that K ext ≥
0, i.e. U ext ( ˆ V ) reduces the volume fluctuations[33]. Chosing the reference volume V ext equal to the aver-age volume V of the isobaric system at imposed P allows,for symmetry reasons, to work at constant mean pressureirrespective of the strength of the external potential [34].The volume fluctuations may then be characterized us-ing the dimensionless parameter λ ≡ K ext / ( K + K ext ).(Since K ext ≥
0, we have 0 ≤ λ ≤ λ →
0, while NVT-statistics should become relevant inthe opposite limit for K ext → ∞ and λ →
1. We shallmonitor various properties, such as the Rowlinson for-mula F Row | λ , as a function of λ and x = P id /K . Wedemonstrate, e.g., the simple lever rule F Row | λ = K [ λ + (1 − λ ) f ( x )] (9)which generalizes eq. (3) and eq. (7) to arbitrary λ . Outline.
In sect. II we consider theoretically vari-ous correlation functions of normal stress contributionsin different λ -ensembles. We begin by summarizing insect. II A the transformation relations for fluctuations be-tween NVT- and NPT-ensembles. Equation (7) is derivedin sect. II D. We turn then in sect. II E to the transfor-mation relations for general λ and demonstrate eq. (9)in sect. II F. Our Monte Carlo (MC) and molecular dy-namics (MD) simulations of several simple coarse-grainedmodels in d = 1, 2 and 3 dimensions are described insect. III. Our theoretical predictions are then checkednumerically in sect. IV. Several well-known but scatteredtheoretical statements are gathered in the appendix. II. THEORETICAL CONSIDERATIONSA. Fluctuations in NVT- and NPT-ensembles
As discussed in the literature [2, 11, 26], a simple av-erage A = h ˆ A i of an observable A does not depend onthe chosen ensemble, at least not if the system is largeenough. (We do thus not indicate normally in whichensemble the average has been taken.) However, a cor-relation function h δ ˆ Aδ ˆ B i of two observables A and B depends on whether V or P are imposed. As shown firstby Lebowitz et al. in 1967 [26], one verifies that [27, 35] D δ ˆ Aδ ˆ B E(cid:12)(cid:12)(cid:12) = D δ ˆ Aδ ˆ B E(cid:12)(cid:12)(cid:12) − KβV ∂A∂P ∂B∂P (10)to leading order. We note that the left hand-side ofeq. (10) must vanish if at least one of the observablesis a function of ˆ V . In this case we have D δ ˆ Aδ ˆ B E(cid:12)(cid:12)(cid:12) = KβV ∂A∂P ∂B∂P . (11)One verifies for ˆ A = ˆ B = ˆ V that eq. (11) is consistentwith eq. (2). For ˆ A = ˆ V and ˆ B = ˆ P one obtains imme-diately the well-known relation [2] − β D δ ˆ V δ ˆ P E(cid:12)(cid:12)(cid:12) = 1 . (12)Similarly, one obtains for ˆ A = ˆ B = ˆ V n that D δ ( ˆ V n ) E(cid:12)(cid:12)(cid:12) = KβV ∂ h ˆ V n i ∂P ! ≈ βK n V n − (13)where the steepest-descent approximation D ˆ V n E /n ≈ D ˆ V E = V (14)for simple averages has been made and V /K = − ∂V /∂P is used again. For the fluctuations of the inverse vol-ume 1 / ˆ V the latter result ( n = −
1) may be rewrittencompactly using eq. (2) as V D δ (1 / ˆ V ) E(cid:12)(cid:12)(cid:12) D δ ˆ V E(cid:12)(cid:12)(cid:12) = 1 (15)where we have changed ≈ to the equal sign for largesystems. That eq. (13) and eq. (15) become exact for V → ∞ can be also seen by using that the distributionof ˆ V in the NPT-ensemble is Gaussian. With ˆ A = ˆ V n and ˆ B = ˆ P one gets similarly D δ ( ˆ V n ) δ ˆ P E(cid:12)(cid:12)(cid:12) = KβV ∂ h ˆ V n i ∂P ≈ − nV n − k B T (16)to leading order for V → ∞ using the same approxima-tion as above. With eq. (12) this gives for n = − − V D δ (1 / ˆ V ) δ ˆ P E(cid:12)(cid:12)(cid:12) D δ ˆ V δ ˆ P E(cid:12)(cid:12)(cid:12) = 1 . (17)The cumulants eqs. (15,17) have been used in the com-putational part of our work to check the precision of thebarostat and to verify whether our configurations are suf-ficiently large for the investigated state point. B. Transformation of pressure auto-correlations
Returning to eq. (10), this implies for ˆ A = ˆ B = ˆ P thetransformation of the pressure fluctuations βV D δ ˆ P E(cid:12)(cid:12)(cid:12) = βV D δ ˆ P E(cid:12)(cid:12)(cid:12) − K, (18)i.e. K may be obtained by measuring the pressure fluc-tuations in both ensembles. As we shall show in para-graph II C, the numerically more convenient Rowlinsonexpression F Row | for K can be derived directly from eq. (18) [13]. In the following we use the more concise no-tation η F ≡ βV h δ ˆ P i for the pressure fluctuations. η F | is also called the “affine dilatational elasticity” η A [13] forreasons which will become obvious in sect. II C. (See alsoappendix A.) Since K > η A ≡ η F | > η F | . Depending on the disorder, η F | is, however, not a negligible contribution. To seethis let us remind that K can also be determined fromthe ( ˆ V , ˆ P )-data measured in an NPT-ensemble using thelinear regression relation K = F reg | ≡ − V D δ ˆ V δ ˆ P E(cid:12)(cid:12)(cid:12) / D δ ˆ V E(cid:12)(cid:12)(cid:12) . (19)Please note that using eq. (12) this reduces to eq. (2).Associated to F reg | is the dimensionless regression coef-ficient [36] c reg | ≡ − D δ ˆ V δ ˆ P E(cid:12)(cid:12)(cid:12) / r D δ ˆ V E(cid:12)(cid:12)(cid:12) D δ ˆ P E(cid:12)(cid:12)(cid:12) (20)which can be also further simplified using eq. (12). In-terestingly, using eq. (2) and eq. (18) one sees that c reg | = p K/η A = q − η F | /η A , (21)i.e. the regression coefficient obtained at constant pres-sure determines the pressure fluctuations at constant vol-ume. Only if the measured ( ˆ V , ˆ P ) are perfectly corre-lated, i.e. c reg | = 1, this implies K = η A and η F | = 0.In fact, for all non-trivial systems one always has c reg | < , thus K < η A and η F | > , (22)i.e. the affine dilatational elasticity η A sets only an upperbound to the compression modulus and a theory whichonly contains the affine response must overpredict K . C. Rowlinson’s formula rederived
MC-gauge.
There is a considerable freedom for defin-ing the instantaneous value of the pressure ˆ P = ˆ P id + ˆ P ex as long as its average P = P id + P ex does not change [11].It is convenient for the subsequent derivations and forour MC simulations (and not in conflict with the alsopresented MD simulations) to define the instantaneousideal pressure ˆ P id by [11]ˆ P id = k B T N/ ˆ V (MC-gauge). (23)Within this “MC-gauge” the thermal momentum fluctu-ations are assumed to be integrated out. This leads tothe usual prefactorˆ V N = exp h − β ( − k B T N log( ˆ V )) i (24)of the remaining partition function. The effective Hamil-tonian H s ( ˆ V ) of a state s of the system thus reads H s ( ˆ V ) = − k B T N log( ˆ V ) + U s ( ˆ V ) + const (25)where the first term on the right hand-side refers to theintegrated out momenta and the second to the total ex-cess potential energy U s ( ˆ V ) expressed as a function ofthe instantaneous volume as shown in appendix A. Non-affine contribution.
An immediate consequenceof the MC-gauge is, of course, that the fluctuations ofˆ P id vanish for the NVT-ensemble and that, hence, η F | = η F , ex | ≡ βV D δ ˆ P E(cid:12)(cid:12)(cid:12) (26)with the instantaneous excess pressure ˆ P ex being com-puted using Kirkwood’s expression, eq. (6). Accordingto eq. (21) the correlation coefficient c reg is thus a func-tion of η F , ex | /η A and vice versa . Affine (Born) contribution.
The task is now to com-pute the pressure fluctuation in the NPT-ensemble. Wenote first for the NPT-weight of a configuration at vol-ume ˆ V that h e − β ( H s ( ˆ V )+ P ˆ V ) i ′ = β ( ˆ P − P ) e − β ( H s ( ˆ V )+ P ˆ V ) (27)where ˆ P ≡ −H ′ s ( ˆ V ) ( ˆ P -gauge) (28) defines the instantaneous total pressure ˆ P = ˆ P id + ˆ P ex .Note that this definition is consistent with the MC-gauge,eq. (23), for ˆ P id and eq. (25). As shown in appendix B,ˆ P ex ≡ − U ′ s ( ˆ V ) is also consistent with the Kirkwood ex-cess pressure, eq. (6), for pair potentials. Using eq. (27)the second moment h ( ˆ P − P ) i can be readily obtainedby integration by parts. This yields η A ≡ η F | = − V D ˆ P ′ ( ˆ V ) E(cid:12)(cid:12)(cid:12) = V D H ′′ s ( ˆ V ) E(cid:12)(cid:12)(cid:12) (29)where apriori the average is understood to be taken overall states s of the system and all volumes ˆ V at imposed P .(The boundary terms for the integration by parts over ˆ V can be neglected for sufficiently large systems since theNPT-weight for ˆ V gets strongly peaked around V .) Itis of importance that the fluctuation η F | has thus beenreduced to a simple average . This allows its computationmore conveniently by NVT-ensemble simulations. Usingeq. (25) it follows further that η A = η A , id + η A , ex with η A , id ≡ V D − k B T N log( ˆ V ) ′′ E(cid:12)(cid:12)(cid:12) ≈ P id (30) η A , ex ≡ V D U ′′ s ( ˆ V ) E(cid:12)(cid:12)(cid:12) (31)where we have used for the ideal contribution that for suf-ficiently large systems V h / ˆ V i ≈ h / ˆ V i ≈ /V . As seenfrom the affine excess energy discussed in appendix A, itfollows for pair potential interactions that η A , ex = η Born + P ex , (32)i.e. both coefficients η A , ex and η Born are equivalent. Westress again that P id , P ex , η Born and η A , ex are simpleaverages and can thus be evaluated readily in both en-sembles using eq. (14). Substituting eqs. (26,31,32) intothe transform eq. (18), this finally confirms eq. (3). D. Correlations at constant P We focus now on the fluctuations of the pressure con-tributions in the NPT-ensemble. According to eq. (32)we have P + η Born = η A ≡ η F | . If the Rowlinson func-tional F Row is measured at imposed P this implies F Row | = η F | − η F , ex | = βV D δ ˆ P E(cid:12)(cid:12)(cid:12) + 2 βV D δ ˆ P id δ ˆ P ex E(cid:12)(cid:12)(cid:12) . (33)As a next step we demonstrate the relations η F , id | ≡ βV D δ ˆ P E(cid:12)(cid:12)(cid:12) = K x (34) η F , mix | ≡ βV D δ ˆ P id δ ˆ P ex E(cid:12)(cid:12)(cid:12) = K x (1 − x ) (35)(with x = P id /K being again the scaling variable)from which eq. (7) is then obtained by substitution intoeq. (33). Remembering eq. (23), eq. (34) is obtained fromthe fluctuations of the inverse volume, eq. (13). The rela-tion eq. (35) describing the coupling of ideal and excesspressure is implied by eq. (16) for n = − P = ˆ P id + ˆ P ex and eq. (34). We emphasize that eq. (7) oreq. (33) do not completely vanish for finite T , i.e. finite x , as does the corresponding stress fluctuation expressionfor the shear modulus G at imposed shear stress τ [13].Having thus demonstrated eq. (34) and eq. (35) andusing the already stated relation for the total pressure,eq. (29), one confirms finally for the fluctuations of theexcess pressure ˆ P ex that η F , ex | ≡ βV D δ ˆ P E(cid:12)(cid:12)(cid:12) = η A , ex − K x (1 − x ) (36)with x = P id /K being the reduced ideal pressure. E. Fluctuations in different λ -ensembles Introduction.
The transformation relations for sim-ple averages and fluctuations between the standard con-jugated ensembles have been given in the 1960s byLebowitz, Percus and Verlet [26]. Focusing on the volume V as the only relevant extensive variable and the conju-gated (reduced) pressure βP of the system as the onlyintensive variable [27, 35] we rederive now their resultsfor generalized λ -ensembles. We illustrate several pointsmade by the histograms presented in fig. 1 for simple 1Dnets of harmonic springs as described in sect. III below. Histogram.
Let us assume a normalized distribution p ( ˆ V ) with a sharp and symmetric maximum at the meanvolume V . Examples for such histograms are given inpanel (a) of fig. 1 for several values λ = K ext / ( K ext + K ).The standard NPT-ensemble corresponds to the value λ = 0. Note that the volume V of the NPT-ensemble at P = 0 is taken as the reference volume V ext of the exter-nal spring potential. Due to this choice neither the meanvolume V nor the pressure P do change with increasing -0.01 0.00 0.01 x = V/V-1 -1 p ( x ) k ext = λ =0 (NPT)k ext =0.0001, λ =0.091k ext =0.001, λ =0.499k ext =0.01, λ =0.917k ext =0.1, λ =0.992 nets of springs:d=1, P=0, δ k=0, T=0.01 ^ (a) -0.01 0.00 0.01 x = V/V-1 -0.010.000.01 δ P ( x ) = P ( x ) - P ( x = ) k ext = λ =0k ext =0.0001, λ =0.091k ext =0.001, λ =0.499k ext =0.01, λ =0.917k ext =0.1, λ =0.992 λ < δ V > d=1, δ k=0,P=0eq.(48) ^ ^ (b) c o ll a p s o f a ll d a t a w i t h sa m e K e q . ( ) FIG. 1: Simple 1D nets without noise ( δk = 0) at tempera-ture T = 0 .
01 and pressure P = 0 (mean volume V = 10 ,number density ρ = 1) for different values λ of the exter-nal spring as indicated: (a) Normalized histograms p ( x ) ofthe rescaled instantaneous volume x = ˆ V /V −
1. The dis-tributions are Gaussian as shown by the line indicated forthe standard NPT-ensemble ( λ = 0). (b) Pre-averaged pres-sure δP ( ˆ V ) traced as a function of the rescaled instantaneousvolume x . The bold line shows eq. (37) for the compressionmodulus K = 1 of the system. Inset: The mean-squared vol-ume fluctuations δV λ decrease linearly with λ according toeq. (44) or eq. (48) as shown by the solid line. λ , i.e. all ensembles correspond to the same thermody-namic state. All distributions are Gaussian, becomingsharper with increasing λ . (The limit λ → δ -function.) The width of the distributionaround its maximum being characterized by the param-eter λ , we introduce the notation δV λ ≡ h δ ˆ V i for themean-squared width. The width δV λ of simple 1D netsis presented in the inset of panel (b) of fig. 1. Observables.
We write ˆ A ( ˆ V ) for the expectation valueof an observable A at a given state s of the (full) phasespace of a system at volume ˆ V . Let us first focus onthe average A ( ˆ V ) of ˆ A ( ˆ V ) taken over all states s . As anexample we show in panel (b) of fig. 1 the mean pressure P ( ˆ V ) averaged over all systems s of a λ -ensemble foundat a given instantaneous volume ˆ V . As one expects, P ( ˆ V ) is seen to decrease linearly to leading order, δP ( ˆ V ) ≡ P ( ˆ V ) − P ( V ) ≈ − KV δ ˆ V , δ ˆ V ≡ ˆ V − V, (37)as stressed by the bold line. Importantly, P ( ˆ V ) does not depend on λ . (As one expects, the statistics deterioratesfor ( ˆ V − V ) /δV λ ≫ A ( ˆ V ) around V yields A ( ˆ V ) = a + a δ ˆ V + 12 a δ ˆ V + . . . (38)where a coefficient a n denotes the n th derivative of A ( ˆ V )with respect to ˆ V taken at ˆ V = V . Simple averages.
The average taken over all properlyweighted volumes of the investigated λ -ensemble thus be-comes to leading order D A ( ˆ V ) E = Z d ˆ V p ( ˆ V ) A ( ˆ V ) ≈ a + a δV λ (39)with a = A ( V ). We have used here that p ( ˆ V ) is sym-metric around V which implies that the linear term ineq. (38) must drop out. The difference of a simple av-erage taken at λ < λ = 1) is then given by the second derivative a timesthe mean-squared width δV λ which increases at most lin-early with V for λ = 0. Since a does not depend onthe ensemble, the correction must decay at least as fast as for the NPT-ensemble discussed in the literature [26].Hence, simple averages computed for any ensemble with0 ≤ λ ≤ Averaged fluctuations.
As a next step let us considerthe correlation function h δAδB i = h AB i − h A ih B i of theobservables A ( ˆ V ) and B ( ˆ V ). Using again the symmetryof the distribution p ( ˆ V ) around V it is seen that h δAδB i = h a b i − h a ih b i + a b δV λ (40)where the coefficients a n and b n stand for the respectivederivatives of A ( ˆ V ) and B ( ˆ V ) taken at ˆ V = V used forthe Taylor expansion around V . The underlined termvanishes since a = A ( V ) and b = B ( V ) are constants.Note also that h δAδB i → λ →
1. This is due tothe fact that we have correlated here the pre-averaged observables A ( ˆ V ) and B ( ˆ V ) instead of the expectationvalues ˆ A and ˆ B which depend not only on the volumeˆ V but also on the state s of the system. Replacing thus A → ˆ A , B → ˆ B , a n → ˆ a n and b n → ˆ b n we have thusin addition to average over all possible states and theunderlined term in eq. (40) remains thus finite h δ ˆ Aδ ˆ B i = h δ ˆ a δ ˆ b i + h ˆ a ˆ b i δV λ . (41)Since h ˆ a ˆ b i is a simple average, this result can be refor-mulated using a more natural notation as D δ ˆ Aδ ˆ B E(cid:12)(cid:12)(cid:12) λ = D δ ˆ Aδ ˆ B E(cid:12)(cid:12)(cid:12) + ∂A∂V ∂B∂V δV λ (42)where | λ indicates that the average is taken over properlyweighted volumes in a general λ -ensemble and | indicatesthe NVT-average ( λ = 1) taken at the maximum of thedistribution p ( ˆ V ). Using eq. (42) for the NPT-limit ( λ =0) this is seen to be identical to eq. (2.11) given in ref. [26].Substracting this reference from the general λ case yields D δ ˆ Aδ ˆ B E(cid:12)(cid:12)(cid:12) λ = D δ ˆ Aδ ˆ B E(cid:12)(cid:12)(cid:12) − (cid:0) δV − δV λ (cid:1) ∂A∂V ∂B∂V . (43)Equation (43) can be further simplified using δV − δV λ = VβK − Vβ ( K + K ext ) = V λβK (44)for the difference of the volume fluctuations in both en-sembles. We have used here eq. (2) for the NPT-ensembleand the fact that for general λ the external spring is par-allel to the system, i.e. the effective modulus must be thesum of the system modulus K and spring modulus K ext .That eq. (44) holds is confirmed by the data presented inthe inset of panel (b) of fig. 1. Following ref. [26] we alsorewrite the derivatives with respect to the volume V asderivatives with respect to the pressure P of the system ∂A∂V ∂B∂V = (cid:18) ∂P∂V (cid:19) ∂A∂P ∂B∂P = (cid:18) KV (cid:19) ∂A∂P ∂B∂P (45)where we have used ∂P/∂V = − K/V in the last step.Using the above three equations this yields finally h δ ˆ Aδ ˆ B i (cid:12)(cid:12)(cid:12) λ = h δ ˆ Aδ ˆ B i (cid:12)(cid:12)(cid:12) − KλβV ∂A∂P ∂B∂P , (46)which compares the correlations in a general λ -ensemble(0 ≤ λ ≤
1) with the correlations in an NPT-ensemble( λ = 0). Note that for λ = 1 this is consistent with theoriginal transformation, eq. (10), derived in ref. [26]. F. Correlations for generalized λ -ensembles Using eq. (46) we restate first several correlations givenabove where at least one of the observables ˆ A and ˆ B is a function of the instantaneous volume ˆ V . Since for λ < λ = 0 or λ = 1 are used.For ˆ A = ˆ B = ˆ V this yields, e.g., h δ ˆ V i (cid:12)(cid:12)(cid:12) λ = h δ ˆ V i (cid:12)(cid:12)(cid:12) − KλβV (cid:18) ∂V∂P (cid:19) (47)= δV (1 − λ ) (48)= k B T VK + K ext (49)restating thus eq. (44). More generally, one confirms forˆ A = ˆ B = ˆ V n that h δ ( ˆ V n ) i (cid:12)(cid:12)(cid:12) λ = n V n − δV (1 − λ ) . (50) The latter result is consistent with eq. (15) which is thusshown to hold for all λ <
1. For ˆ A = ˆ V and ˆ B = ˆ P it isseen that − β D δ ˆ V δ ˆ P E(cid:12)(cid:12)(cid:12) λ = − β D δ ˆ V δ ˆ P E(cid:12)(cid:12)(cid:12) + KλV ∂V∂P = 1 − λ (51)generalizing thus eq. (12). More generally, one sees forˆ A = ˆ V n and ˆ B = ˆ P that − β D δ ( ˆ V n ) δ ˆ P E(cid:12)(cid:12)(cid:12) λ = nV n − (1 − λ ) (52)as expected from eq. (16). One confirms using eq. (52)that eq. (17) must hold for all λ <
1. Interestingly, thealready mentioned linear regression formula F reg is seenusing eq. (48) and eq. (51) to become F reg | λ ≡ − V D δ ˆ V δ ˆ P E(cid:12)(cid:12)(cid:12) λ / D δ ˆ V E(cid:12)(cid:12)(cid:12) λ = k B T VδV = K (53)independent of the ensemble used. As an alternative tothe strain fluctuation relation eq. (49), this allows thusthe determination of K for all λ <
1. The dimension-less correlation coefficient c reg associated to F reg dependshowever on λ c reg | λ = s Kη A s − λ − λK/η A (54)where we have used eq. (55) demonstrated below. Notethat K/η A = 1 implies (as before) c reg | λ = 1 for all λ <
1. For
K/η A < p K/η A at λ = 0 tozero for λ = 1. As one would expect, this shows that themore the volume fluctuations are suppressed by the exter-nal constraint, the more ˆ V and ˆ P must decorrelate. Forthe transformation of the total pressure auto-correlation η F ≡ βV h δP i , eq. (46) simply implies that η F | λ = η F | − Kλ = η A − KK ext K + K ext . (55)Since η A = P id + η A , ex is a simple average computablein any ensemble, eq. (55) may be also used for the de-termination of K . Using again the MC-gauge and thatˆ P = ˆ P id + ˆ P ex , it is seen from eq. (55) that the Rowlinsonfunctional F Row | λ ≡ P + η Born − η F , ex | λ should become F Row | λ = Kλ + η F , id | λ + 2 η F , mix | λ (56)generalizing thus the Rowlinson formula, eq. (3), for λ =1 (remembering that η F , id and η F , mix must vanish in theNVT-ensemble) and eq. (33) for λ = 0. Considering asa next step the fluctuations of the ideal pressure η F , id inthe MC-gauge, eq. (50) implies that η F , id | λ ≡ βV D δ ˆ P E(cid:12)(cid:12)(cid:12) λ = K x (1 − λ ) . (57)Since ˆ P = P id + P ex the latter result together with eq. (52)allows to generalize eq. (58) for the correlations betweenideal and excess pressure contributions. This shows that η F , mix | λ ≡ βV D δ ˆ P id ˆ P ex E(cid:12)(cid:12)(cid:12) λ = K x (1 − x ) (1 − λ ) . (58)Our central result eq. (9) announced in the Introduc-tion is thus simply obtained by substituting eq. (57) andeq. (58) into eq. (56). Finally, using eq. (55) togetherwith eq. (57) and eq. (58) one verifies for the fluctuationsof the excess pressure fluctuations that η F , ex | λ = η A − F Row | λ (59)as obvious from eq. (9). The latter result reduces tothe Rowlinson expression eq. (3) for λ = 1 and using η A = η A , ex + Kx to eq. (36) for λ = 0. III. SOME ALGORITHMIC DETAILS
Introduction.
In order to check our predictions wesampled by MC and MD simulation [11, 12] variousmodel systems for solids and glass-forming liquids in d = 1, 2 and 3 dimensions. Periodic boundary conditionsare used and all systems are first kept at constant pres-sure P using standard barostats [11] as specified below.After equilibrating and sampling in the NPT-ensemble,the volume fluctuations are suppressed either by impos-ing V = ˆ V or by means of a finite spring potential. Var-ious simple averages, such as the Born-Lam´e coefficient η Born , and fluctuations, such as the excess pressure fluc-tuation η F , ex or the Rowlinson functional F Row , are com-pared for all available λ . As expected, all simple averagesare found within numerical accuracy to be identical. Asdiscussed in sect. IV, fluctuations are found to transformfollowing the predictions based on eq. (10) and eq. (46).Specifically, we have verified that eqs. (15,17,51) holdto high precision for all λ . About N = 10 particlesare typically used. Lennard-Jones (LJ) units are usedthroughout this work and k B is set to unity. One-dimensional spring model.
As already seen infig. 1, the bulk of the presented numerical results hasbeen obtained by MC simulation of permanent nets ofideal harmonic springs in strictly d = 1 dimension [37].We use a potential energy U = 12 X l k l ( x l − R l ) (60)with x l being the distance between the connected parti-cles. The reference length R l of the springs is assumedto be constant, R l = R = 1, and the spring constants k l are taken randomly from a uniform distribution of half-width δk centered around a mean value also set to unity.We note for later reference that this implies1 h /k l i = 2 δk log((1 + δk ) / (1 − δk )) (61)as indicated by the bold line in fig. 2. Only simple net-works are presented here where two particles i − i along the chain are connected by one spring l = i , i.e. allforces f l along the chain are on average identical. Glass-forming liquids.
Our two-dimensional (2D) sys-tems are polydisperse Lennard-Jones (pLJ) beads as de-scribed in ref. [13]. The reported three-dimensional (3D)systems refer to binary Kob-Andersen (KA) mixtures [38]sampled by means of MD simulations taking advantageof the LAMMPS implementation [39]. Starting from theliquid limit well above the glass transition temperature T g , both system classes have been quenched [13] deepinto the glassy state at very low temperatures T ≪ T g asmay be seen from fig. 3 [16]. Barostats.
The results reported for d = 1 and d = 2have been obtained using local MC moves for the parti-cles and global MC moves for the barostat [13]. As de-scribed elsewhere [11, 13], an attempted volume change δ ˆ V = ˆ V new − ˆ V old is accepted if ξ ≤ exp( − βδG ) with ξ denoting a uniformly distributed random variable with0 ≤ ξ < δG = δE + δ ˆ V P − k B T N log( ˆ V new / ˆ V old ) . (62)The first term δE stands for the energy difference asso-ciated with the affine displacements of all particles, thesecond term δ ˆ V P imposes the pressure P and the log-arithmic contribution corresponds to the change of thetranslational entropy, i.e. the change of the integratedout momentum contribution discussed above, eq. (24).While a broad range of pressures has been sampled forthe 1D nets, the 2D pLJ beads have been kept at only onepressure, P = 2, for which a glass-transition temperature T g ≈ .
26 has been determined [13].For more general λ -ensembles the increment δU ext ( ˆ V )of the external spring defined in eq. (8) is simply added to δG . We assume throughout this work that V ext ≡ V forthe reference volume of the external spring with V beingthe mean volume of an NPT-ensemble of a given pressure P . Due to the symmetry of the fluctuations around V for the system and the external spring, this is sufficientto keep this pressure constant for all λ [34].For our MD simulations of the KA model we have usedthe Nos´e-Hoover barostat (“fix npt command”) providedby the LAMMPS code [39]. Following Kob and Andersen[38] a constant pressure P = 1 has been imposed for alltemperatures. This choice corresponds to glass-transitiontemperature T g ≈ . IV. COMPUTATIONAL RESULTS
As already noted above, various properties have beencompared for different boundary conditions as character-ized by the parameter λ while keeping the system at thesame state point. We focus first on the comparison ofNPT- ( λ = 0) and NVT-ensembles ( λ = 1) before weturn in sect. IV E to the more general λ -ensembles. δ k K F vol | F Row | η F | - η F | P id2 / η F,id | η F | η F | net d=1: N=10 , T=0.01, P=0 c u s p - s i ngu l a r i t y affine contribution n o n - a f f i n ec o n t r i b u t i o n e q . ( ) e q . ( ) eq.(66) FIG. 2: Compression modulus K for 1D nets at temperature T = 0 .
01 and pressure P = 0 computed using the rescaled vol-ume fluctuations F vol | (filled spheres), the Rowlinson stressfluctuation formula F Row | (crosses), the difference betweenthe total pressure fluctuations in both ensembles (squares)and the fluctuations of the inverse volume P / η F , id | (largespheres). The compression modulus decreases strongly with δk (solid line). Also indicated are the “affine” contribution η F | = 1+2 P id to the compression modulus, eq. (66), and the“non-affine” correction η F | which is seen to increase with δk according to eq. (67) as indicated by the dash-dotted line. -3 -2 -1 T K F vol | F Row | η F | - η F | P id2 / η F,id | net d=1: δ k=0, N=10 , P=0, ρ≅ , P=2KA d=3: N=6912, P=1 KA d a t a s h i ft e d FIG. 3: Compression modulus K vs. temperature T compar-ing again the values obtained using F vol | , F Row | , eq. (18)and eq. (34). The upper data refer to the results obtainedfor the KA model in 3D [13, 38] which have been shifted up-wards (factor 4) for clarity, the data in the middle to systemsof glass-forming 2D pLJ beads at P = 2 [13] and the lowerdata to simple 1D nets of harmonic springs at P = 0 withidentical spring constants ( δk = 0). A. Compression modulus K Comparing different methods.
As shown in figs. 2-4, the compression modulus K may be determined ei-ther using the volume fluctuations in the NPT-ensemble,eq. (2), or using Rowlinson’s stress fluctuation formula,eq. (3), for the NVT-ensemble (crosses). The same val-ues of K are obtained from the transform of the pressure -1.0 -0.5 0.0 0.5 1.0 P F vol | F Row | η F | - η F | K V ref /Vk=1, δ k=0 (no noise)1D net for T=0.01 r e f e r e n ce o ft e n u se d e q . ( ) FIG. 4: Compression modulus K vs. pressure P for 1D netswithout noise ( δk = 0) at T = 0 .
01. The compression modu-lus decreases with P simply since the volume V of the systemlinearly decreases as V = V ref (1 − P/K ref ) with K ref beingthe compression modulus at P = 0 where V = V ref . If K isrescaled with V ref /V (spheres) this yields a pressure indepen-dent material constant K ref . fluctuations, eq. (18), and from the ideal pressure fluctu-ations η F , id | , eq. (34), as indicated by the large sphereswhich thus confirms both relations. Harmonic spring networks.
As seen in fig. 2, the com-pression modulus of the 1D nets decreases strongly with δk . To understand the scaling indicated by the bold linelet us consider a chain of harmonic strings. Since the av-erage force acting on each spring is given by the imposedpressure P = − f ext /A , this implies an average length h x l i = R − P A/k l for a spring constant k l and an av-erage volume V = A P l h x l i . A pressure increment δP thus leads to a volume change δV = − AδP P l /k l . Forthe compression modulus K = − V δP/δV at pressure P this yields K = K ref − P with K ref = ρ ref / h /k l i (63)being the compression modulus of the unstressed refer-ence system at P = 0 and ρ ref = 1 /AR the correspond-ing density. The compression modulus K = K ref at zeropressure is thus inversely proportional to the average in-verse spring constant. For our uniformally distributedspring constants eq. (61) thus implies that K must van-ish as a cusp-singularity for δk → η F | to K (furtherdiscussed in sect. IV B) and the “non-affine” contribution η F | which is seen to increase with δk . The decrease of K is thus due to the increase of the non-affine contribution. Temperature dependence.
Figure 3 shows the temper-ature dependence of the compression modulus for ourthree model systems. As one would expect, K decreaseswith T for the bead systems kept at same pressure reach-ing the ideal gas compressibility K = K id = P id = P forlarge temperatures. As predicted by eq. (63) for strictlyharmonic spring chains, the compression modulus K isfound T -independent for all 1D networks and this irre-spective on the values δk and P sampled (not shown). Pressure dependence.
The pressure dependence of thecompression modulus for 1D nets is investigated in fig. 4.The dashed line indicating the linear behavior predictedby eq. (63) for chains of harmonic springs perfectly fitsthe measured data points. A comment is in order here:Since the springs are permanently fixed and ideal, the intrinsic mechanical properties of the systems do notchange with the external load. In fact, the second deriva-tive of the free energy F ( T, V ) with respect to the volume V becomes constant for these highly idealized systems KV ≡ ∂ F ( T, V ) ∂V = const = K ref V ref (64)with V ref being the reference volume of the system atzero pressure. The difference between the thermody-namic compression modulus K and the constant K ref simply arises since in all thermodynamic relations, suchas eq. (1), the volume V of the current state is takento make the modulus system-size independent (inten-sive) and not a reference volume at a certain pressure P . Since for an idealized elastic body as our 1D nets V = V ref (1 − P/K ref ), this implies the linear relation K = K ref − P (dashed line). Note that for solids the ap-plied pressures are for once small compared to the moduliand the systems can be taken linearly and without (or atleast with negligible) plastic rearrangement from a zero-pressure reference to a finite P . It is thus possible todefine elastic material properties such as K ref which are independent of the applied external stresses. If possible,this is a conceptionally and practically very useful sepa-ration of conditions and effects. However, it is obviouslypointless to take the zero-pressure volume (or the vol-ume of any other pressure) as a reference for a colloidalor polymeric system from which other pressure states aredescribed by a linear extrapolation. Please note that be-ing in principle not very deep (being a matter of defi-nitions), this issue creates significant confusion betweenpeople from different communities. See appendix C foraddition comments on this issue. B. Affine compressibility contribution
The inset of fig. 5 presents for our three model sys-tems the correlation function η F ≡ βV h δ ˆ P i of the totalpressure in the NPT-ensemble measured as a functionof temperature while keeping the mean pressure P con-stant. As shown in sect. II, the “affine dilatational elas-ticity” η A ≡ η F | yields the leading contribution to thecompression modulus K . We note first that η A increaseswith decreasing temperature for our two glass-formingliquids in d = 2 and d = 3 dimensions (due to the in-creasing repulsion of the LJ beads) levelling off belowthe glass-transition temperature T g of each model. Atcontrast, η A is seen to increase monontonously with tem-perature. This qualitatively different behavior needs tobe explained. We remind first that η A can be expressed -3 -2 -1 P id (T) η A + P δ k=0,P=0 δ k=0,P=0.1 δ k=0,P=0.5 δ k=0,P=0.9 δ k=0.5,P=0 δ k=0.9,P=0 δ k=0.999,P=0 -3 -2 -1 T η A
1D net2D pLJ3D KA + P i d D n e t s FIG. 5: Affine dilatational elasticity η A determined from thepressure fluctuations η F | in the NPT-ensemble. Inset: η A asa function of T for 1D nets ( δk = 0, P = 0), 2D pLJ beads( P = 2) and 3D KA binary mixtures ( P = 1). Main panel:Data collapse for η A + P as a function of P id ( T ) for 1D netsfor a broad range of parameters as indicated. The bold linecorresponds to the expected behavior for harmonic springs. by the simple average, eq. (29), which can be evaluatedin any ensemble. For the pair potentials we focus on herethis yields η A = P ex + η Born = P id + η A , ex with η Born beingthe Born-Lam´e coefficient, eq. (4), as we have explicitlychecked for all models [16]. Using eq. (A13) one sees thatfor harmonic springs in d = 1 dimensions η A = P id + 1 V X l k l (cid:10) x l (cid:11) (65)for a given quenched realization of spring constants k l .Since in the NPT-ensemble each spring is uncorrelated,we can use that h x l i − h x l i = k B T /k l for the thermalfluctuation of every spring l . This implies that η A =2 P id + P l k l h x l i /V . Since the average length of a springat constant pressure P is given by h x l i = R − P A/k l andsince the mean volume V is the sum of averaged lengths h x l i times the surface A , it follows using eq. (63) for thecompression modulus K that η A = 2 P id + K ref − P + K ref h k l i h /k l i − − h /k l i P A/R (66)with K ref = ρ ref R / h /k l i and ρ ref = 1 /AR characteriz-ing the zero-pressure reference. For systems of identicalsprings this reduces to η A = 2 P id + K ref − P . With K ref = 1 this is shown by the bold line in the main panelof fig. 5. For polydisperse systems at zero pressure wehave instead η A = 2 P id + ρ ref R h k l i as already shown bythe dashed line in fig. 2. Within the units used this cor-responds numerically also to η A = 1 + 2 P id . Hence, all1D nets without noise or without applied pressure shouldcollapse on the same master curve if η A + P is plotted asa function of the ideal pressure P id ( T ). This is confirmedby the data presented in the main panel of fig. 5.0 -3 -2 -1 P id (T) -3 -2 -1 η F | = η F , ex | δ k=0,P=0 δ k=0,P=0.1 δ k=0,P=0.5 δ k=0,P=0.9 δ k=0.1,P=0 δ k=0.25,P=0 δ k=0.5,P=0 δ k=0.9,P=0 δ k=0.999,P=0 -3 -2 -1 η F | / η F | -1 c r e g | d=2 pLJd=3 KA id eq.(21) D n e t s FIG. 6: Characterization of non-affine elastic response. Mainpanel: Excess pressure fluctuation η F , ex | as a function of P id ( T ) for 1D nets for different δk and P . Systems with no( δk = 0) or little polydispersity and all high-temperature sys-tems collapse on the dashed line indicating η F , ex | = 2 P id .Inset: Dimensionless regression coefficient c reg | for all threemodel systems as a function of η F | / η F | = η F , ex | /η A con-firming eq. (21) as indicated by the bold line. C. Non-affine compressibility contribution
The non-affine deviations from the affine Born con-tribution η A to the compression modulus K are givenaccording to the Rowlinson formula, eq. (3), by thefluctuation η F , ex | of the excess pressure in the NVT-ensemble. Focusing on one low temperature we have al-ready seen η F , ex | in fig. 2 as a function of δk . Plotted asa function of the ideal pressure P id ( T ), the main panelof fig. 6 presents η F , ex | for various δk and P as indi-cated. Since η F , ex | = η A − K , it follows from eq. (63)and eq. (66) characterizing, respectively, the compressionmodulus and the affine dilatational elasticity, that for 1Dnets of harmonic springs η F , ex | = 2 P id + 1 − / h /k l i − P h /k l i (67)where we have used that the mean spring constant h k l i ,the reference length R of the springs and the surface A are all arbitrarily set to unity. Assuming P = 0, eq. (67)is traced in fig. 2 as a function of δk (dash-dotted line).Since the second term vanishes for identical springs, thisrelation implies η F , ex | = 2 P id for all temperatures andpressures as shown in the main panel of fig. 6 by thedash-dotted power-law slope. This is confirmed by thedata presented for δk = 0. As one expects, this alsogives the high-temperature limit for all systems. Notethat η F , ex | only vanishes in the low- T limit of strictlyidentical springs. Since the second term in eq. (67) is afinite constant for δk >
0, the corresponding data mustthus level off for T →
0. Using a simple example wehave thus confirmed the more general finding by Lutsko -3 -2 -1 x = P id /K -20-15-10-50 η F,ex | - data δ k=0,P=0 δ k=0,P=0.1 δ k=0,P=0.5 δ k=0,P=0.9 δ k=0.1,P=0 δ k=0.25,P=0 δ k=0.5,P=0 δ k=0.9,P=0 δ k=0.999,P=0d=2 pLJ P=2d=3 KA P=1 -3 -2 -1 x = P id /K -20-15-10-50 y = η F , m i x | / P i d y = - x low-T limit y = - x high-T limit FIG. 7: Stress fluctuations in the NPT-ensemble. Largespheres refer to 2D pLJ beads for P = 2, squares to 3D KAbeads for P = 1, all other symbols to 1D nets for different δk and P as indicated. Main panel: y = η F , mix | /P id as afunction of the reduced ideal pressure x = P id /K confirmingeq. (35). Inset: Similar scaling for the reduced correlationfunction y = ( η A , ex − η F , ex | ) /P id using the same symbols asin the main panel. The bold line represents eq. (36). [5] that the non-affine contributions to the elastic modulineed not necessarily vanish in the low- T limit. Inter-estingly, due to the P -dependent denominator the low- T non-affine contribution becomes more and more relevantfor larger pressures. It becomes thus increasingly diffi-cult to reach the asymptotic high- T limit (not shown)[41]. As stressed in sect. II B, the correlation coefficient c reg of the instantaneous volumes and pressures in theNPT-ensemble is set by the reduced non-affine contribu-tion η F | / η F | = η F , ex | /η A . This relation is verified inthe inset of fig. 6 presenting for all three studied modelsystem a perfect collapse of the data on the master curve(bold solid line) predicted by eq. (21). Hence, the coef-ficient c reg may be seen as a dimensionless and properlynormalized “order parameter” characterizing the relativeimportance of the non-affine displacements. This alsoshows that it is equivalent to specify for a system either( K, η A ), ( η A , η F , ex | ) or ( K, c reg ). D. Correlations in the NPT-ensemble
As shown in fig. 7, we have checked the predictedcorrelations between the ideal and the excess pressurefluctuations η F , mix | and the auto-correlations of the ex-cess pressure η F , ex | . To make all investigated mod-els comparable the reduced correlation functions y = η F , mix | /P id (main panel) and y = ( η A , ex − η F , ex | ) /P id (inset) are traced as a function of the reduced ideal pres-sure x = P id /K with K as determined independentlyabove. Please note that compared to our theoretical pre-dictions we have scaled the vertical axes of the data notwith K , but using the ideal pressure P id , i.e. the predic-tions have been divided by x . This was done for presen-1tational reasons since thus the low- T asymptote becomesa constant. A perfect data collapse on the (rescaled) pre-diction y = 1 − x (bold line) is observed for all systems.As shown in ref. [29] a similar data collapse is achievedfor the Rowlinson formula if y = F Row | /P id is plotted asa function of x . (Since a related scaling plot is presentedin sect. IV E, we do not repeat this figure here.) Interest-ingly, the latter scaling does not depend on the MC-gaugewhich has been used above to simplify the derivation ofeq. (7). Please note that for our glass-forming liquidsit is not possible (for the imposed pressures) to increase x beyond unity ( K ≥ P id ) and the deviations from thepredicted plateau y = f ( x ) /x ≈ x ≪ E. Generalized λ -ensembles Having characterized the behavior in the NPT-ensemble ( λ = 0) and the NVT-ensemble ( λ = 1) weturn now to the general λ -ensembles. As described at theend of sect. III, these Gaussian ensembles [31, 32] maybe realized by attaching an external harmonic spring,eq. (8), parallel to the system. As we have already seenin fig. 1, this allows to gradually reduce the volume fluc-tuations while keeping constant the average pressure P of the system, i.e. the thermodynamic (average) state ofthe system remains unchanged [34]. Expressed as a func-tion of λ , the mean-squared volume fluctuations decreaselinearly in agreement with eq. (48). Confirming eq. (50),a similar scaling is also found for other moments of thevolume (not shown).Focusing on simple 1D nets the verification of severalpredicted pressure and volume-pressure correlation func-tions is presented in fig. 8. The data presented in panel(a) has been obtained for systems with identical springconstants ( δk = 0) at a low temperature T = 0 .
01. Vol-ume and pressure fluctuations are thus highly correlatedin the NPT-ensemble, i.e. c reg | ≈
1, and the fluctua-tions of the excess pressure η F , ex | in the NVT-ensembleis consequently small. The second system presentedin panel (b) corresponds to a polydispersity δk = 0 . T = 1 for which c reg | ≈ .
44 and η F , ex | /η A ≈ .
8. Note that the vertical axes are madedimensionless by rescaling the data using either the com-pression modulus K or the affine dilatational elasticity η A obtained for λ = 0. We stress first that the com-pression modulus K may be computed irrespective of λ by linear regression (spheres) of the measured ( ˆ V , ˆ P ) aspredicted by eq. (53). The corresponding dimensionlesscorrelation coefficient c reg | λ indicated by the squares iffound to decrease monotonously from its maximum at λ = 0 to zero in the NVT-ensemble. In agreement witheq. (54) the decay is more sudden for values K/η A closeto unity, i.e. if the non-affine contribution is negligibleas in panel (a), and becomes more and more gradual forsmaller K/η A as seen in panel (b). The linear decay ofthe total pressure correlation function η F | λ demonstrates λ F reg | λ /Kc reg | λ η F | λ / η A F Row | λ /Kd=1, δ k=0,T=0.01,P=0: f =0.0199, η A =1.02 (a) e q . ( ) e q . ( ) eq.(53) e q . ( ) ^ λ F reg | λ /Kc reg | λ η F | λ / η A F Row | λ /K η F,ex | λ / η A δ k=0.9,T=0.01,P=0: K=0.6, f =0.5, η A =3.0 (b) e q . ( ) e q . ( ) eq.(53) e q . ( ) e q . ( ) ^ FIG. 8: Verification of various predictions for generalized λ -ensembles for simple 1D nets: (a) identical springs ( δk = 0)at temperature T = 0 .
01 and pressure P = 0 where K = 1, η A = 1 . c reg | ≈ f = 0 . x = P id /K =0 .
01 and (b) polydisperse springs with δk = 0 . T = 1 . P = 0 where K = 0 . η A = 3 . c reg | = 0 .
44 and f = 0 .
53 for x = P id /K = 1 .
68. The lines compare the datawith the respective prediction: the thin solid lines show thecompression modulus obtained using the regression relation,eq. (53), the dash-dotted lines the transformation relation forthe pressure fluctuations η F | λ , eq. (55), the bold solid linesthe key claim, eq. (9), for the Rowlinson functional F Row | λ . that the predicted transformation relation of the totalpressure fluctuations, eq. (55), indeed holds as shown bythe dash-dotted lines. Note that in the limit λ → η F | λ → η F , ex | due to the MC-gauge used forthe instantaneous ideal pressure. Since for the first sys-tem in panel (a) the thermal noise is very week, we have η F | λ ≈ η F , ex | λ for all λ (not shown). Therefore, η F , ex | λ is only presented in the second panel. Our key predictioneq. (9) for the Rowlinson functional is confirmed by thecrosses presenting F Row | λ /K vs. λ .The scaling of F Row | λ is further investigated in themain panel of fig. 9 for a broad range of 1D nets (withparameters as indicated in the figure) and 2D pLJ beads(spheres) for one state point ( T = 0 . P = 2). Ac-cording to eq. (7,) F Row | λ /K − f ( x ) → -3 -2 -1 λ (1-x) -3 -2 -1 F R o w | λ / K - f ( x ) δ k=0,P=0,T=0.01 δ k=0,P=0.9,T=0.001 δ k=0.25,P=0,T=0.1 δ k=0.5,P=0,T=3 δ k=0.9,P=0,T=1 δ k=0.999,P=0,T=1d=2 pLJ P=2,T=0.001 -2 -1 λ x(x-1) -2 -1 η F , m i x / K + x ( x - ) x=P id /K > 1 e q . ( ) e q . ( ) NPT-limit NVT-limitNoise !
FIG. 9: Scaling with λ for various compression moduli K andreduced ideal pressures x = P id /K . Main panel: RescaledRowlinson functional F Row | λ for different polydispersity δk ,pressure P and temperature T as indicated. The spheres in-dicate 2D pLJ beads, all other symbols data obtained for 1Dnets. Note that the horizontal scaling variable λ (1 − f ( x )) = λ (1 − x ) ≥
0. The bold line corresponds to eq. (9). Inset:Rescaled data for η F , mix | λ as a function of λx ( x − x > as scaling variable λ (1 − f ( x )) = λ (1 − x ) for the hor-izontal axis, all data must collapse on the bisection lineaccording to eq. (7). As may be seen from the double-logarithmic plot, this is indeed the case for a broad rangeof data sets. The inset of fig. 9 presents a similar scalingplot for the function η F , mix | λ characterizing the correla-tions of ideal pressure and excess pressure contributionsfor 1D nets. The vertical axis corresponds again to theprediction of the NPT-limit, eq. (35). The horizontalaxis λx ( x −
1) is chosen such that according to eq. (58)all rescaled η F , mix | λ have to fall on the bisection line.(Since the horizontal axis is logarithmic, only data with x > V. CONCLUSION
We have revisited in this paper theoretically and nu-merically various correlations of the normal pressure ˆ P and its contributions ˆ P id and ˆ P ex for isotropic solidsand fluids using simple coarse-grained models [27], e.g.,strictly 1D networks of permanently fixed springs or theKA model for binary mixtures in three dimensions [38].Making more precise several statements made in the ap-pendix of ref. [13] and extending the brief communicationref. [29], we have compared fluctuations in generalized λ -ensembles where the volume fluctuations are tuned bymeans of an external harmonic spring potential allowingto switch gradually between the standard NPT- ( λ = 0)and NVT-ensembles ( λ = 1). We have stressed that thewidely used stress fluctuation formula, eq. (3), for the compression modulus K in the NVT-ensemble may beobtained directly • without the affine volume rescaling trick for theNVT-ensemble (reminded in appendix B) used firstfor liquids by Rowlinson [3] and, • more importantly, without assuming a referenceposition for the particles and a microscopic dis-placement field which is only possible for solids [4]using the general thermodynamic transformation rulesbetween conjugated ensembles [26] and assuming the sys-tems to be sufficiently large ( V → ∞ ) for the given tem-perature T and compression modulus K [28]. The directthermodynamic derivation can readily be adapted to theshear modulus G in isotropic systems [13] and to the moregeneral elastic moduli characterizing anisotropic solids asreminded in appendix C [28].The Rowlinson stress fluctuation functional F Row [11]has been computed deliberately in the unusual NPT-ensemble ( λ = 0) to make manifest the general transformeq. (18) at the heart of the stress fluctuation formalism.We have demonstrated that F Row | = P id (2 − P id /K ), i.e. F Row | vanishes in the low-temperature limit. More gen-erally, we have investigated F Row | λ and other correlationfunctions as a function of the parameter λ characterizingthe volume fluctuations. As announced in the Introduc-tion, eq. (9), the Rowlinson functional is found to interpo-late linearly between the classical ensemble limits (fig. 8).Note that the specification of the compression modulus K and the dimensionless regression coefficient c reg | im-plies η A and η F , ex | (and vice versa ) and together withthe ideal pressure P id (which implies x = P id /K ) thisallows the complete description of all the discussed cor-relation functions at different λ . Our theoretical and nu-merically results, especially eq. (9), may allow to read-ily calibrate (correctness, convergence and precision) thevarious barostats commonly used [11, 39]. In the nearfuture we plan (i) to also consider systems close to afirst-order (e.g., solid to liquid) phase transition gener-alizing the work by Hetherington [31] and (ii) to extendour approach to negative values of λ (following in thatref. [7]) increasing artificially the — then not necessar-ily Gaussian — fluctuations of the extensive variable andmaking thus the system increasingly unstable. Acknowledgments
H.X. thanks the CNRS and the IRTG Soft Matter forsupporting her sabbatical stay in Strasbourg, P.P., C.G.and J.H. the IRTG Soft Matter and F.W. the DAAD forfunding. We are indebted to A. Blumen (Freiburg) forhelpful discussions.3
Appendix A: Affine interaction energy
Affine displacement assumption.
We consider herethe change of the conservative interaction energy U s ( ˆ V )of a configuration s under an imposed dilatational strain assuming that all particles respond to the macroscopicconstraint by an affine microscopic displacement. It isconvenient to introduce a dimensionless scalar character-izing the relative volume change ǫ ≡ δ ˆ V ( ǫ ) / ˆ V (0) ≡ ˆ V ( ǫ ) / ˆ V (0) − V (0) being the instantaneous volume of the unper-turbed reference simulation box ( ǫ = 0). Obviously, forcanonical ensembles the instantaneous volume ˆ V can bereplaced by V . The affinity assumption implies, e.g., x (0) ⇒ x ( ǫ ) = x (0)(1 + ǫ ) /d (A2)for the x -coordinate of particle positions or relative par-ticle distances where the argument (0) denotes again thereference system. Equation (A2) implies that the squareddistance r between two particles transforms as r (0) ⇒ r ( ǫ ) = r (0)(1 + ǫ ) /d . (A3)It follows that dr ( ǫ ) dǫ → d r (0) , (A4) d r ( ǫ ) dǫ → − d ) d r (0) (A5)where we have taken finally for both derivatives the limit ǫ →
0. For the first two derivatives of a general function f ( r ( ǫ )) with respect to ǫ this implies ∂f ( r ( ǫ )) ∂ǫ → d rf ′ ( r ) (A6) ∂ f ( r ( ǫ )) ∂ǫ → d (cid:0) r f ′′ ( r ) + rf ′ ( r ) (cid:1) − d rf ′ ( r ) (A7)for ǫ → Pair potential assumption.
Let us assume now in ad-dition that the interaction energy is given by U s ( ǫ ) = X l u ( r l ( ǫ )) (A8)where the index l labels the interactions between the par-ticles i and j with i < j . We note en passant that for thefirst derivative of the interaction energy U s with respectto the volume eq. (A6) implies − U ′ s ( ˆ V ) = − U ′ s ( ǫ ) /V (0) = − d ˆ V X l r l u ′ ( r l ) (A9) where we have taken ǫ → P ex . As definedin eq. (31), the coefficient η A , ex /V measures the secondderivative of the interaction energy U s ( ˆ V ) of a configura-tion s with respect to an affine dilatational strain. Usingeq.(A7) it follows for the instananeous value ˆ η A , ex thatˆ η A , ex /V ≡ U ′′ s ( ˆ V ) = 1ˆ V (0) U ′′ s ( ǫ )= 1 d ˆ V X l r l u ′′ ( r l ) + r l u ′ ( r l ) − d ˆ V X l r l u ′ ( r l ) (A10)where we have taken again ǫ → Simple averages.
Only the affinity assumption andthe pair potential choice have been used up to now. Themean value η A , ex = h ˆ η A , ex i is then obtained by takingthe thermal average over all configurations s and over allvolumes ˆ V depending on the ensemble. Since all con-tributions to eq. (A10) correspond to simple averages ,thermostatistics tells us that one can replace ˆ V for suffi-ciently large systems by its mean value V , eq. (14). Oneconfirms that the affine dilatational elasticity becomes η A , ex = η Born + P ex with (A11) η Born = 1 d V *X l r l u ′′ ( r l ) + r l u ′ ( r l ) + (A12)as already stated in eq. (4). Note that it is inconsistent toneglect the explicit excess pressure contribution to η A , ex in eq. (A11) but to keep the underlined contribution to η Born which amounts to − P ex /d . The sum of both terms P ex ( d − /d only vanishes in d = 1 dimensions where η A , ex = 1 V *X l x l u ′′ ( x l ) + (A13)as used in sect. IV B. Appendix B: Volume rescaling trick revisited
Introduction.
The stress fluctuation formula, eq. (3),for the compression modulus K of isotropic systems hasbeen derived in sec. II C essentially using the generaltransformation relation eq. (10) between conjugated en-sembles and evaluating the pressure fluctuations in theNPT-ensemble by integration by parts, eq. (29). Histori-cally, eq. (3) has been first derived properly by Rowlinson[3] who computed using the volume rescaling trick thesecond derivative, eq. (1), of the free energy F ( T, V ) = − k B T log[ Z ( V )]. (The NVT-ensemble is thus used and,hence, V = ˆ V in the following.) Since the total system4Hamiltonian may be written as the sum of a kinetic en-ergy and a potential energy, the partition function Z ( V )factorizes in an ideal contribution Z id ( V ) and an excesscontribution Z ex ( V ) on which we focus below. We re-mind [2] that using Z id ( V ) ∼ V N one readily confirms P id = k B T ρ and K id = P id for the ideal contributionsto, respectively, the total pressure P = P id + P ex and thetotal compression modulus K = K id + K ex . Mapping of strained and unstrained configurations.
Since both the perturbed as the unpertubed partitionfunction is sum over all possible particle configurations,one can always compare the contribution of a configu-ration s of the strained system with the contribution ofa configuration of the reference being obtained by theaffine rescaling of all coordinates according to eq. (A2).Using the same notations as in appendix A the interac-tion energy U s ( ǫ ) of the strained system can be expressedin terms of the coordinates (state) of the unperturbedsystem and the explicit metric parameter ǫ . The ex-cess contribution F ex ( ǫ ) = − k B T Z ex ( ǫ ) is thus obtainedfrom Z ex ( ǫ ) = P s exp [ − βU s ( ǫ )] where constant prefac-tors have been omitted and where the sum is taken overall possible configurations s . General conservative potential.
We note for the firsttwo derivatives of the free energy[ln( Z ex ( ǫ ))] ′ = Z ′ ex ( ǫ ) Z ex ( ǫ ) (B1)[ln( Z ex ( ǫ ))] ′′ = Z ′′ ex ( ǫ ) Z ex ( ǫ ) − (cid:18) Z ′ ex ( ǫ ) Z ex ( ǫ ) (cid:19) (B2)and of the partition function Z ′ ex ( ǫ ) = − X s βU ′ s ( ǫ ) e − βU s ( ǫ ) (B3) Z ′′ ex ( ǫ ) = X s ( βU ′ s ( ǫ )) e − βU s ( ǫ ) − X s ( βU ′′ s ( ǫ )) e − βU s ( ǫ ) (B4)with a prime denoting again a derivative with respect tothe indicated argument. Using P ex = − ∂F ex ( V ) /∂V andeq. (B3) and taking the limit ǫ → P ex = D ˆ P ex E with ˆ P ex ≡ − V U ′ s ( ǫ ) | ǫ =0 (B5)which defines the instantaneous excess pressure. (The av-erage taken uses the weights of the unperturbed system.)The excess pressure thus measures the average change ofthe total interaction energy U s ( ǫ ) taken at ǫ = 0. Theexcess compression modulus K ex is obtained using in ad-dition eq. (B2) and eq. (B4) and taking finally the ǫ → βV K ex = h βU ′′ s ( ǫ ) i| ǫ =0 − (cid:16)(cid:10) ( βU ′ s ( ǫ )) (cid:11) − h βU ′ s ( ǫ ) i (cid:17)(cid:12)(cid:12)(cid:12) ǫ =0 . (B6) Being a simple average the first term on the right hand-side corresponds to η A , ex as defined in eq. (31). Usingeq. (B5) the second term is seen to yield the pressurefluctuation contribution η F , ex | . Summarizing all termswe have thus shown that K = K id + K ex = P id + η A , ex − η F , ex | (B7)in agreement with eq. (3) if the Born-Lam´e coefficient isdefined more generally as η Born ≡ η A , ex − P ex . Pair potential choice.
Up to now we have stated P ex and K ex for a general interaction potentiel U s ( ǫ ). Assum-ing the interactions to be described by a pairwise additivepotential, eq. (A8), the Kirkwood relation, eq. (6), is con-firmed using eq. (B5) and eq. (A9). As already noted atthe end of appendix A, one confirms using eq. (A10) that η Born = η A , ex − P ex agrees with eq. (4). Appendix C: General stress fluctuation formalism
Introduction.
We remind here the general stress fluc-tuation formalism derived by Squire, Hold and Hoover [4]and show that Rowlinson’s formula, eq. (3), is a specialcase obtained by symmetry considerations. For conve-nience we introduce the two linear projection operators T [ A αβ ] ≡ d d X α,β =1 A αβ δ αβ (C1) T [ A αβγδ ] ≡ d d X α,β,γ,δ =1 A αβγδ δ αβ δ γδ (C2)with A αβ and A αβγδ being, respectively, second- andforth-rang tensors and δ αβ the Kronecker symbol [36].Greek letters are used for the spatial coordinates α, β, γ, δ = 1 , . . . , d . The following identities are read-ily verified T [ A αβ B γδ ] = T [ A αβ ] T [ B αβ ] (C3) T [ A αβγδ ] = 1 for A αβγδ = δ αβ δ γδ (C4) T [ A αβγδ ] = 1 /d for A αβγδ = δ αγ δ βδ (C5) T [ A αβγδ ] = 1 /d for A αβγδ = n α n β n γ n δ (C6)with n α , . . . being the spatial components of a normal-ized vector, i.e. P α ( n α ) = 1. Thermodynamics and symmetry.
Generalizing thedefinition of the pressure P given in eq. (1), the stresstensor σ αβ may be defined as the first derivative of thefree energy per volume with respect to the linear strain ǫ αβ [1] characterizing the macroscopic deformation of thesystem. Note that in general σ αβ does not vanish at theinvestigated state point, i.e. the systems may be pre-stressed at the reference strain ǫ αβ = 0. The latter pointis crucial for essentially all soft matter systems whichonly assemble because a finite density and/or stress isapplied. It is less important for the classical crystalline5solids where the elastic moduli are normally huge com-pared to the imposed stresses. We denote by δσ αβ theincrement of the stress tensor σ αβ under an increment δǫ αβ generalizing the dilatational strain ǫ used in ap-pendix B. Assuming an infinitessimal strain increment,Hooke’s law reads quite generally [1] δσ αβ = d X γ,δ =1 E αβγδ δǫ γδ (C7)where the elastic moduli E αβγδ stand for the secondderivative ∂ F/∂ǫ αβ ∂ǫ γδ of the free energy per volume atthe given thermodynamic state [1]. Let us now assume apure dilatational strain without shear, i.e. δǫ αβ = ǫ δ αβ .As may be seen from eq. (4.6) of ref. [1], this implies δσ αβ = dKǫδ αβ . Hence, P α δσ αα = d Kǫ . Or, usingeq. (C7) one sees that P α δσ αα = P αβ E ααββ ǫ . Com-paring both expressions, this shows that the compressionmodulus is given by K = T [ E αβγδ ] using the linear pro-jection operator defined above. Stress fluctuation relations.
As described in the liter-ature [4, 5, 12, 20, 24], the elasticity tensor E αβγδ cannumerically be computed from the sum E αβγδ = B αβγδ + C αβγδ K + C αβγδ Born + C αβγδ F , id + C αβγδ F , ex (C8)with contributions as specified below. The compressionmodulus is obtained by applying the linear operator T to each term and by summing up the contributions. Notethat some authors only refer to the sum C αβγδ K + C αβγδ Born + C αβγδ F , id + C αβγδ F , ex as the elasticity tensor [4, 5, 12]. Initial stress contribution.
The (often not included)first term B αβγδ in eq. (C8) may we written as [12] B αβγδ ≡ − σ αβ δ γδ + σ αγ δ βδ + σ αδ δ βγ (C9)= P ( δ αβ δ γδ − δ αγ δ βδ − δ αδ δ βγ ) . (C10)Following Birch [42] we have assumed in the second stepthat the system is isotropically stressed, i.e. σ αβ = − P δ αβ . Since B αβγδ becomes negligible for small to-tal pressures P ≈
0, this term is often not computed.Consistency implies then, however, to set P id = − P ex inthe remaining terms of eq. (C8). Returning to eq. (C10)one sees using eq. (C4) and eq. (C5) that T [ B αβγδ ] = P (cid:18) − d (cid:19) = P − P id d − P ex d . (C11)There is thus no contribution to K from the Birch termfor d = 2 dimensions and a contribution − P for d =1. The leading term P in the second step of eq. (C11)corresponds to the total pressure indicated in eq. (3). Ideal pressure contributions.
The second and theforth term in eq. (C8) correspond to the kinetic (ideal)contributions C αβγδ K ≡ P id ( δ αγ δ βδ + δ αδ δ βγ ) , (C12) C αβγδ F , id ≡ − P id ( δ αγ δ βδ + δ αδ δ βγ ) . (C13) Using again eq. (C4) and eq. (C5) this yields T [ C αβγδ K ] =4 P id /d and T [ C αβγδ F , id ] = − P id /d . Together with theideal pressure contribution from eq. (C11) all the kineticterms sum up (independently of the dimension) to P id inagreement with eq. (3). It is worthwhile to stress that ifthe Birch term B αβγδ is neglected for P id = − P ex , theideal pressure contributions sum up incorrectly as is thecase in ref. [4]. In fact, the traditional separation of the(rather trivial) different kinetic terms stated in eq. (C10)is from the numerical point of view unfortunate. Moreimportantly, this separation is theoretically misleading:The total system Hamiltonian being the sum of the ki-netic energy and the conservative interaction potential,all kinetic contributions to the elastic moduli can (andshould) be factorized out from the start. This leads im-mediately to the obvious total kinetic contribution C αβγδ id = P id δ αβ δ γβ (C14)replacing the aforementioned three terms. Affine Born contribution.
Generalizing the Born co-efficient η Born in eq. (3), the third term C αβγδ Born in eq. (C8)corresponds to the (second derivative) of the energychange assuming the affine displacement of all particles(as in appendix A for a pure dilatational strain). Assum-ing pair interactions it is given by [43] C αβγδ Born ≡ V X l D(cid:0) r l u ′′ ( r l ) − r l u ′ ( r l ) (cid:1) n αl n βl n γl n δl E (C15)where n αl , . . . stand for the normalized vector between thepair of interacting particles labeled by l . Using eq. (C6)one confirms that T [ C αβγδ Born ] = 1 d V *X l r l u ′′ ( r l ) − r l u ′ ( r l ) + = η Born + 2 P ex d (C16)where we have used the Kirkwood formula, eq. (6), in thelast step. Note that the underlined term cancels exactlythe underlined last term in eq. (C11). This is equivalentto the statement made after eq. (A12) in appendix A. Stress fluctuations.
The last two terms indicated ineq. (C8) stand for the correlation function C αβγδ F ≡ C αβγδ F , id + C αβγδ F , ex ≡ − βV (cid:10) δ ˆ σ αβ δ ˆ σ γδ (cid:11) (C17)with δ ˆ σ αβ denoting the fluctuation of the total stress ten-sor σ αβ = σ αβ id + σ αβ ex . Using the factorization of ideal andexcess stress fluctuations in the canonical ensemble atconstant volume V (as shown in sect. II D this is differentin NPT-ensembles), both contributions can be separated.The ideal stress fluctuation contribution has already beenaccounted for above. We remind that the instantaneousexcess pressure ˆ P ex is given by the trace ˆ P ex = T [ ˆ P αβ ex ]over the excess pressure tensor ˆ P αβ ex being the negative6of the instantaneous excess stress tensor. Using eq. (C3)one thus obtains T [ C αβγδ F , ex ] = − βV D T [ δ ˆ P αβ ex ] E (C18)which is identical to the pressure fluctuation term ineq. (3). We have thus confirmed that the general stress fluctuation relation, eq. (C8), is consistent with the Rowl-inson formula for the compression modulus K . [1] L. D. Landau and E. M. Lifshitz, Theory of Elasticity (Pergamon Press, 1959).[2] H. B. Callen,
Thermodynamics and an Introduction toThermostatistics (Wiley, New York, 1985).[3] J. S. Rowlinson,
Liquids and liquid mixtures (Butter-worths Scientific Publications, London, 1959).[4] D. R. Squire, A. C. Holt, and W. G. Hoover, Physica ,388 (1969).[5] J. F. Lutsko, J. Appl. Phys , 2991 (1989).[6] J. P. Wittmer, A. Tanguy, J.-L. Barrat, and L. Lewis,Europhys. Lett. , 423 (2002).[7] K. van Workum and J. de Pablo, Phys. Rev. E , 011505(2003).[8] J.-L. Barrat, J.-N. Roux, J.-P. Hansen, and M. L. Klein,Europhys. Lett. , 707 (1988).[9] G. Papakonstantopoulos, R. Riggleman, J. L. Barrat,and J. J. de Pablo, Phys. Rev. E , 041502 (2008).[10] N. Schulmann, H. Xu, H. Meyer, P. Poli´nska,J. Baschnagel, and J. P. Wittmer, Eur. Phys. J. E , 93(2012).[11] M. Allen and D. Tildesley, Computer Simulation of Liq-uids (Oxford University Press, Oxford, 1994).[12] D. Frenkel and B. Smit,
Understanding Molecular Sim-ulation – From Algorithms to Applications (AcademicPress, San Diego, 2002), 2nd edition.[13] J. P. Wittmer, H. Xu, P. Poli´nska, F. Weysser, andJ. Baschnagel, J. Chem. Phys. , 12A533 (2013).[14] M. Born and K. Huang,
Dynamical Theory of CrystalLattices (Clarendon Press, Oxford, 1954).[15] H. Xu, J. Wittmer, P. Poli´nska, and J. Baschnagel, Phys.Rev. E , 046705 (2012).[16] If a truncated potential is used as for the two glass-forming liquids discussed, some care is needed for thecomputation of η Born . Being a moment of the second po-tential derivative, η Born needs to be corrected using aweighted histogram evaluated at the cutoff as describedin ref. [15]. This correction is a simple average and thesame value is obtained for any λ . It becomes relevant ifeq. (9) is probed for small λ .[17] A. Tanguy, J. P. Wittmer, F. Leonforte, and J.-L. Barrat,Phys. Rev. B , 174205 (2002).[18] A. Tanguy, F. Leonforte, J. P. Wittmer, and J.-L. Barrat,Applied Surface Sience , 282 (2004).[19] F. L´eonforte, R. Boissi`ere, A. Tanguy, J. P. Wittmer, andJ.-L. Barrat, Phys. Rev. B , 224206 (2005).[20] J.-L. Barrat, in Computer Simulations in CondensedMatter Systems: From Materials to Chemical Biology ,edited by M. Ferrario, G. Ciccotti, and K. Binder(Springer, Berlin and Heidelberg, 2006), vol. 704, pp.287—307.[21] A. Zaccone and E. Terentjev, Phys. Rev. Lett. ,178002 (2013). [22] C. Maloney and A. Lemaˆıtre, Phys. Rev. Lett. , 195501(2004).[23] K. Yoshimoto, T. Jain, K. van Workum, P. Nealey, andJ. de Pablo, Phys. Rev. Lett. , 175501 (2004).[24] B. Schnell, H. Meyer, C. Fond, J. Wittmer, andJ. Baschnagel, Eur. Phys. J. E , 97 (2011).[25] F. L´eonforte, A. Tanguy, J. P. Wittmer, and J.-L. Barrat,Phys. Rev. Lett. , 055501 (2006).[26] J. L. Lebowitz, J. K. Percus, and L. Verlet, Phys. Rev. , 250 (1967).[27] We only consider classical systems. The generalization toquantum systems may involve delicate problems.[28] The general stress fluctuation formalism may be appliedclose to a first order phase transition only as long asthe transformation relation between conjugated ensem-bles [26] remains valid. Necessary conditions are that theelastic modulus of interest remains positive definite andthat the fluctuations of the extensive variable X are sym-metric around the main maximum of the distribution p ( ˆ X ). The stress fluctuation formalism becomes incor-rect in general close to a second order phase transition.[29] J. P. Wittmer, H. Xu, P. Poli´nska, F. Weysser, andJ. Baschnagel, J. Chem. Phys. , 191101 (2013).[30] The spring constant k ext has the dimension energy pervolume which implies, as one expects, [ K ext ] = en-ergy/volume for the associated compression modulus.[31] J. Hetherington, J. Low Temp. Phys. , 145 (1987).[32] M. Costeniuc, R. Ellis, H. Touchette, and B. Turkington,Phys. Rev. E , 026105 (2006).[33] A similar external spring potential has been introducedin ref. [7] using a negative spring constant in order toreduce the effective modulus of the total system.[34] The focus of Hetheringtons’s Gaussian ensemble [31], asof related generalizations [32], 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. More im-portantly, Hetherington’s additional weight factor doesin general correspond to a change of the mean inten-sive variable. This is why we have used the Gaussian,eq. (8), centered at the mean volume V = V ext and notjust U ext ( ˆ V ) ∝ ˆ V which would alter the pressure P .[35] With X being the extensive variable, the intensivevariable may be either defined as the derivative I = ∂U ( X ) /∂X of the inner energy U or as the derivative J = ∂S ( X ) /∂X of the entropy S [2]. It is the seconddefinition which is used in ref. [26]. Note that J = − βI for all extensive variables X other than U [2]. In our casewe have X = V , I = − P and J = βP . [36] M. Abramowitz and I. A. Stegun, Handbook of Mathe-matical Functions (Dover, New York, 1964).[37] Ergodicity problems are irrelevant since these systemshave been sampled by means of MC simulation.[38] W. Kob and H. C. Andersen, Phys. Rev. E , 4134(1995).[39] S. J. Plimpton, J. Comp. Phys. , 1 (1995).[40] Qualitatively, this is similar to the predicted [21] andnumerically observed [8, 13] cusp-like singularity G ≈ p − T /T g of the shear modulus G in colloidal and poly-mer glasses at the glass transition temperature T g due to the increase of the non-affine displacements.[41] It is possible to collapse the data by plotting the ratio of η F , ex | and the second term in eq. (67) as function of theratio of P id and a crossover ideal pressure P ∗ id .[42] F. Birch, J. App. Phys. , 279 (1938).[43] Equation (C15) assumes implicitly that the excess stressis computed according the Kirkwood stress expressiongeneralizing eq. (6). Consistency requires that the excessstress fluctuation contribution C αβγδ F , ex is computed usingthe same definition for the instantaneous excess stress. a)(b) i=0i=0 i=1 w a ll w a ll f ext i l=i i−1polydisperse springs constant forceexternal w a ll A = A = l=1 stiff (c) w a ll A = particle position i X X=V springexternal U i=N l=Nl=N