Shear modulus of simulated glass-forming model systems: Effects of boundary condition, temperature and sampling time

The shear modulus G of two glass-forming colloidal model systems in d=3 and d=2 dimensions is investigated by means of, respectively, molecular dynamics and Monte Carlo simulations. Comparing ensembles where either the shear strain gamma or the conjugated (mean) shear stress tau are imposed, we compute G from the respective stress and strain fluctuations as a function of temperature T while keeping a constant normal pressure P. The choice of the ensemble is seen to be highly relevant for the shear stress fluctuations mu_F(T) which at constant tau decay monotonously with T following the affine shear elasticity mu_A(T), i.e. a simple two-point correlation function. At variance, non-monotonous behavior with a maximum at the glass transition temperature T_g is demonstrated for mu_F(T) at constant gamma. The increase of G below T_g is reasonably fitted for both models by a continuous cusp singularity, G(T) is proportional to (1-T/T_g)^(1/2), in qualitative agreement with some recent replica calculations. It is argued, however, that longer sampling times may lead to a sharper transition. The additive jump discontinuity predicted by mode-coupling theory and other replica calculations thus cannot ultimately be ruled out.

The shear modulus G of two glass-forming colloidal model systems in d = 3 and d = 2 dimensions is investigated by means of, respectively, molecular dynamics and Monte Carlo simulations. Comparing ensembles where either the shear strain γ or the conjugated (mean) shear stress τ are imposed, we compute G from the respective stress and strain fluctuations as a function of temperature T while keeping a constant normal pressure P . The choice of the ensemble is seen to be highly relevant for the shear stress fluctuations µF(T ) which at constant τ decay monotonously with T following the affine shear elasticity µA(T ), i.e. a simple two-point correlation function. At variance, non-monotonous behavior with a maximum at the glass transition temperature Tg is demonstrated for µF(T ) at constant γ. The increase of G below Tg is reasonably fitted for both models by a continuous cusp singularity, G(T ) ∝ (1 − T /Tg) 1/2 , in qualitative agreement with some recent replica calculations. It is argued, however, that longer sampling times may lead to a sharper transition. The additive jump discontinuity predicted by mode-coupling theory and other replica calculations thus cannot ultimately be ruled out. Stress fluctuation formalism. Among the fundamental properties of any solid material are its elastic constants [1][2][3]. They describe the macroscopic response to weak external deformations and encode information about the potential landscape of the system. In their seminal work [4] Squire, Holt and Hoover derived an expression for the isothermal elastic constants extending the classical Born theory [1] to canonical ensembles of imposed particle number N , volume V and finite (mean) temperature T . They found a correction term C αβγδ F to the Born expression for the elastic constants which involves the mean-square fluctuations of the stress [4][5][6][7] C αβγδ with β = 1/k B T being the inverse temperature, k B Boltzmann's constant,σ αβ the instantaneous value of a component of the stress tensor σ αβ = σ αβ and Greek letters denoting the spatial coordinates in d dimension. Perhaps due to the fact that the demonstration of the "stress fluctuation formalism" in Ref. [4] assumed explicitly a well-defined reference position and a displacement field for the particles, it has not been appreciated that this approach is consistent (at least if the initial pressure is properly taken into account) with the well-known pressure fluctuation formula for the compression modulus K of isotropic fluids [6][7][8][9] which has been derived several * Electronic address: joachim.wittmer@ics-cnrs.unistra.fr decades earlier, presumably for the first time in the late 1940s by Rowlinson [8]. In any case, the stress fluctuation formalism provides a convenient route to calculating the elastic properties in computer simulations which has been widely used in the past [7,[10][11][12][13][14][15][16][17][18][19][20][21]. It has been generalized to systems with nonzero initial stress [10], hard-sphere interactions [12] and arbitrary continuous potentials [5], or to the calculation of local mechanical properties [20].
Low-temperature limit and non-affine displacements. In particular from Ref. [5] it has become clear that the stress fluctuation correction to the leading Born term does not necessarily vanish in the zero-temperature limit. This is due to the fact that when a system is subjected to a homogeneous deformation, the ensuing particle displacements need not follow the macroscopic strain affinely [13][14][15][16]18]. The stress fluctuation term quantifies the extent of these non-affine displacements, whereas the Born term reflects the affine part of the particle displacements. How important the non-affine motions are, depends on the system under consideration [18]. While the elastic properties of crystals with one atom per unit cell are given by the Born term only, stress fluctuations are significant for crystals with more complex unit cells [5]. They become particularly pronounced for polymer like soft materials [21] and amorphous solids [7, 11, 13-18, 20, 22]. This is revealed by recent simulation studies of various glass formers [11], like Lennard-Jones (LJ) mixtures [13][14][15][16], polymer glasses [7,20] or silica melts [22]. Many of these simulation studies concentrate on the mechanical properties deep in the glassy state, exploring for instance correlations between the non-affine displacement field and vibrational anomalies of the glass ("Boson peak") [16,22], the mechanical heterogeneity at the nanoscale [20,23] or the onset of molecular plasticity in the regime where the macroscopic deformation is still elastic [24].
Shear stress and shear modulus. Surprisingly, the temperature dependence of the elastic constants of glassy materials does not appear to have been investigated extensively using the stress fluctuation formalism [7,11,19]. Following the pioneering work by Barrat et al. [11], we present in this work such a study where we focus on the behavior of the shear modulus G of two isotropic colloidal glass-formers in d = 3 and d = 2 dimensions sampled, respectively, by means of molecular dynamics (MD) and Monte Carlo (MC) simulations [6,9,25]. Since we are interested in isotropic systems one can avoid the inconvenient tensorial notation of the full stress fluctuation formalism as shown in Fig. 1. Focusing on the shear in the xy-plane we use τ ≡ σ xy for the mean shear stress andτ ≡σ xy for its instantaneous value. As an important contribution to the shear modulus we shall monitor the shear stress fluctuations as a function of temperature T . Quite generally, the shear modulus G may then be computed according to the stress fluctuation formula [7,11] where µ A stands for the so-called "affine shear elasticity" corresponding to an assumed affine response to an external shear strain γ. For pairwise interaction potentials it can be shown that i.e. µ A comprises the so-called "Born-Lamé coefficient" µ B , a simple moment of the first and the second derivative of the pair potential with respect to the distance of the interacting particles [1, 7,18], and the excess contribution P ex to the total pressure P . These notations will be properly defined in Sec. II B below where a simple derivation will be given. Since Eq. (3) is not the only operational definition of the shear modulus discussed in this work, an index has been added reminding that the shear stressτ is used for the determination of G.
Physical motivation and systems of interest. If an infinitesimal shear stress τ is applied at time t = 0, an isotropic solid exhibits, quite generally, an elastic response quantified in strength by a finite static (tindependent) shear modulus (G > 0) as shown by the top curve in Fig. 1(b), whereas a liquid has a vanishing modulus (G = 0) and flows on long time scales as shown by the bottom curve. The shear modulus G can thus serve as an "order parameter" distinguishing liquid and solid states [26][27][28][29][30][31]. Upon melting, the shear modulus of crystalline solids is known to vanish discontinuously with T . A natural question which arises is that of the behavior of G(T ) for amorphous solids and glasses in the vicinity of the glass transition temperature T g . Interestingly, qualitative different theoretical predictions have been put forward by mode-coupling theory (MCT) [31,32] and some versions of replica theory [30] predicting a discontinuous jump at the glass transition whereas other versions of replica theory [28,29] suggest a continuous transition in the static limit of large sampling times t. (See Ref. [29] for a topical recent discussion of the replica approach.) The numerical characterization of G(T ) for colloidal glass-former is thus of high current interest. Since the stress fluctuation formalism can readily be added to the standard simulation codes of colloidal glass-formers [11,18,33] at imposed particle number N , box volume V , box shape (γ = 0) and temperature T , called NVγT-ensembles, this suggests the use of Eq. (3) to tackle this issue. Note that the latter stress fluctuation relation should also be useful for wide range of related systems, beyond the scope of the present work, including the glass transition of polymeric liquids [7], colloidal gels [34] and self-assembled networks [26] as observed experimentally for hyperbranched polymer chains with sticky end-groups [35] or bridged networks of telechelic polymers in water-oil emulsions [36,37].
Validity through a solid-liquid transition. This begs the delicate question of whether the general stress fluctuation formalism and especially Eq. (3) for the shear modulus only holds for elastic solids with a well-defined displacement field, as suggested in the early literature [4][5][6], or if the approach may actually be also applied through a solid-liquid transition up to the liquid state [8,21] with µ F → µ A and thus G = G τ τ → 0 as it should for a liquid. By reworking the theoretical arguments leading to Eq. (3) it will be seen that this is in principal the case, at least if standard thermodynamics can be assumed to apply. The latter point is indeed not self-evident for the colloidal glasses we are interested in where (i) some degrees of freedom are frozen on the time scale probed and (ii) the measured response to an imposed infinitesimal shear τ may depend on the sampling time t, especially for temperatures T around T g as shown in Fig. 1(b). To be on the safe side, it thus appears to be appropriate to test numerically the applicability of the stress fluctuation formula, Eq. (3), using the corresponding strain fluctuation relations which are more fundamental on experimental grounds.
Strain fluctuation relations. Panel (a) of Fig. 1 shows a schematic setup for probing experimentally the shear modulus G of isotropic systems confined between two rigid walls. We assume that the bottom wall is fixed and the particle number N , the (mean) normal pressure P and the (mean) temperature T are kept constant. In linear response either a shear strain γ = tan(α) or a (mean) shear stress τ may be imposed. The thin solid line indicates the original unstrained reference system with γ = 0 and τ = 0, the dashed line the deformed state. From both the mechanical and the thermodynamical point of view, the shear modulus is defined as where the derivative may also be taken at τ = 0. In practice, one obtains G by fitting a linear slope to the (γ, τ )-data at vanishing strain. Obviously, the noise level is normally large for small strains. From the numerical point of view it would thus be more appropriate to work at imposed mean stress τ = 0, to let the strain freely fluctuate and to sample the instantaneous values (γ,τ ). The shear modulus G is then simply obtained by linear regression where the index indicates that both the instantaneous valuesγ andτ are used. (Note that imposing the strain γ also fixes the instantaneous strainγ which thus cannot fluctuate.) The associated dimensionless correlation coefficient allows to determine the quality of the fit. We do not know whether Eq. (6) and Eq. (7) have actually been used in a real experiment, however since no difficulty arises to work in a computer simulation at constant τ and to record (γ,τ ), we use this linear regression as our key operational definition of G which is used to judge the correctness of the stress fluctuation relation, Eq. (3), computed at constant γ. To avoid additional physics at the walls, we use of course periodic boundary conditions [9] (the central image being sketched) as shown by the panels on the right-hand side of Fig. 1. As discussed in Sec. II C, it is possible to determine G γτ at constant volume (NVτ T-ensemble) as shown by panel (c) or constant normal pressure (NPτ T-ensemble) as shown by panel (e). The latter ensemble has the advantage that one can easily impose in addition the same normal pressure P for all temperatures T as it is justified for experimental reasons. If thermodynamics holds, it is known (Sec. II A) that should hold [38]. Using Eq. (6) this implies a second operational definition which has the advantage that only the instantaneous shear strainγ has to be recorded. The main technical point put forward in this paper is that all operational definitions yield similar results even for temperatures where a strong dependence on the sampling time t is observed. In other words, Eq. (10) states that this time dependence applies similarly (albeit perhaps not exactly) to the various moments computed for the three operational definitions of G in different ensembles.
Outline. The paper is organized as follows: Several theoretical issues are regrouped in Sec. II. We begin the discussion in Sec. II A by reminding a few useful general thermodynamic relations. A simple derivation of the stress fluctuation formula, Eq. (3), for the shear modulus G τ τ in constant-γ ensembles is given in Sec. II B. (The analogous stress fluctuation formula for the compression modulus K in constant-V ensembles [8] is rederived in Appendix A 2.) Shear stress fluctuations µ F in ensembles with different boundary conditions are discussed in Sec. II C. For constant-τ ensembles µ F = µ F | τ will be shown to reduce to the affine shear elasticity µ A , i.e. the functional G τ τ must vanish for all temperatures. As pointed out recently [39], the truncation of the interaction potentials used in computational studies implies an impulsive correction for the affine shear elasticity µ A . This is summarized in Sec. II D. The numerical models used in this study are presented in Sec. III. We turn then in Sec. IV to our main numerical results. Starting with the high-T liquid limit, Sec. IV A shows that all our operational definitions correctly yield a vanishing shear modulus. Moving to the opposite low-T limit we present in Sec. IV B the elastic response for a topologically fixed network of harmonic springs constructed from the "dynamical matrix" of our 2D model glass system [13,14]. As shown in Sec. IV C, the same behavior is found qualitatively for low-T glassy bead systems. The temperature dependence of stress fluctuations and shear moduli in different ensembles is then discussed in Sec. IV D. As stated above, one central goal of this work is to determine the behavior of the shear modulus around T g . We conclude the paper in Sec. V. Thermodynamic relations and numerical findings related to the compression modulus K are regrouped in Appendix A.

