L Klochko

J Baschnagel

J P Wittmer

H Meyer

O Benzerara

A N Semenov
email: al.ni.semenov@gmail.com

Theory of length-scale dependent relaxation moduli and stress fluctuations in glass-forming and viscoelastic liquids

The spacio-temporal correlations of the local stress tensor in supercooled liquids are studied both theoretically and by MD simulations of a two-dimensional (2D) polydisperse Lennard-Jones system. Asymptotically exact theoretical equations defining the dynamical structure factor and all components of the stress correlation tensor for low wave-vector q are presented in terms of the generalized (q-dependent) shear and longitudinal relaxation moduli, G(q,t) and K(q,t). We developed a rigorous approach (valid for low q) to calculate K(q,t) in terms of certain bulk correlation functions (for q = 0), the static structure factor S(q) and thermal conductivity κ. The proposed approach takes into account both the thermostatting effect and the effect of polydispersity. The theoretical results for the (q,t)-dependent stress correlation functions are compared with our simulation data, and an excellent agreement is found for q b 0.5 (with b, the mean particle diameter) both above and below the glass transition without any fitting parameters. Our data are consistent with recently predicted (both theoretically and by simulations) long-range correlations of the shear stress quenched in heterogeneous glassy structures.

I. INTRODUCTION

Mechanical stress governs the flow and deformation dynamics in liquids or amorphous solids (glasses). 1,2 The stress correlation functions therefore provide important information on dynamical properties of amorphous materials. [START_REF] Hansen | Theory of Simple Liquids[END_REF] In particular, it is well-known that time-dependent correlations of the shear stress σ (t) are related to the shear relaxation modulus, G(t), [START_REF] Forster | Hydrodynamic Fluctuations, Broken Symmetry, and Correlation Fluctuations[END_REF][START_REF] Kriuchevskyi | [END_REF][6][7][8][9][10] which is one of the central rheological functions of a material:

σ (t + t ′ )σ (t ′ ) = T V G(t) + const ( 1
)
where V is the system volume and T is its temperature in energy units (T = k B T abs where T abs is the standard absolute temperature), and 'const ' is a constant [START_REF] Kriuchevskyi | [END_REF]6,9,10,12 that vanishes in an equilibrium liquid, but may take a non-zero value for systems with finite (quasi-) static shear modulus. A similar relation between the correlation function, C p (t), of instantaneous pressure p(t), and the bulk compression modulus K 0 (t) was recently discussed 10 (see also Refs. 7, 8, 16, and 17). K 0 (t)

defines the pressure response after a small uniform compression of the system volume (V →

V (1 -ε), ε ≪ 1) at t = 0:

∆p(t) = εK 0 (t), t > 0 (2)
In addition to K 0 (t), there are two other response functions related to compression and useful for many purposes 10 : (i) The adiabatic compression modulus K A (t) characterizes a thermally isolated system whose energy and entropy are constant at t > 0 (no heat exchange with its environment).

