G George

L Klochko

A N Semenov

J Baschnagel

J P Wittmer

Ensemble fluctuations matter for variances of macroscopic variables

Ensemble fluctuations matter for variances of macroscopic variables

Introduction

Let us consider a stochastic dynamical variable x(t), a generalized coordinate characterizing a large physical system, like certain density fields averaged over the system volume. Ensembles of discrete time series x = {x i = x(t i ), i = 1, . . . , I} are sampled with the data sequence taken at equally spaced times t i = iδt from t 1 = δt up to the "sampling time" ∆t = Iδt. 1 We focus on the ensemble average v and the standard deviation δv of the (empirical) variance

2 v[x] ≡ 1 I I i=1 x 2 i - 1 I 2 I i,j=1
x i x j .

Extending recent work on stress fluctuations [START_REF] Schnell | [END_REF][3][4][5][6][7][8][9][10][11][12][13][14][15] we want to give a systematic and uncluttered overview of three general points of relevance for a large variety of problems in condensed matter [START_REF] Ferry | Viscoelastic properties of polymers[END_REF][START_REF] Graessley | Polymeric Liquids & Networks: Dynamics and Rheology[END_REF][START_REF] Doi | The Theory of Polymer Dynamics[END_REF][START_REF] Rubinstein | Polymer Physics[END_REF][START_REF] Hansen | Theory of simple liquids[END_REF][START_REF] Götze | Complex Dynamics of Glass-Forming Liquids: A Mode-Coupling Theory[END_REF][START_REF] Chaikin | Principles of condensed matter physics[END_REF], material modeling [START_REF] Tadmor | Continuum Mechanics and Thermodynamics[END_REF][START_REF] Tadmor | Modeling Materials[END_REF] and in computational physics [START_REF] Allen | Computer Simulation of Liquids[END_REF][START_REF] Landau | A Guide to Monte Carlo Simulations in Statistical Physics[END_REF]. One important motivation is that many physical quantities can be obtained by equilibrium molecular dynamics (MD) or Monte Carlo (MC) simulations [START_REF] Allen | Computer Simulation of Liquids[END_REF][START_REF] Landau | A Guide to Monte Carlo Simulations in Statistical Physics[END_REF] using fluctuation relations [START_REF] Lebowitz | [END_REF]. Studying how the respective variances v and their standard deviations δv evolve with the computational feasible length ∆t of the production runs of the simulations is thus of particular interest.

a joachim.wittmer@ics-cnrs.unistra.fr 1 We frequently switch between a discrete and a continuous representation i ↔ t, I ↔ ∆t, hi ↔ h(t), ci ↔ c(t), . . . 2 The empirical variance is defined here without the usual "Bessel correction" [START_REF] Press | Numerical Recipes in FORTRAN: the art of scientific computing[END_REF]. Equation [START_REF] Press | Numerical Recipes in FORTRAN: the art of scientific computing[END_REF] is the formal definition of v[x] which coincides with the genuine variance of x(t) only in the limit ∆t ∝ I → ∞.

We assume here that x(t) is a stationary stochastic process respecting the time-translational invariance [START_REF] Van Kampen | Stochastic processes in physics and chemistry[END_REF]. Our first point is that the expectation value v for sampling times ∆t smaller then the terminal relaxation time τ is not necessarily a ∆t-independent constant as often tacitly assumed. 3 This is seen (Sec. 2.2) from the "stationarity relation" [5,7,10,[12][13][14][15]

v = 2 I 2 I-1 i=1 (I -i) h i with h i-j = (x i -x j ) 2 /2 (2)
being the autocorrelation function (ACF) characterizing the mean-square displacements of the data entries x i . 4Hence, v generally depends on I or ∆t and this is especially relevant if the ACF h i = h(t i ) increases strongly for t ≈ ∆t.

Our second and most central point concerns the standard deviation δv of v[x]. It has been observed for shearstress fluctuations [10,[12][13][14][15] that δv may become rather large and of the order of the mean value v if h(t) varies strongly for t ≈ ∆t, i.e. the mean behavior standard experimental or theoretical work focuses on [START_REF] Ferry | Viscoelastic properties of polymers[END_REF][START_REF] Doi | The Theory of Polymer Dynamics[END_REF][START_REF] Tadmor | Continuum Mechanics and Thermodynamics[END_REF][START_REF] Götze | Complex Dynamics of Glass-Forming Liquids: A Mode-Coupling Theory[END_REF] gets masked by strong fluctuations. Reworking Ref. [15] this can be simply understood assuming a stationary Gaussian stochastic process which implies that

δv = δv G [h] with (3)
δv 2 G [h] ≡ 1 2I 4 I i,j,k,l=1
g 2 ijkl and

g ijkl ≡ (h i-j + h k-l ) -(h i-l + h j-k ).
as shown in Sec. 3.3. By analyzing the functional δv G [h] it will be seen (Secs. 3.4 and 3.5) that while δv(∆t) must remain small for h(t ∼ ∆t) ≈ constant, δv(∆t) becomes generally large if ∆t is similar to the characteristic time of an efficient relaxation pathway corresponding to a strong change of h(t) for t ≈ ∆t.

Our third key point emphasizes one limitation of Eq. ( 3) which hinges on the ergodicity of the stochastic process. If the system is (strictly or in practice) non-ergodic, i.e. if independently created trajectories c are restricted to different meta-basins of the generalized phase space, this implies as shown in Sec. 4 that δv(∆t) → ∆ ne = constant for ∆t τ ne τ b (4) with τ b being the typical relaxation time of the metabasins, τ ne a crossover time defined below and ∆ ne the static standard deviation of the quenched variances v c of the configurations c. In this limit δv(∆t) must thus differ from δv G (∆t) ∝ 1/ √ ∆t for ∆t τ b . However, as argued in Sec. 5, in the common case where the observables x(t) average over many, more or less decoupled microstates, the quenched v c become similar with increasing system size and, hence, δv → δv G [h] in the macroscopic limit even for non-ergodic systems.

Various relations and issues discussed theoretically in Secs. 2-5 are illustrated for different coarse-grained model systems in Sec. 6 and Appendix D. The paper concludes in Sec. 7 with a summary and an outlook to future work. Numerically more convenient reformulations of Eq. ( 3) are given in Appendix A. The definitions of the instantaneous shear stress and the corresponding Born-Lamé coefficient are reminded in Appendix B. The three coarse-grained models simulated are presented in Appendix C together with some technical details related to the data processing (Appendix C.5).

Stationary stochastic processes 2.1 Introduction

Having measured and stored the I entries x i of a time series x = {x i , i = 1, . . . , I} various functionals O[x] may be computed, e.g., the moments m αβ

and correspondingly for h g s [x] with

h s,i = 1 2 (x i+s -x i ) 2 = x 2 i + x 2 i+s 2 -c s,i . (6)
Obviously, c g 0 [x] = m 21 [x] and h g 0 [x] = 0. Averages over a given time series are called "t-averages". Since the functionals O[x] are obtained in general from correlated data entries, ensemble averages . . . of fluctuation-type functionals may depend on the sampling time ∆t. This is not the case for "simple averages" [START_REF] Allen | Computer Simulation of Liquids[END_REF]9,10] for which the ensemble average over independent trajectories and the taverage commute. For instance, we have

m α1 = 1 I I i=1 x α i = 1 I I i=1
x α i ∝ ∆t 0 (7) since the ensemble average x α i is ∆t-independent. Interestingly, the commutation of both averaging-operators is not possible for m αβ with β = 1. An argument ∆t often marks below a property being not a simple average.

Stationarity

We suppose that the time series is taken from a stationary stochastic process whose joint probability distribution does not change when shifted in time [START_REF] Van Kampen | Stochastic processes in physics and chemistry[END_REF]. Correlation functions such as x i x j thus only depend on the difference s = |i -j| of the discrete indices i and j. We thus define

c s = c g s [x] and h s = c 0 -c s = h g s [x] (8)
with 0 ≤ s < I in terms of c g s [x] and h g s [x] defined in Sec. 2.1. Note that both c s and h s are simple averages, i.e. they do not depend on ∆t [7,10]. Note also that c 0 = m 21 = m 21 [x] and h 0 = 0. See Sec. 3.2 for a subtle point related to the fluctuations δc s and δh s .

Due to the assumed stationarity, the ensemble average v = v[x] of Eq. ( 1) becomes [START_REF] Landau | A Guide to Monte Carlo Simulations in Statistical Physics[END_REF][START_REF] Allen | Computer Simulation of Liquids[END_REF]5,7,[10][11][12][13][14][15]

v(∆t) = 1 I I i=1 x 2 i - 1 I 2 I i=1 x 2 i + 2 I 2 I-1 k=1 (I -k) x k+1 x 1 = c 0 (1 -I -1 ) - 2 I 2 I-1 k=1 (I -k)c k = 2 I 2 I-1 i=1 (I -i) h i (9)
as already stated in the Introduction, Eq. ( 2). Note that in the last step it was used that h s = c 0 -c s and 2

I 2 I-1 k=1 (I -k) = 1 -1/I. ( 10
)
In statistical mechanics Eq. ( 9) is closely related to the equivalence of the Green-Kubo and the Einstein relations for transport coefficients [START_REF] Hansen | Theory of simple liquids[END_REF][START_REF] Landau | A Guide to Monte Carlo Simulations in Statistical Physics[END_REF][START_REF] Allen | Computer Simulation of Liquids[END_REF]13,14]. 5 Albeit the mentioned ∆t-dependence is well known [START_REF] Allen | Computer Simulation of Liquids[END_REF][START_REF] Landau | A Guide to Monte Carlo Simulations in Statistical Physics[END_REF] it is emphasized here for systematic reasons and since ∆t-effects for such fluctuations are rarely checked [START_REF] Schnell | [END_REF]4]. We also remind [10][11][12][13][14][15] that in the continuum limit for large I = ∆t/δt, Eq. ( 9) reads

v(∆t) = 2 ∆t 2 ∆t 0 dt (∆t -t) h(t) (11)
with h(t) being the continuum limit of h s . This result may be restated equivalently using the inverse relation

h(t) = [v(t)t 2 /2
] with a prime denoting a derivative with respect to time [13,15]. Using that m 21 is a simple average Eq. ( 11) implies that

m 12 (∆t) = m 21 - 2 ∆t 2 ∆t 0 dt (∆t -t) h(t). ( 12
)
The ensemble averages v(∆t) and m 12 (∆t) thus depend in general on the sampling time ∆t. However, the ∆tdependence disappears, it h(t) becomes constant. For instance, this is the case, if h(t) plateaus in an intermediate, sufficiently large, time window, i.e. h(t) ≈ h p = c(0) -c p with h p and c p being constants. We then have