A. Reminder of thermodynamic relations
We assume in the following that standard thermodynamics and thermostatistics [40] can be applied and remind a few useful relations. Let us denote an extensive variable by X and its instantaneous value byX, the conjugated intensive variable by I and its instantaneous value byÎ [41]. With F (T, X) being the free energy at temperature T and imposed X, the intensive variable I and the associated modulus M are given by Alternatively, one may consider the Legendre transformation H(T, I) ≡ F (T, X) − XI for which In our case, the extensive variable is X = V γ, the intensive variable I = τ [42] and the modulus M = G [43]. We emphasize that there is an important difference between extensive and intensive variables: If the extensive variable X is fixed,X cannot fluctuate. If, on the other hand, the intensive variable I = Î is imposed, not onlŷ X but alsoÎ may fluctuate [44]. In our case this means, thatγ ≡ γ for the NVγT-ensemble, while the shear stress fluctuations µ F defined above do clearly not vanish in the NVτ T-ensemble shown in Fig. 1(d). Assuming I to be imposed, it is well known that [40] δX(βδÎ) with β = 1/k B T being the inverse temperature. For our variables X = V γ and I = τ , Eq. (15) implies Eq. (8) from which it is seen that Eq. (6) is consistent with Eq. (9). The latter relation is also obtained directly using Eq. (16).
As discussed in textbooks [9,40] a simple average A = Â of some observable A does not depend of the chosen ensemble, at least not if the system is large enough (N, V → ∞). A correlation function δÂδB of two observables A and B may differ, however, depending on whether X or I are imposed. As demonstrated by Lebowitz, Percus and Verlet in 1967 [45] one verifies that .
and that both fluctuations must become similar in the limit where the modulus M becomes small. For the shear stress fluctuations µ F , Eq. (18) corresponds to the transformation and one thus expects µ F to be larger for NVτ T-systems than for NVγT-systems, while in the liquid limit where G → 0 the imposed boundary conditions should not matter. We shall verify Eq. (20) numerically below. A glance at the stress fluctuation formula, Eq. (3), suggests the identity Assuming Eq. (8) and Eq. (21), it follows for the correlation coefficient c γτ defined above that If the measured (γ,τ ) were indeed perfectly correlated, i.e. c γτ = 1, this implies G = µ A and thus µ F | V γ = 0. In fact, for all situations studied here we always have i.e. the affine shear elasticity sets only an upper bound to the shear modulus and a theory which only contains the affine strain response must overpredict G. We show now directly that Eq. (3) and, hence, Eqs. (21,22,23) indeed hold and this under quite general conditions.

B. Shear modulus at imposed shear strain
Introduction. As shown in Fig. 1(a), we demonstrate Eq. (3) from the free energy change associated with an imposed arbitrarily small pure shear strain γ in the xyplane assuming a constant particle number N , a constant volume V and a constant temperature T (NVγTensemble). It is supposed that the total system Hamiltonian may be written as the sum of a kinetic energy and a potential energy. Since for a plain shear strain at constant volume the ideal free energy contribution does not change, i.e. is irrelevant for τ and G, we may focus on the excess free energy contribution due to the conservative interaction energy of the particles. The (excess) partition function Z ex (0) of the unperturbed system at γ = 0 is the Boltzmann-weighted sum over all states s of the system which are accessible within the measurement time t. The argument (0) is a short-hand for the unstrained reference. Following the derivation of the compression modulus K by Rowlinson [8] the partition function of the sheared system is supposed to be the sum over the same states s, but with a different metric [46] corresponding to the macroscopic strain which changes the total interaction energy U s (γ) of state s and, hence, the weight of the sheared configuration for the averages computed. This is the central hypothesis made. Interestingly, it is not necessary to specify explicitly the states of the unperturbed or perturbed system, e.g., it is irrelevant whether the particles are distinguishable or not or whether they have a well-defined reference position for defining a displacement field. Shear stress and modulus for general potential. Assuming Eq. (24) to hold we compute now the mean shear stress τ and the shear modulus G for a general interaction potential U s (γ). We note for later convenience that for the derivatives of the free energy and for the derivatives of the excess partition function where a prime denotes the derivative of a function f (x) with respect to its argument x. Using Eq. (28) and taking finally the limit γ → 0 one verifies for the shear stress that defining the instantaneous shear stress [42]. Note that the average taken is defined as using the weights of the unperturbed system. The shear stress thus measures the average change of the total interaction energy U s (γ) taken at γ = 0. The shear modulus G is obtained using in addition Eq. (27) and Eq. (29) and taking finally the γ → 0 limit. Confirming thus the stress fluctuation formula Eq. (3) stated in the Introduction, this yields with µ A being the "affine shear elasticity" already mentioned above and µ F as defined in Eq.
(2). The comparison of Eq. (32) and Eq. (20) confirms Eq. (21). Comment. The affine shear elasticity µ A corresponds to the change (second derivative) of the total energy which would be obtained if one actually strains affinely in a computer simulation a given state s without allowing the particles to relax their position. As shown for athermal (T → 0) amorphous bodies [13,14,18], the positions of the particles of such a strained configuration will of course in general change slightly to minimize the interaction energy relaxing thus the elastic moduli. This is also of relevance for thermalized solids where the nonaffine displacements of the particles are driven by the minimization of the free energy. It is for this reason that the shear-stress fluctuation term µ F ≥ 0 must occur in Eq. (32) correcting the overprediction of the shear modulus G by the affine shear elasticity µ A , Eq. (23). This point has been overlooked in the early literature [1] and only appreciated much later [4,5,8,13,14,47] as discussed in Barrat's review [18]. Interestingly, as has been shown by Lutsko [5], µ F and other similarly defined stress fluctuations become temperature independent in the harmonic ground state approximation for T → 0. Probing the stress fluctuations in a low-temperature simulation allows thus to determine the elastic moduli of athermal solids.
Please note that up to now we have not used explicitly the coordinate transformations (metric change) associated with an affine shear strain. This is needed to obtain operational definitions for the instantaneous shear stresŝ τ (and thus for τ and µ F ) and the affine shear elasticity.
Coordinate transformation. As shown in Fig. 1(a), we assume a pure shear strain which transforms the xcoordinate of a particle position or the distance between two particles as leaving all other coordinates unchanged [46]. The squared distance r 2 between two particles thus transforms as r 2 (0) ⇒ r 2 (γ) = r 2 (0) + 2γx(0)y(0) + γ 2 y 2 (0) (35) where we need to keep the γ 2 -term for calculating the shear modulus below. We note for later reference that this implies dr 2 (γ) dγ = 2x(0)y(0) + 2γy 2 (0) (36) d 2 r 2 (γ) dγ 2 = 2y 2 (0) (37) for, respectively, the first and the second derivative with respect to γ. Let us now consider an arbitrary function f (r(γ)). Using Eq. (36) it is readily seen that its first derivative with respect to γ may be written as In the second step we have introduced the components n x = x(0)/r(0) and n y = y(0)/r(0) of the normalized distance vector between both particles. More generally, we denote by n α the α-component of a normalized vector with α = x, y, . . . Stating only the small-γ limit for the second derivative of f (r(γ)) with respect to γ we note finally that lim γ→0 ∂ 2 f (r(γ)) ∂γ 2 = r 2 f ′′ (r) − rf ′ (r) n 2 x n 2 y + rf ′ (r) n y n y where we have dropped the argument (0) for the distance r of the unperturbed reference system.