(ii) The isothermal compression modulus K T (t) is relevant to the ideal case when the system temperature is kept perfectly constant. These two functions characterize the system as such, whereas the pressure correlation function C p (t) depends also on the thermostatting conditions, and so it is non-universal in this respect (see note [START_REF] K T (t | are universal functions by their definition implying a perfect temperature adjustment to keep constant either T or the entropy S[END_REF] ).

In particular, such dependence on the thermostatting mechanism applies to most of molecular dynamics (MD) simulation studies. 10 This effect is akin to the well-known fact that static correlators (like C p (0)) generally depend on the statistical ensemble. 18 It is possible, however, to obtain the universal (thermostat-independent) material functions, K A (t) and K T (t), based on cross-and auto-correlation functions of pressure and temperature. Some recipes on how to do it are described in Ref. 10 (see also Ref. 11).

Complex and glass-forming liquids are known to be highly heterogeneous 19 , so the spaceresolved correlation functions are required to characterize their structure. Accordingly, in the present paper we study the stress correlation and response functions depending on both wavevector (q) and time (t). In particular, we consider all components of the tensorial correlation function C αβ α ′ β ′ (q,t) (defined in eq. 69) for the q-dependent stress, σ αβ (q,t). The latter is the Fourier transform of the stress tensor field σ αβ (r,t):

σ αβ (q,t) = (1/V ) σ αβ (r,t) exp(-iq • r)d d r (3)
where α, β refer to Cartesian components and d is the space dimension. To simplify the analysis we consider here only strictly two-dimensional systems (2D, d = 2), although many results obtained below are more general. The stress correlation functions are first predicted theoretically (see section III B), and then the results are compared with our MD simulation data (section III C) on a 2-dimensional polydisperse Lennard-Jones (pLJ) system described in Appendix A.

The theory of time-space resolved stress correlations was developed in several recent studies using either the Zwanzig-Mori projection operator formalism 14,20,21 or the fluctuation-dissipation theorem (FDT) 6 . Both approaches considered isothermal systems (formally assuming infinite thermal conductivity) and eventually specify the theories in the hydrodynamic limit of small q (in 2D in Ref. 6, and in arbitrary dimension d in Refs. 20 and 21). The results for the stress correlation tensor obtained within these theories are in harmony with each other. In particular, both theories predict a long-range power law tail (1/r d ) for long-time shear stress correlations, C xyxy (r,t), for supercooled viscoelastic liquids and amorphous systems below the glass-transition temperature.

This prediction is in qualitative agreement with simulation data on stress correlations in 2D 22 and 3D binary LJ models 23 . Recent theoretical, experimental and simulation studies of static strain correlations in 2D colloidal glass-formers reveal an Eshelby-strain pattern with a characteristic 1/r 2 long-range decay near the glass transition. [START_REF] Illing | [END_REF] A similar 1/r 3 decay for strain correlations has been also predicted and tested against simulations and experiments on 3D Brownian hard-sphere colloids 15 . However, to the best of our knowledge, a full quantitative comparison of theoretical and simulation results on q,t-resolved stress correlations has not been performed yet. Here we present such a comparison in 2D. In doing so, we also extend the previous theory 6 by including the effects due to heat transport and particle-size polydispersity, both of which are important for a quantitative analysis of the stress correlation tensor for the simulated glass-forming system.

In the next section we introduce the generalized q,t-dependent longitudinal modulus K(q,t)

and argue that for q = 0 it is essentially universal (independent of the ensemble/thermostat) in contrast to K(0,t) and various correlation functions of volume-averaged quantities (section II C).

A theoretical approach allowing to obtain K(q,t) based on standard time-dependent correlation functions (at q = 0) is presented as well. We then describe a recipe for efficient numerical calculation of K(q,t) taking into account both heat conduction and polydispersity of a molecular or colloidal system (see sections II D, II E, II F), and finally obtain the longitudinal response function for the (2D pLJ) model Lennard-Jones system defined in Appendix A. In the 3rd section it is first explained how to obtain the dynamical structure factor S(q,t) based on K(q,t); the theoretical approach is then verified by comparison with the simulation data for the 2D pLJ model. Furthermore, in section III B we present the basic relations (mostly coming from the FDT) defining the stress-correlation tensor C αβ α ′ β ′ (q,t) in terms of the generalized shear (G(q,t)) and longitudinal (K(q,t)) moduli. All independent components of the stress-correlation tensor are then calculated numerically as described therein. The theoretical and simulation results for C αβ α ′ β ′ (q,t)

are compared in section III C followed by a further discussion of the generalized elastic moduli in section IV revealing long-range spacio-temporal stress correlations. The main results are summarized in section V. Finally, Appendix B describes the theory for the temperature autocorrelation function used to determine the thermal conductivity κ.

II. THE Q-DEPENDENT RELAXATION MODULI

A. Uniform deformations, q = 0

In the present study we are mostly interested in two relaxation moduli, G(q,t) and K(q,t)

corresponding to shear and longitudinal deformation, respectively. [START_REF] Forster | Hydrodynamic Fluctuations, Broken Symmetry, and Correlation Fluctuations[END_REF][START_REF] Balucani | Dynamics of the liquid state[END_REF][START_REF] Semenov | [END_REF]27 Their definition for q = 0 is well-known: G(0,t) ≡ G(t) is the classical shear relaxation modulus related to the (ensembleaveraged) shear stress response, ∆σ xy (t), upon a canonical-affine xy-shear deformation γ when the coordinates (r) and velocities (v) of all particles are instantly changed at t = 0 as

r y → r y + γr x , v x → v x -γv y (4)
(for shear along y-axis with gradient in x-direction):

∆σ xy (t) = G(t)γ (5)
for γ → 0. G(t) is therefore the linear response function (quadratic and higher-order terms in γ and other perturbation magnitudes are neglected here and in what follows).

The longitudinal modulus K(0,t) ≡ K L (t) is defined in a similar way as a response of longitudinal pressure, p xx = -σ xx , to a small uniaxial compression ε = -ε xx of the system at t = 0 (while the transverse area, i.e. the cross-section perpendicular to the deformation axis, is kept constant),

r x → r x (1 -ε) , v x → v x (1 + ε) (6)
The ensemble-averaged longitudinal pressure increment is then

∆p xx (t) = K L (t)ε (7)
It is well-known that in isotropic systems the two moduli fully define the stress response to any affine deformation [START_REF] Landau | Theory of Elasticity[END_REF][START_REF] Tadmor | Modeling Materials[END_REF] . In particular, eq. 7, can be supplemented with

∆p yy (t) = M(t)ε
where

M(t) = K L (t) -2G(t)
is the transverse modulus. The bulk compression modulus K 0 (t)

(cf. Ref. 10) is therefore related to K L : K L (t) = K 0 (t) + (2 -2/d)G(t) (8)
where d is the space dimension.

The moduli G(t) and M(t) can be thought of as time-dependent generalizations of the elastic coefficients of isotropic solids [START_REF] Landau | Theory of Elasticity[END_REF][START_REF] Tadmor | Modeling Materials[END_REF] : G(t) corresponds to the shear modulus µ and M(t) to the Lamé coefficient λ . With this correspondence the longitudinal modulus

K L (t) = M(t) + 2G(t)
agrees with the known relation from the linear elasticity (K L = λ + 2µ) and the bulk compression modulus K 0 (t) with that for the classical static bulk modulus (K 0 = λ + 2µ/d).

B. Non-uniform deformations, q = 0 In this section we formulate a microscopic definition of the generalized q-dependent relaxation moduli (which is alternative to the usually adopted way to define them based on phenomenological hydrodynamic equations [START_REF] Hansen | Theory of Simple Liquids[END_REF][START_REF] Balucani | Dynamics of the liquid state[END_REF] ).

It is well-known that in the general case the deformation rate is defined by the flow velocity field v(r,t). [START_REF] Bird | Dynamics of Polymeric Liquids: Fluid Mechanics[END_REF] The microscopic definition of its Fourier transform, v(q,t)

= V -1 v(r,t) exp -iq • r d d r, is: v(q,t) = J(q,t)/ρ (9)
where ρ = ∑ a m a /V is mean mass density, J is momentum density,

J(q,t) = V -1 ∑ a m a v a (t) exp -iq • r a (t) (10)
a = 1, . . ., N is the label of a particle, and m a , v a and r a are the mass, velocity and position of the a-th particle, respectively. The deformation rate is related to the tensor of velocity gradients

γαβ = ∂ v β (r,t) ∂ r α = ∑ q iq α v β (q,t) exp iq • r (11)
where the sum runs over all q-modes defined by the system size, and α, β are Cartesian components. Here we assume that the system boundaries are fixed (undeformed), 13 and that is why there is no deformation at q = 0.

Let us consider a single deformation mode with q = 0. Without loss of generality we can assume that q is parallel to the x-axis. Then γxy = γe iqx (with γ = iqv y ) is the generalized shear deformation rate, and εxx = -εe iqx with ε = -iqv x correspond to the generalized longitudinal deformation. An instant shear deformation mode with amplitude γ is emerged at t = 0 if γ = γδ (t) is imposed by means of an external force field; this deformation is equivalent to imposing v y = -iγδ (t)/q, where δ is the Dirac's δ -function. Similarly, an instant external longitudinal deformation with amplitude ε corresponds to v x = iεδ (t)/q. The resultant particle displacement field is

u(r) = i q (ε x -γ ŷ) e iqx ( 12
)
where x is the unit vector along the x-axis and ŷ along the y-axis. To render the whole deformation canonical the particle velocities should be changed accordingly (cf. eqs. 4, 6):

v → v + xe iqx (εv x -γv y ) (13)
Note that at q = 0 the above equation agrees with eqs. 4, 6.

The generalized moduli can then be defined via the ensemble-averaged stress responses to the specified instant shear and longitudinal deformations at a wave-vector q:

∆σ xy (q,t) = G(q,t)γ, ∆p xx (q,t) = K(q,t)ε (14)
The above equations are valid provided that no further deformation occurs at t > 0:

J(q,t) ≡ 0 at t > 0 (15)
However, the very stress generated by the initial deformation at t = 0 would lead to a flow giving rise to further deformations of the system. This flow is governed by the fundamental momentum

equation 3,6 ∂ J α (q,t) ∂t = iq β σ αβ (q,t) (16)
Therefore, in the general case the stress response comes as a linear superposition of contributions from all the previous small deformations. So, for example, the first eq. 14 should be replaced by [START_REF] Bird | Dynamics of Polymeric Liquids: Fluid Mechanics[END_REF] ∆σ xy (q,t) = G(q,t)γ(0)

+ t 0 + G(q,t -t ′ )dγ(t ′ ) (17)
where dγ(t ′ )/dt ′ is the deformation rate defined by v(q,t). To avoid this complication, one can apply a perturbative external force F(q,t) to the system (adding it to the rhs of eq. 16) in order to prevent any deformation at t > 0 (that is, to render eq. 15 valid at all times after the initial deformation). In real space this is equivalent to application of an external force field with

F a (t) = m a A(t)e iq•r a (18)
being the force on the a-th particle. Here A(t) is proportional to F(q,t). (Indeed, recalling the definition of the Fourier transform we get F(q,t) = 1 V ∑ a F a (t)e -iq•r a = ρA(t)). Obviously, the forces F a (t) provide a coherent acceleration of all particles (with amplitude A(t)). 26 Since physical variables like r a (t), v a (t) are necessarily real, changing q to -q always leads to complex conjugation of q-dependent variables. Hence, for example, σ αβ (-q,t) = σ * αβ (q,t) (cf. eq. 3). Therefrom we get G * (q,t) = G(-q,t), K * (q,t) = K(-q,t). Besides, the isotropy of the system also ensures that these response functions depend only on |q|, hence they must be real.

C. Low-q regime for G(q,t) and K(q,t)

Having defined the generalized relaxation moduli, we now turn to the problem of their calculation. G(q,t) at q = 0 is a central rheological function, and there are well-known techniques to obtain it either experimentally [START_REF] Bird | Dynamics of Polymeric Liquids: Fluid Mechanics[END_REF] or based on simulation data. [START_REF] Forster | Hydrodynamic Fluctuations, Broken Symmetry, and Correlation Fluctuations[END_REF][START_REF] Kriuchevskyi | [END_REF][7][8][9][10][START_REF] Allen | Computer Simulation of Liquids[END_REF] A commonly employed efficient way to calculate G(t) is to use eq. 1. As for the q-dependence of G, we stick to the idea that G(q,t) is uniformly continuous at q → 0 3,6,20 , hence the following approximation is reasonable:

G(q,t) ≃ G(t) for λ ≫ b ( 19
)
where λ = 2π/q is the wavelength, and b is the molecular size (i.e., the mean particle diameter for the pLJ system we consider). Moreover, we are confident that the q-dependence of G is regular, so the deviation of G(q,t) from G(t) is proportional to q 2 . For molecular systems with finite interaction range this can be rigorously proved for t = 0 (details to be given elsewhere). For t > 0 a weak q-dependence of G(q,t) at low q's comes from the following physical argument: let's divide the system in cells of size l, b ≪ l ≪ λ ; within each cell the deformation is almost uniform, so the stress response is nearly given by G(t); any effect of stress relaxation in a distant cell cannot be transferred to a given cell because any deformation (at t > 0) is forbidden by definition of G(q,t), so there is no mechanism by which a structural relaxation of a given cell can be coupled with that in a distant cell (the locality principle).

It may seem that the above argument is also applicable to the longitudinal modulus K(q,t), which would imply the relation K(q,t) ≃ K(q = 0,t) = K L (t). However, this approximation does not hold. On the one hand, it is well known that fluctuations of volume-averaged (i.e., q = 0) variables, like longitudinal stress σ xx (t) or pressure p(t), depend on the statistical ensemble 27,32 .

On the other hand, the related moduli K L (t) and K 0 (t) also depend on whether the dynamics is energy-conserving or thermostatted and thus, in the latter case, on the nature and parameters of the thermostat used in the simulation. 10 This is so because a compression of the system leads to a temperature rise whose evolution is controlled by the thermostat, and thus entails a thermostatdependent pressure response. 10 By contrast, the longitudinal modulus K(q,t) for q > 0 is universal: it does not depend on the statistical ensemble or the presence of a thermostat. The main reason for this is that fluctuations of physical fields at q > 0 are not coupled with volume-averaged (i.e., q = 0) fluctuations of the mean T which are controlled by the thermostat. Therefore, lim q→0 K(q,t) is universal, while K(q = 0,t) = K L (t) is not. As a result, K(q,t) can be discontinuous at q = 0. [START_REF]We do not use subscript 'L' in K(q,t) since a q-dependent compression can be only longitudinal[END_REF] Noteworthily, the effects considered just above are not relevant for G(t) since the shear stress σ xy is not coupled with any scalar field (in particular, temperature fluctuations) for symmetry reasons.

Hence, it does not matter for G(t) whether the dynamics is energy-conserving or thermostatted.

To sum up, G(q,t) just weakly (negligibly in the regime q b ≪ 1) depends on q because the local shear stress response to a prescribed shear deformation depends on the local structure only. By contrast, the longitudinal stress response to a small inhomogeneous extension ε (at a small wavevector q) can significantly depend on q since such a response (associated with K(q,t)) depends not only on the deformation ε, but also on the deformation-induced temperature (T ) increment which is (essentially) a slow conserved variable whose relaxation dynamics is non-local. As a result K(q,t) gets non-local dynamically, i.e., at t > 0, due to the stress-temperature coupling (by contrast K(q,t = 0) can be approximated by its q → 0 limit just like G(q,t)).

D. Calculation of K(q,t) for monodisperse systems

In the previous section we showed that a relation like eq. 19 does not work for the longitudinal modulus K(q,t). Is it possible at all to find K(q,t) based on correlation functions at q = 0? The answer is positive, but the problem is not exactly trivial as clarified below.

Recently 10 we showed how to obtain two universal bulk compression relaxation moduli, K A (t)

for adiabatic compression and K T (t) for isothermic compression, based on correlation data from thermostatted simulations. We also described how to calculate two other universal functions: the time dependent isochoric heat capacity per particle, c v (t), and the thermal pressure p T (t). 10,11 The latter function is defined via the isochoric pressure response (at constant volume), ∆p(t), to an instant small temperature rise, ∆T , at t = 0:

∆p(t) = p T (t)∆T (20)
Its static limit is a thermodynamic derivative: p T (∞) = (∂ p/∂ T ) V . As shown below, K(q,t) at low q can be expressed in terms of these functions. The result also depends on the thermal conductivity κ of the system. If κ = 0, there is no heat transfer, so K(q,t) must be nearly equal to the adiabatic longitudinal modulus K AL (t) (here we apply the same locality principle as for G(q,t)):

K(q,t) ≃ K AL (t), λ ≡ 2π/q ≫ b ( 21
)
where

K AL (t) = K A (t) + (2 -2/d)G(t) (22)
by virtue of eq. 8.

In the general case (finite κ) the modulus K(q,t) is lower. To obtain it we note that in addition to the adiabatic deformation of a cell applied at t = 0 (with compression strain ε = -ε xx ), some heat h(t) can be transferred to it later on (at t > 0) by thermal conduction. The heat itself leads to a pressure increment

∆p (h) (t) = t -∞ dh(t ′ ) dt ′ χ(t -t ′ )dt ′ (23)
where χ(t) is the relevant response function, and h(t) is the heat (per particle) injected into the cell by the time t (the heat injected during a time-period between t and t + ∆t therefore equals h(t + ∆t)h(t)). Therefore, the total longitudinal pressure increment is

∆p xx (t) = K AL (t)ε + ∆p (h) (t) (24)
where the first term in the rhs is due to adiabatic xx-compression, and the last term is the stress increment due to thermal conductivity leading to the energy increment h(t). Here and below we assume for simplicity that the wave-vector q is oriented along the x-axis, so the x-dependence of pressure increment ∆p xx (and other variables: ∆T , h, ε) is given by the factor e iqx which is systematically omitted. (Note that the heat-induced stress tensor change is ∆σ

(h)
αβ (t) = -∆p (h) (t)δ αβ since energy is a scalar variable.)

In order to find χ(t) we note that the heat h(t) leads to a temperature perturbation, ∆T (h) (t), such that

h(t) = t -∞ c v (t -t ′ ) d dt ′ ∆T (h) (t ′ )dt ′ (25)
where c v (t) is time-dependent isochoric heat capacity per particle. Eqs. 23, 25 can be simplified by using a modified Laplace transform 10,11 referred to in what follows as 's-transform' and defined for a function X (t) as

X (t) → X (s) = s ∞ 0 X (t)e -st dt (26)
Eqs. 23, 25, 20 then read in terms of s-transforms:

∆p (h) (s) = h(s)χ(s) (27)
h(s) = c v (s)∆T (h) (s) (28)
∆p (h) (s) = p T (s)∆T (h) (s) (29)
The above equations lead to

χ(s) = p T (s)/c v (s) (30)
and (cf. eq. 24)

∆p xx (s) = K AL (s)ε + h(s)p T (s)/c v (s) (31)
It therefore remains to find h(s). It is defined by the heat transport equation

n 0 ∂ h(t) ∂t = κ∇ 2 T = -κq 2 ∆T ( 32
)
where n 0 = N/V is the mean concentration, and the Laplacian ∇ 2 is replaced with -q 2 in view of the harmonic dependence of all perturbations given by the factor e iq•r . Here

∆T (t) = ∆T (i) (t) + ∆T (h) (t) (33)
is the total temperature increment due to two effects: (i) initial instant adiabatic compression, ε = -ε xx , at t = 0; (ii) the heat conduction. Applying the FDT to an energy-conserving canonical ensemble we get:

∆T (i) (t) = ε T n 0 χ(t) (34)
(The latter equation can be obtained by using eqs. (B2), (B9) of Ref. 10 and noting that ∆T in eq.

(B2) is actually ∆T (i) , and that the 2nd term in the 1st eq. (B9) represents ∆p (h) with h = ε v T .)

Physically, eq. 34 comes from the fact that ∆p (h) as predicted by the FDT is proportional to

C pT (t) = N T 2 δ p(t + t ′ )δ T (t ′ ) (35)
while ∆T (i) must be proportional to C T p (t) which is equal to C pT (t). [START_REF]mean ensemble average supplemented with the averaging over t ′ , and δ X ≡ X -X for any X = T (t)[END_REF] In the static limit, t → ∞, eq. 34 simply comes from eq. 27 and the thermodynamic relation:

T (∂ p/∂ E) V = -(∂ T /∂V ) S ,
where E, S are total energy and entropy, respectively.

As a result, using eqs. 28, 34 we get:

∆T (s) = ε T n 0 p T (s) + h(s) /c v (s) (36)
Next, we use the above equation, eq. 31, and eq. 32 whose s-transform reads sh(s) = -(κ/n 0 )q 2 ∆T (s), to obtain ∆p xx (t) = K(q,t)ε with

K(q, s) = K AL (s) - T n 0 p T (s) 2 c v (s) 1 + sn 0 c v (s)/ κq 2 -1 (37)
We anticipate that this equation is valid at low q's with relative accuracy of ∼ ( b/λ ) 2 (cf. eq. 19

for G(q,t)). Still the sign ' =' is used here and below to simplify equations (whose approximate nature must be kept in mind). For κ → 0 the 2nd term in the above equation vanishes, so the modulus K(q,t) becomes purely adiabatic and eq. 21 is recovered. By contrast, for κ → ∞ eq. 37 becomes

K(q, s) = K AL (s) - T n 0 p T (s) 2 c v (s) ≡ K T L (s) (38)
In this case K(q,t) = K T L (t), where

