Skip to main content
Log in

Ensemble fluctuations matter for variances of macroscopic variables

  • Regular Article - Flowing Matter
  • Published:
The European Physical Journal E Aims and scope Submit manuscript

Abstract

Extending recent work on stress fluctuations in complex fluids and amorphous solids we describe in general terms the ensemble average \(v(\Delta t)\) and the standard deviation \(\delta v(\Delta t)\) of the variance \(v[{\mathbf {x}}]\) of time series \({\mathbf {x}}\) of a stochastic process x(t) measured over a finite sampling time \(\Delta t\). Assuming a stationary, Gaussian and ergodic process, \(\delta v\) is given by a functional \(\delta v_{\mathrm {G}}[h]\) of the autocorrelation function h(t). \(\delta v(\Delta t)\) is shown to become large and similar to \(v(\Delta t)\) if \(\Delta t\) corresponds to a fast relaxation process. Albeit \(\delta v = \delta v_{\mathrm {G}}[h]\) does not hold in general for non-ergodic systems, the deviations for common systems with many microstates are merely finite-size corrections. Various issues are illustrated for shear-stress fluctuations in simple coarse-grained model systems.

Graphic abstract

This is a preview of subscription content, log in via an institution to check access.

Access this article

Price excludes VAT (USA)
Tax calculation will be finalised during checkout.

Instant access to the full article PDF.

Institutional subscriptions

Fig. 1
Fig. 2
Fig. 3
Fig. 4
Fig. 5
Fig. 6
Fig. 7
Fig. 8
Fig. 9
Fig. 10
Fig. 11
Fig. 12

Similar content being viewed by others