Pair interaction potentials.
We assume now that the interaction energy U s (γ) of a configuration s is given by a pair interaction potential u(r) between the particles [6,9] with r(γ) being the distance between two particles, i.e.
where the index l labels the interaction between the particles i and j with i < j. Using the general result Eq. (30) for the shear stress it follows using Eq. (38) that where we have dropped the argument (0) in the γ → 0 limit. We have thus rederived for the shear stress (α = x, β = y) the well-known Kirkwood expression for the general excess stress tensor [9] σ ex The affine shear elasticity µ A is obtained using Eq. (33) and Eq. (40). This yields with the first contribution µ B being the so-called "Born- It corresponds to the first term in Eq. (40). The second contribution σ ex yy stands for the normal excess stress in the y-direction, i.e. the α = β = y component of the excess stress tensor.
Isotropic systems. In this work we focus on isotropic materials under isotropic external loads. The normal (excess) stresses σ ex αα for all directions α are thus identical, i.e.
with P ex ≡ P − P id being the excess pressure, P the total pressure, P id = k B T ρ the ideal gas pressure and d the spatial dimension. This implies that the affine shear elasticity is given by µ A = µ B − P ex as already stated in the Introduction, Eq. (4). We note that for isotropic d-dimensional systems it can be shown [21,39] that the Born-Lamé coefficient can be simplified as The d-dependent prefactor stems from the assumed isotropy of the system and the mathematical formula [48] (n α n β ) (δ αβ being the Kronecker symbol [48]) for the components of a unit vector in d dimensions pointing into arbitrary directions. By comparing the excess pressure P ex , Eq. (46), with the underlined second contribution to the Born coefficient in Eq. (47) this implies that It is thus inconsistent to neglect the explicit excess pressure in Eq. (4) while keeping the second term of the Born-Lamé coefficient µ B . This approximation is only justified if the excess pressure P ex is negligible (and not the total pressure P ).

C. Shear stress fluctuations in different ensembles
Different ensembles. Being simple averages neither P ex , µ A or µ B does depend for sufficiently large systems on the chosen ensemble, i.e. irrespective of whether γ or τ is imposed one expects to obtain the same values. As already reminded in Sec. II A, this is different for fluctuations in general [9]. For instance, it is obviously pointless to use the shear-strain fluctuation formulae (20), one actually expects to observe G τ τ = 0 for all temperatures T if our thermodynamic reasoning holds through the glass transition up to the liquid state. We will test numerically this non-trivial prediction in Sec. IV.
Up to now we have only considered for simplicity of the presentation NVγT-and NVτ T-ensembles at constant volume V . Since on general grounds pure deviatoric and dilatational (volumetric) strains are decoupled (both strains commute), one expects the observable G τ τ also to be applicable in the NPγT-ensemble shown in Fig. 1(d) and the observables G γτ and G γγ to be applicable in the NPτ T-ensemble as sketched in Fig. 1(e). In the latter ensemble one also expects to find G τ τ = 0 and µ F | τ = µ A . This will also be tested below.
Direct derivation of µ F for imposed τ . We demonstrate now directly Eq. (21) for the shear stress fluctuations using the general mathematical identity with x being an unconstrained coordinate, A(x) some property, f (x) ≡ −u ′ (x) a "force" with respect to some "energy" u(x) and the average being Boltzmann weighted, i.e. . . . ∝ dx . . . e −βu(x) . (It is assumed in Eq. (50) that u(x) decays sufficiently fast at the boundaries.) Following work by Zwanzig [49] we express first the instantaneous shear stressτ given above by the Kirkwood formula, Eq. (42), by the alternative virial repre- The second term, sometimes called the "inner virial", stands for a sum over all particles i with y i being the y-coordinate of the particle and f x,i the x-coordinate of the force acting on the particle. This contribution indeed vanishes on average as can be seen using Eq. (50), which shows that τ = τ as it must. The stress fluctuations may thus be expressed as where we have used the integration by part, Eq. (51). This step requires that all particle positions are unconstrained and independent (generalized) coordinates. Note that this is possible at imposed τ , but cannot hold at fixed γ. Assuming then pair interactions and using Eq. (42) one can reformulate Eq. (53) within a few lines. Without further approximation this yields where the sum runs now over all interactions l between particles i < j. x l and y l refer to components of the distance vector of the interaction l and f x,l to the xcomponent of the central force f (r l ) = −u ′ (r l ) between a particle pair at a distance r l . Note that in Eq. (54) the stress fluctuations stemming from different interactions are decoupled, i.e. µ F | τ characterizes the self-or twopoint contributions of directly interacting beads. Since f x (r) = −u ′ (r)x/r and, hence, this is identical to the affine shear elasticity µ A derived at the end of Sec. II B. We have thus confirmed Eq. (21). Interestingly, the argument can be turned around: (1) The stress fluctuation contribution term µ F | τ = µ A is obtained by the simple derivation given in this paragraph.
As already stressed, our thermodynamic derivation of Eq. (3) is rather general, not using in particular a welldefined displacement field for the particle positions.

D. Impulsive corrections for truncated potentials
Truncation. It is common practice in computational condensed matter physics [6,9] to truncate a pair interaction potential u(r), with r being the distance between two particles i and j, at a conveniently chosen cutoff r c . This allows to reduce the number of interactions to be computed and energy or force calculations become O(N )processes. However, the truncation also introduces technical difficulties, e.g., instabilities in the numerical solution of differential equations as investigated especially for the MD method [9,50]. Without restricting much in practice the generality of our results, we assume below that • the pair potential u(r) is short-ranged, i.e. that it decays within a few particle diameters, • it scales as u(r) ≡ u(s) with the reduced dimensionless distance s = r/σ l where σ l characterizes the length scale of the interaction l and • the same reduced cutoff s c = r c /σ l is set for all interactions l.
For instance, for monodisperse particles with constant diameter σ, as for the standard Lennard-Jones (LJ) potential [9], the scaling variable becomes s = r/σ and the reduced cutoff s c = r c /σ. The effect of introducing s c is to replace u(s) by the truncated potential with H(s) being the Heaviside function [48]. Even if Eq. (57) is taken by definition as the new system Hamiltonian, it is well known that impulsive corrections at the cutoff have to be taken into account in general for the excess pressure P ex and other moments of the first derivatives of the potential [6]. This is seen by considering the derivative of the truncated potential Shifting of truncated potential. These corrections with respect to first derivatives can be avoided by considering a properly shifted potential [6]  Correction to the Born coefficient. As shown in Ref. [39], the standard shifting of a truncated potential is, however, insufficient in general for properties involving second (and higher) derivatives of the potential. This is particulary the case for the Born coefficient µ B , Eq. (47), which is required to compute the affine shear elasticity µ A , Eq. (4). The second derivative of the truncated and shifted potential being this implies that the Born coefficient reads µ B =μ B − ∆µ B withμ B being the bare Born coefficient and ∆µ B the impulsive correction which needs to be taken into account. The latter correction may be readily computed numerically from the configuration ensemble using [39] ∆µ being a weighted radial pair distribution function which is related to the standard radial pair distribution function g(r) [51].
Mixtures and polydisperse systems. Below we shall consider model systems for mixtures and polydisperse systems where σ l may differ for each interaction l. For such mixed potentials u(s), u t (s) and u s (s) and their derivatives take in principal an explicit index l, i.e. one should write u l (s), u t,l (s), u s,l (s) and so on. (This is often not stated explicitly to keep a concise notation.) For example one might wish to consider the generalization of the monodisperse LJ potential, Eq. (56), to a mixture or polydisperse system with where ǫ l and σ l are fixed for each interaction l. In practice, each particle i may be characterized by an energy scale E i and a "diameter" D i . The interaction parameters ǫ l (E i , E j ) and σ l (D i , D j ) are then given in terms of specified functions of these properties [14]. The extensively studied Kob-Andersen (KA) model for binary mixtures of beads of type A and B [52], is a particular case of Eq. (62) with fixed interaction ranges σ AA , σ BB and σ AB and energy parameters ǫ AA , ǫ BB and ǫ AB characterizing, respectively, AA-, BB-and AB-contacts. The expression Eq. (61) remains valid for such explicitly l-dependent potentials.
Shear modulus and compression modulus. Since withG τ τ being the uncorrected bare shear modulus. As we shall see in Sec. IV A, the correction ∆µ B is of importance for the precise determination of G τ τ close to the glass transition where the shear modulus must vanish.
An impulsive correction has also to be taken into account for the compression modulus K = K pp computed from the normal pressure fluctuations in a constant volume ensemble as discussed in Appendix A. Using the symmetry of isotropic systems, Eq. (A24), one verifies [39] withK pp being the uncorrected compression modulus.