v(∆t) ≈ h p = c(0) -c p = constant, m 12 (∆t) ≈ m 21 -h p = c p = constant. ( 13
)
Equation ( 13) also holds, if c(t) tends to a constant for times much longer than the terminal relaxation time τ of the system. Then, c p in Eq. ( 13) is replaced by the longtime limit c ∞ = lim t→∞ c(t) = lim ∆t→∞ m 12 (∆t).

Linear response and generalized modulus

The functions h(t) and c(t) can be related to the linear response to an external perturbation conjugate to x(t).

Let R(t) denote the linear response function of the system to a weak external perturbation that is instantaneously switched on at t = 0 and held constant for t > 0. By virtue of the fluctuation-dissipation theorem one can show that [START_REF] Hansen | Theory of simple liquids[END_REF][START_REF] Doi | The Theory of Polymer Dynamics[END_REF]5,10]

R(t) = R A -h(t) = (R A -c(0)) + c(t) (14)
with R A = R(0) being a constant characterizing the initial response of the system after the external perturbation is 5 This may be better seen from the continuum representation of Eq. ( 16) rewritten as

t 0 ds R(s) = d dt [(M (t)t 2 /2].
If the left-hand side converges to a constant η for t → ∞ this implies M (t)t 2 → 2ηt. are given by Eq. ( 14), where we have additionally set c∞ = 0, and M (∆t) and v(∆t) by means of Eqs. ( 16) and ( 18). The two solid horizontal lines mark the intermediate pseudo-plateau for τ1 t τ2. v(∆t) and M (∆t) are seen to converge much more slowly to the respective plateau values than the corresponding response functions h(t) and R(t).

applied. 6 For elastic properties this constant is given by a Born-Lamé affine modulus (Appendix B) [10][11][12][13][14][15]. R(t) is a simple average just as h(t) and c(t). Note that

R ∞ ≡ lim t→∞ R(t) = R A -c(0) + c ∞ ( 15
)
may in general be finite. We rewrite now Eq. ( 11) in terms of R(t) as

M (∆t) ≡ R A -v(∆t) = 2 ∆t 2 ∆t 0 dt (∆t -t)R(t) (16)
with M (∆t) being the "generalized dynamical modulus" [10]. Although this modulus does in general depend on ∆t, it becomes constant

M (∆t) → R ∞ for ∆t/τ → ∞. Being a second integral over R(t) = [M (t)t 2 /2] , M (∆t)
is a smoothing function statistically better behaved than R(t) and containing in general information about both the reversibly stored energy and the dissipation processes.

Generalized Maxwell model

Response functions are often fitted using the generalized Maxwell model

R(t) = R ∞ + pmax p=1 H p exp(-t/τ p ) (17)
with H p and τ p being, respectively, the amplitude and the relaxation time of the mode p [START_REF] Ferry | Viscoelastic properties of polymers[END_REF][START_REF] Rubinstein | Polymer Physics[END_REF]. Note that R(0) = R ∞ + pmax p=1 H p . (Commonly, one considers logarithmic time scales for R(u) with u ∝ log(t) and the modes are distributed logarithmically in time [START_REF] Ferry | Viscoelastic properties of polymers[END_REF][START_REF] Provencher | [END_REF].) Using Eq. ( 16) this implies [10]

M (∆t) = R ∞ + pmax p=1 H p g Debye (∆t/τ p ) (18)
with g Debye (x) = 2 [exp(-x) -1 + x] /x 2 being the Debye function well known in polymer science [START_REF] Doi | The Theory of Polymer Dynamics[END_REF][START_REF] Rubinstein | Polymer Physics[END_REF]. Figure 1 presents both R(t) and M (∆t) for a generalized Maxwell model with two modes with τ 1 τ 2 . The upper solid horizontal line indicates an intermediate pseudo-plateau, Eq. (13). Note also that h

(t) ≈ v(∆t) ≈ R A -R ∞ for t ≈ ∆t τ 2 .
Since v(∆t) and M (∆t) are second integrals over h(t) and R(t), they converge less rapidly to the respective intermediate or terminal plateau values. (h(t) being a monotonically increasing function implies h(t) > v(t) and R(t) < M (t).) As shown by Fig. 1, the determination of a plateau value by means of Eq. (11) or Eq. ( 16) may thus be tedious [4,5,10].

3 Ergodic Gaussian processes

Gaussian variables

Let us consider a Gaussian variable y of variance σ 2 . Since (y -y ) 4 = 3σ 4 we have

z 2 -z 2 = 2σ 4 for z = (y -y ) 2 , (19)
i.e. the variance of the variance z of y is twice the squared variance of y. We assume now that the time series x is a Gaussian process [START_REF] Van Kampen | Stochastic processes in physics and chemistry[END_REF]. (The main physical reason why this assumption holds for many systems is discussed in Sec. 5.) The mean m 11 [x] is thus a Gaussian variable and Eq. ( 19) holds for y = m 11 [x]. Assuming that y = m 11 = 0 by symmetry or by shifting of the data and using that

m α1 [x] β = m αβ [x] this implies [15] δm 2 12 = m 14 -m 2 12 = 2(δm 2 11 ) 2 = 2m 2 12 . (20)

δc and δh for Gaussian processes

Let us next discuss the typical fluctuations of the ACFs c s and h s defined in Sec. 2.2. There are two meaningful ways to define the variances. One characterizes the fluctuations of c g s [x] and h g s [x] by means of

δc g s 2 (I) = c g s [x] 2 -c g s [x] 2 ( 21
)
δh g s 2 (I) = h g s [x] 2 -h g s [x] 2 . ( 22
)
This allows to get the variances and the error bars for the numerical most accurate way to compute c s and h s . The trouble with this definition is that, since the gliding averages are performed first and since the data entries x i are correlated in time, Eq. ( 21) and Eq. ( 22) depend on these correlations in an intricate way. 7 This may mask the fact that the data have a Gaussian distribution. A second way to characterize the fluctuations is to measure in a first step c s,i and h s,i (cf. Sec. 2.1), to take then the ensemble averages

δc 2 s,i = c 2 s,i -c s,i 2 and δh 2 s,i = h 2 s,i -h s,i 2 (23)
and only as the last step (last loop) to take the arithmetic average over all I -s possible indices i, i.e.

δc 2 s = 1 I -s I-s i=1 δc 2 s,i , δh 2 s = 1 I -s I-s i=1 δh 2 s,i . (24)
Assuming x to be Gaussian, y = (x i+s -x i )/ √ 2 is a Gaussian variable of zero mean. According to Eq. ( 19) this implies the important relation

δh 2 s = y 4 -y 2 2 = 2 y 2 2 = 2h 2 s . (25)
In a similar way we find: δc 2 s = c 2 0 +c 2 s . For the fluctuations of R s = R A -h s with R A being constant Eq. ( 25) yields in turn δR 2 s = 2h 2 s . The latter relation may even hold if R A is not strictly constant. This is relevant for the Born-Lamé coefficients considered in Sec. 6 and Appendix D.

δv = δv G [h] for Gaussian processes

We turn now to the derivation of Eq. (3) for the variance

δv 2 ≡ v[x] 2 -v[x]
2 . Using Eq. (1) this may be written

δv 2 = T 2 + T 4 -T 3 with (26)
T 2 ≡ δm 2 21 = m 21 [x] 2 -m 21 [x] 2 = 1 I 2 ij x 2 i x 2 j - 1 I 2 ij x 2 i x 2 j T 4 ≡ δm 2 12 = m 12 [x] 2 -m 12 [x] 2 = 1 I 4 ijkl x i x j x k x l - 1 I 4 ijkl x i x j x k x l T 3 ≡ 2 cov(m 21 , m 12 ) ≡ 2 ( m 21 [x]m 12 [x] -m 21 [x] m 12 [x] ) = 2 I 3 ikl x 2 i x k x l - 2 I 3 kl x 2 i x k x l
where the sums run over all I data entries. As we have assumed that the stochastic process is stationary and Gaussian, Wick's theorem must hold [START_REF] Van Kampen | Stochastic processes in physics and chemistry[END_REF][START_REF] Doi | The Theory of Polymer Dynamics[END_REF] x

i x j x k x l = x i x j x k x l + x i x k x j x l + x i x l x j x k . (27)
Setting in addition c i-j = x i x j it is thus readily seen that the three terms in Eq. ( 26) can be rewritten as

T 2 (∆t) = 2 I 2 ij c 2 i-j (28)
T 4 (∆t) = 2 I 4   ij c i-j   2 (29)
T 3 (∆t) = 4 I 3 s,i,j c i-s c j-s . (30)
Note that T 4 = δm 2 12 = 2m 2 12 in agreement with Eq. ( 20). Numerical more convenient reformulations of T 2 , T 4 and T 3 are given in Appendix A. Importantly, Eqs. [START_REF] Landau | A Guide to Monte Carlo Simulations in Statistical Physics[END_REF][START_REF] Van Kampen | Stochastic processes in physics and chemistry[END_REF][START_REF] Provencher | [END_REF]30) are equivalent to the more compact formula [15]

δv 2 G [c] = 1 2I 4 i,j,k,l g 2 ijkl with (31)
g ijkl = (c i-j + c k-l ) -(c i-l + c j-k )
which looks rather similar as Eq. ( 3). That this holds can be verified by straightforward expansion of Eq. (31). Note that the squared terms c 2 i-j + . . . with two different indices contribute to T 2 , the terms c i-j c k-l + . . . with four different indices to T 4 and the terms c i-j c i-l + . . . with three different indices to T 3 .

With a and b being real constants it follows directly from Eq. ( 31) that

δv G [a] = 0 and δv G [b(f -a)] = |b| δv G [f ] ( 32
)
for any function f (t). Specifically, δv

G [c] = δv G [h]
. This demonstrates finally that Eq. ( 3) is equivalent to Eq. ( 31) and, hence, to Eqs. [START_REF] Landau | A Guide to Monte Carlo Simulations in Statistical Physics[END_REF][START_REF] Van Kampen | Stochastic processes in physics and chemistry[END_REF][START_REF] Provencher | [END_REF]30). It may also be useful to replace c(t) by c(t) -c ∞ or -for thermodynamic equilibrium systems -by the linear response function R(t), Eq. ( 14). We discuss now in Sec. 3.4 some general properties of δv G [f ] and in Sec. 3.5 the behavior of δv G [f ] for various test functions f (t) not necessarily being ACFs.

Some general properties of δv G [f ]

Assuming a constant function f (t) = a one obtains from either Eqs. [START_REF] Van Kampen | Stochastic processes in physics and chemistry[END_REF][START_REF] Provencher | [END_REF]30) or using the corresponding continuum relations that