K T L (t) = K T (t) + (2 -2/d)G(t)
is the isothermic longitudinal modulus defined in eq. 38. Note that the difference

K AL (t) -K T L (t) is equal to K A (t) -K T (t) (cf.
eq. (D23) of Ref. 10) as it should be. The general eq. 37 can be rewritten in terms of the isothermic relaxation modulus K T L :

K(q, s) = K T L (s) + T p T (s) 2 / n 0 c v (s) + κq 2 /s (39)
Finally, the relaxation modulus K(q,t) can be obtained by doing an inverse s-transform of eq. 37 or eq. 39 provided that the response functions K AL (s), p T (s), c v (s) are known. A rigorous way to obtain these functions based on correlation data from thermostatted simulations was described recently 10,11 . Noteworthily, the approach to obtain the general relations between time and wavevector dependent response functions and correlation functions developed here and in refs. 10 and 11 is similar to the generalized hydrodynamic theory [START_REF] Götze | [END_REF][START_REF] Kadanoff | [END_REF] . The main differences are that (i) we now adapt the theory rendering it applicable to the thermostatted simulations with imperfect temperature control, so that the universal (thermostat independent) response functions can be calculated based on non-universal dynamical correlation functions, and (ii) we systematically use the s-transform (instead of merely the Laplace transform of time-dependent functions) which renders many dynamical relations between response functions look pretty similar to their static counterparts (see, e.g., eqs. (D23) and (D25) of ref. 10, or eq. 77 below).

In the next section we propose rather simple approximate ways to obtain K(q, s) which largely avoid the need to do direct or inverse Laplace transforms.

E. Approximate relations for K(q,t)

In most thermostatted simulation schemes the isothermicity condition only holds at long enough time-scales, t ≫ τ damp , where τ damp is the characteristic damping time associated with the thermostat employed in the simulations. As a result, a simple linear relation between, say, the dynamical heat capacity, c v (t), and the energy autocorrelation function gets violated, so one has to resort to more complicated equations 10,11 . A simple way out is to set the thermostat parameters so as to decrease the characteristic time τ damp below the time-scales of interest. We adopted this way in simulations of the 2D pLJ system involving the Nosé-Hoover thermostat (see Appendix A):

the chosen thermostat-related Nosé time, τ damp ∼ 0.01, ensured that T ≈ const at t 0.1, which means that all the correlation data (with the sampling time-step of δt = 0.05) are nearly isothermic, maybe except for t = 0, 0.05. 34 Then, all correlations involving temperature fluctuations can be neglected, and it is reasonable to approximate

c v (t) with c v0 (t) = C E (0) -C E (t), 35 p T (t) with p T 0 (t) = C pE (0) -C pE (t), and K T (t) with K 0 (t) = T n 0 [C p (t) -C p (0)] + η A ,
C p (t) = N T 2 δ p(t + t ′ )δ p(t ′ ) (40)
C pE (t) = 1 T 2 δ p(t + t ′ )δ E(t ′ ) , C E (t) = 1 T 2 N δ E(t + t ′ )δ E(t ′ ) (41)
Here E is the total energy of the system, η A = K 0 (0) is its affine compression modulus, δ X (t) = X (t) -X , and X means time and ensemble average of a variable X . Applying these approximations in eq. 39 leads to

K(q, s) ≃ K * (q, s) ≡ K L (s) + T p T 0 (s) 2 / n 0 c v0 (s) + κq 2 /s (42)
where

K L (s) = K 0 (s) + (2 -2/d)G(s).
Using simulation data on the correlation functions for the system described in Appendix A we established that the approximation, eq. 42, works remarkably well: it is asymptotically exact both at short and long t, and at intermediate times its relative error is 0.1%, which is comparable to (or smaller than) the statistical error bar of our data.

Eq. 42 can be used to obtain the longitudinal modulus in the time domain:

K(q,t) ≈ K * (q,t) ≡ K L (t) + ∆K * (q,t) (43)
where the 1st term is

K L (t) = K 0 (t) + (2 -2/d)G(t) ≈ K T L (t),
and the s-transform of the 2nd term is

∆K * (q, s) = T n 0 p T 0 (s) 2 ψ(s) , ψ(s) ≡ c v0 (s) + κ n 0 q 2 /s (44)
It leads to the following equation for ∆K * (q,t):

t 0 ∆K * (q,t -t ′ )dψ(t ′ ) = T n 0 t 0 p T 0 (t -t ′ )dp T 0 (t ′ ) (45)
where

ψ(t) = c v0 (t) + κ n 0 q 2 t (46)
Eq. 45 is similar to eq. (D7) of Ref. 10 and can be solved numerically in the same way. The results for our 2D pLJ model system at q = q min ≡ 2π/L (L is the system box size) and T = 0.4, 0.3, 0.26, 0.24, 0.2 are shown in Fig. 1 (black curves from bottom to top). The reduced thermal conductivity, κ * = κq 2 min /n 0 = (2π) 2 κ/N, was obtained based on temperature autocorrelation data as described in Appendix B: κ * = 0.0426, 0.0474, 0.0466, 0.0471, 0.0490 for T = 0.4, 0.3, 0.26, 0.24, 0.2, respectively. It is apparent that κ increases just slightly as T is decreased.

A simple approximate expression for ∆K * (q,t) can be obtained by assuming that p T 0 and c v0 are constants, then doing analytically the inverse s-transform of eq. 44, and then replacing p T 0 and c v0 with time-dependent functions:

∆K * (q,t) ≈ T n 0 p T 0 (t) 2 c v0 (t) exp - κq 2 t n 0 c v0 (t) (47)
This equation is analogous to eq. (A5) of Ref. 10. The above approximation of ∆K * (q,t) was used to calculate K(q,t) ≈ K L (t) + ∆K * (q,t). The results are also shown in Fig. 1 by red curves. One can observe that the approximation, eq. 47, works very well in a wide T -range (both above and below the glass transition temperature T g ≈ 0.26). [START_REF]The glass-transition temperature T g ≈ 0.26 was determined in prior Monte Carlo (MC) and MD studies[END_REF] The relative error is 1% for T = 0.4, 0.3 and 0.25% for T = 0.26, 0.24, 0.2.

Eqs. 43, 47 were therefore employed to calculate K(q,t) as described in the next sections.

F. Polydispersity effects for K(q,t)

At long times comparable with the total sampling time ∆t max (cf. Appendix A), t ∼ ∆t max (more precisely, at t ≫ τ T , beyond the temperature wave relaxation time τ T = n 0 c vs / κq 2 , where

c vs = c v (t → ∞)
is the static heat capacity per particle) the correction ∆K * vanishes, so K(q,t) ≃ K L (t), and the isothermic static limit of the longitudinal modulus, K T Ls = K T L (t → ∞), can be approximated as

K T Ls ≃ K L (t ∼ ∆t max ) ≃ K T L (t ∼ ∆t max ) (48)
On the other hand, the equilibrium static modulus K e (q) can be obtained from the generalized compressibility equation [START_REF] Balucani | Dynamics of the liquid state[END_REF]27,[START_REF] Klochko | [END_REF]

S(q) = n 0 T /K e (q) (49)
where

S(q) = 1 N ∑ n e iq•r n 2 ( 50
)
is the static structure factor. Using the S(q) data for the 2D pLJ system we found that K e (q) from eq. 49 (= n 0 T /S(q)) is significantly lower than the static modulus defined in eq. 48 even in the liquid regime, ∆t max ≫ τ α , where τ α is the structural α-relaxation time.

What is the reason for such a contradiction? In fact, it is related to the well-known matter that the standard compressibility equation 3 is violated in polydisperse systems, so

lim q→0 S(q) = n 0 T /K T Ls (51)
where the long-time limit of the isothermic longitudinal modulus

K T L (t), K T Ls ≃ lim t→∞ K T L (t),
defines the isothermic compressibility κ T in the liquid regime: κ T ≃ 1/K T Ls . Physically this discrepancy is due a coupling of composition fluctuations (which are absent in monodisperse systems) with fluctuations of the total concentration field n(r) at any (however small) wave-vector q. The mechanism of such coupling was elucidated in Ref. 36. To illustrate it let us consider a bidisperse mixture of small and large particles. A small longitudinal deformation ε of the system would lead to a concentration wave δ n(x) = n 0 εe iqx . Initially such deformation does not change the local composition of the system (the number fraction of small particles) which remains homogeneous, so δ n(x) implies a perturbation of the particle volume fraction δ φ (x) ≈ εe iqx leading to a high free energy cost which is proportional to ε 2 K T Ls with high static modulus K T Ls for a nearly incompressible system. However, an exchange of small and large particles leading to a higher fraction of small particles in denser regions (and higher fraction of large particles in less dense regions) can significantly reduce δ φ leading to a lower energy cost. As a result, the emerging compositional wave gives rise to a significant decrease of K(q,t) at long times from K T Ls to K e (q) = n 0 T /S(q), since S(q) significantly increases with the polydispersity degree. [START_REF] Klochko | [END_REF] An analogous process of compositional relaxation leads to the terminal relaxation stage of the dynamical structure factor S(q,t) in polydisperse systems: [START_REF] Klochko | [END_REF] (see also note [START_REF]The importance of interdiffusion processes for glass-forming systems with size polydispersity was discussed in Ref. 60[END_REF] )

S(q,t) ≃ (S(q) -n 0 T /K T Ls ) exp -q 2 tD s , t ≫ τ T , τ α (52)
where the α-relaxation time τ α normally coincides with the terminal relaxation time for G(t).

The time-scale of this process is defined by the mean inter-diffusion coefficient D inter of the particles, which nearly equals their self-diffusion coefficient D s , D inter ≈ D s , in weakly polydisperse systems. [START_REF] Klochko | [END_REF] The terminal relaxation time of K(q,t) is related to D s as well, but there is an important difference: while the terminal relaxation of S(q,t), eq. 52, is mostly driven by the ideal-gas entropy of mixing of polydisperse particles, the K(q,t) relaxation process is driven by the much higher elastic free energy of the longitudinal deformation which is proportional to K T Ls . More precisely, according to the FDT [START_REF] Klochko | [END_REF] , the function S(q,t) defines the concentration response δ n q (t)

to a weak external field U (r,t) = U q e iq•r applied to the system at t ≥ 0:

δ n q (t) = n 0 T U q [S(q) -S(q,t)] (53)
(ensemble-averaging as assumed in the l.h.s. of the above equation). Importantly, the definition of K(q,t) implies that an analogous field (denoted as U * q in this case) is applied to the system to keep the longitudinal deformation (at wave-vector q) constant (i.e., independent of time), so that the concentration perturbation δ n * q (t) = n def is constant as well. The difference between the two cases is that U q in eq. 53 is constant (at t ≥ 0), while U * q is time-dependent. A meaningful comparison between the two cases demands that the terminal δ n q (t → ∞) in eq. 53 must be equal to n def (to ensure that all particle components are eventually distributed in the same way in both cases). Therefore we get

U q = n def T n 0 S(q) (54)
Let us consider the very beginning of the particle interdiffusion process (τ T , τ α ≪ t ≪ 1/(q 2 D s )).

The interdiffusion rate is then driven by the external field U q (in the case of eq. 53) and by U * q in the case of longitudinal relaxation modulus, hence

D * /D s = U * q /U q (55)
where D * is the effective interdiffusion constant relevant for K(q,t). The amplitude U * q can be found by minimization of the effective free energy

1 2 K T Ls δ n * q /n 0 2 -U * q δ n * q leading to U * q = K T Ls n def /n 2 0 ( 56
)
Using eqs. 54, 55, 56 we get the effective diffusion constant for the terminal relaxation of the longitudinal modulus:

D * = D s K T Ls /K e (q) (57)
The terminal compositional relaxation stage of K(q,t) is therefore given by the following equation (an alternative derivation of eq. 58 is presented in section IV):

K(q,t) ≃ K T Ls + (K e (q) -K T Ls ) 1 -exp -q 2 tD s K T Ls /K e (q) , t ≫ τ T , τ α (58)
This equation assumes that the self-diffusion is Fickian. In the general case (in particular, at T T g ) it is the mean-square particle displacement (MSD) that matters, so tD s must be replaced with

h 0 (t) = 1 2d r a (t + t ′ ) -r a (t ′ ) 2 ( 59
)
where averaging over all particles (index a) is assumed in the rhs Combining thus modified eq. 58 with eqs. 47, 43 we get the final approximate result:

K(q,t) ≈ K L (t) + T n 0 p T 0 (t) 2 c v0 (t) exp - κq 2 t n 0 c v0 (t) + + (K e (q) -K T Ls ) 1 -exp -q 2 h 0 (t)K T Ls /K e (q) (60)
This equation is used in the next section to predict different q,t-resolved correlation functions.

The functions entering the rhs of eq. 60 have been considered above: The longitudinal modulus K L (t) is defined in eq. 8, the bulk shear modulus G(t) is obtained based on the stress-fluctuation relation; 10 the bulk compression modulus K 0 (t), the thermal pressure p T 0 (t) and the heat capacity c v0 (t) are obtain based on pressure and energy auto-and cross-correlation functions as described in ref. 10. h 0 (t) comes from the particle MSD; K e (q) is defined by the standard structure factor S(q) using eq. 49. Finally the heat conductivity κ is obtained as described in Appendix B. The eventual time-dependence of K(q,t) at T = 0.4 and 0.3 is illustrated in Fig. 2.

III. THEORY FOR q,t-RESOLVED TENSORIAL STRESS CORRELATION FUNCTIONS A. The dynamical structure factor Before turning to stress correlations we chose to verify the basic result of the previous section, eq. 60. To this end let us consider the dynamical structure factor

S(q,t) = 1 N ∑ a,a ′ exp iq • r a (t + t ′ ) -r a ′ (t ′ ) (61)
where a, a ′ run over all particles in the system, and an averaging over t ′ is assumed. For macroscopically isotropic systems only the magnitude of q is relevant, [START_REF]Macroscopic isotropy of the studied systems was demonstrated[END_REF] so S(q,t) can be averaged over orientations of q. The linear response theory leads to the following rigorous relation [START_REF] Balucani | Dynamics of the liquid state[END_REF]27 between K(q,t) and S(q,t):

S(q) -S(q, s) = T ms 2 /q 2 + K(q, s)/n 0 -1 (62)
where S(q, s) is the modified Laplace transform (defined in eq. 26) of S(q,t). The above equation shows that ∂ S(q,t) ∂t = 0 and ∂ 2 S(q,t) ∂t 2 = -T m q 2 at t = 0. The relation, eq. 62, comes from the FDT (or the Zwanzig-Mori projection operator formalism) [START_REF] Hansen | Theory of Simple Liquids[END_REF][START_REF] Forster | Hydrodynamic Fluctuations, Broken Symmetry, and Correlation Fluctuations[END_REF][START_REF] Götze | Complex Dynamics of Glass-Forming Liquids, a mode-coupling theory[END_REF] ; it is general and exact provided that all particles have the same mass. It allows to predict S(q,t) based on the known K(q,t). The latter was calculated using eq. 60 and the classical timedependent correlation functions for q = 0 involved in it (like K L (t), c v0 (t), h 0 (t)) obtained using simulation data for the 2D pLJ system (Appendix A). A multi-exponential approximation was then applied to the function K(q,t) at a given q (cf. point 6 of section IV for further discussion):

K(q,t) = ∑ i K i e -γ i t (63)
Here i = 0, 1, . . ., i m , with non-negative γ i and real K i , γ 0 = 0 and K 0 > 0. K 0 = K e (q) ≡ n 0 T /S(q)

in the liquid regime (T > T g ) when a full relaxation of the system can be achieved. (The number of modes in eq. 63 and in the analogous eq. 76 was around 25.) At low temperatures, T T g , K 0 should be treated as a quasi-static longitudinal modulus at t ∼ ∆t max , K 0 = K T Ls (here ∆t max is the total sampling time, see Appendix A). The s-transform of K therefore is

K(q, s) = ∑ i K i s s + γ i ( 64
)
The general expressions for S(q,t) and S(q, s) can be then deduced from eqs. 62, 64:

S(q, s) = Y 0 + ∑ k Y k s s -r k , S(q,t) = Y 0 + ∑ k Y k exp(r k t) (65)
where k = 1, 2, . . ., i m + 2, r k are the roots of the equation smn 0 /q 2 + K(q, s)/s = 0 (hence, r k are roots of a polynomial of order i m + 2) and Y 0 , Y k are constants. Furthermore, using again eq. 62

and eq. 65 we get:

S(q, s) = Y 0 + ∑ k Y k + ∑ k Y k r k s -r k , Y 0 = S(q) -∑ k Y k (66)
where

Y k = n 0 T r 2 k ∑ i K i (r k + γ i ) 2 - mn 0 q 2 -1 (67)
Hence, finally the dynamical structure factor becomes

S(q,t) = S(q) + ∑ k Y k [exp(r k t) -1] (68)
The theoretical predictions based on eq. 68 are compared in Fig. 3 with the simulation data for S(q,t) obtained directly using its definition, eq. 61 (the amplitudes K i involved in eq. 67 have been obtained by minimization of the mean-square deviation between the prescribed K(q,t) and its multi-exponential approximation, eq. 63). A good agreement between the theory and simulation data for S(q,t) is obvious, which means that the basic approximation, eq. 60, works well both above and below T g . 44

B. Stress correlation functions

The time-space resolved stress correlation tensor is defined as

C αβ α ′ β ′ (q,t) = V T δ σ αβ q,t + t ′ ) δ σ α ′ β ′ -q,t ′ ) (69)
where δ σ αβ q,t = σ αβ q,t -σ αβ q,t ′ (note that σ αβ = 0 for q = 0). In two dimensions the full correlation tensor C αβ α ′ β ′ of a macroscopically isotropic system is completely defined by four independent components 6 :

C 1212 ≡ C s , C 1111 ≡ C , C 1122 ≡ C ⊥ , and C 2222 ≡ C 2 , where
axis '1' is chosen to be parallel to q, while axis '2' is perpendicular to it (cf. eq. 103). As shown elsewhere [START_REF] Hansen | Theory of Simple Liquids[END_REF][START_REF] Forster | Hydrodynamic Fluctuations, Broken Symmetry, and Correlation Fluctuations[END_REF][START_REF] Balucani | Dynamics of the liquid state[END_REF][START_REF] Götze | Complex Dynamics of Glass-Forming Liquids, a mode-coupling theory[END_REF][START_REF] Götze | [END_REF] the first two correlation functions (C s and C ) are directly related to the generalized moduli considered above:

C s (q, s) = ρs 2 G(q, s) ρs 2 + q 2 G(q, s) (70) C (q, s) = ρs 2 K(q, s) ρs 2 + q 2 K(q, s) (71)
where (as before) ρ = mn 0 is the mean mass density, s is the Laplace parameter conjugate to time, and the functions C s (q, s), etc. are the s-transforms defined in eq. 26. These equations come from the FDT (the linear response theory) or from the Zwanzig-Mori projection operator formalism; they are exact as such in the thermodynamic limit. An analogous equation for C ⊥ (q, s) comes from the FDT as well 6 :

C ⊥ (q, s) = ρs 2 M(q, s) ρs 2 + q 2 K(q, s) (72)
Here M(q, s) is the s-transform of the transverse modulus M(q,t) defined by the normal stress response in the direction perpendicular to the elongation axis

M(q,t) = lim ε 11 →0 ∆σ 22 (q,t) ε 11 ( 73
)
where ε 11 is the elongation strain applied at t = 0. At q = 0 this modulus is directly related to the shear and longitudinal moduli:

M(0,t) = K(0,t) -2G(0,t) = K L (t) -2G(t)
due to the global isotropy of the system. This relation stays approximately valid also at low q = 0 since a shear response is local (cf. section II C), while the relevant non-local (conserved) fields (like energy density and composition) that affect the normal pressure are scalar and hence provide equal contributions to M(q,t) and K(q,t), so M(q,t) ≃ K(q,t) -2G(q,t), 2π/q ≫ b (74

)
This approximation is used in the next section.

C 2 (q, s) ≃ K(q, s) - q 2 M(q, s) 2 ρs 2 + q 2 K(q, s) (75)
At this point it is worth stressing again that in the general case (when the longitudinal stress relaxation is coupled with 'conserved' scalar fields like energy and composition), it is important to take into account the q-dependence of the moduli. Eq. 75 is approximate: it is valid asymptotically at low q's (i.e., in the hydrodynamic regime λ = 2π/q ≫ b).

The above four relations for the stress correlation functions are verified below based on our simulation data for the 2D pLJ system.

C. Comparison with simulation data

The stress correlation functions can be predicted using the relations of the previous section involving two generalized relaxation moduli, G(q,t) and K(q,t). These relaxation functions are approximately defined in eqs. 19, 60. The function C (q,t) can be obtained based on the known K(q,t) in exactly the same way as the dynamical structure factor, cf. section III A. C s (q,t) was obtained analogously starting with the multi-exponential approximation of G(q,t), cf. eq. 63:

G(q,t) ≃ G(t) = ∑ i G i e -γ i t (76)
with real G i and positive γ i , leading eventually to a similar representation of C s (q,t), but generally with complex amplitudes and rates, cf. eq. 68. We checked that different approximations of G(t)

(different choices of the relaxation rate set {γ i }) lead to virtually identical results for C s (q,t).

The remaining two relations, eqs. 72, 75, can be used to obtain C ⊥ (q,t) and C 2 (q,t) in a similar fashion. The theoretical predictions for all the four stress correlation functions are compared with our simulation data at different temperatures T and wave-vectors q. The results are shown in is τ slow ∼ 10 4 . By contrast, the total sampling time ∆t max = 10 5 is not enough for its relaxation at T = 0.24. Furthermore, at T = 0.2 (well below the glass transition) any relaxation of C 2 (q,t) is hardly visible for t 500 (cf. Fig. 4(e)). As for all other correlation functions (C , C ⊥ , C s ), they always relax during a time τ atn 500 which is nearly the same for all the 3 functions at T ≤ T g .

Apparently τ atn is the attenuation time for the relevant sound waves (both shear and longitudinal, cf. points 4 and 8 of section IV for further discussion), while τ slow is proportional to the structural α-relaxation time τ α . Overall, for q = q min ≡ 2π/L we observe an excellent agreement between the theory and simulation data for all the stress correlation functions, all temperatures and in the whole sampling time range. A close inspection of the graphs at low T 's reveals some disagreement in a time range around t ∼ 200. This feature may seem surprising given the agreement at both shorter and longer times t. However, there is a simple reason for it: the stress correlation functions based on the simulation data were smoothed at t 100 to accelerate their calculation, and this procedure led to a reduction of the oscillation amplitude as can be observed in Figs. 4 (d),(e) (see red curves there).

The effect of the wave-vector q is demonstrated in Figs. 5, 6, 7. Generally, the oscillation regime becomes shorter and with higher frequency, as expected, for larger q. Beyond this regime the stress correlations vanish except for autocorrelations of the second normal stress, C 2 . At T = 0.3 the positive monotonic long-time tail of C 2 becomes more pronounced at higher q (see Fig. 5 (b),(c)). It appears that the theory slightly underestimates the terminal relaxation time of C 2 at the highest q's (see Fig. 5(c)). Otherwise the agreement between the theory and simulations stays quite good up to q ∼ 8q min corresponding to q b ∼ 0.5. Below the glass transition, at T = 0.24, the agreement is nearly as good for C s , C , C ⊥ (see Fig. 6). However, to get a similar agreement for the autocorrelation function of transverse stress, C 2 , for the highest q's we had to vertically shift the theoretical curve by a constant ∆: C 2 → C 2 + ∆, with ∆ = 1.9 for n q2 ≡ q/q min 2 = 32, and ∆ = 2 for n q2 = 64. Such a shift was also applied at T = 0.20 (see Fig. 7): ∆ = 4.7 for n q2 = 8, ∆ = 1.2 for n q2 = 32 and ∆ = 1.13 for n q2 = 64. The origin of this shift is discussed in point 12 of section IV. Moreover, Fig. 6 indicates that the terminal relaxation time of the simulated C 2 slightly increases at the highest q's (n q2 = 32, 64) at T = 0.24 (see Fig. 6). To sum up, it would be fair to say that the agreement between the theory and simulations is very good for the functions C s , C at all temperatures and wave-vectors with n q2 ≤ 64 (which corresponds to q b 0.5). A small discrepancy (with relative error typically within 5%) can be observed for C 2 (q,t) in the terminal relaxation regime for the highest q's at low T .