III. ALGORITHMIC AND TECHNICAL ISSUES
Introduction. The computational data presented in this work have been obtained using two extremely well studied models of colloidal glass-formers which are described in detail elsewhere [13,14,18,52]. We present first both models and the molecular dynamics (MD) and Monte Carlo (MC) [6,9,25] methods used to compute them in the different ensembles (NVγT, NVτ T, NPγT and NPτ T) compared in this study. We comment then on the quench protocol and locate the glass transition temperature T g for both models via dilatometry, as shown in Fig. 2. This is needed as a prerequisite to our study of the elastic behavior in Sec. IV. As sketched in Fig. 1, periodic boundary conditions are applied in all spatial directions with L x = L y = . . . for the linear dimensions of the d-dimensional simulation box. Time scales are measured in units of the LJ time τ LJ = (mσ 2 /ǫ) 1/2 [9] for our MD simulations (m being the particle mass, σ the reference length scale and ǫ the LJ energy scale) and in units of Monte Carlo Steps (MCS) for our MC simulations. LJ units are used throughout this work (ǫ = σ = m = 1) and Boltzmann's constant k B is also set to unity.
Kob-Andersen model. The Kob-Andersen (KA) model [52] for binary mixtures of LJ beads in d = 3 dimensions has been investigated by means of MD simulation [6] taking advantage of the public domain LAMMPS implementation [33]. We use N = n A + n B = 6912 beads per simulation box and molar fractions c a = n A /n = 0.8 and c b = n B /n = 0.2 for both types of beads A and B. Following Ref. [52] we set σ AA = 1.0σ, σ BB = 0.88σ and σ AB = 0.8σ for the interaction range and ǫ AA = 1.0ǫ, ǫ BB = 0.5ǫ and ǫ AB = 1.5ǫ for the LJ energy scales. Only data for the usual (reduced) cutoff s c = r c /σ l = 2.5 are presented. For the NPγT-ensemble we use the Nosé-Hoover barostat provided by the LAMMPS code ("fix npt command") [33] which is used to impose a constant pressure P = 1 for all temperatures. Starting with an equilibrated configuration at T g = 0.7 well above the glass transition temperature T g = 0.41 (as confirmed in the inset of Fig. 2), the system is slowly quenched in the NPγT-ensemble. The imposed mean temperature T (t) varies linearly as T (t) = T i − Rt with R = 5 × 10 −5 being the constant quench rate. After a tempering time t eq = 10 4 we compute over a sampling time t = 5 · 10 4 various properties for the isobaric system at fixed temperature such as the moduli G τ τ represented by the filled spheres in Fig. 10. Fixing only then the volume (NVγT-ensemble) the system is again tempered (t eq = 10 4 ) and various properties as the ones recorded in Table I are computed using a sampling time t = 5 · 10 4 . Note that in the NVγT-ensemble the equations of motion are integrated using a velocity Verlet algorithm [9] with a time step δt = 0.01 and systems are kept at the imposed temperature T using a Langevin thermostat [9] with friction coefficient γ = 0.5. By redoing for several temperatures around and below T g the tempering and the production runs we have checked that the presented results remain unchanged (within numerical accuracy) and that ageing effects are absolutely negligible. Averages are taken over four independent configurations. As one may expect from the discussion in Sec. II C, very similar results have been found for both ensembles for simple averages, for instance for the affine shear elasticity µ A , and for the shear modulus G τ τ , Eq. (3). Unfortunately, we are yet unable to present data for the KA model data at imposed shear stress τ (NVτ T-, NPτ T-ensembles).
Polydisperse LJ beads. A systematic comparison of constant-γ and constant-τ ensembles has been performed, however, by MC simulation of a specific case of the generalized LJ potential, Eq. (62), where all interaction energies are identical, ǫ l = ǫ, and the interaction range is set by the Lorentz rule σ l = (D i + D j )/2 [51] with D i and D j being the diameters of the interacting particles. Only the strictly two-dimensional (2D) case  I: Some properties for the KA model at normal pressure P = 1 and shear stress τ = 0: the mean density ρ, the mean energy per volume eρ, the impulsive correction term ∆µB (Sec. II D), the (corrected) Born-Lamé coefficient µB, the affine shear elasticity µA = µB − Pex, the shear modulus Gττ obtained in the NVγT-ensemble, the hypervirial ηB and the affine dilatation elasticity ηA = 2P id + ηB + Pex discussed in Appendix A and the compression modulus Kpp obtained for the NVγT-ensemble. The affine elasticities µA and ηA are the natural scales for, respectively, the shear modulus G and the compression modulus K.
(d = 2) is considered. Following Ref. [14], the bead diameters of this polydisperse LJ (pLJ) model are uniformly distributed between 0.8σ and 1.2σ. For the examples reported in Sec. IV we have used N = 10000 beads per box. We only present data for a reduced cutoff s c = 2.0s 0 with s 0 = 2 1/6 being the the minimum of the polydisperse LJ potential, Eq. (62). Various properties obtained for the pLJ model are summarized in Table II.
Local MC moves for the pLJ model. For all ensembles considered here local MC moves are used (albeit not exclusively) where a particle is chosen randomly and a displacement δr, uniformly distributed over a disk of radius δr max , from the current position r of the particle is attempted. The corresponding energy change δE is calculated and the move is accepted using the standard Metropolis criterion [25]. The maxium displacement distance δr max is chosen such that the acceptance rate A remains reasonable, i.e. 0.1 ≤ A ≤ 0.5 [25]. As may be seen from the inset in Fig. 3 where the acceptance rate A(T ) is shown as a function of temperature T , we find that temperatures below T = 0.1 are best sampled using δr max ≈ 0.01, whereas δr max ≈ 0.1 is a good choice for the interesting temperature regime between T = 0.2 and T = 0.3 around the glass transition temperature T g = 0.26. (See the main panel of Fig. 2 for the determination of T g .) Various data sets are presented in Sec. IV for a fixed value δr max = 0.1 which allows a reasonable comparison of the sampling time dependence for temperatures around the glass transition. This value is not necessarily the best choice for tempering the system and for computing static equilibrium properties at the given temperature, especially below T ≈ 0.2. Collective plane-wave MC moves. Local MC jumps in the NVγT-ensemble become obviously inefficient for relaxing large scale structural properties [25] for the lowest temperatures computed for the pLJ model (T ≤ 0.1) as may be seen, e.g., from the finite, albeit very small shear stresses τ which happen to occur naturally in some of the quenched configurations (especially if the quench rate R is too large) and can only with difficulty be relaxed using local jumps alone [53]. We have thus crosschecked and improved the values obtained using only local MC moves by adding global MC move attempts for all particles with longitudinal and (more importantly) transverse plane waves commensurate with the simulation box. The random amplitudes and phases of the waves are assumed to be uniformly distributed. Note that the maximum amplitude a max (q) of each wave is chosen inversely proportional to the length q = |q| of the wavevector in agreement with continuum theory. This allows to keep the same acceptance rate for each global mode if q ≪ 1. As one expects, the plane-wave displacement attempts becomes inappropriate for large q where continuum elasticity breaks down [13,14,18] and, hence, the acceptance rate A(q) gets too large to be efficient [54]. Note that the maximum amplitude a max (q) has to be chosen much smaller for the longitudinal waves as for the transverse ones, i.e. since the systems are essentially incompressible, the transverse modes naturally are more effective for quenching and equilibrating the systems. These global plane-wave moves have been added for the tempering of all presented configurations below T = 0.3 with t eq = 10 7 and for production runs over sampling times t = 10 7 for static properties for T < 0.2. Details concerning the collective plane-wave MC moves must be given elsewhere. MC sampling of different ensembles. In this paper we shall be concerned more with the consequences of affine displacements associated with collective MC moves changing the overall box volume V and shear strain γ. Using again a Metropolis criterion for accepting a suggested change of V or γ, these pure dilatational and pure shear strain fluctuations are used to impose a mean normal pressure P = 2 and a mean shear stress τ = 0. For a shear strain altering MC move one first chooses randomly a strain fluctuation δγ with |δγ| ≤ δγ max . Assuming an imposed shear stress τ the suggested δγ is accepted if with ξ being a uniformly distributed random variable with 0 ≤ ξ < 1. The energy change δE associated with the affine strain may be calculated using Eq. (35). Since τ = 0 is supposed in this study, the last term in the exponential drops out. (As noted in Sec. II B, the volume remains unchanged in a pure shear strain and, hence, no translational entropy contribution appears.) The maximum shear strain step δγ max is fixed such that at the given temperature the acceptance rate remains reasonable [55]. As described in more detail in Ref. [9], one similarly chooses a relative volume change δV /V with V being the current volume [55]. The volume altering move is accepted at an imposed pressure P if where δE may be computed using Eq. (A16). The logarithmic contribution corresponds to the change of the translational entropy. Using Eq. (65) after each MCS for local MC moves (and global plane-wave MC moves if these are added) realizes the NVτ T-ensemble shown in Fig. 1(c), using Eq. (66) the NPγT-ensemble shown in Fig. 1(d) and using in turn both strain fluctuations the NPτ T-ensemble in Fig. 1(e). We remind that NVγT-and NPγT-ensembles are used to determine G τ τ and NVτ Tand NPτ T-ensembles to obtain G γτ and G γγ . As shown in Fig. 2, we quench the configurations starting from T i = 0.5 in the NPγT-and the NPτ T-ensemble using a constant quench rate R. Note that the smallest quench rate R = 10 −8 for one configuration in the NPγT-shown did require alone a run over six months using one core of an Intel Xeon E5410 processor. Interestingly, similar results are obtained with the NPτ T-ensemble only using a rate R = 10 −7 . (Due to the additional computation for the shear strain fluctuations this is, however, only a factor five faster.) The configurations created by the NPγT-quenches are used (after tempering) for the sampling in the NVγT-and NPγT-ensembles, the configurations obtained by the NPτ T-quenches for the sampling in the NVτ T-and NPτ T-ensembles. Only two independent configurations have been sampled following the full protocol for each temperature and each of the four ensembles. Instead of increasing further the number of configurations (which will be done in the future) we have focused in the present preliminary study on long tempering times and production runs over t = 10 7 MCS which have been redone several times for a few temperatures (summing up to total runs of ≈ 5 · 10 7 MCS for T = 0.2 and T = 0.25) to check for equilibration problems and to verify that ageing effects can be ignored. We firmly believe that the breaking of the time translational invariance is not a relevant issue for the strong sampling time dependence reported below for the shear moduli. Glass transition temperature. The inset of Fig. 2 presents for the KA model how the glass transition temperature T g may be obtained using standard dilatometry where the rescaled density log(1/ρ) is traced as a function of temperature T as in Ref. [7] for a similar polymer glass problem. This reveals a linear low-T (solid line) and a linear high-T regime. The intercept of both lines confirms the well-known value T g ≈ 0.41 from the literature [52]. The main panel presents the unscaled number density ρ for pLJ systems at constant normal pressure P = 2 as a function of temperature T for different quench rates R in MCS. Only local MC moves with a fixed maximum particle displacement δr max = 0.01 for all temperatures T are used for the examples given. The crosses refer to a quench in the NPτ T-ensemble where the box shape is allowed to fluctuate at τ = 0. All other data refer to the NPγT-ensemble with γ = 0. The solid line and the dashed line are linear fits to the best density estimates  indicated by the stars obtained for, respectively, low and high temperatures. By matching both lines one determines T g ≈ 0.26. A similar value is given if log(1/ρ) is plotted as a function of T . The main panel of Fig. 3 shows the interaction energy per particle e of the pLJ model as a function of temperature for the same quench protocols as in Fig. 2. The solid line and the dashed line are linear fits to the energy for, respectively, low and high temperatures. This confirms T g ≈ 0.26 for the pLJ model [56].