2T 2 = 2T 4 = T 3 = 4a 2 , (33)
i.e. δv 2 G = T 2 + T 4 -T 3 must vanish in agreement with Eq. (32). This is of relevance for very short sampling times ∆t where f (t) ≈ f (0) = f 0 or if f (t) has an intermediate plateau extending over several orders of magnitude. The summand g 2 ijkl in Eq. ( 31) must remain small, if f (t) is not rigorously, but only nearly constant. The typical summand g 2 can be estimated by the typical slope on logarithmic time scales [15]

g(∆t) ≈ f (∆t) -f (∆t/2) ≈ df (t)/d log(t))| t≈∆t . (34)
One thus expects

δv 2 G [f ] ≈ ijkl g 2 ijkl /I 4 ≈ g(∆t) 2 . ( 35
)
For instance, f (t) may decrease for t τ as f (t) ≈ b exp(-(t/τ ) β ) + f ∞ with constants β > 0. Equations (34) and (35)

lead then to δv G [f ] ≈ |b|(∆t/τ ) β for ∆t τ. (36)
In the opposite limit of very large ∆t τ , the leading scaling dependence is obtained by replacing in Eqs. (66-68) the upper integration bounds by τ and f (t) by a ≈ f (τ ) -f ∞ using Eq. (32). This implies

T 2 ≈ a 2 τ /∆t, T 4 ≈ T 3 ≈ a 2 (τ /∆t) 2 . ( 37
)
In other words, δv 2 G is dominated for ∆t/τ 1 by T 2 = δm 2 21 , i.e. δv G ∝ 1/ √ ∆t as expected for ∆t/τ uncorrelated subintervals. Adding heuristically the short and the long time behavior, Eq. ( 35) and Eq. ( 37), yields the phenomenological approximation [15]

δv 2 G [f ] ≈ g(∆t) 2 + (f (τ ) -f ∞ ) 2 (τ /∆t) (38)
which is useful for processes with one main dominant relaxation process.

δv G [f ]

for test functions f (t)

Introduction

To illustrate some properties of the non-linear functional δv G [f ] we discuss now several test functions f (t). Not all presented f belong to the space of legitimate ACFs c or h of stationary stochastic processes. We remind [START_REF] Hansen | Theory of simple liquids[END_REF][START_REF] Götze | Complex Dynamics of Glass-Forming Liquids: A Mode-Coupling Theory[END_REF] that a legitimate ACF may not change too strongly (especially not discontinuously) and must not violate the Wiener-Khinchin theorem on the power spectrum of the signal stating that the Fourier transform (FT) of c(t) is given by the squared FT of x(t) [START_REF] Hansen | Theory of simple liquids[END_REF][START_REF] Götze | Complex Dynamics of Glass-Forming Liquids: A Mode-Coupling Theory[END_REF][START_REF] Allen | Computer Simulation of Liquids[END_REF]. A general (necessary and sufficient) criterion for a function f (t) to be a legitimate ACF is thus [START_REF] Hansen | Theory of simple liquids[END_REF][START_REF] Götze | Complex Dynamics of Glass-Forming Liquids: A Mode-Coupling Theory[END_REF]

8 f (ω) ≡ ∞ 0 f (t) cos(ωt)dt ≥ 0 for any real ω. ( 39
)
This ensures that f (0) ≥ |f (t)| ≥ 0 and f (ω = 0) = ∞ 0 dtf (t) ≥ 0. Taking advantage of the affine transform Eq. (32) we often set without loss of generality f (0) = 1 and f (t) → 0 for t → ∞. If there is only one characteristic time it is also set to unity. 41) with H2 = 0.5. Also given is f (t) for τ2 = 10000 (solid line with circles). δvG[f ] becomes bimodal with increasing τ2/τ1 with a minimum slightly below τ2 and a second separate maximum at ≈ 5τ2.

Maxwell model

One of the few cases where δv G [f ] can be calculated analytically is the Maxwell model (Debye decay) f (t) = exp(-t). This model is especially of relevance for the selfassembled network systems considered below in Sec. 6.2. Since f (ω) = 1/(1 + ω 2 ) > 0 for all ω, f (t) is a legitimate ACF as expected. Note first that v(∆t) = 1 -g Debye (∆t) with g Debye (∆t) being the Debye function introduced in Sec. 2.2, Eq. ( 18). The three contributions T 2 , T 4 and

T 3 to δv 2 G [f ] = T 2 + T 4 -T 3 are T 2 = 2g Debye (2∆t), T 4 = 2g Debye (∆t) 2 , (40)
T 3 = 4 ∆t 3 × -e -2∆t + (2∆t + 8)e -∆t + 4∆t -7 .
Since g Debye (x) ≈ 2/x for large x we have δv G ≈ 2/∆t for large ∆t. The analytical solution for the Maxwell model is indicated by a bold solid line in Fig. 2. This exact result may be used for testing the numerical determination of δv G [f ] by means of Eqs. (66,67,68).

Two-step relaxation

In view of the presented simulations it is useful to discuss an example for systems with two relaxation processes similar to Fig. 1. Of interest is the limit where f (t) develops an intermediate plateau f (t) ≈ f p for τ 1 t τ 2 with τ 1 corresponding to a fast, local process and τ 2 to a slow, collective relaxation. One expects δv G (∆t) to become bimodal with a first maximum around τ 1 followed by a 1/ √ ∆t-decay and a second maximum around τ 2 followed by a second 1/ √ ∆t-decay. The minimum between both maxima should systematically become deeper with increasing plateau width. Figure 2 presents numerically obtained δv G [f ]-data for

f (t) = H 1 exp(-t/τ 1 ) + H 2 exp(-t/τ 2 ) (41)
with H 1 = τ 1 = 1 and H 2 = 0.5 for the amplitude of the second mode. As for all generalized Maxwell models

f (ω) = pmax p=1 H p τ p 1 + (ωτ p ) 2 > 0, (42)
i.e. Eq. ( 41) is a legitimate ACF. We scan τ 2 over several orders of magnitude as indicated in the figure. We indicate f (t) for the longest second relaxation time, τ 2 = 10000, at the top of the figure (solid line with circles). For large τ 2 /τ 1 one observes for δv G (∆t) two well separated maxima of same shape but different amplitudes ∝ H p . Note that the ratio of the two dashed horizontal lines is H 1 /H 2 = 2. The decay from both maxima is given by δv G ≈ H p 2τ p /∆t.

Stretched and compressed exponentials

Another natural generalization of the one-mode Maxwell model (β = 1) is seen in Fig. 3 where we present δv

G [f ] for f (t) = f β (t) ≡ exp(-t β ). f β (t)
is a "stretched" exponential for β < 1 and a "compressed" exponential for β > 1. It can be readily checked numerically that Eq. ( 39) only holds for β ≤ 2 but not for larger exponents β which do not correspond to ACF of stationary stochastic processes. To see this let us just mention two cases. Since f (ω) ∝ exp(-ω 2 /4) for β = 2, Eq. ( 39) holds for the Gaussian model and it thus also does for even more gently decreasing (less compressed) functions with β < 2. On the other hand f β (t) becomes for β → ∞ equivalent to the cusp singularity be readily calculated analytically and this exact formula is used in Fig. 3.) The cusp singularity is not a legitimate ACF since f = sin(ω)/ω may be negative, i.e. Eq. ( 39) does not hold. As may be seen from the main panel, all δv G [f β ] have a maximum between ∆t ≈ 4 (large β) and ∆t ≈ 10 (small β). As expected from f β (t) → f cusp (t) for β → ∞, it is seen that δv G (β) becomes increasingly similar to the standard deviation of the cusp model (bold solid line), i.e. the peaks become systematically higher, sharper and more lopsided with increasing β. The power-law slopes β (thin solid lines) observed for ∆t 1 are expected from Eq. (36). All models decrease as δv G ≈ a/∆t for large ∆t in agreement with Eq. (37). The amplitude a of this ultimate decay is the largest for the cusp model (a = 4) and the smallest for the Maxwell model (a = 2).

