Communication: Pressure fluctuations in isotropic solids and fluids

Among the fundamental properties of any equilibrium system are its elastic moduli characterizing the fluctuations of its extensive and/or conjugated intensive variables.^1–7 The isothermal compression modulus K of an isotropic solid or fluid may thus be obtained in the NPT-ensemble at imposed particle number N, pressure P, and temperature T from the fluctuations

$δV̂=V̂−V$ of the instantaneous volume

$\hat{V}$

$V̂$ around its mean value

$V = \langle \hat{V}\rangle$

$V=⟨V̂⟩$ according to the strain fluctuation relation² 

with k_B being Boltzmann's constant. Equivalently, K may be obtained in a canonical NVT-ensemble using Rowlinson's stress fluctuation relation^4,8,9

with β = 1/k_BT being the inverse temperature,

$P̂ ex$ the instantaneous excess pressure contribution, and η_B a Born-Lamé coefficient⁹ which for pairwise additive potentials becomes a simple sum of moments of derivatives of the potential with respect to the particle distance.⁸ In this Communication, we emphasize that the stress fluctuation representation of the elastic moduli^4–6 may be obtained directly from the well-known transformation rules between conjugated ensembles.¹⁰ Focusing on the compression modulus K this is made manifest by computing Rowlinson's expression K_row deliberately for NPT-ensembles where the volume is allowed to freely fluctuate. We show that

with P_id being the ideal pressure contribution. We demonstrate first Eq. (3) by considering theoretically the fluctuations of the instantaneous normal pressure

$P̂$ and its ideal and excess contributions

$\hat{P}_\mathrm{id}$

$P̂ id$ and

$\hat{P}_\mathrm{ex}$

$P̂ ex$ in both conjugated ensembles. These correlations are then checked numerically by means of Monte Carlo (MC) simulation of simple coarse-grained model systems.

As discussed in the literature,^2,8,10 a simple average

$A=⟨Â⟩$ of an observable

${\cal A}$

$A$ does not depend on the chosen ensemble, at least not if the system is large enough (V → ∞). A correlation function

$\langle \delta \hat{A}\delta \hat{B}\rangle$

$⟨δÂδB̂⟩$ of two observables

${\cal A}$

$A$ and

${\cal B}$

$B$ may differ, however, depending on whether V or P are imposed. As shown by Lebowitz, Percus, and Verlet,¹⁰ one verifies that

where K = −V∂P/∂V² has been used. For

$Â=B̂=P̂$⁠, this implies the transformation

i.e., the compression modulus K may be obtained from the difference of the pressure fluctuations in both ensembles. Interestingly, the numerically more convenient Rowlinson expression K_row for NVT-ensembles can be derived from Eq. (5)⁹ without using a microscopic displacement field (only possible for solids)⁵ and avoiding the volume rescaling trick used originally for liquids.⁴ 

There is a considerable freedom for defining the instantaneous pressure

$P̂=P̂ id +P̂ ex$ as long as its average P = P_id + P_ex does not change.⁸ It is convenient for the subsequent derivations and the presented MC simulations to define the instantaneous ideal pressure

$\hat{P}_\mathrm{id}$

$P̂ id$ by

and the instantaneous excess pressure

$P̂ ex$ by the Kirkwood expression.^8,9 Within this “MC-gauge” the thermal momentum fluctuations are assumed to be integrated out and the (effective) Hamiltonian

${\cal H}_s$

$Hs$ of a state s of the system may be written

with

$Us$ being the total excess potential energy.

In the following, the concise notation

$η≡βV⟨δP̂2⟩$ is used. An immediate consequence of the MC-gauge is, of course, that the fluctuations of

$\hat{P}_\mathrm{id}$

$P̂ id$ vanish for the NVT-ensemble and that, hence,

Since K > 0 for a stable system, Eq. (5) implies η|_P > η|_V. Depending on the disorder, η|_V is, however, not a negligible contribution as assumed (implicitly) by Born.¹ For solids it measures the effect of non-affinedisplacements under an imposed macroscopic strain.^6,7,9

The second moment of any intensive variable computed in an ensemble, where its mean value is imposed, is obtained readily by integration by parts. Using Eq. (7), this shows that

where a prime denotes a derivative with respect to the indicated variable. Albeit the indicated averages are taken over all states s and all volumes

$V̂$ at imposed P, being simple averagesthey can also be evaluated for sufficiently large systems in the NVT-ensemble yielding identical results. Denoting the last term in Eq. (9) by η_{A, ex}, one can show that it is equivalent for pair interaction potentials to the already mentioned Born-Lamé coefficient:⁹ η_{A, ex} = η_B + P_ex. Substituting Eqs. (8) and (9) into the Legendre transform Eq. (5), this confirms Eq. (2).

We focus now on stress fluctuations in the NPT-ensemble. By comparing with Eq. (5), one sees that if the Rowlinson formula K_row is applied at imposed P, this must yield