IV. DISCUSSION

1. As we pointed out in section II C, while the bulk longitudinal modulus K L (t) depends on the thermostatting method, the generalized longitudinal modulus K(q,t) is universal, that is independent of a coupling of the system with a thermostat. Such universality also holds for the bulk adiabatic and isothermic moduli, K AL (t) and K T L (t), which can be obtained using the theoretical approach proposed in Ref. 10. For low q the generalized modulus K(q,t) of a monodisperse system can be expressed in terms of K AL (t) or K T L (t) and other bulk response functions, cf. eqs. 37, 39.

Noteworthily, the isothermic longitudinal modulus K T L (t) can be defined via the normal pressure (p xx ) response to a small compression in the x-direction, ε = -ε xx (more precisely, a canonic affine transformation of coordinates and velocities of all particles defined in eq. 6) at constant transverse area (of the cross-section perpendicular to the deformation axis) and constant temperature: ∆p xx (t) = K T L (t)ε. The adiabatic longitudinal modulus K AL (t) is defined in a similar way:

K AL (t) = lim ε→0 ∆p xx (t)/ε
, where it is the total energy E that is kept constant after the canonicaffine uniaxial compression at t = 0. In a similar fashion, c pL (t) (cf. Appendix B) is the heat capacity at constant longitudinal normal pressure p xx and constant transverse area (no deformation in the plane perpendicular to the x-axis). c pL (t) is different from the standard isobaric heat capacity c p (t); it is related to c v (t):

c pL (s) = c v (s)/ 1 - T n 0 p T (s) 2 c v (s)K AL (s) = c v (s) K AL (s) K T L (s) (77)
The 2. The generalized longitudinal modulus K(q,t) for monodisperse systems is defined in eqs. 37, 39. It is remarkable that these equations can be obtained rather easily by first assuming that all response functions (K AL , p T , c v ) are constants and therefore can be expressed via the thermodynamic derivatives. Then, in the final relation each constant must be replaced by a function of the Laplace parameter s, i.e., by the modified Laplace transform (cf. eq. 26) of the corresponding response function. We emphasize that this trick, which was already hinted at 10,[START_REF] Götze | [END_REF] , works precisely with the modified Laplace transforms, and generally not in the time domain.

3. It is noteworthy that the heat conductivity κ (involved in eqs. 37, 39) was obtained by fitting the simulation data for the autocorrelation function of temperature, C T (q,t), with the theoretical equation coming from the generalized hydrodynamics (see Appendix B). It turns out that κ is nearly constant (independent of q) at low q, λ = 2π/q ≫ b. We found however that κ somewhat decreases with q for q b 0.5. The frequency dependence of the conductivity κ (or it dependence on the Laplace parameter s) can be totally neglected for low q (λ ≫ b) since the microscopic time ∼ τ v relevant for the heat conduction (here τ v is the time of particle velocity relaxation by collisions and vibrations) is much shorter than the characteristic time of thermal diffusion process, τ T = 1/ D T q 2 , cf. eq. 78.

4. The generalized longitudinal modulus K(q,t) shows a shoulder at t ∼ τ T 10 ÷ 100 corresponding to the transition from adiabatic to isothermic response due to heat conduction (cf.

Figs. 1,2). The characteristic time of this process, τ T = 1/ D T q 2 , is defined by the isochoric thermal diffusivity

D T = κ n 0 c vs (78)
which is relevant here since K(q,t) is a response to a small strain (ε) kept constant at t > 0 is the equilibrium specific heat). All the three parameters involved in the above equation (thermal conductivity κ, particle concentration n 0 and the heat capacity per particle, c v ) 45 just weakly depend on temperature T for the 2D pLJ system we studied. Therefore, the thermal diffusion time τ T is almost independent of T : τ T ∼ 50 for q = q min (cf. Fig. 1). Well above T g ≈ 0.26 the time τ T for q = q min is longer than the structural relaxation τ α as defined by the shear relaxation modulus G(t), however, τ α exceeds τ T at T 0.3 (cf. Fig. 10 in Ref. 10).

Remarkably, the attenuation time τ L of oscillations in the dynamical structure factor S(q,t)

and the stress-correlation functions (C (q,t), C ⊥ (q,t), C 2 (q,t)) shows a similar q-dependence as τ T (see eq. 84). However, the attenuation time τ L is longer than τ T and the ratio τ L /τ T ∼ 2 ÷ 5 increases at low T (cf. Figs. 3,4). A comparison of eqs. 84 and 78 shows that this effect is due to a decrease of the nondimensional ratio η L / (κm) on cooling below T g , where η L is the relevant (for attenuation) longitudinal viscosity (cf. the text around eq. 84).

5. The terminal stage of relaxation for the dynamical structure factor S(q,t) of a polydisperse system is defined in eq. 52; it is related to composition fluctuations (cf. Ref. 36). Eq. 52 is based on the assumption that the free energy cost of a relevant ('soft') composition fluctuation mainly comes from the ideal-gas entropy, hence the relevant interdiffusion constant is close to the selfdiffusivity D s of the particles. This assumption is generally valid for systems with sufficiently weak polydispersity. It was verified for the 2D pLJ system we studied. [START_REF] Klochko | [END_REF] Furthermore, the terminal relaxation of the longitudinal modulus K(q,t), eq. 58, can be derived from eq. 52 based on the general eq. 62: For t ≫ τ L the term ms 2 /q 2 in the rhs of eq. 62 can be neglected (since the attenuation time of oscillations in K(q,t) is much longer than the oscillation period), hence we get K(q, s) ≃ n 0 T S(q) -S(q, s)

, S(q, s) ≃ (S(q)n 0 T /K T Ls ) s s + q 2 D s

The above equations lead to eq. 58 upon doing the inverse s-Laplace transform.

Finally, the validity of eq. 52 at t 100, and the relevance of composition fluctuations is demonstrated in Fig. 10, where we plot the dynamical structure factor, S(q,t), based on our simulation data, its long-time approximation of eq. 52 and the dynamical composition structure factor, S ψ (q,t), again directly obtained by simulations. Here S ψ (q,t) is defined analogously to eq. 61 but with the prefactor κ a κ a ′ before the exponential, where κ a = b 2 a /b 2 -1 is positive for large and negative for small particles. The function S ψ (q,t) is thus related to the correlation function of the composition order parameter ψ(r,t) = ∑ a κ a δ (rr a ) which equals to zero on the average. Fig. 10 clearly shows that the long-time relaxation of S(q,t) is an interdiffusion process following the decay of composition fluctuations quantified by S ψ (q,t). 6. To facilitate direct and inverse Laplace transformations of the stress-correlation functions (cf. eqs. 70, 71) involving the generalized moduli, G(q,t), K(q,t), we employed multi-exponential approximations of their time-dependencies (see eqs. 76, 63). This approach formally corresponds to the generalized Maxwell model [START_REF] Bird | Dynamics of Polymeric Liquids: Fluid Mechanics[END_REF] for the shear relaxation modulus G(t). Noteworthily, the Maxwell model implicitly assumes that inertial effects are negligible for the structural dynamics.

In this case the generalized Brownian dynamics (involving stochastic moves in the configurational space) are applicable, and the detailed balance dictates that all amplitudes in eq. 76 must be positive, so the relaxation function G(q,t) must monotonically decrease in time. Such relaxation behavior is exhibited, in particular, with the MC dynamics. A similar behavior is expected for K(q,t). The molecular dynamics (MD) employed in the present study does account for inertial effects (which are normally important at short t), hence negative amplitudes (G i or K i ) must be allowed as well. Note that in this general case the relaxation moduli can still be uniformly approximated with multi-exponentials to an arbitrarily high precision as follows from the Weierstrass-Bernstein approximation theorem.

7. Using eqs. 70, 71, 72, 74, 75 one can easily get the stress correlations at t = 0 (at no timeshift):

C s (q, 0) = G A (q), C (q, 0) = K AL (q) (79) C ⊥ (q, 0) ≃ K AL (q) -2G A (q), C 2 (q, 0) ≃ K AL (q) (80)
where G A (q) = G(q, 0) and K AL (q) = K(q, 0) are q-dependent generalizations of affine (bulk) shear modulus µ A and longitudinal modulus η AL ≡ η A + (2 -2/d)µ A , and µ A = G(0), η AL = K AL (0). [START_REF]The affine moduli µ A and η A are discussed[END_REF] The moduli G A (q) and K AL (q) define the immediate stress response to the generalized shear (γ > 0, ε = 0) and longitudinal (γ = 0, ε < 0) deformations, eqs. 12, 13. These deformations are not affine for q = 0, but they are always adiabatic. Note that eqs. 79 are exact, while eqs. 80 are expected to be approximately valid for 2π/q ≫ b. Their precision can be assessed using simulation data for the stress correlation functions at t = 0. Note that the adiabatic moduli K AL (q) = K(q, 0) and

K AL (t) = K(q → 0,t) are different functions, but K AL (q → 0) = K AL (t → 0).
8. Turning to the stress correlation function at t > 0, we observe that initial shear (σ 21 ) and longitudinal (σ 11 ) stresses generate body forces f on the system elements:

f 2 = iqσ 12 e iqx , f 1 = iqσ 11 e iqx (81)
(here we assume for simplicity that the wave-vector q and therefore axis 1 are parallel to the xaxis). These forces in turn generate shear and longitudinal sound waves. Their frequencies ω s , ω L are defined by the relevant shear, c s , and longitudinal, c L , sound velocities:

ω s = qc s , ω L = qc L ; c s ≃ G ′ (q, ω s )/ρ, c L ≃ K ′ (q, ω L )/ρ (82)
where G ′ (q, ω) is the storage shear modulus equal to the real part of G(q, s = iω):

G(q, s = iω) = G ′ (q, ω) + iG ′′ (q, ω) (83)
The longitudinal storage modulus K ′ (q, ω) is related to K(q, s) (the s-transform of K(q,t)) in a similar way. The shear oscillations with frequency ω s are reflected in the time-dependence of C s (q,t), while the longitudinal frequency ω L is relevant for C (q,t), but also for the dynamical structure factor S(q,t) (cf. Fig. 3) and auto-and cross-correlations of normal stresses, C 2 (q,t), C ⊥ (q,t) (cf. Fig. 4). With the data for the pLJ system at hand (cf. Ref. 10 and Fig. 1) we estimate c s ∼ 2 ÷ 4, c L ∼ 8 ÷ 9 using eqs. 82 in the most interesting temperature regime T = 0.3 ÷ 0.2 (around the glass transition at T g ≈ 0.26). These estimates correspond to the oscillation periods for q = q min , 2π ω s ∼ 40 ÷ 25, 2π ω L ∼ 12, in agreement with the data shown in Fig. 4. The attenuation times of these oscillations (τ s for C s (q,t), and τ L for other functions) are defined by effective viscosities η s , η L :

1/τ s ≃ η s 2ρ q 2 , 1/τ L ≃ η L 2ρ q 2 (84)
where the relevant viscosities are related to the dynamical loss moduli, G ′′ (q, ω) and K ′′ (q, ω) [START_REF]smallest root s of equation ρs 2 /q 2 + G(q, -s) = 0, cf. the denominator in the rhs of[END_REF] η s = G ′′ (q, ω s )/ω s , η L = K ′′ (q, ω L )/ω L (85)

Note that the viscosities defined in the above equation are dynamical in nature; they are generally different from the classical viscosities obtained in the low-frequency limit. The correlation functions presented in Fig. 4 show comparable sound attenuation times, τ s ∼ τ L , in the regime 0.2 ≤ T ≤ 0.3; the times increase from τ L ∼ 100 to τ L ∼ 300 on cooling. This means that the relevant effective viscosities are also comparable and decrease as the system is cooled below the glass transition. The dynamic viscosities η s , η L therefore show an opposite temperature dependence as compared to the steady shear viscosity which is proportional to the α-relaxation time.

Interestingly, the low-frequency (s → 0) limit of the longitudinal compression modulus K(q, s), obtained using eqs. 43, 44, is:

K hydro-limit (q, s) ≃ K T Ls + s ζ + 2 1 - 1 d η + K T Ls (γ -1) s + D T q 2 , ( 86
)
where K T Ls (= η A -TC p (0)/n 0 ) is the isothermal compression modulus, γ = c ps /c vs with c ps and c vs being respectively the equilibrium specific heat at constant pressure and at constant volume, D T = κ/ (n 0 c vs ), the thermal diffusivity, and ζ and η the bulk and shear viscosity, respectively.