f cusp (t) ≡ H(t)-H(t-1). (δv G [f cusp ] can
The inset of Fig. 3 shows the ratio δv G /v| max taken at the maximum of δv G (∆t) for a broad range of the exponent β. This shows a monotonic increase with β approaching from below the ratio ≈ 1.21 of the cusp model (bold horizontal line). The ratio is ≈ 0.55 for the Maxwell model and ≈ 0.82 for the Gaussian (vertical arrow). Importantly, the standard deviations thus become of the same order as the average behavior for the most rapidly decaying legitimate ACFs with β ≤ 2.

Logarithmic creep

Logarithmically slow ACFs are expected for hopping processes in systems with a broad distribution of barriers and are generally observed in glass-forming fluids [START_REF] Ferry | Viscoelastic properties of polymers[END_REF][START_REF] Hansen | Theory of simple liquids[END_REF][START_REF] Götze | Complex Dynamics of Glass-Forming Liquids: A Mode-Coupling Theory[END_REF]. The general scaling relation Eq. (35) suggests As shown by the squares in Fig. 4 we use

δv G [f ] ≈ 1.55|b| if f (t) ≈ a -b ln(t) (43
f (t) = H 1 e -t/τ1 + a -b ln(t) (1 -e -t/τ2 ) e -t/τ3 (44)
with H 1 = 10 and τ 1 = 1 for the Maxwell model added to mimic the typical microscopic relaxation and τ 2 = 1 and τ 3 = 10 10 setting, respectively, the lower and the upper cutoff of the logarithmic creep. (τ 3 is irrelevant for the presented ∆t-range and the constant a is arbitrary.) The strong Maxwell mode dominates δv G (∆t) below ∆t ≈ 10 3 . Interestingly, as marked by the left arrow deviations from the 1/ √ ∆t-decay (dash-dotted line) expected for the Mawell mode are already observed at ∆t ≈ 10 2 . Only after a broad crossover regime (about three decades) the plateau (solid horizontal line) expected from Eq. ( 43) is reached. The latter model demonstrates how a rather small additional logarithmic creep may lead to strong deviations from an expected 1/ √ ∆t-decay.

Non-ergodic stochastic processes

Our key relation Eq. ( 3) and its various reformulations may obviously fail if one of the stated or implicit assumptions does not apply for the particular ensemble of time series. For instance, strong non-Gaussian contributions may be present in a specific time or frequency range leading to the failure of Wick's theorem, Eq. ( 27). We want in this subsection to address an important assumption not yet explicitly stated. In fact it was assumed that the stochastic process under consideration is ergodic, i.e. all independently created trajectories, called here "configurations", are able to explore given enough time the complete (generalized) phase space. The averages which appear in Wick's theorem, can thus be either obtained by averaging over independent configurations c or by averaging ofter subsets of one extremely long trajectory. To see that this condition matters let us consider a strictly non-ergodic system where the configurations c are trapped in subspaces of the total phase space (since the terminal relaxation time τ of the system diverges).

∆ 2 ne ≡ var(v c ) = 1 N c Nc c=1 v 2 c - 1 N c Nc c=1 v c 2 (45)
being the variance of the N c quenched variances v c = lim ∆t→∞ v[x c ] of the independent configurations. Obviously, ∆ ne vanishes for identical v c . This holds indeed for ergodic systems for ∆t τ (with the finite τ replacing τ b ), but in general not for non-ergodic systems.

On the other hand, for small ∆t the non-ergodicity constraint should not matter much and one expects δv 2 ≈ δv 2 G . Interpolating between both ∆t-limits a useful approximation for non-ergodic systems may be written as

δv 2 (∆t) ≈ δv 2 G (∆t) + ∆ 2 ne for τ b ∆t τ (46)
motivated by the idea that δv 2 is the sum of two variances describing the independent fluctuations within each configuration and between the different configurations. Moreover, Eq. ( 46) suggests the operational definition

δv G (∆t ! = τ ne ) = ∆ ne (47)
identifying τ ne as the crossover time between both limits. Quite generally, τ ne τ b . 9 A rigorous justification of the above interpolation formula Eq. ( 46) will be given elsewhere [30]. We only outline here the general idea. To understand the discrepancy between δv and δv G for (strictly) non-ergodic systems it is necessary to introduce an extended ensemble of time series x ck where for each of the N c independent configurations c one samples N k time series k [31]. 10 Obviously, the time series k of the same configuration c are correlated (being all confined in the same basin) and k-averaged expectation values and variances may then depend on the configuration c. It thus becomes relevant in which order caverages over configurations c and k-averages over time series k of a given configuration c are performed. Three variances of v[x ck ] must be distinguished: the total variance δv 2 tot = δv 2 int +δv 2 ext and its contributions δv 2 int , the typical internal variance within the meta-basins, and δv 2 ext , characterizing the dispersion between the basins. The present paper focuses on the total standard deviation δv tot . (The index tot is dropped elsewhere.) Importantly, if the trajectory of each configuration c remains essentially Gaussian, Wick's theorem can be applied to δv int as before. This implies δv int ≈ δv G . Moreover, since δv ext (∆t) ≈ ∆ ne for τ b ∆t τ , this leads to Eq. ( 46). Variances due to independent physical causes are naturally additive. We remind that the variance of the blackbody radiation is the sum of a variance describing the Rayleigh-Jeans part of the spectrum (wave aspect) and 9 For volume-averaged density fields τ b must be to leading order system-size independent while τ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 ∆ttemp τ b and by sampling then N k time intervals ∆t separated by constant spacer time intervals ∆tspac τ b . N k is assumed to be arbitrarily large and the k-averaged properties δvint and δvext do thus neither depend on N k nor the total sampling time ∆tmax = N k (∆t + ∆tspac). The latter point may become a delicate issue if the non-ergodicity constraint (τ → ∞) is not strictly obeyed. of a variance describing the Wien part (discrete particle aspect) [32]. Interestingly, as in the blackbody radiation analogy, the two contributions δv int (internal basin fluctuations) and δv ext (fluctuations between basins) to δv tot have also different statistics. This is manifested by their different system size dependences as will be shown now.

5 System-size effects Due to the central limit theorem [START_REF] Van Kampen | Stochastic processes in physics and chemistry[END_REF] the stochastic process of many systems is to a good approximation Gaussian since the data entries x i are averages over N m 1 microscopic (often unknown or inaccessible) contributions x im . Specifically, we shall consider below the instantaneous shear stress σi = dr σi,r /V being the volume average over the local shear stress σi,r . For such intensive field averages N m corresponds to the number of local volume elements dV computed, i.e. N m ≈ V /dV . Albeit these microscopic contributions x im may be correlated, i.e. they may not all fluctuate independently, the fluctuations of the x i commonly decrease with increasing N m . Since v ∝ 1/N m for uncorrelated variables x im , it is often useful to incorporate this reference in the definition of the data entries by rescaling x i ⇒ √ N m x i . (This is done in Sec. 6.1 by rescaling the stress by √ V .) For perfectly uncorrelated microscopic variables x im subject to a finite quenched random field this leads to

v ∝ h ∝ δv G ∝ N 0 m and (
v c ≡ x 2 c -x 2 c = 1 N m × 1 N m m v cm . ( 49
)
v is then the ensemble average over all configurations c. Using that also the variances v cm of each microstate are decorrelated we obtain in turn

∆ 2 ne ≡ var[v c ] = 1 N 3 m × 1 N m m var[v cm ] . (50)
Since the m-averages (brackets) in Eq. ( 49) and Eq. ( 50) become N m -independent, this implies ∆ 2 ne ≈ v 2 /N m which in turn confirms γ ext = 1/2. Equation ( 48) also holds for fluctuating density fields with a finite V -independent correlation length ξ for sufficiently large systems (V ξ d ). In this case N m in Eq. ( 48) is simply replaced by the number of independent subvolumes V /ξ d . A smaller exponent γ ext < 1/2 is expected for long-range and scale-free spatial correlations. In agreement with Eq. ( 47) and assuming δv G ∝ 1/∆t β with β ≈ 1/2 we have

τ ne ∝ N γext/β m , (51)
i.e. the crossover time increases with N m and ∆ 1/β ne ∝ 1/τ ne . Details will be given elsewhere [30]. The generally important point is here that ∆ ne decreases and τ ne increases with the system size if γ ext > 0 and thus δv → δv G for N m → ∞ according to Eq. ( 46).

6 Shear-stress fluctuations

Introduction

The theoretical results presented above should be useful for the analysis of general time series x of stochastic processes which are essentially stationary and Gaussian. We illustrate this for the shear-stress fluctuations measured numerically for the coarse-grained model systems sketched in Fig. 5. See Appendix C for further details of the model systems and Appendix B for the definition of the instantaneous shear stress σ and the corresponding instantaneous affine shear modulus μA . The stochastic process x(t) is obtained by rescaling σ(t) ⇒ x(t) ≡ βV σ(t) (52) with β = 1/T being the inverse temperature (setting Boltzmann's constant k B to unity) and V the (two-or threedimensional) volume of the system. With this rescaling v[x], Eq. ( 1), characterizes the empirical shear-stress fluctuations of the time series and the expectation value v(∆t) is equivalent to the "shear-stress fluctuation" µ F (∆t) considered in previous publications on the stress-fluctuation formalism for elastic moduli [4,10,[12][13][14][15]33]. For consistency with the theoretical considerations we keep the general notations defined above. Since these notations differ from the ones widely used for sheared elastic bodies [START_REF] Ferry | Viscoelastic properties of polymers[END_REF][START_REF] Rubinstein | Polymer Physics[END_REF][START_REF] Doi | The Theory of Polymer Dynamics[END_REF][START_REF] Hansen | Theory of simple liquids[END_REF]4,[10][11][12][13]15] the following correspondence list may be useful to the reader:

v(∆t) ↔ µ F (∆t) ≡ µ 0 -µ 1 (∆t) (53)
m 21 ↔ µ 0 ≡ βV σ2 ( 54
)
m 12 (∆t) ↔ µ 1 (∆t) ≡ βV σ2 (55)
R A ↔ µ A = μA (56) R(t) ↔ G(t) = µ A -h(t) (57) M (∆t) ↔ µ sf (∆t) = µ A -µ F (∆t) (58) δv ↔ δµ F ( 59
)
The overbars on the right-hand sides denote the average over a given time series. Note that Eq. ( 57) is the fluctuation dissipation relation [START_REF] Doi | The Theory of Polymer Dynamics[END_REF] for the shear-stress relaxation after an infinitesimal change of the shear strain, i.e., R(t) is the "shear relaxation function", Eq. ( 58) is the corresponding relation for the "generalized shear modulus" [4]. We also remind [4][5][6] that the additional assumption R A = m 21 together with the identities c(0) = m 21 and h(t) = c(0) -c(t) imply that

R(t) = c(t) and M (∆t) = m 12 (∆t). ( 60
)
While R A = m 21 holds indeed under liquid equilibrium conditions, this may become incorrect in general [4][5][6]. In order to test Eq. ( 3) we compare the standard deviation δv(∆t), lumping all N c × N k time series together (cf. Appendix C.5), with δv

G [R] = δv G [h] = δv G [c]
obtained by means of Eqs. (63,64,65) using the measured ACFs.

Self-assembled networks

The TSANET model described in Appendix C.2 is from the rheological point of view very similar to patchy colloids [34,35] or "vitrimers" [36]. Rheologically similar selfassembled transient networks may also be formed by hyperbranched polymer chains with sticky end-groups [37] or microemulsions bridged by telechelic polymers [38][39][40].

As shown in Figs. 6 and 7 of Ref. [10], TSANET is a simple Maxwell fluid, i.e. the shear-stress relaxation modulus, computed by means of Eq. (57) or Eq. ( 60), decays exponentially

R(t) ≈ R M exp(-t/τ M (ν)) for t 1 ( 61
)
with R M ≈ 18 being the plateau modulus set by the equilibrium shear modulus for permanent springs (ν = 0) and τ M (ν) ≈ 16/ν the Maxwell time. One may use as scaling variable the reduced sampling time ∆x = ∆t/τ M for 1 ∆t ≤ ∆t max to collapse data obtained for different ∆t and ν [10]. The short-time decay of R(t) from the initially imposed affine stain, R(t = 0) = R A , to the plateau modulus R M is reasonably described by a compressed exponential (β > 1). A useful formula for R(t) for all t and ν is given by the two-mode approximation

R(t) = b exp(-(t/τ A ) β ) + R M exp(-t/τ M ) (62)
with amplitude b ≈ R A ≈ 33, relaxation time τ A ≈ 0.3 and exponent β ≈ 1.5. Since the hopping moves changing the network connectivity obey detailed balance, changing ν leaves all truly static properties unchanged as may be seen from Fig. 4 of Ref. [10]. For this reason all simple averages, Eq. ( 7), such as the average shear stress σ, the average affine shear modulus R A or the moment m 21 , should not depend on ∆t or ∆x. That this is indeed the case can be seen from Fig. 6. Note that σ = 0 by symmetry. R A and m 21 are roughly equal, R A ≈ m 21 ≈ 33, albeit m 21 fluctuates more strongly for small ν. 11 Also presented in Fig. 6 are m 12 , v and M which are all seen to depend on the reduced sampling time. Note that M ≈ m 12 in agreement with Eq. (60). As expected for a Maxwell model according to Eq. ( 18), M (∆x) = R M g Debye (∆x) holds (bold line). The corresponding relation v(∆x) = R A -R M g Debye (∆x) is indicated by the dash-dotted line. Note that g Debye (x) → 1 for x → 0 and g Debye (x) → 2/x for x 1. This implies that in the liquid limit M (∆x) ≈ 2R M /∆x and v(∆x) ≈ R A -2R M /∆x. Importantly, the ∆x-dependence of these properties is not due to aging or equilibration problems but is caused by the finite time needed for the equilibrium fluctuations to explore the phase space reflecting the stress relaxation process R(t).

We turn now to the characterization of the standard deviation δv. As shown in Fig. 7 δv ≈ δv G is found to hold to high precision for all high hopping frequencies with 1 ∆t max /τ M (ν) ∝ ν. For ν = 0.1 both the short-time relaxation time τ A , characterizing the relaxation of the affine strain, and the Maxwell time τ M (ν) = 16/ν, characterizing the reorganization of the network, are relevant. We have also included the shear-stress relaxation function R(t) which is seen to vanish above τ M . Note that v has a shoulder with v ≈ R A -R M for ∆t τ M and that v ≈ R A for larger ∆t, as expected in the liquid limit. In agreement with Sec. 3.5 we observe for this two-step relaxation process two well separated δv-maxima (cf. Fig. 2). The first relaxation process at τ A is described by a compressed exponential with β = 1.5 (bold dashed line), the second relaxation process due to the reorganization of the spring network for larger ∆t by a Maxwell model (β = 1) with relaxation time τ M (bold solid line). In fact, for all not too small ν δv is given by δv G [R] assuming Eq. (62) to hold.

As shown in Fig. 8 this becomes different for small ν due to quenched shear-stress fluctuations. We present here data obtained for a quenched network with switched off hopping moves (ν = 0). Also indicated is the shearstress fluctuation v (diamonds) which is seen to rapidly increase for small ∆t, corresponding to the relaxation of the imposed affine strain, and to level off for ∆t τ A as indicated by the horizontal dashed line. Since in agreement with Eq. ( 13) h(t) or R(t) become also constant in this time regime (not shown) this implies that δv G must decay as 1/ √ ∆t for ∆t τ A . This is confirmed by the δv G [h]-data (triangles) computed from the measured h(t), revealing after a first regime with δv G ∝ ∆t β and β ≈ 1.5 the expected 1/ √ ∆t-decay. As shown by the bold solid line, a reasonably fit of δv G for all ∆t is obtained using Eq. ( 62). (For ν = 0 the second term in Eq. ( 62) is an irrelevant constant.) While δv is identical (within numerical precision) to δv G for short ∆t it deviates for large ∆t where it levels off, δv → ∆ ne ≈ 0.16, as indicated by the bold dashed horizontal line. As discussed in Sec. 4 the leveling-off is expected for a finite dispersion of the v c . The interpolation formula Eq. ( 46) motivated in Sec. 4 gives a reasonable approximation of δv (stars) matching both limits for ∆t τ A . To leading order, δv is thus given by δv G and, hence, by h(t) or R(t) plus an additional constant. As indicated by the tilted arrow in Fig. 8, Eq. ( 46) slightly overpredicts δv for intermediate ∆t. This suggests that the constant ∆ ne should be replaced by the more general standard deviation δv ext (∆t), describing the ∆t-depending dispersion between configurations, approaching monotonically the large-∆t limit ∆ ne with increasing ∆t from below. This minor difference will be discussed elsewhere [30]. Similar and complementary results for glass-forming polydisperse particles and thin free-standing polymer films are reported in Appendix D.

System-size effects for δv and ∆ ne

We have focused up to now on the variation of the sampling time ∆t, the hopping frequency ν or the temperature T (cf. Appendix D) while keeping fixed other parameters such as the total number of beads n. While most properties discussed above as R A , v or h are defined as intensive properties, i.e. as we have checked their mean values do not or extremely weakly depend on n, this is less obvious for their respective standard deviations [10,41,15]. We address here briefly the n-dependence of the standard deviation δv. All presented systems have roughly the same number density ρ of order unity, i.e. n ≈ V .

δv is presented in Fig. 9 for a broad range of n for the TSANET model and the pLJ particles (cf. Appendicies C.3 and D.1). The TSANET data in panel (a) are plotted as a function of the reduced hopping frequency ∆x ∝ ν for our largest sampling time ∆t = ∆t max = 10 5 . (The short-time behavior around τ A is thus irrelevant.) The dashed horizontal lines indicate ∆ ne obtained as in Fig. 8 from the large ∆t-limit of δv for quenched networks. The bold solid line represents δv G for the one-mode Maxwell model. Not shown for clarity are the δv G [h] obtained for the different n which are found to be essentially n-independent and very similar to the Maxwell model. At variance to this δv is only intensive for sufficiently large ∆x where δv ≈ δv G holds, but not in the low-∆x limit where δv → ∆ ne . Panel (b) presents δv as a function of ∆t for the pLJ particles. Also given are δv G [h]-data for n = 10000 (crosses). The dashed horizontal lines indicate the plateau value ∆ ne for each n. ∆ ne systematically decreases with n. δv and δv G thus become increasingly similar according to Eq. (46). We also note that a scaling collapse of δv is achieved for both models by plotting δv/∆ ne (n) as a function of ∆t/τ ne (n) with τ ne (n) determined according to Eq. (47).

Figure 10 summarizes the system-size dependence of ∆ ne for several models. We compare published results [41,15] with new data obtained for the TSANET model (ν = 0), the pLJ particles (T = 0.2) and the free-standing polymer films (T = 0.05). The TSANET data (circles) are fitted by ∆ ne ∝ 1/n γext with γ ext = 1/2 (bold solid line). Fig. 10. ∆ne vs. n for a binary LJ mixture at T = 0.05 (stars) [41], for a 3D polymer glass at T ≈ 0.1 (triangles) [15], for polymer films at T = 0.05 (diamonds), for the pLJ model at T = 0.2 (squares) and for the TSANET model at ν = 0 (circles). ∆ne decreases with a power-law exponent γext ≈ 1/3 for the amorphous glasses (dashed lines) and γext = 1/2 for the TSANET model (bold solid line).

The observed "strong self-averaging" [START_REF] Landau | A Guide to Monte Carlo Simulations in Statistical Physics[END_REF] suggests that independent localized shear-stress fluctuations with a finite correlation length ξ are responsible for ∆ ne in agreement with Eq. ( 48). This finding is at variance to the somewhat smaller exponent γ ext ≈ 1/3 suggested by recent simulation studies of 2D binary LJ mixtures (stars) [41] 12 and of dense 3D polymer glasses (triangles) [15]. This led us to conclude [15] that local elastic (structural) properties in amorphous systems may cause long-range spatial correlations and a diverging correlation length ξ. As shown by the lower dashed line, the exponent γ ext ≈ 1/3 is also compatible with the new data obtained for our extremely well equilibrated pLJ particles. Unfortunately, only two system sizes have been probed for the 3D polymer system [15]. This makes it difficult to assess whether γ ext depends on the spatial dimension or not. Simulations with a broader range of n are currently sampled to corroborate this point and to verify γ ext in two and three dimensions.

Conclusion

We have discussed systematically the ensemble average v(∆t) and the standard deviation δv(∆t) of the variance v[x], Eq. ( 1), of a time series x measured over a sampling time ∆t. Our aim was to give an uncluttered overview of some relations which may be useful in different fields where the stochastic processes are essentially, albeit perhaps not rigorously, both stationary and Gaussian. We have emphasized first in Sec. 2 that for stationary processes v is given by a weighted sum (integral), Eq. ( 2), over the ACF h(t) = c(0) -c(t) (Sec. 2.2). Assuming an ergodic Gaussian process (Sec. 3) δv was shown in Sec. 3.3 to be given by the functional δv G [h], Eq. ( 3). As discussed in Sec. 3.5 the reduced standard deviation δv G /v taken at the maximum of δv G may become of order unity if h(t) changes rapidly (Fig. 3), i.e. the average behavior v gets masked by strong fluctuations. As emphasized in Sec. 4 Eq. ( 3) cannot hold for non-ergodic systems with a finite dispersion of the frozen variances v c of the different independent configurations c since δv must become constant, ∆ ne , for ∆t τ ne while δv G vanishes. However, if the observable x is the sum of many more or less decoupled microscopic variables the quenched v c become similar with increasing system size and, hence, ∆ ne → 0, τ ne → ∞ and δv → δv G even for non-ergodic systems (Sec. 5).

In the computational part of this work (Sec. 6) we have illustrated some of the relations by applying them to the shear-stress fluctuations (Sec. 6.1) in transient selfassembled networks (Sec. 6.2). Similar results are reported in Appendix D for glass-forming polydisperse particles and free-standing polymer films. Albeit it is non-trivial [START_REF] Hansen | Theory of simple liquids[END_REF]42,43,15] whether the shear-stress trajectories of these systems are sufficiently stationary and Gaussian, all examples reveal qualitatively the same behavior:

all systems are (at least effectively) stationary as shown in panel (c) of Fig. 18; δh(t) 2 ≈ 2h(t) 2 holds to high precision (Fig. 14) as expected for Gaussian processes; δv ≈ δv G holds within numerical precision in the ergodic limit for large hopping frequencies ν (Fig. 7) or temperatures T (Figs. 15 and19b); while δv G ∝ 1/ √ ∆t vanishes in the non-ergodic limit, δv becomes constant, δv → ∆ ne > 0, for large ∆t (Figs. 8, 16 and 19a); -∆ ne decreases with system size (Figs. 9 and 10) suggesting that Eq. ( 3) becomes valid for macroscopic albeit non-ergodic elastic bodies.

The presented work was limited to the characterization of fluctuations and relaxation processes of stationary (equilibrium) stochastic processes, i.e. no external perturbation was applied to directly measure the average response function R(t) or the average modulus M and their, respective, standard deviations δR(t) and δM . Our claim that δR/R or δM/M must generally become large (of order unity) for times where R(t) strongly decays and that these ratios are, moreover, system-size independent may in fact be misleading for the out-of-equilibrium responses of real macroscopic materials. From the theoretical point of view it is an interesting question how to generalize the fluctuation-dissipation relations, connecting the average linear out-of-equilibrium response to the average equilibrium relaxation [START_REF] Doi | The Theory of Polymer Dynamics[END_REF][START_REF] Hansen | Theory of simple liquids[END_REF]10], to describe the sampleto-sample fluctuations.

We have briefly discussed in Sec. 6.3 the system-size effects for various properties focusing on δv and the associated finite plateau value ∆ ne for non-ergodic systems. As shown in Fig. 10 two different exponents γ ext characterize the decay of ∆ ne with n for the perfectly equilibrated TSANET model (γ ext = 1/2) and the quenched amorphous glasses (γ ext ≈ 1/3). The clarification of this difference is beyond the scope of the present paper. We shall also give elsewhere [30] a systematic description of the two different types of standard deviations briefly mentioned in Sec. 4 which must be distinguished for the complete characterization of fluctuations of ensembles of non-ergodic systems. This will allow to further discuss the surprising "breakdown of nonlinear elasticity in amorphous solids" [41] -based on the numerically observed divergence with system size for standard deviations associated with higherorder nonlinear analogs of the elastic shear modulusclaimed at variance to the every day experience that sufficiently large amorphous (plastic) bodies are well behaved according to standard continuum mechanics [START_REF] Ferry | Viscoelastic properties of polymers[END_REF][START_REF] Tadmor | Continuum Mechanics and Thermodynamics[END_REF].

A Reformulations of δv G [f ]

Since for large I the sums over two, three or even four indices stated in Sec. 3.3 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

T 2 (∆t) = 2 I 2 If 2 0 + 2 I-1 k=1 (I -k)f 2 k (63) T 4 (∆t) = 2 I 4 If 0 + 2 I-1 k=1 (I -k)f k 2 (64)
with f k = f (t = t k ) for an arbitrary function f (t). Let us define the sum S(s, I) ≡ I i=1 f i-s . Note that S(s, I) 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

T 3 (∆t) = 4 I 3 I s=1 S(s, I) 2 = 8 I 3 I/2 s=1 S(s, I) 2 (65)
becomes also of order 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 = ∆t/δt the three terms further simplify to

T 2 (∆t) = 4 ∆t 2 ∆t 0 dt (∆t -t)f (t) 2 (66)
T 4 (∆t) = 8 ∆t 4 ∆t 0 dt (∆t -t)f (t) 2 (67)
T 3 (∆t) = 8 ∆t 3 ∆t/2 0 dt (η(t) + η(∆t -t)) 2 (68)
using η(t) ≡ t 0 du f (u) for the last contribution.

B Shear stress and Born-Lamé coefficient

Let us consider a small simple shear strain [START_REF] Tadmor | Continuum Mechanics and Thermodynamics[END_REF] γ in the xy-plane as it would be used to measure the shear-stress relaxation function R(t) [START_REF] Allen | Computer Simulation of Liquids[END_REF][5][6][7]. Assuming that all particle positions r follow an imposed "macroscopic" shear in an affine manner according to r x → r x + γ r y the Hamiltonian Ĥ of a given configuration changes to leading order as [6, 10,

) -Ĥ(γ = 0))/V ≈ σγ + 1 2 μA γ 2 for |γ| 1. 11, 13] ( Ĥ(γ
The instantaneous shear stress σ and the instantaneous Born-Lamé coefficient μA are thus defined as σ ≡ Ĥ (γ)/V | γ=0 and (70)

μA ≡ Ĥ (γ)/V | γ=0 = σ (γ)| γ=0 (71)
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. 13 Assuming a pairwise central conservative potential l u(r l ) with r l being the distance between a pair of monomers l, one obtains the excess contributions [6,14]

σ = 1 V l r l u (r l ) n l,x n l,y and (72) μA
= 1 V l r 2 l u (r l ) -r l u (r l ) n 2 l,x n 2 l,y + 1 V l r l u (r l ) (n 2 l,x + n 2 l,y )/2 (73)
with n l = r l /r l being the normalized distance vector. Note that Eq. ( 72) is strictly identical to the corresponding off-diagonal term of the Kirkwood stress tensor [START_REF] Allen | Computer Simulation of Liquids[END_REF]. We have used a symmetric representation for the last term of Eq. ( 73) 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 μA in terms of the pair correlation functions of the bonded and the non-bonded interactions of the particles. Please note that μA depends on the second derivative u (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 (r) of the potential is not continuous. Unfortunately, this is the case at the cutoff r cut of the LJ potentials used for the pLJ beads (Appendix D.1) and for the polymer films (Appendix D.2). The "bare" μA is thus roughly about 0.2 too high for both models and must be corrected [3]. The ensemble average μA is called µ A in previous publications [START_REF] Schnell | [END_REF][3][4][5][6][7][8][9][10][11][12][13][14][15] and R A in the present work.

C Numerical models and technical details C.1 Introduction

Various issues discussed theoretically in Secs. 2-5 are illustrated in Sec. 6 and in Appendix D for the fluctuations of shear stresses in simple coarse-grained model systems. A sketch of the three14 models studied is given in Fig. 5.

Standard MD and MC methods [START_REF] Allen | Computer Simulation of Liquids[END_REF][START_REF] Landau | A Guide to Monte Carlo Simulations in Statistical Physics[END_REF] 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 B , the typical size of the particles (beads) and the particle mass are all set to unity and Lennard-Jones units [START_REF] Allen | Computer Simulation of Liquids[END_REF] are used throughout this work. For all systems we use periodic boundary conditions [START_REF] Allen | Computer Simulation of Liquids[END_REF][START_REF] Landau | A Guide to Monte Carlo Simulations in Statistical Physics[END_REF] and the number density ρ = n/V is close to unity. The salient features of each model and some algorithmic details are given below.

C.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" [44,45] 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 ν. 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 40000 springs and n = 10000 beads, i.e. ρ = n/V = 1. 15 We also report in Sec. 6.3 on data for the same number densities of particles and springs obtained for n = V = L 2 = 100, 300, 1000, 3000, 30000 and 100000. In addition to the MC moves, changing the connectivity matrix of the network, standard MD simulation with a strong Langevin thermostat [START_REF] Allen | Computer Simulation of Liquids[END_REF] is used to move the particles effectively by overdamped motion. 16 We systematically vary ν using a constant total sampling time ∆t max = 10 5 . For each trajectory we store every δt = 0.01 instantaneous properties. Our configurations are first equilibrated in the liquid regime (ν = 1). We then gradually decrease ν over several orders of magnitude. The hopping moves are switched off (ν = 0) for quenched networks. All properties are averaged over N c = 100 configurations.

C.3 Polydisperse Lennard-Jones particles

We present in Appendix D.1 data obtained for polydisperse Lennard-Jones (pLJ) particles in d = 2 dimensions [START_REF] Allen | Computer Simulation of Liquids[END_REF]46,47,3,4] where the interaction range σ = (D i + D j )/2 is set by the Lorentz rule [START_REF] Hansen | Theory of simple liquids[END_REF] with D i and D j being the diameters of the interacting particles i and j. Following Refs. [46,47,3,4,48] the particle diameters are uniformly distributed between 0.8 and 1.2. The standard particle number used is n = 10000. Additional results obtained for n = 100, 200, 500, 1000, 2000 and 5000 are given in Sec. 6.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 ρ(T ) of the older data from Ref. [4] reveal two distinct linear slopes which have been used to define a glass transition temperature T g ≈ 0.26. For details see Ref. [4]. We present here new data [48] where we have used in addition MC swap moves [49] 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 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 As can be seen, much higher densities are achieved at low

C.4 Free-standing polymer films

As sketched in panel (c) of Fig. 5, we study by means of MD simulation [START_REF] Allen | Computer Simulation of Liquids[END_REF] of a bead-spring model [50] freestanding polymer films suspended parallel to the (x, y)plane [14]. All unconnected monomers interact with a truncated and shifted LJ potential while connected monomers are bonded by harmonic springs [50,2, 51, 11-13, 15, 14].

The films contain M = 768 monodisperse chains of length N = 16, i.e. in total n = 12288 monomers, in a periodic box of lateral box size L = 23.5. Entanglement effects are irrelevant for the short chains [START_REF] Ferry | Viscoelastic properties of polymers[END_REF][START_REF] Doi | The Theory of Polymer Dynamics[END_REF][START_REF] Rubinstein | Polymer Physics[END_REF][START_REF] Graessley | Polymeric Liquids & Networks: Dynamics and Rheology[END_REF]. As discussed in detail in Ref. [14], ensembles of N c = 100 independent configurations are generated for a broad range of temperatures. Production runs are performed over ∆t max = 10 5 storing data each δ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 ρ 0 = ρ(z ≈ 0) the film thickness is defined as H ≡ N M/ρ 0 L 2 and the film volume as V = HL 2 . As may be seen from the main panel, H decreases monotonically upon cooling with two linear branches fitting the glass (dashed line) and the liquid (solid line) limits. T g ≈ 0.37 is estimated from the intercept of both asymptotes.

C.5 Data handling

As indicated above we equilibrate for each state point of the considered model N c = 100 independent configurations c. This allows to probe all properties accurately. For each configuration c we compute and store one long trajectory with ∆t max /δt ≈ 10 7 data entries. Since we want to investigate the dependence of various properties on the sampling time ∆t we probe for each ∆t max -trajectory N k equally spaced subintervals k of length ∆t ≤ ∆t max with I = ∆t/δ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 k = ∆t max /∆t, i.e. N k and ∆t are thus coupled and the accuracy is better for small ∆t.

D Shear-stress fluctuations continued D.1 Polydisperse Lennard-Jones particles

We turn now to the shear-stress fluctuations of the pLJ model. As mentioned in Appendix C.3, our configurations have been first tempered and annealed at a constant pressure P = 2 using in addition to the standard local MC moves [START_REF] Landau | A Guide to Monte Carlo Simulations in Statistical Physics[END_REF]4] swap moves exchanging the particle diameters [49]. 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 ∆t max = 10 7 MC steps, we indicate here for the latter case also some temperatures sampled with ∆t max = 10 8 .) As a quick way to characterize the terminal relaxation time τ α (T ) one may set, say, R(t ≈ τ α ) ≈ 0.1. This implies for the swap moves that τ α 10 5 for all temperatures above T frac . For instance, τ α (T = 0.2) ≈ 10 4 is three orders of magnitudes smaller than the production time ∆t max = 10 7 . This suggests that our pLJ samples are well equilibrated and one expects R A ≈ m 21 above T frac . As may be seen from the open symbols, this equilibration does not change much the temperature T g where the relaxation time τ α exceeds the production time ∆t max = 10 7 if the swap moves are switched off. As emphasized by the double-headed arrow, a similar value T g ≈ 0.26 is found as in our previous study [4].

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

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 ≈ 0.25 and T ≈ 0.3.) The standard deviation δh of h was calculated using Eq. [START_REF] Tadmor | Modeling Materials[END_REF]. As discussed in Sec. 3.2, 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 α 2 (t) = δh 2 /2h 2 -1 demonstrates that |α 2 (t)| 0.02 for all times and temperatures (not shown). Moreover, since h ≈ v for all T and h = v = R A above T ≈ 0.3, this implies the observed non-monotonic behavior of δh(T ). We have also checked that δR ≈ δh (not shown) as one expects [10,15] since R A barely fluctuates. These findings confirm that Gaussian processes are dominant for all T , even below the demixing transition at T frac .

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