Notes

  1. We frequently switch between a discrete and a continuous representation \(i \leftrightarrow t, I \leftrightarrow \Delta t, h_i \leftrightarrow h(t), c_i \leftrightarrow c(t), \ldots \).

  2. The empirical variance is defined here without the usual “Bessel correction” [1]. Equation (1) is the formal definition of \(v[{\mathbf{x}}]\) which coincides with the genuine variance of x(t) only in the limit \(\Delta t \propto I \rightarrow \infty \).

  3. The longest relaxation time \(\tau \) of glass-forming liquids is generally called \(\tau _{\alpha }\) [16, 20, 21].

  4. The ensemble average \(\left\langle \ldots \right\rangle \) may be computed by taking the arithmetic average over \(N_{\mathrm{c}}\) independently prepared and sampled configurations c. For ergodic systems it is equivalent to sample over \(N_{\mathrm{k}}\) sub-intervals of length \(\Delta t\) of a very long trajectory of length \({\Delta t}_{\mathrm{max}} \gg \tau \).

  5. This may be better seen from the continuum representation of Eq. (15) rewritten as

    $$\begin{aligned} \int _0^t {\mathrm{d}} s \ R(s) =\frac{{\mathrm{d}}}{{\mathrm{d}} t} [(M(t)t^2/2]. \end{aligned}$$

    If the left-hand side converges to a constant \(\eta \) for \(t \rightarrow \infty \) this implies \(M(t)t^2 \rightarrow 2 \eta t\).

  6. Equation (14) holds if the perturbation is a “deformation”. In the case of an externally applied “force” it becomes \(R(t)=h(t)\) [18].

  7. The variances increase with s since the number of data used for the gliding average decreases linearly with s.

  8. According to Bochner’s theorem \({\hat{f}}(\omega ) \ge 0\) if and only if f(t) is a positive-definite function, i.e., all eigenvalues of the matrix \(g_{i,j}=f(t_i-t_j)\) are non-negative [21].

  9. For volume-averaged density fields \(\tau _{\mathrm{b}}\) must be to leading order system-size independent while \(\tau _{\mathrm{ne}}\) diverges in the macroscopic limit [30].

  10. The time series k may be obtained by first tempering the configuration c over a time interval \({\Delta t}_{\text {temp}}\gg \tau _{\text {b}}\) and by sampling then \(N_\text {k}\) time intervals \(\Delta t\) separated by constant spacer time intervals \(\Delta t_{{\mathrm{spac}}} \gg {\text {b}}\). \(N_\text {k}\) is assumed to be arbitrarily large and the k-averaged properties \(\delta v_{\text {int}}\) and \(\delta v_{\text {ext}}\) do thus neither depend on \(N_\text {k}\) nor the total sampling time \({\Delta t}_\text {max}= N_\text {k} (\Delta t+{\Delta t}_{\text {spac}})\). The latter point may become a delicate issue if the non-ergodicity constraint (\(\tau \rightarrow \infty \)) is not strictly obeyed.

  11. Ensemble averages within a meta-basin c, \(\left\langle \ldots \right\rangle _c\), compatible with the non-ergodicity constraint(\(\tau \rightarrow \infty \)) are considered: \(\left\langle \ldots \right\rangle _c\) is obtained for \(N_{{\mathrm{k}}} \rightarrow \infty \). Our task is to compute \(\Delta _{{\mathrm{ne}}}^{2} =\text {var}[v_{c}]\) for \(v_c = \left\langle x^2 \right\rangle _c -\left\langle x \right\rangle _c^2\). Substituting \(x = \sum _mx_{m}/N_{{\mathrm{m}}}\), setting \(v_{cm} = \left\langle x_m^2\right\rangle _c - \left\langle x_m \right\rangle _c^2\) and using that the microstates m are decorrelated yields

    $$\begin{aligned} v_c&=\frac{1}{N_m}\times \left( \frac{1}{N_m}\sum _m v_{cm}\right) \\ Var[v_c]&=\frac{1}{N_m^3}\times \left( \frac{1}{N_m}\sum _m \hbox {var} [v_{cm}]\right) \end{aligned}$$

    where the m-averages (brackets) do not depend on \(N_{{\mathrm{m}}}\). This implies \(\Delta _{{\mathrm{ne}}}^2 \approx v^2/N_{{\mathrm{m}}}\) confirming \(\gamma _{{\mathrm{ext}}}=1/2\), Eq. (48).

  12. In agreement with Eq. (47) and assuming \(\delta v_\text {G}\propto 1/\Delta t^{\beta }\) with \(\beta \approx 1/2\) we have \(\tau _{\text {ne}} \propto N_\text {m}^{\gamma _{\text {ext}/\beta }}\), i.e., the crossover time increases with \(N_\text {m}\) and \(\Delta _\text {ne}^{1/\beta }\propto 1/\tau _{\text {ne}}\).

  13. That two properties are equal on average does, of course, not imply that they are identical since other moments (fluctuations) may be different. This is the case for \(R_\text {A}\) and \(m_{21}\) which have different standard deviations \(\delta R_\text {A} \ll \delta m_{21}\) as seen from Fig. 8 of Ref. [10].

  14. The data of the binary LJ mixture scanned from the first panel of Fig. 2 of Ref. [40] corresponds strictly speaking to the standard deviation \(\delta M\) for the shear modulus \(M=R_\text {A}-v\). Since the fluctuations of \(R_\text {A}\) are negligible, however, as shown elsewhere [10, 15], \(\delta M \approx \delta v\).

  15. This is called Monte Carlo gauge in previous work [4, 5]. Note that an ideal gas has a vanishing shear modulus. Thus, the ideal contributions to the stress-fluctuation formula \(M = R_{\mathrm{A}}- v\) rigorously cancel.

  16. A forth example may be found in Ref. [15] where three-dimensional glass-forming polymer melts are presented.

  17. Due to the strong repulsion of the beads and the high number density, the bead distribution is always macroscopically homogeneous and the overall density fluctuations are weak. This has been checked using snapshots and the standard radial pair correlation function g(r) and its Fourier transform S(q) as discussed elsewhere [10].

  18. This allows us to suppress long-range hydrodynamic modes otherwise relevant for two-dimensional systems.

  19. Different lateral box sizes L and film thicknesses H are presented in Ref. [14].

  20. The vertical box size \(L_z\) is chosen sufficiently large (\(L_z \gg H\)) to avoid any interaction in this direction. The instantaneous stress tensor vanishes outside the films. While this implies for all z-planes within the films that the average vertical normal stress must vanish [14], the tangential normal stresses must be finite at the film surface since otherwise the surface tension \(\Gamma \) [25] would vanish and the film be unstable. Note that \(\Gamma \approx 1.7\) at \(T \approx T_{{\mathrm{g}}}\) decreasing weakly with T. Importantly, the surface tension gives rise to a non-negligible contribution to \(R_{{\mathrm{A}}} \) due to the last term of Eq. (B5). The shear modulus M does not vanish at high temperatures (as it must) if this term is not taken into account as in some studies [54].

  21. In agreement with Ref. [15], \(R_{{\mathrm{A}}}(T) \approx m_{21}(T) \propto \rho (T)\) down to \(T \approx T_{{\mathrm{frac}}}\).

  22. Experimental studies on polymer films [50, 51] rather investigate the (shear-strain) creep compliance J(t) [16]. Since the Laplace transforms of J(t) and R(t) are reciprocally related [16], if one function is precisely known, the other can be calculated in principle

  23. Our results are consistent with the continuous decay of the storage modulus \(G^{\prime }(\omega =const,T)\) shown in Fig. 6 of Ref. [53]. Similar continuous behavior has also been reported for the Young modulus [54].

  24. The visible minor differences are due to numerical difficulties related to the finite time step and the inaccurate integration of the strongly oscillatory R(t) at short times.

References

  1. W. Press, S. Teukolsky, W. Vetterling, B. Flannery, Numerical Recipes in FORTRAN: The Art of Scientific Computing (Cambridge University Press, Cambridge, 1992)

    Google Scholar 

  2. B. Schnell, H. Meyer, C. Fond, J.P. Wittmer, J. Baschnagel, Eur. Phys. J. E 34, 97 (2011)

    Article  Google Scholar 

  3. H. Xu, J.P. Wittmer, P. Polińska, J. Baschnagel, Phys. Rev. E 86, 046705 (2012)

    Article  Google Scholar 

  4. J.P. Wittmer, H. Xu, P. Polińska, F. Weysser, J. Baschnagel, J. Chem. Phys. 138, 12A533 (2013)

    Article  Google Scholar 

  5. J.P. Wittmer, H. Xu, J. Baschnagel, Phys. Rev. E 91, 022107 (2015)

    Article  Google Scholar 

  6. J.P. Wittmer, H. Xu, O. Benzerara, J. Baschnagel, Mol. Phys. 113, 2881 (2015)

    Article  Google Scholar 

  7. J.P. Wittmer, I. Kriuchevskyi, J. Baschnagel, H. Xu, Eur. Phys. J. B 88, 242 (2015)

    Article  Google Scholar 

  8. D. Li, H. Xu, J.P. Wittmer, J. Phys. Condens. Matter 28, 045101 (2016)

    Article  Google Scholar 

  9. J.P. Wittmer, H. Xu, J. Baschnagel, Phys. Rev. E 93, 012103 (2016)

    Article  Google Scholar 

  10. J.P. Wittmer, I. Kriuchevskyi, A. Cavallo, H. Xu, J. Baschnagel, Phys. Rev. E 93, 062611 (2016)

    Article  Google Scholar 

  11. I. Kriuchevskyi, J.P. Wittmer, O. Benzerara, H. Meyer, J. Baschnagel, Eur. Phys. J. E 40, 43 (2017)

    Article  Google Scholar 

  12. I. Kriuchevskyi, J.P. Wittmer, H. Meyer, J. Baschnagel, Phys. Rev. Lett. 119, 147802 (2017)

    Article  Google Scholar 

  13. I. Kriuchevskyi, J.P. Wittmer, H. Meyer, O. Benzerara, J. Baschnagel, Phys. Rev. E 97, 012502 (2018)

    Article  Google Scholar 

  14. G. George, I. Kriuchevskyi, H. Meyer, J. Baschnagel, J.P. Wittmer, Phys. Rev. E 98, 062502 (2018)

    Article  Google Scholar 

  15. L. Klochko, J. Baschnagel, J.P. Wittmer, A.N. Semenov, J. Chem. Phys. 151, 054504 (2019)

    Article  Google Scholar 

  16. J.D. Ferry, Viscoelastic Properties of Polymers (Wiley, New York, 1980)

    Google Scholar 

  17. W.W. Graessley, Polymeric Liquids & Networks: Dynamics and Rheology (Garland Science, London, 2008)

    Google Scholar 

  18. M. Doi, S.F. Edwards, The Theory of Polymer Dynamics (Clarendon Press, Oxford, 1986)

    Google Scholar 

  19. M. Rubinstein, R.H. Colby, Polymer Physics (Oxford University Press, Oxford, 2003)

    Google Scholar 

  20. J.P. Hansen, I.R. McDonald, Theory of Simple Liquids, 3rd edn. (Academic Press, New York, 2006)

    Google Scholar 

  21. W. Götze, Complex Dynamics of Glass-Forming Liquids: A Mode-Coupling Theory (Oxford University Press, Oxford, 2009)

    Google Scholar 

  22. P.M. Chaikin, T.C. Lubensky, Principles of Condensed Matter Physics (Cambridge University Press, Cambridge, 1995)

    Book  Google Scholar 

  23. E.B. Tadmor, R.E. Miller, R.S. Elliot, Continuum Mechanics and Thermodynamics (Cambridge University Press, Cambridge, 2012)

    Google Scholar 

  24. E.B. Tadmor, R.E. Miller, Modeling Materials (Cambridge University Press, Cambridge, 2011)

    Book  Google Scholar 

  25. M.P. Allen, D.J. Tildesley, Computer Simulation of Liquids, 2nd edn. (Oxford University Press, Oxford, 2017)

    Book  Google Scholar 

  26. D.P. Landau, K. Binder, A Guide to Monte Carlo Simulations in Statistical Physics (Cambridge University Press, Cambridge, 2000)

    Google Scholar 

  27. J.L. Lebowitz, J.K. Percus, L. Verlet, Phys. Rev. 153, 250 (1967)

    Article  Google Scholar 

  28. N.G. van Kampen, Stochastic Processes in Physics and Chemistry (North-Holland, Amsterdam, 1992)

    Google Scholar 

  29. S. Provencher, Comput. Phys. Commun. 27, 229 (1982)

    Article  Google Scholar 

  30. G. George, L. Klochko, A. Semenov, J. Baschnagel, J.P. Wittmer, under preparation (2020)

  31. H. Yoshino, J. Chem. Phys. 136, 214108 (2012)

    Article  Google Scholar 

  32. A. Einstein, Phys. Zeitschrift 10, 185 (1909)

    Google Scholar 

  33. F. Smallenburg, L. Leibler, F. Sciortino, Phys. Rev. Lett. 111, 188002 (2013)

    Article  Google Scholar 

  34. S. Roldán-Vargas, F. Smallenburg, W. Kob, F. Sciortino, J. Chem. Phys. 139, 244910 (2013)

    Article  Google Scholar 

  35. D. Montarnal, M. Capelot, F. Tournilhac, L. Leibler, Science 334, 965 (2011)

    Article  Google Scholar 

  36. C. Tonhauser, D. Wilms, Y. Korth, H. Frey, C. Friedrich, Macromol. Rapid Commun. 31, 2127 (2010)

    Article  Google Scholar 

  37. A. Zilman, J. Kieffer, F. Molino, G. Porte, S.A. Safran, Phys. Rev. Lett. 91, 015901 (2003)

    Article  Google Scholar 

  38. G. Hed, S. Safran, Eur. Phys. J. E 19, 69 (2006)

    Article  Google Scholar 

  39. V. Testard, J. Oberdisse, C. Ligoure, Macromolecules 41, 7219 (2008)

    Article  Google Scholar 

  40. I. Procaccia, C. Rainone, C.A.B.Z. Shor, M. Singh, Phys. Rev. E 93, 063003 (2016)

    Article  Google Scholar 

  41. E.E. Ferrero, K. Martens, J.-L. Barrat, Phys. Rev. Lett. 113, 248301 (2014)

    Article  Google Scholar 

  42. L. Klochko, J. Baschnagel, J.P. Wittmer, A.N. Semenov, Soft Matter 14, 6835 (2018)

    Article  Google Scholar 

  43. L. Berthier, E. Flenner, H. Jacquin, G. Szamel, Phys. Rev. E 81, 031505 (2010)

  44. L. Berthier, H. Jacquin, F. Zamponi, Phys. Rev. E 84, 051103 (2011)

    Article  Google Scholar 

  45. J.P. Wittmer, A. Tanguy, J.-L. Barrat, L. Lewis, Europhys. Lett. 57, 423 (2002)

    Article  Google Scholar 

  46. A. Tanguy, J.P. Wittmer, F. Leonforte, J.-L. Barrat, Phys. Rev. B 66, 174205 (2002)

    Article  Google Scholar 

  47. A. Ninarello, L. Berthier, D. Coslovich, Phys. Rev. X 7, 021039 (2017)

    Google Scholar 

  48. S.J. Plimpton, J. Comp. Phys. 117, 1 (1995)

    Article  Google Scholar 

  49. S. Frey, F. Weysser, H. Meyer, J. Farago, M. Fuchs, J. Baschnagel, Eur. Phys. J. E 38, 11 (2015)

    Article  Google Scholar 

  50. P.A. O’Connell, G.B. McKenna, Science 307, 1760 (2005)

    Article  Google Scholar 

  51. P.A. O’Connell, S.A. Hutcheson, G.B. McKenna, J. Polym. Sci. Part B Polym. Phys. 46, 1952 (2008)

    Article  Google Scholar 

  52. P. Chapuis, P.C. Montgomery, F. Anstotz, A. Leong-Hoi, C. Gauthier, J. Baschnagel, G. Reiter, G.B. McKenna, A. Rubin, Rev. Sci. Instrum. 88, 093901 (2017)

    Article  Google Scholar 

  53. K. Yoshimoto, T.S. Jain, P.F. Nealey, J.J. de Pablo, J. Chem. Phys. 122, 144712 (2005)

    Article  Google Scholar 

  54. A. Shavit, R.A. Riggleman, Macromolecules 46, 5044 (2013)

    Article  Google Scholar 

Download references

Acknowledgements

We are indebted to O. Benzerara (ICS, Strasbourg) for helpful discussions and acknowledge important computational resources from the computing center of the University of Strasbourg through GENCI/EQUIP@MESO.

Author information

Authors and Affiliations

Authors

Contributions

JB and JPW designed the research project. The presented theory was gathered from different sources and further developed by ANS and JPW. GG (polymer films), LK (pLJ particles) and JPW (TSANET) performed the simulations and the analysis of the data. JPW wrote the manuscript, benefiting from contributions of all authors.

Corresponding author

Correspondence to J. P. Wittmer.

Appendices

Appendix A: Reformulations of \(\delta v_{\mathrm {G}}[f]\)

Since for large I the sums over two, three or even four indices stated in Sect. 2.4.2 rapidly become numerically unfeasible, it is of importance that the three terms \(T_2\), \(T_4\) and \(T_3\) of Eq. (26) can be simplified to single loops [15]. The first two terms simply become

$$\begin{aligned} T_2(\Delta t)= & {} \frac{2}{I^2} \left( I f_0^2 + 2 \sum _{k=1}^{I-1} (I-k) f_k^2 \right) \end{aligned}$$
(A1)
$$\begin{aligned} T_4(\Delta t)= & {} \frac{2}{I^4} \left( I f_0 + 2 \sum _{k=1}^{I-1} (I-k) f_k \right) ^2 \end{aligned}$$
(A2)

with \(f_k = f(t=t_k)\) for an arbitrary function f(t). Let us define the sum \(S(s,I) \equiv \sum _{i=1}^{I} f_{i-s}\). Note that S(sI) may be computed starting from S(0, I) using the recursion relation \(S(s+1,I) = S(s,I) + f_s - f_{I-s}\). Using this the calculation of

$$\begin{aligned} T_3(\Delta t) = \frac{4}{I^3} \sum _{s=1}^{I} S(s,I)^2 = \frac{8}{I^3} \sum _{s=1}^{I/2} S(s,I)^2 \end{aligned}$$
(A3)

becomes also of order \({{\mathcal {O}}}(I)\). Using the symmetry \(S(s,I)=S(I-s+1,I)\) we have assumed in the last step that I is even. In the continuum limit for large \(I=\Delta t/\delta t\) the three terms further simplify to

$$\begin{aligned} T_2(\Delta t)= & {} \frac{4}{\Delta t^2} \int _0^{\Delta t} \mathrm {d}t \ (\Delta t-t) f(t)^2 \end{aligned}$$
(A4)
$$\begin{aligned} T_4(\Delta t)= & {} \frac{8}{\Delta t^4} \left( \int _0^{\Delta t} \mathrm {d}t \ (\Delta t-t) f(t) \right) ^2 \end{aligned}$$
(A5)
$$\begin{aligned} T_3(\Delta t)= & {} \frac{8}{\Delta t^3} \int _0^{\Delta t/2} \mathrm {d}t \left( \eta (t)+\eta (\Delta t-t) \right) ^2 \end{aligned}$$
(A6)

using \(\eta (t) \equiv \int _0^t \mathrm {d}u \ f(u) \) for the last contribution.

Appendix B: Shear stress and Born–Lamé coefficient

Let us consider a small simple shear strain [23] increment \(\gamma \) in the xy-plane as it would be used to measure the shear-stress relaxation function R(t) by means of a direct out-of-equilibrium simulation [5,6,7, 25]. Assuming that all particle positions \({\underline{r}}\) follow an imposed “macroscopic” shear in an affine manner according to \(r_x\rightarrow r_x+ \gamma \ r_y\) the Hamiltonian \({\hat{H}}\) of a given configuration changes to leading order as [6, 10, 11, 13]

$$\begin{aligned} ({\hat{H}}(\gamma )-{\hat{H}}(\gamma =0))/V \approx {\hat{\sigma }}\gamma + \frac{1}{2} {\hat{\mu }}_\mathrm {A}\gamma ^2 \text{ for } |\gamma | \ll 1.\qquad \end{aligned}$$
(B1)

The instantaneous shear stress \({\hat{\sigma }}\) and the instantaneous Born–Lamé coefficient \({\hat{\mu }}_\mathrm {A}\) are thus defined as

$$\begin{aligned} {\hat{\sigma }}\equiv & {} {\hat{H}}^{\prime }(\gamma )/V|_{\gamma =0} \text{ and } \end{aligned}$$
(B2)
$$\begin{aligned} {\hat{\mu }}_\mathrm {A}\equiv & {} {\hat{H}}^{\prime \prime }(\gamma )/V|_{\gamma =0} = {\hat{\sigma }}^{\prime }(\gamma )|_{\gamma =0} \end{aligned}$$
(B3)

where a prime denotes a functional derivative with respect to the affine small strain transform. All properties considered here refer to the excess contributions due to the potential part of the Hamiltonian, i.e., the ideal contributions are assumed to be integrated out.Footnote 15 Assuming a pairwise central conservative potential \(\sum _l u(r_{l})\) with \(r_{l}\) being the distance between a pair of monomers l, one obtains the excess contributions [6, 14]

$$\begin{aligned} {\hat{\sigma }}= & {} \frac{1}{V} \sum _l r_{l}u^{\prime }(r_{l}) \ n_{l,x}n_{l,y} \ \text{ and } \end{aligned}$$
(B4)
$$\begin{aligned} {\hat{\mu }}_\mathrm {A}= & {} \frac{1}{V} \sum _l \left( r_{l}^2 u^{\prime \prime }(r_{l}) - r_{l}u^{\prime }(r_{l}) \right) n_{l,x}^2 n_{l,y}^2 \nonumber \\+ & {} \frac{1}{V} \sum _l r_{l}u^{\prime }(r_{l}) \ (n_{l,x}^2+n_{l,y}^2)/2 \end{aligned}$$
(B5)

with \({\underline{n}}_l= {\underline{r}}_l/r_{l}\) being the normalized distance vector. Note that Eq. (B4) is strictly identical to the corresponding off-diagonal term of the Kirkwood stress tensor [25]. We have used a symmetric representation for the last term of Eq. (B5) exchanging x and y for the affine transform and averaging over the equivalent x and y directions. Note that this last term automatically takes into account the finite normal pressure of the system. Similar relations are obtained for the xz- and the yz-plane. See Refs. [3, 15] for the corresponding expression of the ensemble average of \({\hat{\mu }}_\mathrm {A}\) in terms of the pair correlation functions of the bonded and the non-bonded interactions of the particles. Please note that \({\hat{\mu }}_\mathrm {A}\) depends on the second derivative \(u^{\prime \prime }(r)\) of the pair potential. As emphasized elsewhere [3], impulsive corrections need to be taken into account due to this term if the first derivative \(u^{\prime }(r)\) of the potential is not continuous. Unfortunately, this is the case at the cutoff \(r_\mathrm {cut}\) of the LJ potentials used for the pLJ beads (“Polydisperse Lennard–Jones particles” of appendix) and for the polymer films (“Free-standing polymer films” of appendix). The “bare” \({\hat{\mu }}_\mathrm {A}\) is thus roughly about 0.2 too high for both models and must be corrected [3]. The ensemble average \(\left\langle {\hat{\mu }}_\mathrm {A}\right\rangle \) is called \(\mu _\mathrm {A}\) in previous publications [2,3,4,5,6,7,8,9,10,11,12,13,14,15] and \(R_\mathrm {A}\) in the present main text.

Appendix C: Numerical models and technical details

1.1 Introduction

Various relations and issues discussed theoretically in Sect. 2 are illustrated in Sect. 3 and in “Appendix D” for the fluctuations of shear stresses in simple coarse-grained model systems. A sketch of the threeFootnote 16 models studied is given in Fig. 5. Standard MD and MC methods [25, 26] are used and in some cases combined. The number of particles n and the temperature T are imposed while the system volume V is allowed to fluctuate in some cases. Boltzmann’s constant \(k_\mathrm {B}\), the typical size of the particles (beads) and the particle mass are all set to unity and Lennard–Jones units [25] are used throughout this work. For all systems we use periodic boundary conditions [25, 26] and the number density \(\rho = n/V\) is close to unity. The salient features of each model and some algorithmic details are given below.

1.2 Transient self-assembled networks

As explained in detail in Ref. [10], we use a simple model for transient self-assembled networks (TSANET) in \(d=2\) dimensions where repulsive “harmonic spheres” [43, 44] are reversibly bridged by ideal springs. As shown in panel (a) of Fig. 5, it is assumed that the springs break and recombine locally with an MC hopping frequency \(\nu \). The particles are monodisperse and the temperature T is set to unity. The body of our numerical results has been obtained using periodic simulation boxes of linear size \(L=100\) containing 40,000 springs and \(n=10{,}000\) beads, i.e., \(\rho = n/V = 1\).Footnote 17 We also report in Sect. 3.3 on data for the same number densities of particles and springs obtained for \(n=V=L^2 = 100\), 300, 1000, 3000, 30, 000 and 100, 000. In addition to the MC moves, changing the connectivity matrix of the network, standard MD simulation with a strong Langevin thermostat [25]Footnote 18 is used to move the particles effectively by overdamped motion. We systematically vary \(\nu \) using a constant total sampling time \({\Delta t}_{\mathrm {max}}= 10^5\). For each trajectory we store every \(\delta t= 0.01\) instantaneous properties. Our configurations are first equilibrated in the liquid regime (\(\nu = 1\)). We then gradually decrease \(\nu \) over several orders of magnitude. The hopping moves are switched off (\(\nu = 0\)) for quenched networks. All properties are averaged over \(N_\mathrm {c}=100\) configurations.

1.3 Polydisperse Lennard–Jones particles

We present in “Polydisperse Lennard–Jones particles” of appendix data obtained for polydisperse Lennard–Jones (pLJ) particles in \(d=2\) dimensions [3, 4, 25, 45, 46] where the interaction range \(\sigma = (D_i + D_j)/2\) is set by the Lorentz rule [20] with \(D_i\) and \(D_j\) being the diameters of the interacting particles i and j. Following Refs. [3, 4, 46, 47] the particle diameters are uniformly distributed between 0.8 and 1.2. The standard particle number used is \(n=10{,}000\). Additional results obtained for \(n=100\), 200, 500, 1000, 2000 and 5000 are given in Sect. 3.3. See panel (b) of Fig. 5 for a snapshot of a configuration at \(T=0.2\). By means of an MC barostat [4] we impose an overall pressure \(P=2\) to equilibrate the systems at a given temperature T.

As shown in Fig. 11, the number density \(\rho (T)\) of the older data from Ref. [4] reveal two distinct linear slopes which have been used to define a glass transition temperature \(T_{\mathrm {g}}\approx 0.26\). For details see Ref. [4]. We present here new data [48] where we have used in addition MC swap moves [47] exchanging pairs of particles of different sizes. Time is measured in terms of the local MC steps per particle as in Ref. [4]. We temper now each of the \(N_\mathrm {c}=100\) independent configurations over \(10^7\) MC steps at each temperature with switched on barostat and local and swap MC moves. We then fix the volume and temper over again \(10^7\) MC steps with local and swap moves and over additional \(10^7\) MC steps only with local MC moves. Production runs only using local MC moves are then performed over \({\Delta t}_{\mathrm {max}}=10^7\). For temperatures between \(T=0.21\) and \(T=0.25\) additional production runs have been done for \(N_\mathrm {c}=20\) over \({\Delta t}_{\mathrm {max}}=10^8\) MC steps. The data are normally sampled in intervals of \(\delta t=10\) MC steps. The number densities obtained from these new simulations are indicated in Fig. 11 (circles). As can be seen, much higher densities are achieved at low temperatures and \(\rho (T)\) follows now the high-temperature slope (solid line) well below \(T_{\mathrm {g}}\). Below \(T_{\mathrm {frac}}\approx 0.15\) the fractionation (demixing) of particles of different size is observed, i.e., similar sized particles get closer. As shown in “Polydisperse Lennard–Jones particles” of appendix, the glass transition temperature \(T_{\mathrm {g}}\approx 0.26\), determined originally from the two density slopes, remains useful since it indicates roughly the temperature below which only using local MC moves the terminal relaxation time \(\tau _{\alpha }(T)\) obtained from the relaxation of the shear-stress fluctuations exceeds \({\Delta t}_{\mathrm {max}}\).

1.4 Free-standing polymer films

As sketched in panel (c) of Fig. 5, we study by means of MD simulation [25] of a bead-spring model [48] free-standing polymer films suspended parallel to the (xy)-plane [14]. All unconnected monomers interact with a truncated and shifted LJ potential while connected monomers are bonded by harmonic springs [2, 11,12,13,14,15, 48, 49]. The films contain \(M=768\) monodisperse chains of length \(N=16\), i.e., in total \(n=12{,}288\) monomers, in a periodic box of lateral box size \(L=23.5\).Footnote 19\(^{,}\)Footnote 20 Entanglement effects are irrelevant for the short chains [16,17,18,19]. As discussed in detail in Ref. [14], ensembles of \(N_\mathrm {c}=100\) independent configurations are generated for a broad range of temperatures. Production runs are performed over \({\Delta t}_{\mathrm {max}}=10^5\) storing data each \(\delta t=0.05\). A central parameter for the description of our films is the film thickness H. As shown in Fig. 12, using a Gibbs dividing surface construction and measuring the midplane density \(\rho _0 = \rho (z \approx 0)\) the film thickness is defined as \(H \equiv N M /\rho _0 L^2\) and the film volume as \(V = H L^2\). As may be seen from the main panel, H decreases monotonically upon cooling with two linear branches fitting reasonably the glass (dashed line) and the liquid (solid line) limits. \(T_{\mathrm {g}}\approx 0.37\) is estimated from the intercept of both asymptotes.

1.5 Data handling

As indicated above we equilibrate for each state point of the considered model \(N_\mathrm {c}= 100\) independent configurations c. This allows to probe all properties accurately. For each configuration c we compute and store one long trajectory with \({\Delta t}_{\mathrm {max}}/\delta t\approx 10^7\) data entries. Since we want to investigate the dependence of various properties on the sampling time \(\Delta t\) we probe for each \({\Delta t}_{\mathrm {max}}\)-trajectory \(N_\mathrm {k}\) equally spaced subintervals k of length \(\Delta t\le {\Delta t}_{\mathrm {max}}\) with \(I=\Delta t/\delta t\) entries. It is inessential for all properties discussed in the present work whether these subintervals do partially overlap or do not. Since overlapping subintervals probe similar information it is, however, numerically not efficient to pack them too densely. We use below \(N_\mathrm {k}= {\Delta t}_{\mathrm {max}}/\Delta t\), i.e., \(N_\mathrm {k}\) and \(\Delta t\) are thus coupled and the accuracy is better for small \(\Delta t\).

Appendix D: Shear-stress fluctuations continued

1.1 Polydisperse Lennard–Jones particles

We turn now to the shear-stress fluctuations of the pLJ model. As mentioned in “Appendix D”, our configurations have been first tempered and annealed at a constant pressure \(P=2\) using in addition to the standard local MC moves [4, 26] swap moves exchanging the particle diameters [47]. That this changes dramatically the stress relaxation and thus the equilibration of the configurations can be seen in Fig. 13 comparing the shear-stress relaxation function R(t) for both methods at different temperatures. It is seen that the case with swap moves (filled symbols) decays orders of magnitude faster than the standard method only using local moves. (While all the data presented elsewhere for the pLJ model refers to productions runs with \({\Delta t}_{\mathrm {max}}=10^7\) MC steps, we indicate here for the latter case also some temperatures sampled with \({\Delta t}_{\mathrm {max}}=10^8\).) As a quick way to characterize the terminal relaxation time \(\tau _{\alpha }(T)\) one may set, say, \(R(t \approx \tau _{\alpha }) \approx 0.1\). This implies for the swap moves that \(\tau _{\alpha }\ll 10^5\) for all temperatures above \(T_{\mathrm {frac}}\). For instance, \(\tau _{\alpha }(T=0.2) \approx 10^4\) is three orders of magnitudes smaller than the production time \({\Delta t}_{\mathrm {max}}=10^7\). This suggests that our pLJ samples are well equilibrated and one expects especially \(R_\mathrm {A}= m_{21}\) to hold for temperatures above \(T_{\mathrm {frac}}\). As may be seen from the open symbols, this equilibration does not change much the temperature \(T_{\mathrm {g}}\) where the relaxation time \(\tau _{\alpha }\) exceeds the production time \({\Delta t}_{\mathrm {max}}=10^7\) if the swap moves are switched off. As emphasized by the double-headed arrow, a similar value \(T_{\mathrm {g}}\approx 0.26\) is found as in our previous study [4].

Fig. 13
figure 13

Shear-stress relaxation function R(t) for pLJ particles for several temperatures T comparing the case with additional swap moves (filled symbols) with the standard method. The thin solid lines indicate the phenomenological stretched exponential \(R(t) \approx 35 \exp (-(t/\tau _{\alpha })^{1/3})\) with \(\tau _{\alpha }\propto 1/T^{5/2}\) fitting the swapped systems for \(T > T_{\mathrm {frac}}\). As shown by the double-headed arrow, \(\tau _{\alpha }(T)\), obtained with local moves only, becomes similar to \({\Delta t}_{\mathrm {max}}=10^7\) slightly below \(T_{\mathrm {g}}\approx 0.26\)

Figure 14 presents various properties of the pLJ model for one constant \(t=\Delta t=10^5\) as functions of the temperature T. As expected for a well-equilibrated liquid, \(R_\mathrm {A}\approx m_{21}\) holds at least down to \(T_{\mathrm {frac}}\). (This was not the case for the older data in Ref. [4].) Moreover, they decrease, more or less, linearly with increasing T.Footnote 21\(R_\mathrm {A}\) is the upper bound of the stress-fluctuation \(v(\Delta t)\) for any \(\Delta t\) since \(M(\Delta t) \equiv R_\mathrm {A}-v(\Delta t) \ge 0\). Note that v(T) increases first monotonically with increasing temperatures until it merges continuously at \(T \approx 0.3\) with \(R_\mathrm {A}\) decreasing then together with \(R_\mathrm {A}\). This implies that after a monotonic decay at low T, the shear modulus M(T) must also vanish continuously at \(T \approx 0.3\) for \(\Delta t=10^5\) (not shown). The relaxation function R, being related to M through the general relation Eq. (15), thus vanishes continuously at about the same temperature.

Fig. 14
figure 14

Various properties obtained for the pLJ model taken at \(t=\Delta t=10^5\) versus T: affine shear modulus \(R_\mathrm {A}\), moment \(m_{21}\), shear-stress fluctuation v and its standard deviation \(\delta v\), mean-square displacement h and its standard deviation \(\delta h\) and shear-stress relaxation function \(R=R_\mathrm {A}-h\). Note that \(R_\mathrm {A}\approx m_{21}\) and \(h \approx v\) for all T and that \(\delta h \approx 2^{1/2} h\) (bold solid line) holds in agreement with Eq. (25)

The mean-square displacement h of the instantaneous shear-stresses (small, filled diamonds), Eq. (8), is more or less identical to v for all T. (Closer inspection reveals that v and h slightly differ between \(T \approx 0.25\) and \(T \approx 0.3\).) The standard deviation \(\delta h\) of h was calculated using Eq. (24). As discussed in Sect. 2.4.1, Eq. (25) must hold for Gaussian stochastic processes. As shown by the bold solid line, this is nicely confirmed to high precision. A log-linear plot of the non-Gaussianity parameter \(\alpha _2(t) = \delta h^2/2h^2-1\) demonstrates that \(|\alpha _2(t)| \ll 0.02\) for all times and temperatures (not shown). Moreover, since \(h \approx v\) for all T and \(h = v = R_\mathrm {A}\) above \(T \approx 0.3\), this implies the observed non-monotonic behavior of \(\delta h(T)\). We have also checked that \(\delta R \approx \delta h\) (not shown) as one expects [10, 15] since \(R_\mathrm {A}\) barely fluctuates. These findings confirm that Gaussian processes are dominant for all T, even below the demixing transition at \(T_{\mathrm {frac}}\).

Also indicated at the bottom of Fig. 14 is the standard deviation \(\delta v\) of v. Figure 15 compares \(\delta v\) (open symbols) and \(\delta v_{\mathrm {G}}\) (lines) as functions of T for a broad range of \(\Delta t\). Naturally, our data get less accurate for \(\Delta t\rightarrow {\Delta t}_{\mathrm {max}}\). The most important point is here that within numerical precision \(\delta v \approx \delta v_{\mathrm {G}}\) for all \(\Delta t\) for large T even slightly below \(T_{\mathrm {g}}\). With increasing \(\Delta t\) the maximum of \(\delta v(T) \approx \delta v_{\mathrm {G}}(T)\) becomes sharper and shifts to lower T. Interestingly, \(\delta v \approx \delta v_{\mathrm {G}}\) even holds for \(T \ll T_{\mathrm {g}}\), but only for \(\Delta t\ll 10^5\). However, while \(\delta v_{\mathrm {G}}\) decreases strongly upon cooling for larger \(\Delta t\) this is not observed for \(\delta v\) being only weakly T-dependent. Thus, Eq. (3) does not hold in this limit.

Fig. 15
figure 15

\(\delta v\) (open symbols) and \(\delta v_{\mathrm {G}}\) (lines) as functions of T for a broad range of \(\Delta t\). While \(\delta v \approx \delta v_{\mathrm {G}}\) holds for large T and small \(\Delta t\), \(\delta v\) is seen to become constant for low T and high \(\Delta t\). An increasingly sharp maximum appears around \(T_{\mathrm {g}}\) shifting with increasing \(\Delta t\) to lower T

Focusing on one low temperature (\(T=0.2\)) Fig. 16 compares v with \(\delta v\) and \(\delta v_{\mathrm {G}}\) as functions of \(\Delta t\). In agreement with Eq. (11), \(v(\Delta t)\) increases monotonically becoming rapidly constant. Since consistently with Eq. (13) also h(t) and R(t) become constant (Fig. 13), this implies \(\delta v_{\mathrm {G}}(\Delta t) \propto 1/\sqrt{\Delta t}\) (bold solid line). At variance to the decay of \(\delta v_{\mathrm {G}}\), \(\delta v \rightarrow \Delta _{\mathrm {ne}}\approx 0.25\) for \(\Delta t\gg \tau _{\mathrm {ne}}\) (dashed horizontal line) similar to the behavior for quenched TSANET systems (Fig. 8). As in Sect. 3.2 Eq. (46) yields a reasonable approximation of \(\delta v\) albeit the shifted data (stars) are slightly above \(\delta v\) at \(\Delta t\approx \tau _{\mathrm {ne}}\).

Fig. 16
figure 16

Shear-stress fluctuation \(\delta v\) and the corresponding standard deviations \(\delta v\) and \(\delta v_{\mathrm {G}}\) versus \(\Delta t\) at \(T=0.2\) which is well below \(T_{\mathrm {g}}\). \(v(\Delta t)\) is seen to become constant above \(\Delta t\approx 10^3\) where \(v \approx h(t) \approx 17.1\). As indicated by the bold solid line, \(\delta v_{\mathrm {G}}\propto 1/\sqrt{\Delta t}\) as expected in this regime. \(\delta v\) levels off for \(\Delta t\gg \tau _{\mathrm {ne}}\approx 200{,}000\)

Fig. 17
figure 17

\(R(t) = R_\mathrm {A}-h(t)\) for free-standing films using half-logarithmic coordinates for a broad range of temperatures T. R(t) increases continuously with decreasing T. Logarithmic creep behavior (thin solid lines) with \(R(t) \approx a - b \ln (t)\) is found for \(T \ll T_{\mathrm {g}}\) and above the glass transition (\(T \approx 0.45\))

1.2 Free-standing polymer films

The rheological properties of free-standing polymer films may be characterized in a computer experiment by means of the shear-stress relaxation function R(t) as shown in Fig. 17.Footnote 22 Albeit we average over \(N_\mathrm {c}=100\) independent configurations it was necessary for the clarity of the presentation to use gliding averages, Eq. (8), i.e., the statistics becomes worse for \(t \rightarrow {\Delta t}_{\mathrm {max}}=10^5\). Interestingly, it is clearly seen that R(t) increases continuously with decreasing T [14]Footnote 23 in perfect agreement with all published experimental [50,51,52] and computational [53, 54] studies. The presented R(t)-data is used below to obtain \(\delta v_{\mathrm {G}}[R]\).

Fig. 18
figure 18

Contributions to \(M=R_\mathrm {A}-v\): a Temperature dependence for \(\Delta t=10^4\). b Double-logarithmic representation of \(m_{21}/R_\mathrm {A}-1\) versus T. c \(\Delta t\)-effects for M and its contributions for \(T=0.3\). Only \(R_\mathrm {A}\) and \(m_{21}\) are strictly \(\Delta t\)-independent. The dashed lines have been obtained using Eq. (11), the solid line using Eq. (15) taking advantage of the directly measured shear-stress relaxation function R(t)

Figure 18 presents various contributions to the generalized shear modulus \(M=R_\mathrm {A}-v=(R_\mathrm {A}-m_{21})+m_{12}\). At variance to Fig. 14 a half-logarithmic representation is used in panel (a) to better show the decay of \(M \approx m_{12}\) for temperatures above the glass transition where \(R_\mathrm {A}\approx m_{21}\). Note that \(R_\mathrm {A}< m_{21}\) below \(T_{\mathrm {g}}\). As may be seen in panel (b) frozen-in out-of-equilibrium stresses are observed upon cooling below \(T_{\mathrm {g}}\) as made manifest by the dramatic increase of the dimensionless parameter \(y=m_{21}/R_\mathrm {A}-1\). The prefactor \(\beta =1/T\) of \(m_{21}\), Eq. (51), implies due to the frozen stresses the observed \(y(T) \propto 1/T\) for \(T \ll T_{\mathrm {g}}\). (Similar behavior has been reported for three-dimensional polymer bulks [13].) Since there is currently no algorithm comparable to the swap algorithm [47] allowing to equilibrate glass-forming polymer melts and films as for the pLJ model, our films are clearly not at thermal equilibrium below \(T_{\mathrm {g}}\), at least not in the sense of an equilibrium liquid. This does not mean that the stochastic process is not effectively stationary. This is addressed in panel (c) where we test for one temperature below \(T_{\mathrm {g}}\) (where \(R_\mathrm {A}< m_{21}\)) that Eqs. (11) and (15) hold for v (dashed lines) and M (solid line). To test these relations R(t) was integrated numerically.Footnote 24 Albeit time-translational invariance apparently holds, this does not mean that no aging occurs (for \(T < T_{\mathrm {g}}\)) but just that this is irrelevant for the timescales considered here.

We address in panel (a) of Fig. 19 the \(\Delta t\)-dependence of \(\delta v\) and \(\delta v_{\mathrm {G}}\) (open symbols) for different temperatures as indicated. \(\delta v \approx \delta v_{\mathrm {G}}\) holds again for T above and around \(T_{\mathrm {g}}\). Note also that \(\delta v \approx \delta v_{\mathrm {G}}\propto 1/\sqrt{\Delta t}\) (bold solid line) for all \(\Delta t\le {\Delta t}_{\mathrm {max}}\) for the largest temperature \(T=0.55\). Interestingly, T-dependent shoulders appear for \(T \approx 0.4\). This is a consequence of the creep-like decay of R(t) in this regime which is approximately fitted by \(R(t) \approx a - b \ln (t)\) as indicated in Fig. 17. According to Fig. 4, one expects a shoulder with \(\delta v_{\mathrm {G}}\approx 1.55 |b|\). That this holds is seen by the upper thin horizontal line using \(b=1.4\) for \(T=0.4\). In the low-T limit we see again that \(\delta v\) becomes constant, \(\delta v \approx \Delta _{\mathrm {ne}}\) for \(\Delta t\gg \tau _{\mathrm {ne}}\approx 800\), while \(\delta v_{\mathrm {G}}\) continues to decrease with \(\Delta t\). The small deviations of \(\delta v_{\mathrm {G}}\) from the \(1/\sqrt{\Delta t}\)-asymptote are caused by the fact that R(t) does not become rigorously constant, \(R(t) \rightarrow R_p > 0\), even for our lowest temperatures. As shown by the upper thin line in Fig. 17, R(t) is fitted by a logarithmic creep. The amplitude b of this low-temperature creep are, however, too small to lead to a clear-cut shoulder for \(\delta v_{\mathrm {G}}\). As shown by the lower thin solid line with \(b = 0.11\) for \(T=0.05\) at least one decade longer production runs are required to make the low-T creep manifest for \(\delta v_{\mathrm {G}}\). A complementary representation is given in penal (b) of Fig. 19 focusing on the temperature dependence of \(\delta v\) and \(\delta v_{\mathrm {G}}\). A very similar behavior as in Fig. 15 for pLJ particles is seen. Most importantly, \(\delta v \approx \delta v_{\mathrm {G}}\) for large and intermediate T. The maximum of \(\delta v \approx \delta v_{\mathrm {G}}\) slightly below \(T_{\mathrm {g}}\) becomes systematically larger with increasing \(\Delta t\) and shifts to lower temperatures. Also it is confirmed that Eq. (3) holds for sufficiently small \(\Delta t\) for all T. As expected, Eq. (3) breaks down for \(\Delta t\gg \tau _{\mathrm {ne}}(T) \approx 800\) in the non-ergodic limit (lowest temperatures). While \(\delta v_{\mathrm {G}}\), measuring the fluctuations within each configuration, is seen to decrease upon cooling, \(\delta v\) becomes constant, \(\delta v \rightarrow \Delta _{\mathrm {ne}}\) (dashed horizontal line), due to the finite dispersion of the quenched \(v_c\) of the independent configurations.

Fig. 19
figure 19

Comparison of \(\delta v\) (filled symbols) and \(\delta v_{\mathrm {G}}[R]\) (open symbols) for free-standing films: a \(\Delta t\)-dependence of \(\delta v_{\mathrm {G}}[R]\) (open symbols) for all indicated T and of \(\delta v\) (filled symbols) for \(T=0.55\), 0.4, 0.2 and 0.05. The \(1/\sqrt{\Delta t}\)-decay for small \(\Delta t\) is shown by the bold solid line, the plateau value \(\Delta _{\mathrm {ne}}\approx 1\) of \(\delta v\) for small T and \(\Delta t\gg \tau _{\mathrm {ne}}\approx 800\) by the bold dashed line and \(\delta v_{\mathrm {G}}\approx 1.55|b|\) (\(b=1.4\) for \(T=0.4\) and \(b=0.11\) for \(T=0.05\)) expected for logarithmic creep by thin horizontal lines. b Temperature dependence for different \(\Delta t\) as indicated

Rights and permissions

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

George, G., Klochko, L., Semenov, A.N. et al. Ensemble fluctuations matter for variances of macroscopic variables. Eur. Phys. J. E 44, 13 (2021). https://doi.org/10.1140/epje/s10189-020-00004-7

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1140/epje/s10189-020-00004-7

Navigation