The viscosities are defined by

ζ = T n 0 ∞ 0 C p (t)dt, η = ∞ 0 G(t)dt (87)
Noteworthily, on the one hand Eq. 86 is consistent with the results obtained for d = 3 and κ → ∞ in Ref. 27 (see text below eq. ( 12) there). On the other hand, for finite κ the term in square brackets of Eq. 86 agrees with Eq. (3.104b) in Götze's book [START_REF] Götze | Complex Dynamics of Glass-Forming Liquids, a mode-coupling theory[END_REF] .

9. It is obvious from Figs. 4 -7 that all stress correlation functions except C 2 (q,t) rapidly vanish beyond the attenuation time τ s or τ L . By contrast, C 2 (q,t) shows a long-time monotonically decaying tail. What is the origin of this difference? Firstly, we observe that the 3 correlation functions, C s , C , C ⊥ involve shear or longitudinal stresses, σ 12 = σ 21 and σ 11 , which generate a body force (cf. eq. 81). Therefore their long-time presence (for t ≫ τ s , τ L ) would contradict mechanical equilibrium of the system, hence these stresses must rapidly vanish also in the glassy regime, T < T g . As for the function C 2 (q,t), it involves only the transverse normal stress σ 22 which does not generate any fluid motion (the body force associated with this component of stress field is identically zero since q is normal to axis 2), so σ 22 (q,t) may stay unchanged for a long time. Generally, the long-time tail of σ 22 (q,t) is defined by the structural processes reflected also in the relaxation moduli G(q,t) ≃ G(t) and K(q,t); it may persist virtually forever in amorphous systems at low T . We therefore arrive at the concept of a quenched (or frozen) stress, which is equal to the mean stress σαβ (q) (=time-average of σ αβ (q,t)): The stress tensor in a glassy state (i.e., a metabasin of the free energy landscape) can be therefore written as δ σ αβ (q,t) = σαβ (q,t) + σαβ (q) (88)

where σαβ (q,t) is the fluctuation part whose time-average is zero. Eq. 88 obviously implies that

lim t→∞ C αβ α ′ β ′ (q,t) = V T σαβ (q) σα ′ β ′ (-q) (89)
From the above discussion it is clear that the only non-zero component of σαβ (q) is σ22 (q) (recall that we consider 2D systems in the present paper). Therefore, C 2 tends to a finite limit at long time,

C 2 (q,t) → V T σ22 (q) 2 ≡ θ (q), t → ∞ (90)
in the amorphous solid regime well below T g (where the relaxation time τ α can be considered as infinite). The ... brackets here imply also averaging over different orientations of q. Such behaviour can be indeed observed in Figs. 4567.

10. As discussed above the stress correlation functions C (q,t), C ⊥ (q,t), C s (q,t) rapidly vanish for t τ L . The behavior of the transverse normal stress autocorrelation function C 2 (q,t) is different at T T g : it shows a slow monotonic decay at t 2τ L (the time τ L is defined in eq. 84). This monotonic regime can be analyzed using eq. 75 leading to the following approximate result:

C 2 (q,t) ≈ 4G(t) 1 - G(t) K(q,t) , t 2τ L (91)
We verified that the above equation does a fair job at t 2τ L and q b 0.5 and agrees with the results of the full numerical analysis based on eq. 75 (cf. section III C) within an accuracy of ∼ 1%.

Moreover, we found that K(q,t) in eq. 91 can be replaced with just K T L (t) since the terminal decay of K(q,t) (by the interdiffusion process) occurs at times longer than the relaxation time τ α of G(t)

for the studied range of wavevector q. Hence (with accuracy of ∼ 5%)

C 2 (q,t) ≈ 4G(t) 1 - G(t) K T L (t) , t 2τ L ( 92
)
This equation shows that the relaxation time of C 2 (q,t) is comparable with τ α .

11. The FDT relations, eqs. 70, 71, 72, and the relation 75, are based on the linear response theory, and are valid for equilibrium systems. The latter condition is not a problem in the liquid state, but amorphous systems (supercooled liquids below T g ) are normally out of the thermodynamic equilibrium. Are the quoted equations useful in this case? The answer is yes. The reason is that while amorphous systems are distributed over the glassy states in a non-equilibrium (and often unknown) way, they can still be well-equilibrated within each metabasin (glassy state). 9 In this case the dynamical (FDT-like) relations stay valid provided that a (generally unknown) timeindependent constant is added in the rhs [START_REF] Kriuchevskyi | [END_REF]9 (this constant term accounts for an unknown relation between the partition function in a glassy state and in its deformed counterpart). Eq. 89 ensures that the correlation functions C s , C , C ⊥ must tend to 0 at long times, so the additional constant must vanish for the corresponding FDT eqs. 70, 71, 72, which therefore stay unchanged also below T g (note that the t → ∞ limit of the theoretical rhs in any of these equations corresponds to s → 0, so the rhs must tend to 0 in this limit). It is only eq. 75 that should be amended as

C 2 (q, s) ≃ K(q, s) - q 2 M(q, s) 2 ρs 2 + q 2 K(q, s) + ∆(q) (93)
to be generally valid also in the glassy regime.

Remarkably, the need for the ∆-correction for C 2 below T g was demonstrated by the data shown in Figs. 6,7 (with a vertical shift of the theoretical C 2 in some cases). It is tempting therefore to attribute this correction simply to a non-equilibrium character of our systems below T g . We, however, do not support this explanation because we are confident that our swap-equilibrated systems must be close to the global thermodynamic equilibrium in a wide T -range including T ≥ 0.2 as argued in Ref. 10. If so, what could possibly be the reason for the ∆-shift effect? We see two possibilities: (i) it could be due to a systematic discrepancy between theory and simulations at higher q's (recall that the theory is valid up to a relative correction of O(q b) 2 ), (ii) or it is simply due to stochastic variations of θ (q) (cf. eq. 90) within the ensemble. Below we show that the lion's share of ∆ seems to be due to the second effect. To find the shift ∆(q) we apply (the inverse s-transform of) eq. 93 in the most relevant long-time regime, t > 2τ L :

∆(q) = θ (q) -θ th (q) (94)
where θ (q) is the long-time plateau of C 2 (q,t) corresponding to the low-s limit of C 2 (q, s) (cf.

the lhs of eq. 93), and θ th (q) is the theoretical prediction for θ (q) corresponding to the rhs of eq. 75 at low s. More precisely,

θ (q) ≡ V T σ22 (q) 2 = 2 (∆t max ) 2 ∆t max 0 C 2 (q,t) (∆t max -t) dt (95)
(cf. eq. 90) where ... brackets here mean averaging over many (m) independent systems and N q independent orientations of q (with the same q ). Similar relations involving the volume-and time-averaged shear stress have been employed in Refs. 5, 9, and 12. The mean stress σ is defined

as 51 σ22 (q) = 1 ∆t max ∆t max 0 σ 22 (q,t)dt (96)
The working definition of θ th (q) is given by the rhs of eq. 95 with the only change: the simulation- based C 2 (q,t) there should be replaced with the theoretical C 2th (q,t) defined as the inverse stransform of the rhs of eq. 75. Our analysis shows (see the previous point) that for t > 2τ L the function C 2th (q,t) can be approximated with eq. 92, so we get in analogy with eq. 95: 56

θ th (q) ≃ 2 (∆t max ) 2 ∆t max 0 4G(t) 1 - G(t) K T L (t) (∆t max -t) dt (97)
Thus, the theoretical θ th (q) is actually predicted to be nearly independent of q for small q, q b ≪ 1.

Obviously θ th is closely related to µ s f (cf. eq. 18 in Ref. 10), and, in any case, θ th is independent of any quenched stress, either at q = 0 or at q > 0 (recall that the latter property is inherent in both G(t) and K T L (t) as defined in eqs. ( 2) and ( 16) of Ref. 10). It is therefore likely that the statistical fluctuations of θ (q) and θ th within the ensemble are independent. As a result we can write (cf.

eq. 94):

var(∆(q)) ≃ var(θ (q)) + var(θ th )

for any chosen q and T . Here var(θ (q)) is the variance of θ (q) averaged over the ensemble of m independent systems and N q orientations of q. Assuming Gaussian statistics of σ22 (q) (cf. ref. 9), it is easy to show that the relative standard deviation of σ22 (q) 2 is just 1, which means that the relative standard deviation of θ (q) (which is averaged over mN q independent values) is δ θ (q)/θ (q) ≈ 1/ mN q (99)