Focusing on one low temperature (T = 0.2) Fig. 16 compares v with δv and δv G as functions of ∆t. In agreement with Eq. ( 11), v(∆t) increases monotonically becoming rapidly constant. Since consistently with Eq. ( 13) Fig. 15. δv (open symbols) and δvG (lines) as functions of T for a broad range of ∆t. While δv ≈ δvG holds for large T and small ∆t, δv is seen to become constant for low T and high ∆t. An increasingly sharp maximum appears around Tg shifting with increasing ∆t to lower T . also h(t) and R(t) become constant (Fig. 13), this implies δv G (∆t) ∝ 1/ √ ∆t (bold solid line). At variance to the decay of δv G , δv → ∆ ne ≈ 0.25 for ∆t τ ne (dashed horizontal line) similar to the behavior for quenched TSANET systems (Fig. 8). As in Sec. 6.2 Eq. ( 46) yields a reasonable approximation of δv albeit the shifted data (stars) are slightly above δv at ∆t ≈ τ ne .

D.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. 18 Albeit we average over N 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 → ∆t max = 10 5 . Interestingly, it is clearly seen that R(t) increases continuously with decreasing T [14]. in perfect agreement with all published experimental [52][53][54] and computational [55,56] studies. The presented R(t)-data is used below to obtain δv G [R]. 54), implies due to the frozen stresses the observed y(T ) ∝ 1/T for T T g . (Similar behavior has been reported for three-dimensional polymer bulks [13].) Since there is currently no algorithm comparable to the swap algorithm [49] allowing to equilibrate glass-forming polymer melts and films as for the pLJ model, our films are clearly not at thermal equilibrium below T 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 g (where R A < m 21 ) that Eq. ( 11) and Eq. ( 16) hold for v (dashed lines) and M (solid line). To test these relations R(t) was integrated numerically. 19 Albeit time-translational invariance apparently holds, this does not mean that no aging occurs (for the Laplace transforms of J(t) and R(t) are reciprocally related [START_REF] Ferry | Viscoelastic properties of polymers[END_REF], if one function is precisely known, the other can be calculated in principle. 19 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. The dashed lines have been obtained using Eq. ( 11), the solid line using Eq. ( 16) taking advantage of the directly measured shear-stress relaxation function R(t).

T < T g ) but just that this is irrelevant for the timescales considered here.