Interestingly, Eq. (10) does not completely vanish for finite T as does the corresponding stress fluctuation expression for the shear modulus G at imposed shear stress τ.⁹ As a next step, we demonstrate the relations

from which Eq. (3) is then directly obtained by substitution into Eq. (10). Returning to the general transformation relation Eq. (4), we first note that the l.h.s. must vanish if at least one of the observables is a function of

$V̂$⁠. With

$\hat{A}= \hat{B}= 1/\hat{V}$

$Â=B̂=1/V̂$⁠, it follows that

making the steepest-descent approximation

$⟨1/V̂⟩≈1/V$ for simple averages and using finally V/K = −∂V/∂P.¹¹ Remembering Eq. (6), this implies Eq. (11). With

$\hat{A}= \hat{P}_\mathrm{id}\break= k_{\rm B}TN/\hat{V}$

$Â=P̂ id =kBTN/V̂$ and

$\hat{B}= \hat{P}$

$B̂=P̂$⁠, one obtains similarly

to leading order for V → ∞. This relation implies finally the claimed correlation between ideal and excess pressure fluctuations, Eq. (12), using

$P̂=P̂ id +P̂ ex$ and the already demonstrated Eq. (11).¹² Please note that in Ref. 9, ideal and excess pressure fluctuations have incorrectly been assumed to be uncorrelated.

The numerical results reported here to check our predictions have been obtained by MC simulation of (i) one-dimensional (1D) nets with permanent cross-links and (ii) two-dimensional (2D) glass-forming liquids. Periodic boundary conditions are used and the pressure P is first imposed using a standard MC barostat.^8,9 After equilibrating and sampling in the NPT-ensemble, the volume is fixed,

$V=V̂$⁠, and various simple means and fluctuations are obtained in the NVT-ensemble at the same state point.¹³ For the 1D nets, we assume ideal harmonic springs, U = ∑_lk_l(x_l − R_l)²/2, with x_l being the distance between the connected particles, the reference length R_l of the springs being set to unity and the spring constants k_l being taken randomly from a uniform distribution of half-width δk centered around a mean value also set to unity. Only simple networks are presented here where two particles i − 1 and i along the chain are connected by one spring l = i, i.e., at zero temperatures all forces f_l along the chain become identical. This implies K ∼ 1/⟨1/k_l⟩. The compression modulus decreases thus strongly with δk as indicated by the bold line in the inset of Fig. 1. Our 2D systems are polydisperse Lennard-Jones (pLJ) beads⁷ kept at a constant pressure P = 2 as described in Ref. 9.

As shown in Fig. 1, the compression modulus K may be determined using the volume fluctuations in the NPT-ensemble, Eq. (1), or using Rowlinson's stress fluctuation formula, Eq. (2), for the NVT-ensemble. The same values of K are obtained from the Legendre transform for the pressure fluctuations, Eq. (5), and from the ideal pressure fluctuations η_id|_P, Eq. (11), which thus confirms both relations. As seen in the inset of Fig. 1, the compression modulus of the 1D nets decreases with δk. Also indicated is the “affine” contribution η|_P to K, measuring the mean spring constant ⟨k_l⟩ = 1, and the “non-affine” contribution η|_V which is seen to increase with δk. The decrease of K is thus due to the increase of the non-affine contribution.

As shown in the inset of Fig. 2, we have also checked the correlations between the ideal and the excess pressure fluctuations η_mix|_P. To make both models comparable, the reduced correlation function y = η_mix|_P/P_id is traced as a function of the reduced ideal pressure x = P_id/K with K as determined independently above. A perfect data collapse on the prediction y = 1 − x (bold line) is observed for all systems. The main panel of Fig. 2 shows finally the scaling of the Rowlinson formula computed in the NPT-ensemble. As before a scaling collapse of the data is achieved by plotting y = K_row|_P/P_id vs. x. The bold line indicates our key prediction, Eq. (3). Interestingly, the latter result does not depend on the MC-gauge which has been used above to simplify the derivation of Eq. (3). Please note that it is not possible to increase x beyond unity for our liquid systems (K ⩾ P_id) and the deviations from the low-temperature plateau y = 2 are thus necessarily small. The additional P_id/K correction has thus been overlooked in our previous publication.⁹ 

Emphasizing the underlying Legendre transform, Eq. (5), of the stress fluctuation formalism, we have investigated here the well-known Rowlinson stress fluctuation expression K_row for the compression modulus, Eq. (2), using deliberately the NPT-ensemble. Correcting several statements made in Ref. 9, it has been demonstrated theoretically and numerically that Eq. (3) holds. The latter result, as the other correlation relations indicated in the paper, may allow to readily calibrate (correctness, convergence, and precision) various barostats commonly used.⁸ 

P.P. thanks the Région Alsace and the IRTG Soft Matter and F.W. the DAAD for funding. We are indebted to A. Blumen (Freiburg) and A. Johner (Strasbourg) for helpful discussions.