Shear-strain and shear-stress fluctuations in generalized
Gaussian ensemble simulations of isotropic elastic networks
Shear-strain and shear-stress correlations in isotropic elastic bodies are investigated both theoretically and numerically at either imposed mean shear-stress () or shear-strain () and for more general values of a dimensionless parameter characterizing the generalized Gaussian ensemble. It allows to tune the strain fluctuations with being the inverse temperature, the volume, the instantaneous strain and the equilibrium shear modulus. Focusing on spring networks in two dimensions we show, e.g., for the stress fluctuations ( being the instantaneous stress) that with being the affine shear-elasticity. For the stress autocorrelation function this result is then seen (assuming a sufficiently slow shear-stress barostat) to generalize to with being the shear-stress relaxation modulus.
A “simple average” of an observable Allen and Tildesley (1994) does not depend on which thermodynamic ensemble it is measured in, at least not if the system is sufficiently large and if each ensemble samples indeed the same thermodynamic state point Lebowitz et al. (1967); Callen (1985); Allen and Tildesley (1994); Frenkel and Smit (2002); Landau and Binder (2000); Thijssen (1999). An example for such a simple average is the affine shear-elasticity characterizing the energy change of an affinely shear-strained elastic body Born and Huang (1954); Wittmer et al. (2013a, 2015a, 2015b) as properly defined below in sect. III. As one may verify numerically Wittmer et al. (2013a), it is irrelevant whether one computes in the -ensemble at constant particle number , volume , shear strain and temperature or in the conjugated -ensemble where the strain is allowed to fluctuate subjected to the constraint that the internal mean shear stress is imposed by the external applied shear stress Callen (1985). In contrast to this, the fluctuation of two observables and may depend on whether an extensive variable or its conjugated intensive variable is imposed Lebowitz et al. (1967); Callen (1985); Allen and Tildesley (1994). 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 rescaled shear strain and the conjugated intensive variable the shear stress . Using the Lebowitz-Percus-Verlet (LPV) transform Lebowitz et al. (1967) it is seen for the (rescaled) mean-squared fluctuation of the instantaneous shear stress that Wittmer et al. (2013a, 2015a)
with being the equilibrium shear modulus. For later convenience the -ensemble is indicated by “” and the -ensemble by “”. Generalizing eq. (1) in the time-domain it has been shown for the stress-stress correlation function that Wittmer et al. (2015a, b)
This assumes that the shear-barostat needed to sample the -ensemble must act very slowly on the time-scales used to probe . This may be realized equivalently by averaging over an ensemble of configurations with quenched strains distributed according to the -ensemble Wittmer et al. (2015a). It is due to the ergodicity breaking generated by the averaging over the quenched ensemble that stays finite for , while must vanish for large times (as one commonly expects for all correlation functions in ergodic systems Hansen and McDonald (2006)). By integration by parts it is seen that is equivalent to the experimentally important shear relaxation modulus Wittmer et al. (2015a). The transform eq. (2) thus implies , i.e. the relaxation modulus may be measured using the equilibrium modulus and the correlation function both determined in the -ensemble Wittmer et al. (2015a).
Generalized Gaussian ensemble.
As sketched in panel (b) of fig. 1, it is straightforward to interpolate between the -ensemble and the -ensemble by imposing an external field
with being the spring constant of the external harmonic spring. The standard -ensemble at is recovered by setting . Note that our approach is conceptually similar to the so-called “Gaussian ensemble” proposed some years ago by Hetherington and others Hetherington (1987); Costeniuc et al. (2006) generalizing the Boltzmann weight of the canonical ensemble by an exponential factor of the instantaneous energy . A similar external spring potential has been also used in the recent “elastic bath” approach by Workum and de Pablo van Workum and de Pablo (2003). Choosing the reference strain equal to the mean strain of the -system at vanishing shear stress allows to work at constant zero mean shear stress irrespective of the strength of the external potential foo (a). All the ensembles considered correspond thus to the same thermodynamic state, i.e. all first derivatives of the energy or the free energy Callen (1985) and all simple averages are identical. As sketched in panel (c) of fig. 1, is not necessarily positive, reducing the strain fluctuations, but may even become negative, which thus amplifies the fluctuations. It has been argued van Workum and de Pablo (2003) that this may allow a more convenient determination of the elastic modulus. Since the external spring is parallel to the system, the combined system and external device have an effective spring constant . Defining the dimensionless parameter the strain fluctuations of the combined system are thus given by
i.e. must vanish linearly with . -ensemble statistics is expected for , while -statistics should become relevant in the opposite limit for , i.e. . The system must become unstable in the limit , i.e. .
The aim of the present study is to generalize the relations for static fluctuations, eq. (1), and dynamical correlation functions, eq. (2), to the more general transformations between Gaussian ensembles characterized by the continuous parameter . In this way we want to make manifest that these transformation relations are generated by the continuous change of the constraint imposed on the fluctuations of the extensive variable. Focusing on spring networks in two dimensions we show, e.g., for the stress-stress fluctuations in different ensembles that
with being given by the affine shear-elasticity mentioned above. Assuming a very slowly acting shear-barostat, which is irrelevant for the system evolution for short times , the above result is then seen to generalize in the time-domain for the stress-stress correlation function
Since , eq. (6) allows the determination of the relaxation modulus from the equilibrium modulus and the stress-stress correlation function for any and . The upper time limit is seen to be set by the diffusion time of the instantaneous strain over the typical strain fluctuations for . As a consequence from eq. (5) and eq. (6) it appears that it is the equilibrium modulus of the system, which generates both transforms
with and for . Similar relations, allowing also to determine the creep compliance Witten and Pincus (2004); Rubinstein and Colby (2003), are obtained for strain-strain and strain-stress correlations functions.
The stated key relations eq. (5) and eq. (6) are justified theoretically in sect. II. The static fluctuations are discussed in sect. II.1 before we address the dynamical correlation functions in sects. II.2, II.3 and II.4 and the macroscopic linear response in sect. II.5. Algorithmic details are given in sect. III where the specific model system considered is introduced in sect. III.1. The canonical affine plane shear transformations used in this work are specified in sect. III.2, the instantaneous shear stress and the instantaneous affine shear-elasticity in sect. III.3. The zero-temperature ground state properties of the two-dimensional elastic network are summarized in sect. III.4. Some technical details related to the finite-temperature simulation of the Gaussian -ensemble using a Metropolis Monte Carlo (MC) scheme as a function of the maximum strain displacement rate , the second key operational parameter of this study, are given in sect. III.5. Section IV presents the numerical results obtained by molecular dynamics (MD) simulations at one finite (albeit small) temperature and one relatively high value of the friction constant of the Langevin thermostat used Allen and Tildesley (1994). The relevant static properties are described in sect. IV.1 before we turn to the dynamical strain-strain (sect. IV.2), strain-stress (sect. IV.3) and stress-stress (sect. IV.4) correlation functions. Our work is summarized in sect. V.
Ii Theoretical considerations
ii.1 Static fluctuations
We begin by restating the LPV transform in a convenient form assuming that the relevant extensive variable is the (rescaled) shear strain and the conjugated intensive variable the shear stress . Following Lebowitz et al. Lebowitz et al. (1967); Allen and Tildesley (1994), one verifies (see also sect. II.A of ref. Wittmer et al. (2013a)) that to leading order
for the transformation of with and from the -ensemble () to the -ensemble (). The more general transformation between arbitrary -ensembles can be found by reworking the saddle-point approximation of ref. Lebowitz et al. (1967) taking into account that the fluctuations around the peak of the distribution of the extensive variable is now not set by the modulus of the system but by the effective modulus . How this may be done can be seen in sect. 2.5 of ref. Wittmer et al. (2013b) for the volume as extensive variable and the (negative) pressure as intensive variable. Rewriting the latter result to the present case we get the generalized LPV transform
which reduces for to eq. (8).
Let us check the ensemble dependence of the rescaled strain-strain fluctuations foo (b). The generalized Gaussian ensemble corresponds to replacing the shear modulus of the system by the effective modulus . As already stated above, eq. (4), this leads to foo (c)
One verifies readily that the postulated LPV transform for general , eq. (9), is consistent with this result. To see this one sets and uses the fact that the strain fluctuations must vanish in the -ensemble, i.e. that the first term on the right hand-side of eq. (9) must vanish.
This result can be also obtained directly by replacing in the definition of the fluctuation by and using then the strain-strain relation, eq. (10).
We turn to the most important stress-stress fluctuation foo (b). We set now in the LPV transform. Since the stress fluctuations do not vanish in the -ensemble the corresponding term must now be included. This yields foo (c)
Interestingly, the contribution may be rewritten using the notation for the strain-average of a property with being the normalized strain-distribution for the considered -ensemble. The total mean stress of the ensemble is thus given by with being the average shear-stress of all configurations of shear-strain . Using that is a Gaussian and that with being the maximum of the distribution, it is then seen that
Using eq. (12) this implies in turn
as one expects to lowest order from a standard saddle-point approximation.
Stress-fluctuation formula for shear modulus.
Substracting the transform for from eq. (12) confirms the key result eq. (5). As shown in ref. Wittmer et al. (2015b), can be reduced by integration by parts to the affine shear-elasticity . Since the latter expression is a simple average, i.e. the same value is obtained in any ensemble, eq. (5) may be further simplified as foo (c)
It is thus sufficient to compute the material constants and in any ensemble to obtain the stress-stress fluctuations as a function of . Note that the special case for corresponds to the well-known stress-fluctuation formula Frenkel and Smit (2002); Squire et al. (1969); Mizuno et al. (2013); Wittmer et al. (2013a)
used in numerous numerical studies to compute the modulus conveniently in the -ensemble Barrat et al. (1988); Wittmer et al. (2002); Tanguy et al. (2002); Wittmer et al. (2013a, b); Mizuno et al. (2013); Flenner and Szamel (2015); Wittmer et al. (2015a, b).
ii.2 Dynamics: Definitions and general relations
We define now several dynamical observables of interest, which will be investigated numerically in sect. IV, and remind some well-known general relations Hansen and McDonald (2006); Doi and Edwards (1986); Allen and Tildesley (1994). Time translational symmetry and time reversal symmetry are assumed. We use and for either the instantaneous shear strain or the shear stress , and for their thermodynamic averages and and for their time-dependent fluctuations.
where the time translational invariance is used in the first step and the time reversal symmetry in the second step. Using again both symmetries one verifies that Hansen and McDonald (2006); Doi and Edwards (1986)
This implies that for large times foo (d). Albeit and thus contain the same information it will be sometimes better for theoretical or numerical reasons to focus on either or . Since for it will be natural, e.g., to consider the double-logarithmic representation of to clarify the power-law scaling of the correlations at short times.
Finite-time dependence of fluctuation estimation.
As defined in sect. II.1, is a thermodynamic ensemble average, hence, a time-independent property. In practice, may, however, often be estimated using a time series with a finite number of more or less correlated entries Allen and Tildesley (1994); Landau and Binder (2000). It is supposed here that these series are sampled with a constant time interval of length , i.e. corresponds to a time . With denoting such a finite average, one samples the -dependent observable foo (b)
with standing for an additional ensemble average over different time series of subsequent data points. We have for and for in an ergodic system. Interestingly, the detailed time-dependence of can be obtained from a weighted integral of the correlation function Landau and Binder (2000); Wittmer et al. (2015a). To see this we note first the identity
which can be verified by straightforward algebra. (See sect. 2.4 of ref. Doi and Edwards (1986).) Using time translational invariance this implies
where the weight stems from the finite length of the trajectory used. We change now to continuous time variables, , and replace the discrete sum by the time integral
Using eq. (20) this yields the general relation
which can be used to obtain if the correlation function is known.
Relaxation time .
Interestingly, the time integral in eq. (25) must become constant in all cases considered below where vanishes ultimately for . (This applies if a finite shear-barostat is switched on.) Defining the characteristic relaxation time
this leads to
ii.3 Dynamics: Additional model assumptions
Let us suppose that the MSD is diffusive for short times, i.e. for with being a diffusion constant. Matching this short time regime with for gives the possibility to operationally define as the crossover time
Different short-time dynamics may suggest of course a different operational definition. Let us further assume an exponentially decaying correlation function with . Using eq. (20) this is seen to be consistent with a diffusive short-time regime and eq. (28). It follows by integration from eq. (25) that
using the Debye function well-known in polymer science Doi and Edwards (1986). For large reduced times this leads to , i.e. by comparison with eq. (27) we have for an exponentially decaying correlation function.
As we shall see in sect. IV.4, it may also occur that the correlation function is more or less constant between a local (barostat independent time) up to a (barostat dependent) time , i.e. essentially with being a constant and the Heaviside function. It follows from eq. (26) that
i.e. to leading order.
ii.4 Dynamics: Ensemble effects
We have just stated various general relations applying to all ensembles. These relations do, however, not allow to relate the dynamical correlations between different . Since some ensembles are more readily computed than others, it would be useful to have a transformation relation such as the LPV transform for static fluctuations considered in sect. II.1. We emphasize that the correlation functions depend in general on the ensemble, i.e. the parameter , and on the dynamics of the shear-barostat. Since in the end we want to describe the intrinsic relaxation dynamics of the system, it is sufficient to focus on the limit where the barostat becomes very slow such that it becomes essentially irrelevant for the system evolution below an upper time . We consider the limit where this time is be much larger than any intrinsic relaxation time of the system. Since we want still to consider meaningful thermal averages with respect to the chosen ensemble, we need either trajectories of a time interval much larger than to allow the full sampling of the phase space or, equivalently, we need to average over independent start configurations representing the ensemble. Under these two assumptions useful transformation relations can be formulated by reworking the generalization of the LPV transform for dynamical correlations replacing the system modulus by the effective modulus . Being beyond the scope of this paper, we proceed by postulating the central scaling relation for the MSD and argue briefly why this relation is natural. We discuss then in turn the strain-strain correlations , the strain-stress correlations and the stress-stress correlations for different . Alternative, more direct ways to derive the relations are mentioned.
We postulate that under the two assumptions made above the MSD does not depend on the ensemble
i.e. the MSD behaves as a simple mean. This scaling postulate is justified by two facts. Firstly, the barostat is (by construction) too weak to change for the evolution of the system. Secondly, that the starting points of the trajectory at are distributed according to the considered ensemble must become an irrelevant higher order effect (vanishing rapidly with the system volume ), since the MSD probes displacements with respect to the starting points, not their absolute values. The fundamental scaling, eq. (31), implies using eq. (20) that
where only the first term on the right hand-side depends on . This leads us finally to the general transformation relation between correlation functions
where stands for the -ensemble with The two static contributions on the right hand-side can be further simplified using results from sect. II.1. This demonstrates the linear relation .
This result is directly obtained by setting in the definition of . We emphasize that it is also consistent with the LPV transform, eq. (9), as one verifies by setting , , and and using that the strain fluctuation term for must vanish.
This result is also obtained from the LPV transform setting and and using again that . The postulated scaling eq. (31) has thus led again to a reasonable result.
Interestingly, as one may see from the -ensemble limit considered in Wittmer et al. (2015a, b), the stress-stress MSD is not expected to simply vanish for as in the two previous cases. (The instantaneous stress fluctuates even at a fixed strain .) However, eq. (33) still holds leading to foo (c)
where the static term on the right hand-side has been simplified using eq. (5). This confirms the key relation eq. (6) announced in the Introduction. Note that the correlation function in the -ensemble must vanish beyond some local time scale which does depend on the network considered but, of course, not on the shear-barostat. For a sufficiently slow barostat the correlation function thus becomes constant
where we have dropped . Using eq. (32) this leads to the remarkable relation
which must hold for all . We note finally that while for the MSD for , holds for all times for according to stress-fluctuation formula, eq. (16).
ii.5 Dynamics: Macroscopic linear response
The experimentally important macroscopic linear response measured by the creep compliance and the shear relaxation modulus Doi and Edwards (1986); Witten and Pincus (2004); Rubinstein and Colby (2003) may be obtained conveniently in an equilibrium simulation at a given using some of the correlation functions discussed above. Please note that being material properties of the given state point (experimentally obtained using a simple average, not a fluctuation) both response functions and do, of course, not depend on .
Let us first consider the creep compliance for . It is assumed that for the system is at thermal equilibrium and the internal mean stress equals the applied external stress of the -ensemble. After imposing at a small increment , the creep compliance measures the ensuing average strain increment . We note en passant that the average internal shear stress does neither immediately reach the new equilibrium value but shows a similar time dependence as the strain. Reworking the arguments put forward by Doi and Edwards, see eq. (3.67) of ref. Doi and Edwards (1986), it is seen that
i.e. the creep compliance is most readily computed using the strain-strain MSD in the -ensemble. As described in sect. III.5, we shall change the shear-strain using a MC shear-barostat which corresponds to a perfectly viscous dynamics. One thus expects for the strain-strain correlations. Using for and eq. (20) this suggest
in agreement with the Kelvin-Voigt model Rubinstein and Colby (2003) representing a purely viscous damper and a purely elastic spring connected in parallel.
Shear relaxation modulus.
The shear relaxation modulus may be obtained from the stress increment for after a small step strain with has been imposed at time . It is well known that the components of the Fourier transformed relaxation modulus , the storage modulus and the loss modulus , are directly measurable in an oscillatory shear strain experiment Rubinstein and Colby (2003); Wittmer et al. (2015b). As seen, e.g., by eq. (32) in ref. Wittmer et al. (2015b), it can be demonstrated by integration by parts that
i.e. the barostat should be irrelevant on the time scales considered. Using the transformation relation between different ensembles, eq. (36), the relaxation modulus may equivalently be obtained for other according to
which holds for all times since the barostat becomes irrelevant for . Note that eq. (43) may also be derived directly (without using the transformation relation between different ensembles) using Boltzmann’s superposition principle for an arbitrary strain history and the standard fluctuation-dissipation theorem for the after-effect function Hansen and McDonald (2006) as shown elsewhere Wittmer et al. (2015a, b). Two immediate consequences of eq. (43) are (i) that only becomes equivalent to for in the liquid limit where (trivially) and (ii) that the shear modulus is only probed by on time scales where must vanish. In principle, it is thus impossible to obtain the static shear modulus of an elastic body only from as often incorrectly assumed Klix et al. (2012); Flenner and Szamel (2015).
Iii Algorithmic details
iii.1 Model Hamiltonian
To illustrate our key relations we present in sect. IV numerical data obtained using a periodic two-dimensional () network of harmonic springs connecting vertices. The model Hamiltonian is given by the sum of a kinetic energy contribution
with being the velocity of vertex and its (assumed) monodisperse mass, and an excess potential
where denotes the spring constant, the reference length and the length of spring . The sum runs over all springs connecting pairs of beads and with at positions and . The vertex mass and Boltzmann’s constant are set to unity and Lennard-Jones (LJ) units Allen and Tildesley (1994) are assumed.
iii.2 Canonical affine shear transformations
While the box volume is kept constant throughout this work, we shall frequently change the shape of the simulation box. As sketched in panel (a) of fig. 1, we perform plane shear transformations of the instantaneous shear strain with an essentially infinitesimal strain increment . We assume that not only the box shape is changed but that the particle positions (using the principal box convention Allen and Tildesley (1994)) follow the imposed macroscopic constraint in an affine manner according to
with all other coordinates remaining unchanged. Albeit not strictly necessary for the demonstration of our key relations, we assume, moreover, that this shear transformation is also canonical Goldstein et al. (2001); Frenkel and Smit (2002). This implies that the -component of the velocity must transform as Wittmer et al. (2015b)
iii.3 Shear stress and affine shear-elasticity
Let denote the system Hamiltonian of a configuration originally at strained using a canonical affine transformation to compactly written as a function of the strain increment . The instantaneous shear stress and the instantaneous affine shear-elasticity may be defined as the expansion coefficients associated to the energy change
with being the reference, i.e. foo (f)
The derivatives and with respect to may be computed as shown in sect. 2.1 of Wittmer et al. (2015b). Similar relations apply for the corresponding contributions and to and for the contributions and to . Using eq. (44) this implies Wittmer et al. (2015b)
for the ideal contributions to the shear stress and the affine shear-elasticity. Note that the minus sign for the shear stress is due to the minus sign in eq. (47) required for a canonical transformation. For the excess contributions one obtains Wittmer et al. (2015b)
with being the normalized distance vector between the particles and . Interestingly, eq. (53) is strictly identical to the corresponding off-diagonal term of the Kirkwood stress tensor Allen and Tildesley (1994). The last term in eq. (54) automatically takes into account the finite normal pressure of the system.
iii.4 Groundstate characterization
As explained elsewhere Wittmer et al. (2013a, 2015b), the specific elastic network used in this work has been constructed using the dynamical matrix of a quenched polydisperse LJ bead glass, i.e. at low temperatures our network has exactly the same mechanical and vibrational properties as the original discrete particle system. Prior to forming the network the bead system was cooled down to using a constant quenching rate and imposing a normal pressure . The original LJ beads are represented in the snapshot shown in fig. 2 by grey polydisperse circles, the permanent spring network created from the quenched bead system by lines between vertices. The dark (black) lines indicate repulsive forces between the vertices, while the light (red) lines represent tensile forces. The line width is proportional to the tension (repulsion). Note that the force network is strongly inhomogeneous with zones of weak attractive links embedded within a strong repulsive skeleton as already discussed in refs. Wittmer et al. (2002); Tanguy et al. (2002). Only a small subvolume of the network is represented. The total periodic box of linear length contains vertices and springs. The monomer density is close to unity.
Finite shear stress .
Due to the construction of the network the total force acting on each vertex of the reference network must vanish at . As seen in the snapshot, this does not imply that the repulsive and/or tensile forces transmitted along the springs must also vanish. Due to the periodic boundary conditions and the constant-strain constraint () the shear stress does not necessarily vanish. For a pair potential such as eq. (45) the relevant excess contribution of the shear stress is readily computed using the Kirkwook expression, eq. (53). As shown in the main panel of fig. 2 (horizontal arrow), it turns out that for the specific network we use throughout this work we have a finite shear stress .
Affine shear-elasticity .
Let us consider a small affine shear strain according to eq. (46). As show in panel (b) of fig. 2, one may now compute using eq. (53) the shear stress for different (squares). As one expects from the definition of the affine shear-elasticity coefficient , eq. (48), this yields the linear relation with as indicated by the dashed line. The same coefficient is also obtained directly from the unstrained configuration using eq. (54).
Equilibrium shear modulus .
The forces acting on the vertices of an affinely strained network do not vanish in general and the system is normally not at mechanical equilibrium. As described elsewhere Tanguy et al. (2002); Wittmer et al. (2015b), we relax these forces by first applying a steepest descend algorithm, i.e. by imposing displacements proportional to the force, and then by means of the conjugate gradient method Thijssen (1999). The ensuing non-affine displacements of the vertices decrease the energy Wittmer et al. (2015b) and the magnitude of the shear stress as may be seen from the large circles indicated in fig. 2. As indicated by the bold solid line these final stresses scale again linearly as
with and . The affine coefficient has thus been replaced by the much smaller equilibrium shear modulus of the groundstate 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 descend scheme) without any external force applied, the system relaxes to .