A. High-temperature liquid limit
We begin our discussion by focusing on the high temperature liquid limit where all reasonable operational definitions of the shear modulus G must vanish. The results presented in Fig. 4 have been obtained for the pLJ model at mean pressure P = 2, mean shear stress τ = 0 and mean temperature T = 1, i.e. far above the glass transition temperature T g ≈ 0.26. Only local MC monomer displacements with a maximum jump displacement δr max = 0.1 are reported here. Instantaneous properties relevant for the moments are written down every 10 MCS and averaged using standard gliding averages [9], i.e. we compute mean values and fluctuations for a given time interval [t 0 , t 1 = t 0 + t] and average over all possible intervals of length t. The horizontal axis indicates the latter interval length t in units of MCS. The vertical axis is made dimensionless by rescaling the moduli with the (corrected) affine shear elasticity µ A = µ B − P ex ≈ 17 from Table II. (Being a simple average, µ A is found to be identical for all computed ensembles.) In agreement with Eq. (23), the ratio G(t)/µ A for all operational definitions drops rapidly below unity. The dashed line indicates the asymptotic power-law slope −1 [21,39]. This can be understood by noting that at imposed τ the shear strain freely diffuses in the liquid limit, i.e. δγ 2 ∼ t, which implies that G ≈ G γτ ∼ 1/t. For both NVτ T-and NPτ T-ensembles we compare our key definition G γτ , Eq. (6), to the observable G γγ , Eq. (9). In both cases we confirm G γτ (t) ≈ G γγ (t) for larger times as suggested by the thermodynamic relation Eq. (8) which has been directly checked [38].
We have also computed the dimensionless correlation coefficient c γτ (t), Eq. (7), which characterize the correlation of the measured instantaneous shear strainsγ and shear stressesτ . As indicated for the NPτ T-ensemble (stars), c γτ vanishes rapidly with t, i.e.γ andτ are decorrelated. The dash-dotted line indicates the power-law exponent −1/2 expected from Eq. (22) and G γτ (t) ∼ 1/t. Note that by plotting c γτ (t) vs. G γτ (t)/µ A , Eq. (22) has been directly verified for both NVτ T-and NPτ Tensembles (not shown).
The stress fluctuation formula, Eq. (3), is represented in Fig. 4 by squares and diamonds for, respectively, NVγT-and NPγT-ensembles. The bare shear modulus G τ τ without the impulsive correction ∆µ B for the Born coefficient µ B discussed in Sec. II D is given by small filled symbols. The impulsive correction ∆µ B ≈ 0.2 (Table II) computed independently from the weighted radial pair correlation function, Eq. (61), is shown by a horizontal line. The corrected moduli G τ τ =G τ τ − ∆µ B vanish indeed as expected. See Ref. [39] for a systematic numerical characterization of ∆µ B as a function of the reduced cutoff s c . As may also be seen there, a truncation correction is also necessary for the KA model. The ∆µ B for different temperatures may be found in Table I for the KA  model and in Table II for the pLJ model. We assume from now on that the impulsive corrections are properly taken into account without stating this technical point explicitly.

B. Dynamical matrix harmonic network
Introduction. We turn now to the opposite low-T limit focusing specifically on the pLJ model at T = 0.001. Obviously, such deeply quenched colloidal glasses must behave as amorphous solids with a constant shear modulus for large times (albeit not infinite times) as sketched in Fig. 1(b) by the top curve. Before presenting in Sec. IV C the elastic properties of these glasses, we discuss first conceptually simpler substitute systems formed by permanent spring networks. We do this for illustration purposes since questions related to ergodicity and ageing are by construction irrelevant and the thermodynamical relations reminded in Sec. II should hold rigorously.
Network Hamiltonian. Assuming ideal harmonic springs the interaction energy reads with K l being the spring constant and R l the reference length of spring l. The sum runs over all springs l between topologically connected vertices i and j of the network. We assume here a strongly and homogeneously connected network where every vertex i is in contact with many neighbors j. Such a network may be constructed from a pLJ configuration at T = 0.001 keeping the particle positions for the positions of the vertices and replacing each LJ interaction l by a spring of spring constant K l and reference length R l . The simplest choice we have investigated is to set K l ≡ 1 and R l ≡ r l which corresponds to a well-defined ground state of energy E = 0 and pressure P = 0 at T = 0 (not shown). We discuss here instead a slightly more realistic choice of K l and R l which corresponds to the same "dynamical matrix" as the pLJ configuration at T → 0 and P = 2, i.e. the same second derivative of the total interaction potential with respect to the particle positions [13,14]. This choice imposes the setting with u(r) = u(s) being the interaction potential for reduced distances s l = r l /σ l ≤ s c . (Impulsive corrections at s l = s c are neglected here.) Sampling. Taking advantage of the fact that all interacting vertices are known, these networks can be computed using global MC moves with longitudinal and transverse planar waves as described in Sec. III. No local MC jumps for individual vertices have been added. Note that although some K l and R l may even be negative, this does not affect the global or local stability of the network. Presumably due to the fact that our spring constants are rather large, we have not observed either any buckling instability if volume fluctuations are allowed. For networks at the same pressure P , shear stress τ and temperature T as the original pLJ configuration, we obtain, not surprisingly, the same affine shear elasticity µ A = 33.9 as for the reference. The same applies for other simple averages.
Shear modulus. Figure 5 presents various measurements of the shear modulus G for different ensembles as a function of sampling time t. As in Sec. IV A, we use gliding averages and the vertical axis has been rescaled by the affine shear elasticity µ A . As can be seen, we obtain in the long-time limit G/µ A ≈ 0.67 (bold line) irrespective of whether the shear modulus is determined from G γτ or G γγ in the NVτ T-or NPτ T-ensembles or using G τ τ in the NVγT-or NPγT-ensembles at fixed shear strain γ. Note that even in these cases the shear modulus decreases first with time albeit the network is perfectly at thermodynamic equilibrium and no ageing occurs by construction. As shown by the dash-dotted line, G(t)/µ A → 1 for short times, i.e. in this limit the response is affine. As indicated by the stars, we have also computed the correlation coefficient c γτ (t) as a function of time. From c γτ = 1 for small t, it decreases somewhat becoming c γτ ≈ √ 0.67 ≈ 0.82 for large times in agreement with Eq. (22). If G τ τ is computed in the "wrong" NVτ T-or NPτ T-ensembles, it is seen to vanish rapidly with time as shown by the filled squares. This decay is of course expected from Eq. (21) as discussed in Sec. II B and Sec. II C.
Shear-stress fluctuations. To emphasize the latter point we show in Fig. 6 the stress fluctuations µ F (t) for different ensembles. This is also done to have a closer look at the time dependence of various contributions to G τ τ (t). As before, the vertical axis has been rescaled with µ A ≡ 33.9. We also indicate the rescaled values of µ B (t) (thin top line) and µ A (t) = µ B (t) − P ex (t) (dashdotted line) computed by a gliding average over time intervals [t 0 , t 1 = t 0 + t]. Being simple averages both properties do not depend on the ensemble probed (not shown) and become virtually immediately t-independent as can be seen from the indicated horizontal lines. This is quite different for the different shear stress fluctuations µ F presented which increase monotonously from µ F ≈ 0 for short times to their plateau value for long times. In agreement with Eq. (21) we find that µ F (t)| τ → µ A for NVτ T-and NPτ T-ensembles and as predicted by  Eq. (20) we confirm µ F (t)| γ → µ A − G for NVγT-and NPγT-ensembles. The asymptotic limit is reached after about 10 4 MCS in the first case for imposed τ , but only after 10 7 for imposed γ. This is due to the fact that µ F | τ corresponds to the fluctuations of the self-contribution of individual springs as discussed in Sec. II C, while µ F | γ also contains contributions from stress correlations of different springs. It should be rewarding to analyze in the future finite system-size effects and the time dependence of three-dimensional (3D) systems using the dynamical matrix associated to the KA model in the T → 0 limit.