where δ θ (q) = var(θ (q)) (here we also used the natural property that σ22 (q) for different wavevectors q are not correlated). The rhs of the above equation for m = 100 and N q = 2 (valid for 0.2 ≤ T ≤ 0.3 and all q's we considered) is 1/ mN q ≈ 0.071. Furthermore, taking into account that K T L (t) is nearly constant at long t and low temperatures, T ≤ 0.24, with K T L /G 5, and using eq. 97 we estimate the relative standard deviation of θ th as

δ θ th /θ th ∼ 1 √ m δ µ s f µ s f (100)
where the shear modulus µ s f and its standard deviation δ µ s f have been calculated at different temperatures in ref. 10. It is more convenient to consider the ratio X (q) ≡ ∆(q)/θ th = θ (q)/θ th -1 as a variable of interest. Taking into account that both δ θ (q)/θ (q) and δ θ th /θ th are small and that θ (q) ≈ θ th , we get (based on eqs. 98)

δ (∆(q)/θ th ) ≈ (δ θ (q)/θ (q)) 2 + δ θ th /θ th 2 1/2 (101)
Finally, we take into account a systematic difference between θ (q) and θ th leading to a non-zero X sys (q) = X . Thus we get

X 2 ≈ X sys (q) 2 + δ (∆(q)/θ th ) 2 (102)
We then calculated the average X 2 based on the simulation results for q = q min , T = 0.2, 0.21, 0.22, 0.23, 024: this way we get the root-mean-square (RMS) of X , RMS ∆(q)/θ th ≈ 0.092.

On the other hand, the similar RMS average of δ (∆(q)/θ th ), calculated using eqs. 99, 100, 101 is δ (∆(q)/θ th ) ≈ 0.08 which is close to δ θ (q)/θ (q) ≈ 1/ mN q . It means that the systematic term X sys is relatively small, and that the observed ∆ can be largely explained by the inevitable statistical error due to a random distribution of σ22 (q). Considering, in addition, all (T, q) combinations with n q2 = q/q min 2 = 8, 32, 64 and T = 0.24, 0.2 (cf. Figs. 6,7) we get RMS ∆/θ th ≈ 0.086, that is a smaller deviation which is closer to the lower bound of RMS(X ), 1/ mN q ≈ 0.071 (the fact that RMS(X ) ≥ 1/ mN q follows from eqs. 99, 101, 102). 57 These results clearly show that the systematic deviation X sys is not resolved here (it is masked by the stochastic part of X = ∆/θ th ): were X sys mattered, its effect would increase with q, while we detected an opposite trend. Therefore, the vertical shift ∆ seems to be mostly a stochastic effect due a limited size of the ensemble.

12. As mentioned in section III B (see also Refs. 6 and 21) the four stress correlation functions analyzed in the present paper define the full time-space resolved tensorial stress correlation function, C αβ α ′ β ′ (q,t), eq. 69, for 2D systems (cf. eq. ( 24) in Ref. 6):

C αβ α ′ β ′ (q,t) = (C 2 (q,t) -2C s (q,t))δ αβ δ α ′ β ′ + (C ⊥ (q,t) -C 2 (q,t) + 2C s (q,t))× q α q β δ α ′ β ′ + q α ′ q β ′ δ αβ /q 2 +C s (q,t) δ αα ′ δ β β ′ + δ αβ ′ δ β α ′ + C (q,t) +C 2 (q,t) -2C ⊥ (q,t) -4C s (q,t) q α q β q α ′ q β ′ /q 4 (103)
In particular, the autocorrelation function of the shear stress (in the fixed coordinate frame, x, y) is 6,21 C xyxy (q,t) = C s (q,t) + C (q,t) +C 2 (q,t) -2C ⊥ (q,t) -4C s (q,t) q 2 x q 2 y /q 4 (104)

The above eqs. 103 and 104 are exactly valid provided that the system is perfectly isotropic. We have already established 10 that bulk properties of our pLJ systems are isotropic to a high precision at T > T f ≈ 0.16 after statistical averaging over a large ensemble. A similar statement concerning the wave-vector q (or distance r) correlation functions is less trivial because here the finite box size (FBS) and the periodic boundary conditions (PBC) can affect correlation functions at large r or low q. We defer to a separate publication a discussion of these and related effects for the stress correlation functions resolved in real and q-space. Here, instead, we limit the discussion to the results for an extended ensemble including (in addition to the original independent configurations) also the 'rotated configurations' obtained by turning the simulation box as a whole over an angle φ with uniform distribution of φ from 0 to 360 • . Such an extended ensemble is perfectly isotropic by construction (that is, all directions are statistically equivalent there), so eqs. 103 and 104 are mathematically exact for it. Obviously, if some anisotropies (like preferential orientation of vectors connecting neighboring particles along the box sides) were present in the original configurations, they get wiped out by the angular averaging implied by the extended sampling. Note, however, that the functions C s (q,t), C 2 (q,t), etc. are obtained here by averaging the corresponding stress correlators over all possible orientations of the wave-vector q (for a given q ), hence they remain the same for the extended ensemble.

At sufficiently long time, t 2τ L , the only surviving correlation is C 2 (q,t). At T < T g the function C 2 (q,t) shows a long-time slowly decaying tail which is nearly independent of q (cp. Figs. 4(d), 6(a), (b) and (c)). Therefore our direct simulation data point to a non-analytical q-dependence of the shear stress correlation function at long times (and small q) according to C xyxy (q,t) ∝ q 2 x q 2 y /q 4 , which indicates that correlations of shear stress in real space are long-range (see below). Moreover, using the theoretical eqs. 75, 92 we can predict the time-dependence of the shear-stress correlations (for t 2τ L ) based on the bulk relaxation moduli:

C xyxy (q,t) ≈ 4G(t) 1 - G(t) K T L (t) q 2
x q 2 y /q 4 (105)

Doing its inverse Fourier transform to real space (cf. Refs. 6, 20, and 21) we get for r ≪ L (the latter condition is needed to avoid the FBS effects):

C xyxy (r,t) ≈ - 1 π 1 r 2 cos(4θ )G(t) 1 - G(t) K T L (t) , t ≫ τ L ∼ 2ρ η L r 2 ( 106
)
where θ is the polar angle between the distance vector r and the x-axis. It means that for any t, τ v ≪ t τ α , there is a wide range of distances where the 1/r 2 power law, eq. 106, is valid (recall that τ v is the microscopic time for particle velocity relaxation by collisions and vibrations).

This result is in harmony with the theories 6,20,21 predicting the 1/r d power-law stress-correlation tail for isotropic glass-forming liquids. Note that strictly speaking eq. 106 is justified here for an extended ensemble described above.

13. In the present study we neglected a q-dependence of the generalized shear relaxation modulus G(q,t), cf. eq. 19. By contrast, some theoretical and simulation studies [START_REF] Semenov | [END_REF][START_REF] Puscasu | [END_REF][41][42][43] report a rather significant q-dependence of G(q,t). There are two main reasons for this effect: (i) It is present in polymer systems [START_REF] Semenov | [END_REF][START_REF] Puscasu | [END_REF]43 where it is simply due to a large size of macromolecular coil: the size b is then b = R coil √ N m l m , where N m ≫ 1 is the polymerization degree and l m is monomer size. (ii) It was also detected in glass-forming liquids near the glass transition, T g . [41][42][43] More precisely, they considered a q-dependence of the generalized shear viscosity η(q) which is related to G(q,t): η(q) = ∞ 0 G(q,t)dt (cf. eq. 87). Quite naturally they found [41][42][43] that η(q) becomes significantly lower than η(0) (say, by a factor of 2) at sufficiently large q ∼ 1/ξ η , where ξ η is the viscosity-based correlation length. Still, at lower q's, qξ η ≪ 1, the q-dependence of η(q) can be neglected. The simulation studies 41,42 deal with a 3D binary mixture of soft spheres and show that the length ξ η ∼ b m /3 (b m is molecular size) is low well above T g , but it increases to ξ η ∼ 4b m on cooling below T g . Simulations of liquid butane 43 reveal a similar trend with ξ η /b m ∼ 3 ÷ 5.

We emphasize that the systems studied in refs. 41-43 are quite different from our system which is both polydisperse and two-dimensional. To assess the q-dependence of G(q,t) for our system we chose to consider the correlation function C 2 (q,t) well below the glass transition, at the lowest T = 0.2 (where the q-effect is expected to be stronger). In this regime C 2 (q,t) shows a long-time plateau, C 2 (q,t) ≈ θ (q) (cf. eq. 90) which, by virtue of eq. 75, is related to the plateau levels, G s (q) and K s (q) of G(q,t) and K(q,t), respectively: θ (q) ≈ 4G s (q) [1 -G s (q)/K s (q)] (cf. eq. 91).

Therefore, the q-dependence of θ is expected to provide an idea on the behavior of the long-time shear relaxation modulus. Somewhat surprisingly, our simulation data point to θ (q) ≈ 45 which is nearly independent of q in a rather wide q-range, for q b < 2.5 (recall that the data considered in the present paper correspond to a more narrow range, q b < 0.5). Furthermore, we found that G(q,t) at t = 0 also stays virtually independent of q for q b < 0.5 (this conclusion is based on the stress-correlation data and the exact relation, G(q,t) = C s (q,t) at t = 0, which comes from eq. 70). Thus, for the 2D pLJ system considered here we do not observe any indication of a qdependence of G(q,t) in the studied q-range, q b 0.5. This statement is also backed by a very good agreement between the shear stress correlation function C s (q,t) obtained by simulations and its theoretical prediction based on the exact eq. 70 and involving the approximation G(q,t) ≈ G(t) (cf. Figs. 4567).

V. SUMMARY AND CONCLUSIONS

1. In the present paper we analyzed spacio-temporal correlations of mechanical stress and other variables both theoretically and by MD simulations of a two-dimensional polydisperse system of LJ particles (2D pLJ) described in Appendix A. The theory is based on the FDT relations, eqs. 62, 70, 71, 72, 75, between the wave-vector (q) dependent correlation functions and the generalized relaxation moduli, G(q,t) and K(q,t), for shear and longitudinal stress, respectively. We argue that at low q (when the wave-length λ = 2π/q strongly exceeds the particle size b) the modulus G(q,t) can be approximated with the classical bulk shear relaxation modulus G(t), that is, a q-dependence of G(q,t) can be neglected at low q (cf. section II C). By contrast, we show that a similar approximation for the longitudinal modulus, K(q,t) (which was adopted before 6,20,21 ) is not valid in practice since the bulk modulus K(0,t) = K L (t) is not universal: it depends on the thermostat employed in simulations 10 , hence K(q,t) may be discontinuous at q = 0 (cf. section II C).

To overcome this issue a theoretical approach was developed (cf. section II D) which allows to obtain K(q,t) in terms of the bulk modulus K L (t), bulk correlation functions of energy and pressure, the particle MSD, h 0 (t), the static structure factor S(q) and the heat conductivity κ. Our theory is akin to the generalized hydrodynamic approach [START_REF] Götze | [END_REF][START_REF] Kadanoff | [END_REF] .

2. We elucidate that the generalized longitudinal modulus K(q,t) at q → 0 cannot be uniformly approximated by a function of time only, and so, in particular, lim q→0 lim t→∞ K(q,t) < lim t→∞ lim q→0 K(q,t). Moreover, in polydisperse systems the first limit can be much smaller than the second, and none of these limits is related to the isothermic bulk compressibility of the system.

The time dependence of K(q,t) for a low q above T g is illustrated in Fig. 8. It involves a strong decrease of K(q,t) at short times, t τ v (with τ v , the vibration/collision time), due to an initial relaxation of the system structure and particle velocities, a further structural relaxation (at nearly constant energy density) followed by an adiabatic shoulder (quasi-plateau), then a further decrease of K(q,t) due to heat transport with characteristic time τ T = n 0 c vs / κq 2 defined by the thermal conductivity κ and wavevector q (the thermal diffusion stage). The thermal process is followed by the isothermic shoulder with K(q,t) ≈ K T L (t) ≈ K T Ls , and the terminal relaxation by particle interdiffusion leading a significant decrease of K(q,t) from the bulk longitudinal isothermic modulus K T L to the equilibrium level K(q, ∞) = n 0 T /S(q) defined by the static structure factor S(q). The characteristic time of the latter process is τ inter = 1/ D inter q 2 , where the interdiffusion constant D inter is close to the particle self-diffusivity D s (for the 2D pLJ system we studied, see Appendix A).

3. The terminal interdiffusion stage for K(q,t) is clearly reflected in the dynamical structure factor, S(q,t), for T = 0.3 shown in Fig. 3. This figure also demonstrates a very good agreement between the direct simulation data for S(q,t) and the theoretical results based on the predicted K(q,t) and the general relation, eq. 62. The proposed theoretical approach to calculate K(q,t) is thus validated numerically. It is noteworthy that well below T g the structure factor S(q,t) exhibits a long-time plateau, S(q,t) ≈ S(q)n 0 T /K T Ls for t 2τ L (see Fig. 9), which corresponds to the isothermic quasi-plateau of the longitudinal modulus, K(q,t) ≈ K T Ls (this plateau occurs before the interdiffusion stage which is not visible since the sampling time ∆t max ≪ τ inter ). The longtime plateau of S(q,t) is related to the polydispersity of the studied system, as discussed in Ref. 36.

4. The obtained theoretical longitudinal modulus K(q,t) was used to predict all components of the stress correlation tensor, C αβ α ′ β ′ (q,t), based on the general relations, eqs. 70, 71, 72, 75. The results are compared with our simulation data on the 2D pLJ system in Figs. 4, 5, 6, 7. Most of the predicted results are in excellent agreement with the simulation data for q/q min ≤ 8, i.e. q b 0.5 (see Figs. 4, 5, 6). The agreement gets slightly worse well below the glass transition temperature T g for the highest q's we studied.

5. The theoretical approach to calculate the time-space resolved stress correlations 6,20,21 was developed for the monodisperse isothermal systems (with infinite thermal conductivity κ). In the present paper the theory is generalized to the case of a polydisperse system with finite κ. The generalized theory we developed is asymptotically exact for q → 0; at finite q it works within a relative error of O(q b) 2 .

6. As demonstrated in Figs. 4, 5, 6, 7, the basic stress-correlation functions, C s , C , C ⊥ , vanish beyond the sound attenuation time τ L . By contrast, the autocorrelation function of the transverse normal stress, C 2 (q,t), shows a shoulder at t > τ L near and below T g (the shoulder transforms to a plateau well below T g ). As discussed in section IV this feature leads to long-range dynamical correlations of the local shear stress σ xy = σ xy (r) according to the power law, 1/r 2 for 2D systems (where r = |r 1r 2 | is the distance between two points) in agreement with theoretical 6,20,21 and simulation 22,23 results. Note that this power law is justified in the present paper for an extended ensemble described in point 12 of section IV. The long-range correlations characterize the (transiently) quenched stress field in a glassy system. 7. The plateau of C 2 (q,t) persists until t ∼ τ α . Interestingly, our stress correlation data indicate that the terminal relaxation of C 2 (q,t) gets slower at higher q's, q b 0.5, with relaxation time longer than τ α based on G(t). This effect becomes stronger as q increases up to q b ∼ 2.

8. Below the glass transition, T < T g , we detected a small discrepancy between the predicted C 2 (q,t) (cf. the rhs of eq. 75) and the direct simulation data on C 2 . This discrepancy, ∆(q), takes positive or negative values depending on T and q. However, as we demonstrated numerically (cf. point 11 of the Discussion), its relative root-mean-square value (denoted as RMS(∆(q)/θ th )) is largely accounted for by the statistics of the time-averaged quenched stress whose relative dispersion is defined by the size of the statistical ensemble. Therefore, the vertical shift ∆ is mostly a stochastic effect due a limited size of the ensemble, which is constant for the considered temperature range. Other contributions to the RMS(∆(q)/θ th ), including that due to the systematic error of the predicted C 2 (q,t), appear to be small. We anticipate, however, that the latter (systematic) contribution may become significant for large wave-vectors (q/q min ≫ 8).

By virtue of the FDT (i.e., the linear response theory 3 ) the correlation function C T is related to the temperature response ∆T q (t) = T q (t) -T q (0 -) induced by a heat wave generated by the following perturbation of particle velocities applied at t = 0 (cf. eq. 4):

v → v 1 + exp iq • r ε v /d (B5) Namely, ∆T q (t)/ε v = TC T (q,t), t > 0 (B6)
where T = T 0 and ε v → 0. Using the above equation one immediately gets

C T (q, 0 + ) = 2/d (B7)
The time dependence of C T can be obtained for low q's using the generalized hydrodynamic approach (cf. section II D). The longitudinal pressure increment we get (cf. eq. 31):

∆p xx (s) = K AL (s)ε(s) + χ(s)h(s) (B8)
where ε xx = -ε is the longitudinal strain eventually generated by the initial heat wave, and h = h i + h t is the heat per particle including the initial heat h i = ε v T injected in the system at t = 0 (as follows from eq. B5) and the heat h t transferred into a cell by thermal conductivity. Note that here, as before, we assume that the wave-vector q is oriented along the x-axis, and the phase factor e iqx is omitted. The time-dependence of h t is given by eq. 32, where the temperature increment ∆T is (cf. eq. 36):

∆T (s) = T n 0 χ(s)ε(s) + h(s) c v (s) (B9)
Another relation defining the strain ε = ε(t) in terms of pressure p xx is given by two equations corresponding to mass and momentum conservation:

∂ n ∂t = -∇ • j, mn 0 ∂ v x ∂t = - ∂ p xx ∂ x = -iq∆p xx (B10)
where n = n(x,t) = n 0 (1 + ε(t)e iqx ) and j = nv ≃ n 0 v. Using the above equations we find ∆T (s) ≡ ∆T q (s) and finally, on using eq. B6, the s-transform of C T (q,t):

C T (q, s) = c v (s) + κq 2 n 0 s + T n 0 p T (s) 2 K T L (s) + mn 0 s 2 /q 2 -1 (B11)
This equation agrees with eq. ( 18) of Ref. 54 provided that the bulk modulus K B there is identified with the isothermic longitudinal modulus (K T L ). Eq. B11 can be also obtained in a simpler heuristic way following the approach outlined in point 2 of the Discussion.

The time-dependence of C T (q,t) obtained directly from our simulation data typically shows a decay superimposed with decaying oscillations. Such a behavior also follows from eq. B11. Doing inverse s-transform and applying approximations similar to those that led to eq. 47 we obtain

C T (q,t) ≈ 1 c pL (t) (r(t) -1)cos (ω(t)t)exp -q 2 tΓ(t) + exp -q 2 t κ n 0 c pL (t) (B12)
where

c pL (t) ≈ c v (t)r(t), r(t) = K AL (t)/K T L (t), ω(t) ≈ K AL (t) mn 0 1/2 , Γ(t) ≈ 1 2mn 0 t 0 K AL (t ′ ) -K AL (t) dt ′
Here c pL (t) is the heat capacity at constant longitudinal normal pressure p xx and constant transverse area (i.e. at no deformation in the plane perpendicular to the x-axis, cf. eq. 77). Note that eq. B12 at t = 0 agrees with eq. B7 since c v (t = 0) = d/2. 11 Eq. B12 for C T (q,t) was used to obtain the thermal conductivity κ at different temperatures by fitting it to the simulation data (the approximations like c v (t) ≈ c v0 (t), valid because the relevant time is sufficiently long, t ≫ τ damp , were used to this end). The results are quoted in section II E. are insensitive to small T -variations.

FIGURE CAPTIONS

FIG. 1. Longitudinal relaxation modulus K(q min ,t) for T = 0.4, 0.3, 0.26, 0.24, 0.2 calculated using eqs. 43, 45 (black curves from bottom to top) and its approximation, eqs. 47, 43 (red curves). Note that the black curves superimpose almost completely onto the corresponding red curves at t 100. FIG. 2. Longitudinal modulus K(q,t) (black curve), adiabatic modulus K AL (t) (red curve) and isothermic modulus K T L (t) (blue curve). Upper panel: T = 0.4, q = q min ; lower panel: T = 0.3, q = √ 8q min .

FIG. 3. The dynamical structure factor S(q,t) vs. t for T = 0.3, q = q min (upper left panel), T = 0.3, q = √ 8q min (upper right), T = 0.26, q = q min (lower left), T = 0.24, q = √ 8q min (lower right). Simulation results based on eq. 61 are shown with red curves; theoretical results, eqs. 68, 67 -black curves.

FIG. 4. Time dependence of correlation functions for shear stress, C s (q,t) (upper left panel); longitudinal stress, C (q,t) (upper right); normal stress cross-correlations, C ⊥ (q,t) (lower left); second normal stress, C 2 (q,t) (lower right) for q = q min ≡ 2π/L. Black curves -theory, red symbols and curves -simulation data. We show symbols only for short times (t ≤ 6) for clarity.

The temperature T = 0.4 (a), 0.3 (b), 0.26 (c), 0.24 (d), 0.2 (e).

FIG. 5. Same functions (as in Fig. 4) for T = 0.3 and n q2 ≡ q/q min 2 = 8 (a), n q2 = 32 (b), 64 (c).

FIG. 6. Same functions (as in Fig. 4) for T = 0.24 and n q2 ≡ q/q min 2 = 8 (a), showing all relaxation stages as described in the text.

FIG. 9. The dynamical structure factor S(q,t) for T = 0.2, q = q min . Simulation results based on eq. 61 (red curve), theoretical results, eqs. 68, 67 (black curve). The blue horizontal line indicates the long-time isothermic plateau level = S(q)n 0 T /K T Ls .

FIG. 10. The dynamical structure factor S(q,t) based on simulation data (black line), its exponential approximation, eq. 52 (blue curve), and the correlation function of composition fluctuations, 1.33S ψ (q,t), obtained by simulations (red curve). T = 0.4 and q = √ 8q min for all the curves.

The 4th relation involving C 2

 2 was derived in Refs. 6, 20, and 21 for monodisperse isothermic systems (i.e., for the limit κ → ∞). Following the argument of Ref. 6 (based on the concept of structural stress noise) we show that in the general case (including polydisperse systems) the relation between C 2 and relaxation moduli reads

Figs. 4 , 5 , 6 , 7 .

 4567 Figs.4, 5, 6, 7. At T = 0.4, 0.3 (above the glass transition) there is a complete agreement for all correlation data (the black theoretical curve is virtually invisible at t > 6 being completely overlapped by the red data curve, see Fig.4(a),(b)). The shear stress correlations C s show a few decaying oscillations reflecting transient shear acoustic waves. As for the correlations of normal stresses (C , C ⊥ , C 2 ), they show much more pronounced oscillations which decay more slowly. These oscillations correspond to longitudinal sound waves. Their decay time (∼ 100 ÷ 200) increases as the temperature is lowered. Turning to the autocorrelation function C 2 , we note a remarkable feature emerged at T = 0.3: the oscillations are now significantly asymmetric pointing to a monotonic slowly decaying process. This slow process is even more apparent at T = T g = 0.26, where it persists beyond the oscillation relaxation time. The relaxation time of the slow process at T = 0.26

  above equation involve s-transforms of the time-dependent functions; it was obtained in analogy with the argument considered in section II D. Therefore c pL (s)/c v (s) = K AL (s)/K T L (s); this equation is analogous to the more standard relation c p (s)/c v (s) = K A (s)/K T (s) (cf. eqs. (D23), (D25) in Ref. 10).

FIG. 7 .

 7 FIG.6. Same functions (as in Fig.4) for T = 0.24 and n q2 ≡ q/q min 2 = 8 (a), n q2 = 32 (b), 64 (c). The theoretical data for C 2 are vertically shifted by ∆ = 1.9 for n q2 = 32, and ∆ = 2 for n q2 = 64.

  Fig. 1

  Fig.2

  Fig.4 (a)

  Fig.9

Fig. 10 S

 10 Fig. 10

DATA AVAILABILITY

The data that support the findings of this study are available from the corresponding author upon reasonable request.

Conflicts of interest

There are no conflicts of interest to declare.

ACKNOWLEDGMENTS

A grant of computer time at the HPC computing cluster of the University of Strasbourg is gratefully acknowledged. L.K. was supported by a doctoral contract from the University of Strasbourg in the framework of the IRTG "Soft Matter Science".

Appendix A: 2D pLJ system and simulation algorithm

We study a 2-dimensional (2D) system of polydisperse Lennard-Jones (pLJ) particles with the same mass m. 10,11,17,[START_REF] Klochko | [END_REF]48 The particle diameters are uniformly distributed around the mean diameter b, between 0.8 b and 1.2 b. The particle size polydispersity index is b 2 / b2 -1 ≈ 0.013. Each pair of particles (l = (a, a ′ )) interacts with the energy u LJ (s = r/b l ), where the interaction range b l is defined by the Lorentz rule, b l = (b a + b a ′ ) /2, and the LJ potential u LJ (s) = 4ε s -12s -6 was truncated at s cut = 2 7/6 and shifted to avoid discontinuity at s = s cut . All quantities are given in LJ units: m, b, ε and the Boltzmann constant k B are set to 1. We used the standard molecular dynamics (MD) simulations (velocity-Verlet algorithm with MD time-step t MD = 0.005) with periodic boundary conditions as implemented in the LAMMPS code. 49 The system of N = 10 4 particles was first tempered well in the liquid regime, at constant temperature T 0 = 1 and pressure p 0 = 2 using the Nosé-Hoover thermostat and barostat to prepare m = 100 well-equilibrated independent configurations which were then continuously cooled down with rate dT /dt = -10 -5 at p 0 = 2. These slow cooling runs allowed us to estimate the dilatometric glass-transition temperature, T g ≈ 0.26, in agreement with the previous studies 17 .

To better equilibrate the quenched configurations at T ≤ 0.5 we used a hybrid approach involving Monte-Carlo (MC) dynamics. The MC part involves a combination of local particle displacements and non-local particle swaps. [START_REF] Ninarello | [END_REF] The system volume fluctuations (controlled by an MC barostat [START_REF] Frenkel | Understanding Molecular Simulation[END_REF] to impose p 0 = 2) were allowed as well. This way an ensemble of m independent configurations (m = 100 for 0.2 ≤ T ≤ 0.3, m = 50 for T > 0.3, m = 20 for T < 0.2) was tempered over 10 7 MC steps at constant pressure. Then the instantaneous volume was fixed and the system equilibrated further at constant volume over 10 7 MC steps (with local and swap moves), and over the same time with local moves only. This MC approach was successful in equilibrating the pLJ system below T g , down to T f = 0.16 below which particle demixing occurs. 10 Each configuration was then tempered for 2 • 10 5 LJ time units with MD Nosé-Hoover (NH) dynamics at constant pressure p 0 = 2 to equilibrate the velocities, and then over the same time at constant volume.

Finally, the production runs (served to obtain all correlation functions) were performed in the NVT ensemble during the total sampling time ∆t max = 10 5 . The time space between the data entries was δt = 0.2 for q-dependent variables, and δt = 0.05 for all other variables. The thermal mass parameter Q of the NH thermostat was Q = 10T /3 corresponding to a high Nosé frequency ω Q = (2T Nd/Q) 1/2 ∼ 110 and a short thermostat relaxation time (the Nosé time)

The linear dimension of the simulation box is L ∼ 100, the system volume V = L 2 , and the mean concentration n 0 = N/V ∼ 1 in all the cases.

Appendix B: Temperature autocorrelation function and the heat conduction coefficient

The q-dependent temperature autocorrelation function can be defined as

where ... imply the ensemble and time-shift averaging (over t ′ ), and

is a Fourier component of the temperature fluctuation field (the index a runs over all particles).

Here

is the excess kinetic temperature of particle a, T 0 is the mean temperature of the system (which is prescribed by the thermostat), and d is the space dimension. Note that δ T a (t) = 0 and T q (t) = 0, hence δ T q (t) = T q (t) for q = 0. Note also a direct way to calculate C T (q,t) based on simulation data:

which is analogous to the standard definition of the dynamical structure factor S(q,t) (see eq. 61).