We address in panel (a) of Fig. 19 the ∆t-dependence of δv and δv G (open symbols) for different temperatures as indicated. δv ≈ δv G holds again for T above and around T g . Note also that δv ≈ δv G ∝ 1/ √ ∆t (bold solid line) for all ∆t ≤ ∆t max for the largest temperature T = 0.55. Interestingly, T -dependent shoulders appear for T ≈ 0.4. This is a consequence of the creep-like decay of R(t) in this regime which is approximately fitted by R(t) ≈ a -b ln(t) as indicated in Fig. 17. According to Fig. 4, one expects a shoulder with δv G ≈ 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 δv becomes constant, δv ≈ ∆ ne for ∆t τ ne ≈ 800, while δv G continues to decrease with ∆t. The small deviations of δv G from the 1/ √ ∆t-asymptote are caused by the fact that R(t) does not become rigorously constant, R(t) → R p > 0, even for our lowest temperatures. As shown by the upper thin line in Fig. 17 A complementary representation is given in penal (b) of Fig. 19 focusing on the temperature dependence of δv and δv G . A very similar behavior as in Fig. 15 for pLJ particles is seen. Most importantly, δv ≈ δv G for large and intermediate T . The maximum of δv ≈ δv G slightly below T g becomes systematically larger with increasing ∆t and shifts to lower temperatures. Also it is confirmed that Eq. (3) holds for sufficiently small ∆t for all T . As expected, Eq. (3) breaks down for ∆t τ ne (T ) ≈ 800 in the non-ergodic limit (lowest temperatures). While δv G , measuring the fluctuations within each configuration, is seen to decrease upon cooling, δv becomes constant, δv → ∆ ne (dashed horizontal line), due to the finite dispersion of the quenched v c of the independent configurations.

  [x] ≡ ( I i=1 x α i /I) β . As stated in the Introduction we focus in this work on the variance v[x] = m 21 [x] -m 12 [x]. Note that v[x] = 0 for I = 1. It is also useful to consider functionals with a discrete time lag s (with s = 0, . . . , I -1) such as the "gliding average" [25] c g s [x] ≡ 1 I -s I-s i=1 c s,i with c s,i = x i+s x i