C. Low-temperature glass limit
We return now to the colloidal glasses where the topology of the interactions is, of course, not permanently fixed (at least not at finite temperatures) and the shear stresses may thus fluctuate more strongly. As an example, we present in Fig. 7 the shear modulus G(t) as a function of sampling time t for pLJ beads at mean pressure P = 2, mean shear stress τ = 0 and mean temperature T = 0.001, i.e. the same conditions as in the previous subsection. All simple averages, such as the affine shear elasticity µ A = 33.9, are the same for all ensembles albeit different quench protocols have been used. For all examples given we use local MC jumps with δr max = 0.01 which allows to keep a small, but still reasonable acceptance rate A ≈ 0.161. We have also performed runs with additional global MC moves using longitudinal and transverse plane waves. This yields similar data which is shifted horizontally to the left (not shown).
As one would expect, the shear moduli presented in Fig. 7 are similar to the ones shown in Fig. 5 for the permanent spring network. The filled triangles show G τ τ computed in the "wrong" NPτ T-ensemble. As expected from Eq. (20) these values vanish rapidly with time. If on the other hand G γτ , G γγ and G τ τ are computed in their natural ensembles, all these measures of G are similar and approach, as one expects for such a low temperature, the same finite asymptotic value G/µ A ≈ 0.46 indicated by the bold horizontal line. The stars indicate the squared correlation coefficient c 2 γτ (t) which is seen to collapse perfectly on the rescaled modulus G γτ (t)/µ A which shows that Eq. (22) even holds for small times before the thermodynamic plateau has been reached. Similar reduced shear moduli G τ τ /µ A have been also obtained for the KA model in the low-T limit as may be seen from Table I. These results confirm that about half of the affine shear elasticity is released by the shear stress fluctuations in agreement with Refs. [13,14]. Interestingly, this value is slightly smaller as the corresponding value for the permanent dynamical matrix model presented above. Since some, albeit very small, rearrangement of the interactions must occur for the pLJ particles at finite temperature, this is a reasonable finding.
The rescaled shear stress fluctuations µ F (t)/µ A for the same ensembles are presented in Fig. 8. While the two-point correlations µ B (t)/µ A (thin horizontal line) and µ A (t)/µ A (dash-dotted horizontal line) do not dependent on time, one observes again that µ F (t)/µ A increases monotonously from zero to its long-time asymptotic plateau value. In agreement with Eq. (21) we find µ F (t)| τ → µ A for the NVτ T and the NPτ T ensembles, while µ F (t)| γ → µ A − G (dashed horizontal line) for the NVγT and the NPγT ensembles. (Similar behavior has been obtained for the KA model in the NVγT-ensemble for temperatures T ≤ 0.1.) As may be better seen for the data set obtained for T = 0.010 using the NVγTensemble (small filled squares) where µ A ≈ 33.7, the stress fluctuations µ F (t)| γ do not rigorously become con- stant, but increase extremely slowly on the logarithmic time scales presented. There is, hence, even at low temperatures always some sampling time dependence if the interactions are not permanently fixed as for the topologically fixed networks discussed above. However, this effect becomes only sizeable above T = 0.2 for the KA model and above T = 0.1 for the pLJ model. We shall now attempt to characterize this temperature dependence.

D. Scaling with temperature
Stress fluctuations. We turn now to the characterization of the temperature dependence of the shear modulus G and various related properties. We focus first on the best values obtained for the longest simulation runs available. Figure 9 presents the stress fluctuations obtained for the KA model using the NVγT-ensemble (open spheres) and for the pLJ model using NVγT-and the NPτ T-ensembles. The affine shear elasticity µ A obtained for each system is represented by small filled symbols. As shown for both models by a bold solid line for the low-T and by a dashed line for the high-T regime, µ A (T ) decays roughly linearly with T . Interestingly, the slopes match precisely at T = T g for both models. As may be seen from the large triangles for the NPτ T-ensemble, we confirm for all temperatures that µ F | τ ≈ µ A as expected from Eq. (21). For systems without box shape fluctuations, it is seen that µ F | γ becomes non-monotonous with a clear maximum at T ≈ T g . In the liquid regime for T > T g where the boundary conditions do not matter, we obtain µ F | τ ≈ µ F | γ ≈ µ A , i.e. all shear stress fluctuations decay with temperature. Below T g the boundary constraint (γ = 0) becomes more and more relevant with decreasing T which reduces the shear stress fluctuations. According to Eq. (20), µ F | γ must thus decrease with increasing G as the system is further cooled down. That µ F | γ must necessarily be non-monotonous with a maximum at T g is one of the central results of the presented work. Shear modulus as a function of temperature. The temperature dependence of the (unscaled) shear moduli G(T ) is presented in Fig. 10. Data obtained by MD simulation of the KA model (P ≈ 1, T g ≈ 0.41) are indicated for the NVγT-ensemble (open spheres) and the NPγT-ensemble (filled spheres). All other results refer to the best values obtained by MC simulations of the pLJ model (P = 2, T g ≈ 0.26) using both local and global moves and different ensembles as indicated. For the pLJ model the values G τ τ obtained in the NVγT ensemble (squares) are seen to be essentially identical for all T to the moduli G γτ (crosses) and G γγ (triangles) obtained in ensembles with strain fluctuations. As shown by the small filled triangles, G τ τ vanishes in the NPτ T ensemble for all T in agreement with Eq. (21). (We remind that the impulsive correction ∆µ B has to be taken into account if G τ τ is used.) Decreasing the temperature further below T g the shear moduli are seen to increase rapidly for both models. As shown by the solid and the dashed lines both models are well fitted by a cusp-singularity with fit constants g 1 ≈ 23 for the KA model and g 1 ≈ 15 for the pLJ model where we have set in both cases g 0 = 0. Confirming the MD simulations by Barrat et al. [11], this suggests that the transition is very sharp, albeit continuous in qualitative agreement with the predictions from Ref. [28,29]. If we fit the data close to T g with an additional off-set g 0 , a slightly negative value is obtained which is not compatible with MCT [31,32]. However, admittedly both the number of data points close to T g and their precision (due to the small number of independent configurations) are yet not sufficient to rule out completely a positive, albeit (presumably) small off-set.
A slightly different representation of the data is given can only be seen below T ≈ 0.2 and only below T = 0.01 one observes a t-independent plateau over two orders of magnitude. We emphasize that we present here a sampling time effect and not a time correlation function. The same sampling time effect is observed after additional tempering.
in Fig. 11 where the reduced shear modulus y = G(T )/µ A (T ) is plotted as a function of the reduced temperature x = T /T g for both models. Using the independently determined glass transition temperature T g and affine shear elasticity µ A (T ) as scales to make the axes dimensionless, the rescaled data are found to collapse ! This is of course a remarkable and rather unexpected result considering that two different models in two different dimensions have been compared. The bold line indicates the cusp-singularity y = 0.45(1 − x) 1/2 with a prefactor which is compatible to the ones used in Eq. (69) considering the typical values of µ A (T ) in both models. Whether this striking collapse is just due to some lucky coincidence or makes manifest a more general universal scaling is impossible to decide at present. (Please note that the corresponding data for the compression modulus shown in Fig. 15 does not scale.) No rescaling of the vertical axis is necessary for the squared correlation coefficient c 2 γτ shown in the inset for two different ensembles of the pLJ model. The same continuous cusp-like decay (bold line) is found as for the modulus G. This is again consistent with the thermodynamic reasoning put forward at the end of Sec. II A, Eq. (22). Sampling time dependence. A word of caution must be placed here. The results presented in Figs. 9-11 correspond to the longest simulation runs we have at present been able to perform. Obviously, this does not mean that they correspond to the limit for asymptotically long sampling times, if the latter limit exists, or at least an intermediate plateau of respectable period of time. Focusing on the pLJ model in the NPτ T-ensemble we attempt in Fig. 12 to give a tentative characterization of the sampling time dependence. We plot G γτ (t) as a function of t for different temperatures as indicated. To make the data comparable all simulations for T ≥ 0.1 have been performed with local MC jumps of same maximum distance δr max = 0.1. As we have already seen in Fig. 4, the shear modulus decays inversely with time in the liquid regime at high temperatures. This limit is indicated by the thin line. With decreasing T the (unscaled) shear modulus increases and a shoulder develops. Note that even for T = 0.25, i.e. slightly below the glass transition temperature T g = 0.26 estimated via dilatometry (Fig. 2), G γτ (t) decays strongly with time. A more or less t-independent shoulder (on the logarithmic scales used) can only be seen below T ≈ 0.2. We emphasize that we present here a sampling time effect expressing the fact that the configuration space is more and more explored with increasing t and not a time correlation function [9]. Please note also that this time dependence has nothing to do with equilibration problems or ageing effects. Time translational invariance is perfectly obeyed in our simulations as we have checked by rerunning the sampling after additional tempering over at least 10 7 MCS for all temperatures. Similar t-dependencies for different temperatures T have also been observed for G γτ (t) and G τ τ (t) for the pLJ model and for G τ τ (t) for the KA model (not shown).
A different representation of the data for both models is shown in Fig. 13 where the reduced shear modulus y = G(T )/µ A (T ) is plotted against the reduced temperature x = T /T g for different sampling times t as indicated. Similar behavior is found for both models. The shear modulus decreases with t and this the more the larger the temperature T . For smaller temperatures the different data sets appear to approach more readily as one expects from the t-independent shoulder for the pLJ model in Fig. 12. As a consequence the glass transition is seen to become sharper with increasing t. The data approach the continuous cusp-singularity indicated by the bold line. However, it is not clear from the current data if the data converge to an at least intermediately stable behavior. Longer simulation runs are clearly necessary for both models to clarify this issue. Similar data have also been obtained for the pLJ model from the independent determination of G γτ (t) and G γγ (t) in the NVτ T-ensemble and for G τ τ (t) in the NVγT-and NPγT-ensemble.