Fig. 1 .

 1 Fig. 1. Illustration of properties discussed in Sec. 2.4 using a two-modes version of Eq. (17) with amplitudes R∞ = 1, H1 = 6, H2 = 3, i.e. R(0) = RA = R∞ + H1 + H2 = 10, and relaxation times τ1 = 1 and τ2 = 100000. h(t) and c(t) are given by Eq. (14), where we have additionally set c∞ = 0, and M (∆t) and v(∆t) by means of Eqs. (16) and (18). The two solid horizontal lines mark the intermediate pseudo-plateau for τ1 t τ2. v(∆t) and M (∆t) are seen to converge much more slowly to the respective plateau values than the corresponding response functions h(t) and R(t).

Fig. 2 .

 2 Fig. 2. δvG[f ] vs. ∆t for the one-mode Maxwell model with H1 = τ1 = 1 (bold solid line) revealing a maximum at ∆t ≈ 5 and a final decay δvG ≈ 2/∆t. The other data refer to the two-step relaxation model Eq. (41) with H2 = 0.5. Also given is f (t) for τ2 = 10000 (solid line with circles). δvG[f ] becomes bimodal with increasing τ2/τ1 with a minimum slightly below τ2 and a second separate maximum at ≈ 5τ2.

Fig. 3 .

 3 Fig. 3. δvG[f ] vs. ∆t for f (t) = f β (t) ≡ exp(-t β ) with β = 1 corresponding to the one-mode Maxwell model, β = 2 to a Gaussian and β → ∞ to the cusp model fcusp(t) ≡ H(t) -H(t-1). Only exponents β ≤ 2 correspond to legitimate ACFs. Note that δvG ∝ ∆t β for ∆t 1 (thin solid lines for β = 0.3, 0.5 and 1) and δvG ∝ 1/ √ ∆t for ∆t 1. Inset: δvG/v| max vs. β. The vertical arrow marks the ratio ≈ 0.82 for β = 2, the horizontal line the ratio ≈ 1.21 for β → ∞.

Fig. 4 .

 4 Fig. 4. δvG[f ] vs. ∆t for logarithmic creep. The circles indicate Eq. (43) for b = 1, the squares Eq. (44) for H1 = 10, τ1 = τ2 = 1, b = 0.1 and τ3 = 10 10 and the dash-dotted line the 1/ √ ∆tdecay of the short-time Maxwell model. The arrows mark the broad crossover between the Maxwell model and the plateau expected for b = 0.1.

2 c

 2 48) ∆ ne ≡ lim ∆t→∞ δv(∆t) ∝ N -γext m with γ ext = 1/2. Due to Eq. (3) the N m -independence of δv G is implied by the N m -independence of h. To see that γ ext = 1/2 it is sufficient to compute (∆t-independent) ensemble averages . . . c for the different quenched meta-basins c compatible with the non-ergodicity constraint (τ → ∞). (Note that . . . c is obtained for N k → ∞.) Substituting x = m x m /N m , setting v cm = x 2 m c -x m and using that the microstates m are decorrelated yields

Fig. 5 .

 5 Fig. 5. Model systems considered computationally: (a) Transient self-assembled elastic networks (TSANET) created by reversibly breaking and recombining springs with an attempt frequency ν per spring. The spring s thus connects the beads i and j on the left and the beads i and k on the right. Results presented in Sec. 6.2. (b) Monte Carlo (MC) simulations of polydisperse Lennard-Jones (pLJ) particles with larger beads being red, smaller beads darker. Results presented in Appendix D.1. (c) Molecular dynamics (MD) simulations of thin free-standing films of glass-forming polymers. Results presented in Appendix D.2.

Fig. 6 .

 6 Fig. 6. Various properties obtained for the TSANET model vs. the reduced sampling time ∆x = ∆t/τM for ∆t = 10 5 : mean shear stress σ, affine shear modulus RA, contributions m21 and m12 to the shear-stress fluctuation v = m21 -m12 and stress-fluctuation formula M = RA -v for the shear modulus. The prediction for a Maxwell model is indicated by the bold solid line for M and by the dashed-dotted line for v. Note that RA ≈ m21 and M ≈ m12 for all ∆x.

Fig. 7 .

 7 Fig. 7. Comparison of δv(∆t) and δvG(∆t) for a high hopping frequency ν = 0.1 confirming that δv ≈ δvG for all ∆t. Also included are the shear-stress relaxation function R(t) and the shear-stress fluctuation v(∆t). Note that v(∆t) has a shoulder below τM ≈ 160 (lower dashed line) and a plateau with v ≈ RA for ∆t τM (upper dashed line). The thin solid line on the left indicates the power law δvG ∝ ∆t β with β = 1.5 expected from Eq. (36) and Eq. (62), the thin solid line on the right the final 1/ √ ∆t-decay.

Fig. 8 .

 8 Fig. 8. Shear-stress fluctuation v and the standard deviations δv and δvG vs. ∆t for ν = 0. The bold solid line corresponds to δvG obtained assuming a compressed exponential with β = 1.5. While δv ≈ δvG holds for ∆t τne ≈ 4000, δv → ∆ne for larger ∆t. The stars indicate Eq. (46) using ∆ne = 0.16.

Fig. 9 .

 9 Fig. 9. δv for different n as indicated: (a) δv(∆x) for TSANET networks with constant ∆t = 10 5 . The bold solid line indicates δvG for the Maxwell model, Eq. (40), the dashed horizontal lines ∆ne obtained for quenched networks (ν = 0). (b) δv(∆t) for pLJ particles at T = 0.2. δvG obtained from R(t) for n = 10 4 is indicated by crosses, the asymptotic power-law slope 120/ √ ∆t by the thin solid line. The deviations from Eq. (3) vanish with increasing n.

7

 7 MC steps only with local MC moves. Production runs only using local MC moves are then performed over ∆t max = 10 7 . For temperatures between T = 0.21 and T = 0.25 additional production runs have been done for N c = 20 over ∆t max = 10 8 MC steps. The data are normally sampled in intervals of δt = 10 MC steps. The number densities obtained from these new simulations are indicated in Fig. 11 (circles).