A. Summary
The stress fluctuation formalism is a powerful method for computing elastic moduli for canonical ensembles using simulation boxes of constant volume and shape [4][5][6][7]. Focusing on the shear modulus G of two wellknown glass-forming colloidal model systems in d = 3 and d = 2 dimensions, which we have sampled by means of MD and MC simulations, we have addressed the general question of whether the stress fluctuation method may be used through a solid-liquid transition where a well-defined displacement field ceases to be defined. Deliberately setting the experimentally motivated linear regression Eq. (6) as the fundamental definition, the shear modulus G has been computed from the strain and stress fluctuations (assuming a fixed finite measurement time window t) in ensembles where either a shear strain γ = 0 (NVγT-and NPγT-ensembles) or a (mean) shear stress τ = 0 (NVτ T-and NPτ T-ensembles) have been imposed (Fig. 1). Working at constant mean normal pressure we have computed and compared the temperature dependence for various simple averages and fluctuations con-tributing to the shear modulus in different ensembles. For the KA model only the currently available results for NVγT-and NPγT-ensembles have been presented. As has been stressed in Sec. II B and Sec. II C, the stress fluctuation representation G τ τ at fixed shear strain does not rely in principle on a solid-like reference state for a displacement field of tagged particles if thermodynamics can be assumed to hold. It thus holds formally through the glass-transition temperature T g up to the liquid regime (albeit with a trivial value G = 0) as does the better known relation for the compression modulus K (Appendix A) [8]. As emphasized in Sec. II D, impulsive corrections due to the truncation of the pair potentials have to be properly taken into account, especially at high temperatures (Fig. 4), for the precise determination of the affine shear elasticity µ A = µ B − P ex [39]. Confirming Eq. (10), it has been shown for the pLJ model that G τ τ at γ = 0 has the same long-time behavior as the conceptually more direct observables G γτ and G γγ obtained using the strain fluctuations at constant mean shear stress τ = 0. The standard thermodynamic relations comparing different ensembles appear thus to hold also in practice (at least within our numerical precision) through the glass transition for all T and this albeit 1. some degrees of freedom get quenched at low temperatures on the time window t computationally accessible, 2. all operational definitions of G are transient in the sense that they vanish for t → ∞ (Figs. 12 and 13) and 3. it is not self-evident that in the latter limit and for large temperatures γ and τ can still be treated as a pair of conjugated thermostatistical variables.
The latter point compactly epitomized by the thermodynamic relation Eq. (8) has been shown to hold, however, with remarkable precision up to very high temperatures for the pLJ model. As predicted by general Legendre transformation, G τ τ vanishes for all temperatures T , if τ rather than γ is imposed (Figs. 5, 7 and 10). The shear-stress fluctuations µ F | τ are thus given by the affine response under an external load µ A (Fig. 9), i.e. µ F | τ reduces to a simple two-point pair correlation function if pair potentials are considered as discussed in Sec. II C. The same holds above T g for constant γ, since the boundary conditions are irrelevant for the liquid state. This symmetry with respect to the boundary conditions (τ ↔ γ) is broken at the glass transition for a large, but finite t: at constant γ the stress fluctuations reveal a strong non-monotonous behavior with a clear maximum at T g (Fig. 9). The shear modulus is the order parameter characterizing this symmetry breaking. Alternatively, as we have discussed in Sec. II C, G τ τ /µ A ≡ 1 − µ F | γ /µ A may be seen as an order parameter comparing the ratio of the non-affine to the affine shear responses. Since µ A > µ F | γ for T < T g , this implies that the stress fluctuations contain higherorder correlations reducing the free energy of the system.
The increase of G below T g is reasonably fitted for both models by a continuous cusp singularity, Eq. (69), in qualitative agreement with recent replica calculations [28,29]. A jump discontinuity, as suggested by modecoupling theory [31,32] and another replica theory [30], is not compatible with the currently available data shown in Fig. 10 and Fig. 11. However, as shown in Fig. 12 and Fig. 13, our data depend strongly on sampling time t and the computation of larger t could lead to a sharper transition. Thus a jump discontinuity cannot be ruled out completely. (We are not aware of any other investigation of the sampling time effect for the shear modulus.) At present we believe, however, that it would be surprising if our data could be reconciled with a discontinuity at T g . It seems more likely that below T g , where the choice of the boundary conditions does matter as shown, the definition of the glassy shear modulus used in Refs. [30][31][32] do not correspond to the key operational definition G γτ of the shear modulus considered by us. In any case, our work shows that any theory and numerical scheme used to determine the shear modulus using a stress or strain fluctuation relation (in the low-wavevector limit or at a small finite q) must cleary specify and take into account whether the shear stress or the shear strain are macroscopically imposed. Otherwise, such an approach is void.