Fig. 11 .

 11 Fig. 11. Number density ρ of the pLJ model at normal pressure P = 2 as a function of temperature T . The old data (triangles) from Ref.[4], obtained from a continuous cooling process (with rate 10 -7 ) using local MC moves only, reveal two distinct linear slopes which were used to determine a glass transition temperature Tg ≈ 0.26. Using in addition swap moves much higher densities have been achieved (circles).

Fig. 12 .

 12 Fig. 12. Film thickness H and glass transition temperature Tg of free-standing polymer films. Inset: Number density profile ρ(z) for T = 0.5 with z = 0 corresponding to the center of mass of each film. Using the measured midplane density ρ0(T ) (horizontal line) yields H ≡ N M/ρ0L 2[14]. Main panel: H as a function of T . Tg ≈ 0.37 is obtained from the intercept of the linear extrapolations of the glass (dashed line) and liquid (solid line) limits[14].

Fig. 13 .

 13 Fig. 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) ≈ 35 exp(-(t/τα) 1/3 ) with τα ∝ 1/T 5/2 fitting the swapped systems for T > T frac . As shown by the double-headed arrow, τα(T ), obtained with local moves only, becomes similar to ∆tmax = 10 7 slightly below Tg ≈ 0.26.

Fig. 14 .

 14 Fig. 14. Various properties obtained for the pLJ model taken at t = ∆t = 10 5 vs. T : affine shear modulus RA, moment m21, shear-stress fluctuation v and its standard deviation δv, meansquare displacement h and its standard deviation δh and shearstress relaxation function R = RA -h. Note that RA ≈ m21 and h ≈ v for all T and that δh ≈ 2 1/2 h (bold solid line) holds in agreement with Eq. (25).

Fig. 16 .

 16 Fig. 16. Shear-stress fluctuation δv and the corresponding standard deviations δv and δvG vs. ∆t at T = 0.2 which is well below Tg. v(∆t) is seen to become constant above ∆t ≈ 10 3 where v ≈ h(t) ≈ 17.1. As indicated by the bold solid line, δvG ∝ 1/ √ ∆t as expected in this regime. δv levels off for ∆t τne ≈ 200000.

Fig. 17 .

 17 Fig. 17. R(t) = RA -h(t) for free-standing films using halflogarithmic coordinates for a broad range of temperatures T . R(t) increases continuously with decreasing T . Logarithmic creep behavior (thin solid lines) with R(t) ≈ a-b ln(t) is found for T Tg and above the glass transition (T ≈ 0.45).

Figure 18

 18 presents various contributions to the generalized shear modulus M = R A -v = (R 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 ≈ m 12 for temperatures above the glass transition where R A ≈ m 21 . Note that R A < m 21 below T g . As may be seen in panel (b) frozen-in out-of-equilibrium stresses are observed upon cooling below T g as made manifest by the dramatic increase of the dimensionless parameter y = m 21 /R A -1. The prefactor β = 1/T of m 21 , Eq. (

Fig. 18 .

 18 Fig. 18. Contributions to M = RA -v: (a) Temperature dependence for ∆t = 10 4 . (b) Double-logarithmic representation of m21/RA -1 vs. T . (c) ∆t-effects for M and its contributions for T = 0.3. Only RA and m21 are strictly ∆t-independent.The dashed lines have been obtained using Eq. (11), the solid line using Eq. (16) taking advantage of the directly measured shear-stress relaxation function R(t).

  , 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 δv G . As shown

Fig. 19 .

 19 Fig. 19. Comparison of δv (filled symbols) and δvG[R] (open symbols) for free-standing films: (a) ∆t-dependence of δvG[R] (open symbols) for all indicated T and of δv (filled symbols) for T = 0.55, 0.4, 0.2 and 0.05. The 1/ √ ∆t-decay for small ∆t is shown by the bold solid line, the plateau value ∆ne ≈ 1 of δv for small T and ∆t τne ≈ 800 by the bold dashed line and δvG ≈ 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 ∆t as indicated.

  If t and ∆t exceed the typical relaxation time τ b of these basins, h(t) and v(∆t) must become constant. As shown in Sec. 3.4 and Sec. 3.5, δv

G ∝ 1/ √ ∆t for ∆t τ b . At variance to this δv → ∆ ne becomes constant with

The longest relaxation time τ of glass-forming liquids is generally called τα[START_REF] Hansen | Theory of simple liquids[END_REF][START_REF] Götze | Complex Dynamics of Glass-Forming Liquids: A Mode-Coupling Theory[END_REF][START_REF] Ferry | Viscoelastic properties of polymers[END_REF].

The ensemble average . . . may be computed by taking the arithmetic average over Nc independently prepared and sampled configurations c. For ergodic systems it is equivalent to sample over N k sub-intervals of length ∆t of a very long trajectory of length ∆tmax τ .

Equation (14) holds if the perturbation is a "deformation". In the case of an externally applied "force" it becomes R(t) = h(t)[START_REF] Doi | The Theory of Polymer Dynamics[END_REF].

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

According to Bochner's theorem f (ω) ≥ 0 if and only if f (t) is a positive-definite function, i.e. all eigenvalues of the matrix gi,j = f (ti -tj) are non-negative[START_REF] Götze | Complex Dynamics of Glass-Forming Liquids: A Mode-Coupling Theory[END_REF].

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 RA and m21 which have different standard deviations δRA δm21 as seen from Fig.8of Ref.[10].

The data of the binary LJ mixture scanned from the first panel of Fig.2of Ref.[41] corresponds strictly speaking to the standard deviation δM for the shear modulus M = RA -v. Since the fluctuations of RA are negligible, however, as shown elsewhere[10,15], δM ≈ δv.

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 = RA -v rigorously cancel.

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

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].[START_REF] Ferry | Viscoelastic properties of polymers[END_REF] This allows us to suppress long-range hydrodynamic modes otherwise relevant for two-dimensional systems.

Experimental studies on polymer films[52,53] rather investigate the (shear-strain) creep compliance J(t)[START_REF] Ferry | Viscoelastic properties of polymers[END_REF]. Since

Acknowledgments

We are indebted to O. Benzerara (ICS, Strasbourg) for helpful discussions and acknowledge important computational resources from the HPC cluster of the University of Strasbourg.

Author contribution statement

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.


domain: sde info info.info-mo

text:
Jean-Christophe Fabre
email: fabrejc@supagro.inra.fr

Michael Rabotin

David Crevoisier

Aline Libres

Cécile Dagès

Roger Moussa

Xavier Louchart

Philippe Lagacherie

Damien Raclot

Marc Voltz

Michaël Rabotin

Glibmm Gtkmm

OpenFLUID : an open-source software environment for modelling fluxes in landscapes

ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Landscapes are complex systems where many processes interact in time and space. In agro-ecosystems, these processes are mainly physical processes, including hydrological-processes, biological processes and human activities. Modelling such systems requires an interdisciplinary approach, coupling models coming from different disciplines, developed by different teams.

In order to support collaborative works, involving many models coupled in time and space for integrative simulations, an open software modelling platform is a relevant answer.

OpenFLUID is a software framework and an operational platform for integrative modelling and simulation of landscapes functionning. It allows to build spatial simulations based on i) pluggable models which are coupled, ii) a digital representation of the landscapes.

contact@openfluid-project.org

OpenFLUID framework

OpenFLUID simulations can be run either using the command line interface or the graphical user interface.

During simulations, the plugged models compute and exchange variables and events attached to every spatial graph node. The plugged observers can monitor the simulation continuously, for exporting data, controlling values, or for performing any non intrusive operation. The data storage and exchanges over the spatial graph for the model coupling is managed by the platform engine.

Simulation results

atmospheric processes

Landscapes are represented as a hierarchical graph where nodes represent the spatial units (fields, ditches, roads, rivers, ...) and can embed rich information (geometry, physical properties, simulation variables,...). Edges represent the connections (from-to, child-parent) between these spatial units.

Objectives: 1) studying the effect of intrafield variability on the runoff generation and flow pathways, 2) comparing different implementations of the same hydrological process.

Methods: A 1200 m² vineyard field (Roujan, France) is divided into more than 1000 surface units. The coupled model is made of a rainfall-runoff model and three surface transfer models which are interchanged for implementations comparison.

Risk analysis of water contamination by pesticides

Hydrological network topology effects Simplified modelling of traffic on Manhattan streets

OpenFLUID relies on many open-source libraries for the core architecture and functioning (Boost, GLib, libXML2, libCURL), the management of GIS data (GDAL, OGR), the spatial algorithms (GEOS) and the GUI features (GTK).

All of these libraries licenses allow a non-restrictive reuse.

OpenFLUID is licensed under the terms of the GPLv3 license, with a special exception.

This special exception allows to plug and distribute models which are not compliant with the standard GPL license.

OpenFLUID uses a collaborative approach either for its development and the scientific applications. OpenFLUID make the sharing of knowledge and source codes easier.

OpenFLUID modelling Real landscape

Objectives: Using OpenFLUID in a very different context of usual applications, demonstrating its abilities in modelling various types of spatial fluxes. Methods: A simplified coupled model of road traffic transfer and traffic lights switching, applied to a spatial representation of streets in a urban context. The OpenFLUID framework can be embedded into the Python language using the PyOpenFLUID module, or into the GNU R statistical environment using the ROpenFLUID package.

The pluggable models can be easily developed de novo or encapsulate existing models using the dedicaded development environment. They are written in C++ or compatible language (C, Fortran, ...) models repositories store for personal or shared reuse reuse of existing models direct use of models

Models development and sharing coupled model simulation monitoring

The input dataset can be built by hand using text editors or GIS environments, or using specifc tools such as GeoMHYDAS (Lagacherie et al, 2010)

Simulation execution

Spatial graph

This poster is licensed under a Creative Commons license

  OpenFLUID software, documentation and collaborative support are available on the web site: http://www.openfluid-project.org OpenFLUID source code is available on the GitHub open-source hosting platform: http://github.com/OpenFLUID/ Objectives: Studying the impact of hydrological network density on erosion and surface water accumulation. Methods: 8000 simulations using a rainfall-runoff-transfer coupled model, with differents hydrological networks densities and topologies (Levavasseur et al. 2012). low density high density balanced density Spatial processes are represented as one or many models. These models simulate the dynamics of processes by computing data attached to every node of the graph, and can rely on rich information embedded by nodes and edges. The spatial coupling of models is based on controlled data exchanges over the spatial graph.

  Puissalicon catchment (100ha, France) is divided into 190 surface units (fields, roads, banks, buffers) and 40 reach segments (ditches, pipes, rivers). The coupled model is made of hydrological models and fate-and-transfer pesticide models.

@OpenFLUID Input dataset

		model for process A	model for processes B and C	model for process D
	digital landscape
	representation, coupled model configuration, monitoring definition, ...	command line interface		export simulation data
				control
				simulation values
				simulation
		spatial graph	scheduler
		manager
		graphical interface