B. Outlook
Beyond the presented work. Future work should clearly focus on the more detailed description of the sampling time dependence. By extrapolating appropriately for the large-t asymptotic behavior, this may allow to settle the theoretical debate. As emphasized above, one shortcoming of the present work is that the time-scale used in our MC simulations was slightly arbitrary, especially if additional non-local particle moves are used in the low-T limit and if box volume and shape changing strains are included. The comparison of t-dependent properties for different temperatures T becomes thus delicate. Our MC study for the 2D soft particles should thus in any case be recomputed with Langevin thermostat MD dynamics as used for the data of the KA model presented. The t-dependence of the latter model, which is a better glass former than the 2D pLJ model [56], has still to be worked out for the key operational definition G γτ in the NVτ T-or NPτ T-ensembles [57]. We plan also in the near future to reconsider more carefully our previous investigation of the glass transition of polymer melts [7] taking properly into account the impulsive truncation corrections and comparing the shear moduli G γτ ≈ G γγ (NPτ T-ensemble) and G τ τ (NPτ T-ensemble) around T g .
In addition, it should be rewarding to compare our shear moduli with the values obtained from the displacement of particles following the procedure chosen in Ref. [31]. This procedure relies on the idea that the particles fluctuate around a well-defined reference position. In the low-T limit where the interaction network can be replaced by the dynamical matrix, this should yield the same results. However, for larger temperatures where the harmonic approximation must break down, this approach is questionable (both for the analysis of experimental and computational data) and should be compared to the moduli obtained directly from the strain and stress fluctuations of the overall simulation box. A simple scalar model. We are currently working on an extremely simplified scalar model for permanent and transient networks in d = 2 dimensions which may help to clarify the theoretical debate [58]. As sketched in the panel (a) of Fig. 14, in this model beads (filled spheres) are glued on rigid rods which are only allowed to move horizontally. In its most simple implementation this is done by means of global MC moves using transverse plane waves with a wavevector pointing in the y-direction. The glued beads are connected by ideal harmonic springs where the spring constants K l and reference lengths R l are chosen randomly according to fixed, narrow or broad distributions. The stiffness of the contacts between rods may differ (even including negative K l and R l as for the dynamical matrix studied in Sec. IV B). The shear modulus may be computed either from the shear strain fluctuations in the NVτ T-or the shear stress fluctuations in the NVγT-ensemble. The shear stresses can be relaxed by either allowing the particles to move along the rods using longitudinal plane waves in the x-direction and/or by breaking and reconnection of springs. The latter step is currently done by applying a chemical potential for the springs. Physically, the study of such transient networks is motivated by systems of hyperbranched polymer chains with sticky end-groups [35] or telechelic polymers in water-oil emulsions [36,37]. The associated relaxation frequency ω l ∝ exp(−E b /T ) for a given spring l is assumed to be set by an activation barrier E b which may differ for each spring. Depending on the distribution of this activation barrier the shear modulus G(t) thus decays more or less rapidly with sampling time t as sketched in panel (b) of Fig. 14. The shear modulus becomes constant, if the random percolating network is permanent (ω l = 0, E b /T → ∞) or the relaxation negligible (top solid line). The system behaves as a liquid (bottom line) for large (mean) ω l , while a shoulder is seen in the intermediate frequency range. While this is confirmed by preliminary simulations, we do at present not understand the role of the quenched noise injected in the model which is seen to strongly influence the shear stress fluctuations µ F (especially if fluffy layers happen to appear). Also the understanding of finite system-size effects is of course crucial for such a low-dimensional model. (Only systems containing N = 10 4 particles have been considered so far.) As sketched in panel (c) the key question is to describe accurately for asymptotically long simulation runs and asymptotically large box sizes the shape of G(T ) around the glass transition temperature T g ( E b ) where the shear modulus vanishes and µ F | γ should again have a maximum. the extensive variable is here the volume X = V = V and the conjugated intensive variable the negative pressure I = −P = − P whereV andP stand for the instantaneous values of the volume and the pressure. With these notations Eq. (A1) is obviously consistent with Eq. (12) or Eq. (14). In analogy with G γτ for the shear modulus G, the compression modulus K may be determined in an NPγT-or an NPτ T-ensemble using the linear regression relation where the index vp indicates that both the measured instantaneous volumeV and pressureP have been used. The corresponding correlation coefficient is Since from Eq. (15) we have δV δP = −k B T if P is imposed, this implies that K vp should be identical to the better known relation [40] which corresponds to the operational definition G γγ discussed above. The compression modulus K may of course also be determined using the pressure fluctuations in systems of imposed volume V as sketched in Fig. 1(c). As (to our knowledge) first derived by Rowlinson [8], the corresponding stress fluctuation formula reads being the fluctuation of the total normal pressure (corresponding to µ F ) and η A the "affine dilatational elasticity" (corresponding to µ A ). We shall give below a proper operational definition for η A . Here we only note that due to the general Legendre transformation Eq. (18) for intensive variables one knows immediately that This implies by comparison with Eq. (A5) that in analogy to Eq. (21) for the shear, i.e. η A can be directly be determined from the pressure fluctuations in NPγT-and NPτ T-ensembles. As before for the correlation coefficient c γτ , Eq. (22), it follows from Eq. (A7) that Since c vp ≤ 1, the affine dilatational elasticity sets an upper bound to the compression modulus K.

Rowlinson's stress fluctuation formula
Introduction. We reformulate here briefly the derivation of the stress fluctuation relation K pp for the compression modulus K of an isotropic system at imposed volume V by Rowlinson [8]. With V (0) being the volume of the unperturbed simulation box we assume a small relative volume change with all coordinates equally deformed as, e.g.
for the x-coordinate in d-dimensions. (The argument (0) denotes again the reference system.) Using a similar rescaling trick as for the shear modulus G in Sec. II B, the interaction energy U s (ǫ) of a configuration s of the strained system is expressed in terms of the coordinates (state) of the unperturbed system and the explicit metric parameter ǫ. Due to the volume change we have now also to take into account the ideal gas (kinetic) contributions. General excess contributions. We derive first the excess contribution P ex to the total pressure P = P id + P ex and the excess contribution K ex to the total compression modulus K = K id + K ex valid for an arbitrary conservative potential. The general relations Eqs. (26-29) stated above for an imposed shear strain γ still hold with the relative volume change ǫ replacing γ. Using Eq. (28) it follows for the excess pressure that defining the instantaneous excess pressure [42]. In the second step we have taken the limit ǫ → 0 and have dropped the argument (0). The average is again evaluated using the weights of the unperturbed system, Eq. (31). Using Eq. (27) and Eq. (29) one obtains for the compression modulus which is again understood to be taken at ǫ = 0. The first term η A,ex ≡ U ′′ s (ǫ) /V in Eq. (A12) corresponds to the change of the system energy assuming an affine strain transformation for all particle positions. It gives the excess contribution η A,ex to the total affine dilatational elasticity η A = η A,id + η A,ex . The second term corresponds to the excess contribution η F,ex to the total normal pressure fluctuation η F = η F,id + η F,ex mentioned in Sec. A 1. It corrects the overprediction of the affine strain contribution η A,ex [18]. Compression modulus. The main panel of Fig. 15 shows the best long-time values of the compression modulus K for both models. The vertical axis is rescaled by η A . The rescaled stress fluctuation formula, Eq. (A14), for the NVγT-ensemble is given for both models. For the pLJ model we also indicate K vp (crosses) and K vv (large triangles) for the NPτ T-ensemble. No rescaling with η A is necessary for the squared correlation coefficient c 2 vp (stars) also presented. A perfect collapse of all pLJ data is found confirming the equivalence of all operational definitions used. For both models we observe a strong increase of K/η A around x ≈ 1 with decreasing temperature. Interestingly, the rescaled data of both models is similar but not identical. The mentioned increase around x ≈ 1 is clearly more pronounced for the KA model. Note that if the unscaled K is directly plotted against T , this effect becomes even stronger (not shown) consistently with the scaling of η A presented in the inset. We note finally that for very low temperatures we observe for both models that the ratio K/η A = c 2 γτ approaches unity, i.e. the compression modulus is essentially dominated by the affine response. We remind that this is different for the reduced shear modulus for both models which is observed to approach G/µ A ≈ 1/2 in the low-T limit (Fig. 11), i.e. the shear stress fluctuations µ F | γ are more relevant than the pressure fluctuations η F | V in this limit.
Sampling time dependence. For the pLJ model we present in Fig. 16 the sampling time dependence of various operation definitions of K for the liquid limit at temperature T = 1.0. The horizontal axis indicates the sampling time t in units of MCS, the vertical axis is made dimensionless using η A as natural scale. As shown for K vp and K vv obtained in the NPτ T-ensemble and for K pp obtained using Eq. (A14) in the NVγT-ensemble, all operational definitions applied in their natural ensemble converge rapidly within t ≈ 10 3 MCS to the same plateau. (This is a general finding for all temperatures.) The plateau value for T = 1.0 is indicated by the bold solid line. Note that K/η A is smaller unity as expected from the discussion at the end of Sec. A 1.
The small filled symbols refer to K pp computed in the "wrong" NPτ T-ensemble. The small diamonds have been computed using directly Eq. (A14). As indicated by the dashed line, this yields K pp ≈ 2P id for large times. Only if we compute K pp = η A − η F using Eq. (A15), it does vanish properly as shown by the small triangles. This shows numerically that Eq. (A14), albeit perfectly fine for constant-V ensembles, is just not the correct formulation allowing to make manifest the general Legendre transformation Eq. (A6). This is, however, achieved using Eq. (A15). This finding holds rigorously for all temperatures as can be seen from the scaling of the small triangles in the inset of Fig. 15.  [54] In principle, it would be better to change the maximum amplitude of the plane waves adaptively in order to fix an acceptance rate A(q, T ) = const for all wavevectors and temperatures. This has yet not been done for the present study. As shown in Sec. IV B for the sampling of fixed networks constructed from the dynamical matrix corresponding to the colloidal glass, some systems have also been computed successfully only using global moves. It is thus not stictly necessary to include the computationally expensive local MC moves. In any case, the global planewave MC algorithm requires a more systematic theoretical and numerical analysis to improve its computational efficiency.
[55] For the temperature quench of the configurations sufficiently small values for δγmax for a shear strain altering MC move and for δVmax for a volume change have been chosen such that the acceptance rate of both MC steps remains reasonable at the lowest temperatures (T < 0.1).
As for the global plane-wave MC moves, it might better to change these maximum amplitudes adaptively while the quench is performed to keep fixed constant acceptance rates.
[56] Similar values for Tg have also been obtained for both models following the standard diagnostics described in the literature [11,52] of various dynamic properties such as the scaling of the mean-square displacements, the self-diffusion coefficient or the Fourier transformed density autocorrelation function in the low-wavevector limit. From the analysis of these dynamical properties it can be seen that compared to the KA model the pLJ model is a much more gradual glass-former qualitatively similar to the purely repulsive soft-sphere model studied in