Strain correlation functions in isotropic elastic bodies: large wavelength limit for two-dimensional systems

Abstract

Strain correlation functions in two-dimensional isotropic elastic bodies are shown both theoretically (using the general structure of isotropic tensor fields) and numerically (using a glass-forming model system) to depend on the coordinates of the field variable (position vector r in real space or wavevector q in reciprocal space) and thus on the direction of the field vector and the orientation of the coordinate system. Since the fluctuations of the longitudinal and transverse components of the strain field in reciprocal space are known in the long-wavelength limit from the equipartition theorem, all components of the correlation function tensor field are imposed and no additional physical assumptions are needed. An observed dependence on the field vector direction thus cannot be used as an indication for anisotropy or for a plastic rearrangement. This dependence is different for the associated strain response field containing also information on the localized stress perturbation.

---

I. Introduction

A. General background

A tensor field assigns a tensor to each point of the mathematical space, in our case for simplicity a two-dimensional Euclidean vector space with Cartesian coordinates and an orthonormal tensor basis.^1–4 Tensor fields are used in differential geometry,¹ general relativity,^5,6 in the analysis of stress and strain in materials^7–9 and in numerous other applications in science and engineering. Tensor fields are experimentally^10,11 or numerically^4,12–32 probed by means of correlation functions^33–35 of their components and, importantly, these correlation functions are themselves components of tensor fields.⁴ See Appendix A for a brief review. Assuming translational invariance, correlation functions are naturally best analyzed, both for theoretical^12–14,17 and numerical^4,18,33 reasons, in a first step as functions of the wavevector q in reciprocal space. The dependence on the spatial field vector r in real space can then be deduced (cf. Appendix B) in a second step by inverse Fourier transformation (FT). This was done, e.g., in our recent analysis⁴ of the spatial correlations of the (time-averaged) stress tensor fields in amorphous glasses formed by polydisperse Lennard-Jones (pLJ) particles deep in the glass regime (cf. Section IIIA). It can thus be shown that all stress correlation functions (both in reciprocal as in real space) can be described by means of one “Invariant Correlation Function” (ICF) in reciprocal space characterizing the typical ensemble fluctuations of the quenched normal stress components in reciprocal space perpendicular to the wavevector q. Under additional but rather general assumptions⁴ this ICF is given in the large-wavelength limit by a thermodynamic quantity, the equilibrium Young modulus of the system.

B. Investigated case study

As another example of the general procedure we shall investigate in the present work the correlation functions of the instantaneous strain tensor field ε_αβ(r) in real space. These may be readily obtained³³ from the components of the tensor field

c_αβγδ(q) = 〈ε_αβ(q)ε_γδ(−q)〉 (1)

in reciprocal space with being the Fourier transformed strain tensor field components. (The average 〈…〉 will be specified below) An example for the autocorrelation function c₁₂₁₂(r) of the shear strain ε₁₂(r) is given in Fig. 1 for the same two-dimensional model system already used in ref. 4 and 18. Interestingly, the correlation function is seen to strongly depend both on the orientation of the field vector r (panel (a)) and on the rotation angle α of the coordinate system (panel (b)). Since the simulated system can be shown to be perfectly isotropic down to a few particle diameters,^4,18,36–39 these findings beg for an explanation. Expanding on our recent work on stress correlations,^4,17,18 this behavior can be traced back to the fact that correlation functions of tensor fields of isotropic systems must be components of a generic isotropic tensor field (cf. Section IIB). This field is shown below (cf. Section IIID) to be completely described in terms of two ICFs c_L(q) and c_T(q) in reciprocal space (q = |q| being the magnitude of the wavevector). These ICFs characterize the independent fluctuations of the longitudinal and transverse strain components ε_L(q) and ε_T(q). Due to the equipartition theorem of statistical physics c_L(q) and c_T(q) are given by^7,10,11,21in the large-wavelength limit with β = 1/k_BT being the inverse temperature, V the d-dimensional volume of the system and λ and μ two macroscopic Lamé coefficients.^7,8 All strain correlation functions are thus imposed on large scales. In turn this explains without any additional physical input the octupolar pattern† observed in Fig. 1 (cf. Section IVC and Appendix D) and shows that strain correlations in elastic bodies must necessarily be long-ranged. This is different for the closely related but distinct tensorial response field being the tensorial product of correlation functions and the imposed tensorial perturbation. As emphasized in Sections IIE and V, the response field thus contains additional information due to the source term and its symmetry.

Fig. 1 Autocorrelation function c₁₂₁₂(r) of the strain field component ε₁₂(r) obtained from our colloidal glasses in two dimensions: (a) unrotated frame with coordinates (r₁, r₂), (b) frame rotated by an angle α = 30° (rotations marked by “′”). Albeit the system is isotropic, the correlation function is strongly angle dependent, revealing an octupolar symmetry. While each pixel corresponds in (a) and (b) to the same spatial position r, the correlation functions differ by the angle α. c₁₂₁₂(r) is positive (red) along the axes and negative (blue) along the bisection lines of the respective axes.

C. Outline

We begin in Section II with some general theoretical considerations on isotropic tensor fields. Technical points concerning the model system and the data production of tensorial fields on discrete grids are discussed in Section III. This is followed in Section IV by the presentation of our main numerical results. The strain response due to an imposed stress point source is discussed in Section V. A summary and an outlook are given in Section VI. More details may be found in the Appendix both on the theoretical background (cf. Appendices A and D) and on computational issues (cf. Appendices B and C).

II. General considerations

A. Isotropic tensors and tensor fields

Isotropic systems, such as generic isotropic elastic bodies,^7–9 simple and complex fluids,^34,40–42 amorphous metals and glasses,^{23–25,27–32,43} polymer networks and gels,^40,41 foams and emulsions^20,22 or, as a matter of fact, our entire universe⁵ are described at least on some scales by isotropic tensors and isotropic tensor fields (cf. Appendix A2).^1,3,9 It is well known^3,9 that the components of isotropic tensors remain unchanged under an orthogonal coordinate transformation (including rotations and reflections). For instance,for the forth-order elastic modulus tensor of an isotropic body (cf. Appendix C3)^8,9 with “*” marking an arbitrary orthogonal transformation (cf. Appendix A1). This implies (cf. Appendix A4) that E_αβγδ is given by two invariants, e.g., the two Lamé coefficients λ and μ. Importantly, this does not hold for isotropic tensor fields.^3,4,17 For instance, for a forth-order correlation function in reciprocal space the isotropy condition becomeswith q* being the “actively” transformed wavevector (cf. Appendix A2).

B. Structure of isotropic correlation functions

Assuming in addition the system to be achiral and two-dimensional (cf. Appendix A3) it can be shown⁴ that correlation functions of second-order tensor field components must take the following mathematical structurein terms of four ICFs i_n(q), the coordinates q̂_α of the normalized wavevector q̂ and the Kronecker symbol δ_αβ. Legitimate correlation functions of isotropic systems may thus depend on q̂_α and, hence, on the orientation of the wavevector and of the coordinate system. While the isotropy of the system may not be manifested by one correlation function, it is crucial for the structure of the complete set of all correlation functions given by eqn (5). We note finally that it is useful to express the above ICFs in terms of an alternative set of ICFs c_L(q), c_T(q), c_⊥(q) and c_N(q) given bySee Appendix A4 for more details.

C. Planar harmonic basis functions

Instead of using the components q̂_α one may, quite generally, express all isotropic tensor fields in two dimensions in terms of the orthogonal planar harmonic basis functions cos(pθ) and sin(pθ) with q̂₁ = cos(θ) and q̂₂ = sin(θ) and p = 0, 2 and 4. (See Appendix D for more details.) For instance, it follows from eqn (5) thatHence, if the invariant i₄(q) is sufficiently large, c₁₂₁₂(q) must reveal an octupolar pattern. Due to eqn (B18) derived in Appendix B3, this alternative representation is especially useful for performing the inverse FT to real space. This also shows that the corresponding correlation function in real space must have the same mathematical properties.

D. Response to point source

Let us consider the second order tensor field R_αβ(q) obtained by the contractionwith a symmetric but not necessarily isotropic tensor s_αβ using the standard summation convention over repeated indices.^1,3 (For convenience we have introduced the system volume V.) We shall call R_αβ(q) the “response field” (in reciprocal space) and s_αβ the “point source tensor”. In fact, using eqn (B4) and (B6) it is seen that in real space the tensor s_αβ/V corresponds to a “point source” s_αβδ(r) (using Dirac's delta function) and R_αβ(q) becomesusing . We shall say more about the specific linear strain response in real space in Section V but focus here on the generic response in reciprocal space. Being symmetric the source tensor may be diagonalized by a rotation of the coordinate system where s₁₂ = s₂₁ = 0 and s₁₁ and s₂₂ become the two (in general not identical) eigenvalues. Hence,We emphasize that the sum must be taken over all eigenvalues of the source tensor, i.e. two for the presented two-dimensional case. (The failure to sum properly over all tensorial contributions to R_αβ(q) leads to incorrect angular dependences.) Importantly, R_αβ(q) thus contains information over both the system, characterized by the correlation functions, and the imposed source term.

E. Different types of source terms

If we now assume that not only c_αβγδ(q) is an isotropic tensor field but that, moreover, s_αβ is isotropic, i.e. s₁₁ = s₂₂, the product theorem eqn (A7) discussed in Appendix A2 implies that R_αβ(q) must also be an isotropic tensor field. According to eqn (A15) it is given by

R_αβ(q) = k₁(q)δ_αβ + k₂(q)q̂_αq̂_β (11)

in terms of two invariants k₁(q) and k₂(q) which can in turn be expressed in terms of the invariants of c_αβγδ(q) and s_αβ. R_αβ(q) can thus at most be quadrupolar (p = 2). Specifically,

R₁₂(q) = k₂(q)q̂₁q̂₂ ∝ sin(2θ) (12)

which is distinct from c₁₂₁₂(q), cf.eqn (7). Importantly, in many physical situations the source is in fact not isotropic and thus in turn the response field not consistent with eqn (11). We remind that according to a popular model of localized plastic failure by means of “shear transformation zones”^23,44–47 two orthogonal twin force dipoles of opposite signs may be imposed at the origin.‡ This suggests to consider the case s₁₁ = −s₂₂. It follows then from eqn (5) and (10) that

R₁₂(q) ∝ i₄(q)q̂₁q̂₂(q̂₁² − q̂₂²) ∝ sin(4θ). (13)

The (non-isotropic) response field R₁₂(q) thus is in this case octupolar as the correlation field c₁₂₁₂(q), however, shifted by an angle π/8. It is readily seen by inverse FT that the same general behavior applies in real space.

III. Technical issues

A. Algorithm, configurations and frames

We investigate amorphous glasses in two dimensions formed by pLJ particles^{4,18,36–39,48} which are sampled by means of Monte Carlo (MC) simulations.³³ See Appendix C for details (Hamiltonian, units, cooling and equilibration procedure, data production, generation and analysis of tensor fields on discrete grids). We focus on systems containing n = 10 000 and 40 000 particles at a working temperature T = 0.2. This is much lower than the glass transition temperature T_g ≈ 0.26,³⁶i.e. for any computationally feasible production time the systems behave as solid elastic bodies.³⁸N_c = 200 completely independent configurations c are prepared using a mix of local and swap MC hopping moves^38,49 while the presented data are computed using local MC moves only. For each c we store time-series containing N_t = 10 000 “frames” t computed using equidistant time intervals. As described in Appendix C3, the elastic modulus tensor is isotropic and determined by the two Lamé coefficients λ ≈ 38 and μ ≈ 14.^36,38

B. Sampled discrete tensorial fields

As shown in Fig. 2, a discrete square grid is used to store and to manipulate the various fields needed for the microscopic description. The standard lattice constant for the grid in real space is a_grid ≈ 0.2. The displacement field u(r) in real space is determined for each frame t using a standard method^10,11,21 from the displacement vector of each particle using as reference position the time-averaged particle position (cf. Appendix C5). We obtain then from the Fourier transformed displacement field the strain tensor field^8,17in reciprocal space. Using the correlation function theorem for FTs (cf. Appendix B) the strain correlations functions in reciprocal space are given by eqn (1) where the average is taken over all t and c. Both strain and correlation function fields in reciprocal space are defined to be dimensionless (cf. Appendix B1). We emphasize by a prime “′” all tensor field components obtained in a coordinate system rotated by an angle α (with α = 0 being the original unrotated system). Specifically, the correlation functions are obtained using the components and in the rotated frame.

Fig. 2 Two-dimensional (d = 2) square lattice with a_grid being the lattice constant and n_L = L/a_grid the number of grid points in one spatial dimension. The filled circles indicate microcells of the principal box, the open circles some periodic images. The spatial position r of a microcell is either given by the r₁- and r₂-coordinates (in the principal box) or by the distance r = |r| from the origin (large circle) and the angle θ.

C. Natural Rotated Coordinates

All the tensorial fields introduced above depend on the orientation of the coordinate system. Importantly, we consider these properties in a first step in “Natural Rotated Coordinates” (NRC) where for each wavevector q the coordinate system is rotated until the 1-axis coincides with the q-direction. We mark these new tensor field components by “°” to distinguish them from standard rotated tensor field components (marked by primes “′”). Note that for all wavevectors q. Using the components and we obtain (as before) the strain tensor . Importantly,in agreement with eqn (14). We thus only have two independent components of the strain tensor field in NRC. We alternatively write for convenience , , and for the longitudinal and transverse components of the displacement and strain tensor fields. Note thatandi.e. displacement and strain fields in NRC contain essentially the same information.

D. Correlation functions in NRC

The correlation functions may for finite N_c not only depend on q but also on q̂. Consistently with eqn (A18) and ref. 4, 17 and 18 we thus operationally define the ICFs , , and by averaging over all wavevectors with |q| ≈ q. However, implies immediately thatThe two remaining non-trivial ICFs c_L(q) and c_T(q) are called, respectively, the “longitudinal ICF” and the “transverse ICF”. As already noted in the Introduction, according to the equipartition theorem c_L(q) and c_T(q) are given for sufficiently large wavelengths by the Lamé coefficients λ and μ. The stated eqn (2) can be readily obtained from published work^7,10,11,21 using eqn (16) and (17) to substitute the displacement fields u_L(q) and u_T(q) in NRC by the corresponding strain fields ε_L(q) and ε_T(q).

IV. Main numerical results

A. Measured longitudinal and transverse ICFs

We turn now to the numerical results of this work. Fig. 3 focuses on the two non-vanishing correlation functions obtained in reciprocal space and NRC. All correlation functions are rescaled by βV having thus the dimension of an inverse modulus. As can be seen for the two indicated particle numbers n, a data collapse for different system sizes is observed, confirming the expected volume scaling. The inset presents the (not yet spherically averaged) correlation functions and as functions of the wavevector angle θ for one small wavevector with q ≈ 0.1. As expected for isotropic systems, these correlation functions are θ-independent (apart from a small noise contribution due to the finite number of independent configurations). The main panel presents the q̂-averaged longitudinal and transverse ICFs βVc_L(q) and βVc_T(q) as functions of q. The expected large-wavelength limit eqn (2) is indicated in both panels by bold horizontal lines. As can be seen in the main panel, it is well confirmed for q ≪ 1 over at least one order of magnitude where we have used the known values of λ and μ. We remind that eqn (2) has been used in various experimental and numerical studies^10,11,21 to fit λ and μ. c_L(q) and c_T(q) characterize the typical length of the complex random variables ε_L(q) and ε_T(q). Their distributions and correlations will be discussed elsewhere.⁵⁰ We note finally that the increase of the ICFs from the low-q asymptotics visible for q > 1 correlates with the deviation of the total static structure factor S(q) from its low-q plateau (cf.Fig. 6).

Fig. 3 Rescaled correlation functions in NRC and reciprocal space. The bold dashed and solid lines indicate the expected low-q limit eqn (2). Inset: and vs. θ for q ≈ 0.1. Main panel: Semi-logarithmic representation of ICFs βVc_T(q) and βVc_L(q) vs. q.

B. Correlation functions in reciprocal space

While remaining in reciprocal space we consider next coordinate frames which are either unrotated (α = 0) or rotated as in Fig. 1(b) using the same angle α for all q. According to eqn (5) the correlation functions c_αβγδ(q) of isotropic achiral systems in two dimensions depend quite generally on the four ICFs c_L(q), c_T(q), c_N(q) and c_⊥(q). Due to eqn (18) the last two of these ICFs must vanish while c_L(q) and c_T(q) are given by eqn (2). Let us introduce for later convenience the two “creep compliances”This yields in the original coordinatesfor q → 0. The dots mark irrelevant constant contributions.§ See Appendix D for more details. For correlation functions in rotated coordinate systems one merely needs to substitute θ by x = θ − α. These relations are put to a test in Fig. 4 where we focus for clarity on the reduced shear-strain autocorrelation function for n = 40 000. The angular dependences are presented in the main panel for one wavevector in the low-q limit. Focusing on the first term in eqn (20) we have taken off the mean constant average over all θ (corresponding to the dots). Importantly, all data for different α are seen to collapse when plotted as a function of the scaling variable x. Obviously, this simple scaling (without characteristic angle) would not hold for anisotropic systems. To obtain a precise test of the q-dependence of c_αβγδ(q) we project out the angular dependences usingfor p = 2 and p = 4. For convenience the prefactor of the integral is chosen such that P[cos(2θ),2] = P[cos(4θ),4] = 1. The result for the shear-stress autocorrelation function with p = 4 is shown in the inset of Fig. 4. In agreement with eqn (20) the presented data is given by J₁/8 (solid line) for sufficiently small wavevectors. Equivalent results have been obtained for the other correlation functions mentioned above.

Fig. 4 Rescaled correlation function for n = 40 000. Main panel: Angle dependence of vertically shifted f(q) for q ≈ 0.1. Data collapse is observed using the reduced angle x = θ−α. The bold solid line indicates the prediction, eqn (20). Inset: Comparison of P[f,p](q) for p = 4 with the predicted low-q limit J₁/8 (bold solid line).

C. Correlation functions in real space

We turn finally to the correlation functions in real space. As shown in Appendix D, inverse FT implieswith x = θ − α being the difference of the angles θ and α indicated in Fig. 1. The same large-r limit holds also for and for . Moreover,for r ≫ 1, i.e. a dipolar symmetry is expected. A verification of the r-dependence is obtained using again (now in real space) the projection P[f,p](r), cf.eqn (24). Focusing on n = 40 000 several rescaled correlation functions f(r) are presented in Fig. 5. In agreement with eqn (25) the indicated first three cases collapse for p = 4 and r ≳ 20 on J₁/4πr² (bold solid line). This confirms the octupolar symmetry of these correlation functions. Confirming eqn (26) the last case with f(r) = −β(c₁₁₁₁(r) − c₂₂₂₂(r))/2 collapses onto J₂/4πr² (dashed line). p = 2 is used here in agreement with the predicted quadrupolar symmetry of this correlation function. Similar results are obtained for other particle numbers n.

Fig. 5 P[f,p](r) for various correlation functions and modes p for n = 40 000. The bold solid line marks the prediction J₁/4πr² for the first three cases, the dashed line the prediction J₂/4πr² for the last one.

V. Linear response to point stress

A. Time-dependent strain correlations

Correlation functions describe quite generally the linear response to a small imposed perturbation.^7,34,35,42 Being tensorial fields, just like the correlation function fields, the response fields must in general depend on the direction of the field vector and on the orientation of the coordinate system. As already emphasized in Section II4 and Section II5, the response fields contains information of both the system and the imposed source and the source term may either be isotropic or anisotropic. We elaborate here this general point focusing, naturally, on correlation functions of the instantaneous strain field [small epsilon, Greek, circumflex]_αβ(r) = ε_αβ(r, t). Extending very briefly our discussion to the time domain, let us introduce the time-dependent correlation functions

c_αβγδ(q, t) = 〈ε_αβ(q, t)ε_γδ(−q, t = 0)〉 (27)

of the strain fields in reciprocal space with t being the “time lag”.³³ Naturally, this reduces to eqn (1) for t → 0. This definition allows us to take advantage of the general “Fluctuation-dissipation theorem” (FDT) of statistical mechanics as stated, e.g., in ref. 34, 35 and 42. We thus anticipate immediate generalizations of the present study for time-dependent tensorial correlation and response fields which will be discussed elsewhere.

B. Fluctuation-dissipation theorem

Let us switch on at time t = 0 a small perturbationto the Hamiltonian of the system with δσ_αβ(r) being an imposed external stress field. This is equivalent to the application of an appropriate external perturbative force field to each particle. For a general “growth function”⁴² in response to a sudden application of a step field, such as eqn (28), the relevant FDT relations are stated (for scalar fields) by eqn (3.65) and (3.67) of ref. 42. The mean strain increment δε_αβ(r, t) induced by this perturbation is then given in real space by a convolution integral for the time-dependent correlation functions c_αβγδ(r, t) and the stress perturbation δσ_αβ(r). Using eqn (B6) this relation may be written more compactly in reciprocal space as

δε_αβ(q, t) = βV [c_αβγδ(q, t = 0) − c_αβγδ(q, t)]δσ_γδ(q) (29)

where the summation over repeated indices is essential and cannot be omitted. Note that δε_αβ(q, t) = 0 for t ≤ 0 and that, since all c_αβγδ(q, t) are continuous functions of time, the creep response δε_αβ(q, t) must also be continuous, especially at t = 0.⁴⁰ The time-dependent creep is thus determined by the time-dependent correlation functions and the imposed stress perturbation. We are interested here only in the static equilibrium response a long time after the perturbation is switched on. The time-dependent strain correlation functions (computed for the unperturbed Hamiltonian at switched off external perturbation ) must, of course, vanish

c_αβγδ(q, t) → 0 for t → ∞. (30)

Hence, eqn (29) reduces to

δε_αβ(q) = βVc_αβγδ(q)δσ_γδ(q) for t → ∞ (31)

with denoting the long-time creep and standing for the spatial correlation function without time lag, eqn (1), as everywhere else in this paper.

C. Response to point source

Following the discussion at the end of Section II, we investigate now the long-time creep for a point source

δσ_αβ(r) = s_αβδ(r) (32)

with s_αβ being a symmetric 2 × 2 matrix (of dimension “stress × volume = energy”). According to the FDT relation eqn (31) this implies

δε_αβ(q) = βc_αβγδ(q)s_γδ for t → ∞ (33)

and an equivalent relation in real space. As in Section IID it is convenient to diagonalize s_αβ by an appropriate rotation of the coordinate system. The perturbation becomes therefore equivalent to that of two small force dipoles²³ oriented along the eigenvectors. The real-space analogue of eqn (10) is thus given by

δε_αβ(r) = βc_αβ11(r)s₁₁ + βc_αβ22(r)s₂₂. (34)

Using in addition eqn (D13) we may replace for isotropic systems the real space correlation functions by the corresponding invariants ĩ_n(r) in real space. Introducing the scalars s₁ = s_γγ/2 and s₂ = r̂_αs_αβr̂_β/2 this may be written quite generallyTaking now advantage of the specific results for strain correlations presented in Section IVC and in Appendix D the invariants ĩ_n(r) are given by eqn (D14), i.e. we may quite generally express δε_αβ(r) in terms of the two creep compliances J₁ and J₂, cf.eqn (19).

We also note that the term s_αβ in the last line of eqn (35) must be isotropic, i.e. s₁ = s₁₁ = s₂₂ = 2s₂, to obtain an isotropic second-order tensor field in agreement with eqn (A20). As expected from the more general argument given in Section 2, the shear strain increment

becomes quadrupolar in this case.

As already emphasized in Section IIE, the source tensor need not necessarily be isotropic albeit the system is isotropic. To be specific, let us consider the “shear transformation zone” model for localized plastic failure involving two orthogonal force dipoles of opposite signs.²³ Hence, s₁₁ = −s₂₂ and s₁ = 0 and s₂ = s₁₁(r̂₁² − r̂₂²)/2. Using eqn (13) or eqn (35) this yields

for the shear strain response.

As one expects on general grounds, all reference to statistical physics, i.e. the inverse temperature β, drops out in both cases. Moreover, the shear strain response naturally strongly depends on the type of source term applied at the origin: for force dipoles of same sign it is quadrupolar and proportional to J₂ while for dipoles of opposite sign it gets octupolar and proportional to J₁.

VI. Conclusion

A. Summary

Strain correlation functions. The present work has focused on correlation functions of components of strain tensor fields in two-dimensional, isotropic and achiral elastic bodies. This was done theoretically using

• the general mathematical structure of isotropic tensor fields as summarized in Section IIB, cf.eqn (5), and in more detail in Appendix A and

• the well-known equipartition theorem of statistical physics for macroscopic strain fluctuations in reciprocal space, cf.eqn (2).

Numerically we have tested our predictions by means of glass-forming particles deep in the glass regime. This shows that these correlation functions may depend on the coordinates of the field variable (q_α in reciprocal space or r_α in real space) and implies in turn that they depend in general on the direction of the field vector and on the orientation of the coordinate system. Scaling with x = θ − α these angular dependencies are distinct from those of ordinary anisotropic systems. Importantly, correlation functions of strain tensor fields are components of an isotropic forth-order tensor field, eqn (5), being characterized by the two ICFs c_L(q) and c_T(q). With the asymptotic plateau values being given by two Lamé coefficients, eqn (2), all strain correlation functions are determined and all (finite) real-space strain correlations must be long-ranged decaying as 1/r² (cf.Fig. 5). We thus obtain similar results as in our recent study on correlation functions of stress tensor fields.⁴ Note that time-averaged stress fields have been probed in the latter study while correlations of instantaneous strain fields have been considered here. Our numerical findings do agree with other studies of strain correlations^29,30,32 being, however, now traced back to the isotropy of the system and the tensor field nature of the probed correlations. Importantly, we have given here a complete and asymptotically exact description for the correlation functions of strain tensor fields of isotropic elastic bodies. No additional physical assumption is thus needed (for sufficiently small wavevectors).

Response to tensorial point sources. We also discussed the associated linear response fields as defined in general mathematical terms by the tensorial contraction of the correlation function tensor by means of a source tensor and, more physically, by the FDT relation for the strain increment due to an imposed small stress perturbation, cf.eqn (28) and (29). Naturally, the response must by definition contain information from both the correlation functions, characterizing the system, and from the imposed source tensor which may not be isotropic. We have emphasized that the summation over repeated indices must be properly performed, i.e. the response field is not given by one correlation function times a scalar but by the sum over all eigenvalues of the source tensor. For this reason response and correlation fields, albeit closely related, have in general different angular dependences, e.g., the shear strain correlation function c₁₂₁₂(r) in an isotropic system must be octupolar, cf.eqn (25), while the shear strain response δε₁₂(r) may be either quadrupolar for an isotropic source, cf.eqn (36), or octupolar for an anisotropic source corresponding to two force dipoles of opposite signs, cf.eqn (37). Albeit all contributing correlation functions are isotropic the response field is anisotropic in the latter case due to the source. It is thus important to not lump together correlation functions and response fields. Mesoscopic elasto-plastic models^45,46 thus must specify not only the correlation functions but also the source tensors characterizing the local plastic events.

B. Outlook

Our work suggests several natural extensions:

• The general mathematical framework for isotropic tensor fields and the discussed relations and numerical procedure for correlation and response fields naturally generalize to higher spatial dimensions, especially for the three-dimensional case.

• The present work has focused on Euclidean spaces and Cartesian coordinates. A generalization for systems embedded in non-Euclidean spaces, say for glasses on spheres,^51,52 and more general curvilinear coordinate systems^1,5 may be worked out.

• The present work has focused on the large-wavelength limit (q → 0). More generally, one may express the longitudinal and transverse ICFs c_L(q) and c_T(q) for finite q as

in terms of the generalized longitudinal and transverse elastic moduli L(q) and G(q) (with L(q) → λ + 2μ and G(q) → μ for small q).¶ ^17,48,50 It can be shown⁵⁰ that both the isotropicity and the harmonicity of the strain modes assumed in the derivation of eqn (38) are well justified for the present model up to q ≈ 1 while deviations become relevant for larger wavevectors, especially around the main peak of the static structure factor S(q).

• A further generalization of the current work concerns time-dependent correlation functions c_αβγδ(q, t) as defined in eqn (27). These can be again expressed viaeqn (5) in terms of (now time-dependent) longitudinal and transverse ICFs c_L(q, t) and c_T(q, t). These time-dependent ICFs are given in turn by time-dependent creep compliance material functions which can be related to the two time-dependent material functions L(q, t) and G(q, t).¹⁷ Strain correlation functions may thus reveal octupolar pattern whenever the invariant

|i₄(q, t)| = |c_L(q, t) − 4c_T(q, t)| (39)

is sufficiently large. Since i₄(q, t) must become q-independent for small q, a long-range decay with

c_αβγδ(r, t)⋍ 1/r^d (40)

is generally expected for the time-dependent correlation functions in isotropic d-dimensional systems. Using eqn (29) similar long-range relations are predicted for the associated dynamical response fields.

• It may be also of interest to characterize correlations of tensor fields of different order. For instance, the forth-order elastic modulus field E_αβγδ(r)^26,53 may be characterized by a correlation function tensor of order eight.⁵³ Strong angular dependencies are expected based on our formalism. For isotropic systems these correlation functions must again adopt a general mathematical structure in terms of a small finite number of ICFs. Once these ICFs are characterized (theoretically or numerically using NRC) in the low-q limit all correlation functions are again determined.

Author contributions

J. P. W. designed and wrote the project benefitting from contributions of all co-authors.

Conflicts of interest

There are no conflicts to declare.

Appendix A Summary of isotropic tensor fields

1 Background

Isotropic systems are described by “isotropic tensors” and “isotropic tensor fields”. We give here a brief recap of various useful aspects already presented elsewhere.^3,4 Quite generally, a tensor field assigns a tensor to each point of the mathematical space, in our case a d-dimensional Euclidean vector space.³ We denote an element of this vector space by the “spatial position” r (real space) or by the “wavevector” q for the corresponding reciprocal space. The relations for tensor fields are formulated in reciprocal space since this is more convenient both on theoretical and numerical grounds due to the assumed spatial homogeneity (“translational invariance”). The corresponding real space tensor field is finally obtained by inverse FT.

For simplicity we assume Cartesian coordinates with an orthonormal basis {e₁,…, e_d}.^1,3,9 Greek letters α, β, … are used for the indices of the tensor (field) components. A twice repeated index α is summed over the values 1, …, d, e.g., q = q_αe_α with q_α standing for the vector coordinates. This work is chiefly concerned with tensors T^(o) = T_α1…αoe_α1… e_αo of “order” o = 2 and o = 4 and their corresponding tensor fields with components depending either on r or q. The order of a component is given by the number of indices. Note that

for the d^o coordinates in real and reciprocal space (with denoting the FT as discussed in Appendix B).

We consider linear orthogonal coordinate transformations (marked by “*”) with matrix coefficients c_αβ given by the direction cosine .³ We remind that³

under a general orthogonal transform. For a reflection of the 1-axis we thus have, e.g.,i.e. we have sign inversion for an odd number of indices equal to the index of the inverted axis. The field vector remains unchanged by these “passive” transforms albeit its coordinates change.

2 Definitions, properties and construction of general isotropic tensors and tensor fields

Isotropic tensors. Components of an isotropic tensor remain unchanged by any orthogonal coordinate transformation,^3,9i.e.As noted at the end of Section A1 the sign of tensor components change for a reflection at one axis if the number of indices equal to the inverted axis is odd. Consistency with eqn (A4) implies that all tensor components with an odd number of equal indices must vanish, e.g.,

T₁₂ = T₁₁₁₂ = T₁₂₂₂ = 0. (A5)

Isotropic tensor fields. The corresponding isotropy condition for tensor fields is given by³with which reduces to eqn (A4) for q = 0. Please note that the fields on the left handside of eqn (A6) are evaluated with the original coordinates of the vector field variable q while the fields on the right handside are evaluated with the transformed coordinates. Another way to state this is to say that the left hand fields are computed at the original vector q = (q₁,…,q_d) while the right hand fields are computed at the “actively transformed” vector . It is for this reason that finite components with an odd number of equal indices, e.g., T₁₂₂₂(q) ≠ 0, are possible in principle for finite q.

Natural rotated coordinates. Fortunately, there is a convenient coordinate system where the nice symmetry eqn (A5) for isotropic tensors can be also used for isotropic tensor fields. In these “Natural Rotated Coordinates” (NRC) the coordinate system for each wavevector q is rotated until the 1-axis coincides with the q-direction, i.e. with q = |q|. These tensor field components in NRC are marked by “°” to distinguish them from standard rotated tensor fields (marked by primes “′”) where the same rotation is used for all q. If in addition is an even function of its field variable q (as in the case of achiral systems for even order o) it can be shown⁴ that all tensor field components with an odd number of equal indices must vanish.

Product theorem for isotropic tensor fields. Let us consider a tensor field C(q) = A(q) ⊗ B(q) with A(q) and B(q) being two isotropic tensor fields and ⊗ standing either for an outer product, e.g. c_αβγδ(q) = A_αβ(q)B_γδ (q), or an inner product, e.g. c_αβγδ(q) = A_αβγν(q)B_νδ(q). Hence,using in the second step a general property of tensor (field) products, due to eqn (A2), and in the third step eqn (A6) for the fields A(q) and B(q) where q* stands for the “actively” transformed field variable. C(q) is thus also an isotropic tensor field. One may use this theorem to construct isotropic tensor fields from known isotropic tensor fields A(q) and B(q).

3. Summary of assumed symmetries

All second-order tensors in this work are symmetric, T_αβ = T_βα, and the same applies for the corresponding tensor fields in either r- or q-space. This is, e.g., the case for the strain field , cf.eqn (14), or the source tensor s_αβ needed for a response field, cf.eqn (8). We assume for all forth-order tensor fields that

T_αβγδ(q) = T_βαγδ(q) = T_αβδγ(q) (A8)

T_αβγδ(q) = T_γδαβ(q) (A9)

and

T_αβγδ(q) = T_αβγδ(−q). (A10)

Note that eqn (A10) is necessarily valid both for achiral and chiral two-dimensional isotropic systems. Forth-order tensor fields are often constructed by taking outer products⁹ of second-order tensor fields. We consider, e.g., correlation functions 〈T̂_αβ(q)T̂_γδ(−q)〉 with T̂_αβ(q) being an instantaneous second-order tensor field. eqn (A8) then follows from the symmetry of the second-order tensor fields. The evenness of forth-order tensor fields, eqn (A10), is a necessary condition for achiral systems. It implies that T_αβγδ(q) is real if T_αβγδ(r) is real and, moreover, eqn (A9) for correlation functions since

〈T̂_αβ(q)T̂_γδ(−q)〉 = 〈T̂_γδ(q)T̂_αβ(−q)〉. (A11)

As already emphasized, all our systems are assumed to be isotropic, i.e., eqn (A6) must hold for ensemble-averaged tensor fields.

4. General mathematical structure

General structure of tensors. Isotropic tensors of different order are discussed, e.g., in Section 2.5.6 of ref. 9. Due to eqn (A5) all such tensors of odd order must vanish. The finite isotropic tensors of lowest order are thus

T_αβ = k₁δ_αβ, (A12)

T_αβγδ = i₁δ_αβδ_γδ + i₂(δ_αγδ_βδ + δ_αδδ_βγ) (A13)

with k₁, i₁ and i₂ being invariant scalars. Please note that all symmetries stated above hold, especially also eqn (A5). Note that the symmetry eqn (A8) was used for the second relation, eqn (A13). Importantly, this implies that only two coefficients are needed for a forth-order isotropic tensor. As a consequence, the elastic modulus tensor is completely described by two elastic moduli (cf. Section 9.3).

General structure of tensor fields. We restate now the most general isotropic tensor fields for 1 ≤ o ≤ 4 and focusing on two-dimensional systems (d = 2) compatible with the symmetries stated in Section 7.3. With l_n(q), k_n(q), j_n(q) and i_n(q) being invariant scalar functions of q = |q| we have^3,4

T_α(q) = l₁(q)q̂_α (A14)

T_αβ(q) = k₁(q)δ_αβ + k₂(q)q̂_αq̂_β (A15)

T_αβγ(q) = j₁(q)q̂_αδ_βγ + j₂(q)q̂_βδ_αγ + j₃(q)q̂_γδ_αβ + j₄(q)q̂_αq̂_βq̂_γ (A16)

for finite wavevectors q. See ref. 4 for a derivation, generalizations for d > 2 and a discussion of the limit q → 0. Terms due to the invariants k₁(q), i₁(q) and i₂(q) are independent of the coordinate system. All other terms depend on the components q̂_α of the normalized wavevector q̂ and thus on the orientations of the field vector and of the coordinate system.

Alternative representation for forth-order tensor fields.

It is convenient to define the four functionsusing NRC. For an isotropic system these four functions can only depend on the wavenumber q but not on the direction q̂ of the wavevector q. Importantly, all other components are either by eqn (A8) and (A9) identical to these invariants or must vanish for an odd number of equal indices as reminded in Section A2. The d⁴ = 16 components are thus completely determined by the four invariants and this for any q. T_αβγδ(q) is then obtained by the inverse rotation to the original unrotated frame using eqn (A2). It is readily seen thatbeing consistent with eqn (6).

Isotropic tensor fields in real space. We have formulated above all tensor fields in terms of the wavevector q and its components since it is most convenient to start the analysis in reciprocal space. The above results also hold, however, in real space. This implies, e.g., for isotropic (and achiral) fields in two dimensions that

T_αβ(r) = k̃₁(r)δ_αβ + k̃₂(r)r̂_αr̂_β (A20)

for r > 0 with k̃_n(r) and ĩ_n(r) denoting the invariants in real space and r̂_α = r_α/r components of the normalized vector r̂ = r/r. As already stated, eqn (A1), the tensor field components in real and reciprocal space are related by FT. Note that k̃_n(r) and ĩ_n(r) are in general not the FTs of, respectively, k_n(q) and i_n(q). For the important case that the invariants in reciprocal space are q-independent constants it follows quite generally thatfor r > 0. (Additional δ(r)-terms arise at the origin. The constant invariants k₁, i₁ and i₂ only contribute to these terms.) That this holds can be readily shown using relations put forward in Appendix B and Appendix D.

Appendix B: Useful Fourier transformations

1. Continuous Fourier transform

We consider real-valued functions f(r) in d dimensions. As in ref. 4, 18 and 54 the Fourier transform (FT) from “real space” (variable r) to “reciprocal space” (variable q) is defined bywith V being the volume of the system. The inverse FT is then given byNote that f(r) and f(q) have the same dimension. For notational simplicity the function names remain unchanged. We remind the FTswith δ(r) being Dirac's delta function. Let us consider the spatial convolution functionin real space. With and this implies according to the “convolution theorem”⁵⁵We also remind for completeness that the spatial correlation functionof real-valued fields g(r) and h(r) becomes according to the “correlation theorem”⁵⁵

c(q) = g(q)h*(q) = g(q)h(−q) (B8)

with marking the conjugate complex. For auto-correlation functions, i.e. for g(r) = h(r), this simplifies to (“Wiener–Khinchin theorem”)

c(q) = g(q)g(q) = |g(q)|², (B9)

i.e. the Fourier transformed auto-correlation functions are real and ≥0 for all q. Moreover, we shall consider correlation functions c(r), eqn (B7), being even in real space, c(r) = c(−r), and thus also in reciprocal space, c(q) = c(−q) = c*(q), i.e. c(q) is real.

2. Discrete Fourier transform on microcell grid

All fields f(r) are stored on a regular equidistant d-dimensional grid as shown in Fig. 2 for d = 2. Periodic boundary conditions are assumed.³³ The discrete FT and its inverse becomewith and being discrete sums over n_V = n^d_L = V/a^d_grid grid points in, respectively, real or reciprocal space. As shown in Fig. 2 we label the grid points in real and reciprocal space usingTo take advantage of the Fast-Fourier transform (FFT) routines⁵⁵ the number of grid points in each spatial direction n_L = L/a_grid is an integer-power of 2.

3. Fourier transform of planar harmonic functions

As discussed in the main part, all correlation functions in reciprocal space become in the large-wavelength limit independent of the magnitude q of the wavevector q but depend on its angle θ_q (and, more generally, on the angle difference θ_q −α for rotated coordinate frames). As noted in Section IIC, these angular dependencies can be uniquely expressed in terms of the planar harmonic basis functions cos(pθ_q) and sin(pθ_q) with p being an integer. We denote these orthogonal basis functions by b_p(θ_q). We thus need to compute the inverse FTs of f(q) = b_p(θ_q)/V. More specifically, we are interested in modes with p = 2 and p = 4. Additional constant terms (p = 0), such as the ones indicated by dots in eqn (20)–(22), are irrelevant leading merely to δ(r)-contributions at the origin. For d = 2 eqn (B2) becomeswith θ_r being the angle of r̂ = (cos(θ_r), sin(θ_r)). We make now the substitution θ = θ_q − θ_r and use that⁵⁶

cos(pθ + pθ_r) + cos(−pθ + pθ_r) = 2cos(pθ)cos(pθ_r)

sin(pθ + pθ_r) + sin(−pθ + pθ_r) = 2cos(pθ)sin(pθ_r).

We remind that following eqn (9.1.21) of ref. 56 the integer Bessel function J_p(z) may be writtenwhich leads toFor finite r we may rewrite this aswhere we have set . As may be seen from eqn (11.4.16) of ref. 56 the latter integral becomes||

I_p(t) → 1 for t → ∞ and p > −2 (B17)

from which we obtain the final resultNote that f(r) = f(−r) and that f(r) is real for even p. Generalizing the above argument it is seen that f(r) ∝ 1/r^d for higher dimensions d.

Appendix C Computational details

1 Simulation model

We consider systems of polydisperse Lennard-Jones (pLJ) particles in d = 2 dimensions where two particles i and j of diameter D_i and D_j interact by means of a central pair potential^{4,18,36–39,48}being the reduced distance according to the Lorentz rule.³⁴ This potential is truncated and shifted^4,33 with a cutoff s_cut = 2s_min given by the minimum s_min of u(s). Lennard-Jones units³³ are used throughout this study, i.e. ε = 1 and the average particle diameter is set to unity. The diameters are uniformly distributed between 0.8 and 1.2. We also set Boltzmann's constant k_B = 1 and assume that all particles have the same mass m = 1. The last point is irrelevant for the presented Monte Carlo (MC) simulations.³³ Time is measured in units of MC steps (MCS) throughout this work.

2. Parameters and configuration ensembles

We focus on systems with n = 10000 and n = 40 000 particles. We first equilibrate = 200 independent configurations c at a high temperature T = 0.55 in the liquid limit. These configurations are adiabatically cooled down using a combination of local MC moves³³ and swap MC moves exchanging the sizes of pairs of particles.^38,49 In addition, an MC barostat³³ imposes an average normal stress P = 2.^36,38 At the working temperature T = 0.2 we first thoroughly temper over Δτ = 10⁷ all configurations with switched-on local, swap and barostat MC moves and then again over Δτ = 10⁷ with switched-on local and swap moves and switched-off barostat moves. The final production runs are carried out at constant volume V only keeping local MC moves. Under these conditions, T = 0.2 is well below the glass transition temperature T_g ≈ 0.26 determined in previous work.^36,38 Due to the barostat used for the cooling the box volume V = L^d differs slightly between different configurations c while V is identical for all frames t of the time-series of the same configuration c. In all cases the number density is of order unity. For each particle number n and each of the independent configurations c we store ensembles of time series containing = 10 000 instantaneous “frames” t. These are obtained using the equidistant time intervals δτ = 1000 for n = 10 000 and δτ = 100 for the other system sizes.

3. Macroscopic linear elastic properties

The amorphous glasses formed by pLJ particle systems at a pressure P = 2 and a temperature T = 0.2 ≪ T_g ≈ 0.26 are for sampling (production) times Δτ ≤ 10⁷ MCS reversible linear elastic bodies whose plastic rearrangements can be neglected for all practical purposes.^{18,21,36,38,57,58} Moreover, these systems can be shown to be isotropic above distances corresponding to a couple of particle diameters.^18,38 Following eqn (A13) the forth-order elastic modulus tensor for isotropic systems may be written^8,9

E_αβγδ = λδ_αβδ_γδ + μ(δ_αγδ_βδ + δ_αδδ_βγ) (C2)

in terms of the two isothermic Lamé moduli λ and μ. As described in detail elsewhere^{21,36,57–60} we have determined λ and μ either by means of strain fluctuations, e.g., by letting the box volume V fluctuate at imposed pressure P,⁵⁷ or using the stress-fluctuation formalism at fixed volume and shape of the simulation box.^59,60 This shows that λ ≈ 38 and μ ≈ 14. We have verified especially that similar values are obtained for n ≥ 5000 and using different components, say E₁₁₁₁ and E₂₂₂₂ for λ + 2μ, and that the fluctuations of E_αβγδ|_c for independent configurations c become negligible for n ≥ 5000.

4. Discrete fields on square grid

We turn now to the relevant microscopic tensor fields as functions of either the spatial position r (real space) or the wavevector q (reciprocal space). The different fields are stored on equidistant discrete grids as sketched in Fig. 2. The same n_L is used for both spatial directions and for all configurations and frames of a given particle number n. As already mentioned, the box volume V = L² fluctuates slightly between different configurations c (at same n) due to the barostat used for the cooling, tempering and equilibration of the systems. Accordingly, a_grid also differs between different configurations c. These fluctuations become small, however, with increasing system size. If nothing else is mentioned we report data obtained using a lattice constant a_grid ≈ 0.2. As shown in Fig. 6 for the rescaled transverse ICF βVc_T(q) plotted using a double-logarithmic representation, there is no need to further decrease a_grid. Even the very large grid constant a_grid ≈ 3.2 gives, apparently, the correct large-wavelength asymptote βc_T(q) ≈ 1/4μ indicated by the dashed horizontal line.

Fig. 6 Rescaled transverse ICF βVc_T(q) for different grid constants a_grid as indicated. The open symbols have been obtained using eqn (C4), the thin solid line using eqn (C5) and n_L = 2048. Importantly, we obtain the same results in all cases where qu_rms ≪ 1 and qa_grid ≪ 1. Even a rather coarse grid, say for n_L = 64, is sufficient to confirm the expected large-wavelength limit (horizontal dashed line). The total static structure factor S(q) is shown for comparison (solid line). The “dip” of S(q) around q ≈ 4 is caused by the polydispersity of the particles as emphasized elsewhere.⁴⁸S(q) and βVc_T(q), at least for sufficiently small a_grid, have both a strong peak located similarly at q ≈ 6.5 (arrow).

5. Displacement fields

As in previous experimental and numerical studies^10,11,21 the displacement field u(r) is constructed from the instantaneous spatial positions r_a of the particles a with respect to their reference positions r_a. As reference position r_a we have used either the average particle position determined using a long trajectory or the particle position after a rapid quench to T = 0. Having not observed any significant quantitative difference between both methods we only report here data computed using the first one. We thus get first the displacement vector u_a = r_a − r_a for each a. By construction the average displacement vector 〈u_a〉 must vanish. We find

u_rms ≡ 〈u_a²〉 > ^1/2 ≈ 0.13 (C3)

for the root-mean-squared average u_rms (sampled over all particles, frames and configurations). The instantaneous displacement field may then be defined by^10,11,21assuming for the moment an infinitessimal fine grid, i.e. a_grid → 0. (Different prefactors have been used in ref. 10, 11 and 21). We remind that the δ(r)-function has the dimension “1/volume”. By definition u(r) has thus the same dimension “length” as the displacement vector u_a. Following the common definition of the particle flux density,³⁴ the reference position r̃_a in the δ-function may be replaced by the time-dependent position r_a, i.e. the displacement field may alternatively be defined byBoth operational definitions are compared for the transverse ICF βVc_T(q) in Fig. 6 where data obtained using eqn (C4) are indicated by open symbols. In reciprocal space we obtainfor eqn (C4) and similarly for eqn (C5) with r_a replacing r̃_a. Since r_a = r̃_a+ u_a we have to leading order

exp(−iq·r_a) ≈ exp(−iq·r_a)(1 − iq·u_a…) (C7)

for q|u(q)| ≪ 1. Both operational definitions eqn (C4) and (C5) thus become equivalent for qu_rms ≪ 1. Due to the small typical (root-mean-squared) displacement u_rms, eqn (C3), this holds for all sampled q as may be seen from the data presented in Fig. 6. Due to the center-of-mass convention for all particle displacements the volume integral over u(r) must vanish and, equivalently, we have u(q = 0) = 0 in reciprocal space for each instantaneous field. We also remind that the two coordinates of the displacement field in NRC are the longitudinal component u_L(q) ≡ u₁°(q) and the transverse component u_T(q) ≡ u₂°(q).

In numerical practice, the continuous field vector r of eqn (C4) and (C5) corresponds to the discrete point on the grid, eqn (B12), closest (using the minimal image convention) to, respectively, the reference position r_a or the particle position r_a. Strictly speaking, we thus obtain by means of eqn (B10) the FT with respect to their respective closest grid points. (In principle one could directly without approximation compute the displacement field u(q) using eqn (C6) in reciprocal space. Unfortunately, this leads to an additional loop over all n particles for each wavevector q.) The differences between these definitions become neglible for qa_grid ≪ 1. That this holds can be clearly seen from Fig. 6 where we have varied a_grid over more than one order of magnitude.

6. Linear strain fields

Using the displacement field u(r) the linear (“small”) strain tensor field is defined by^8,9Due to eqn (B3) this becomesin reciprocal space as already stated in Section IIIB. (Obviously, ε_αβ = ε_βα for any r in real space and any q in reciprocal space.) Note that both ε_αβ(r) and its FT ε_αβ(q), cf.eqn (B1), are dimensionless fields. Due to our definitions and conventions the macroscopic strain ε_αβ(q = 0) is assumed to vanish. Using eqn (B12) we numerically determine the relevant components of the (symmetric) strain tensor field ε_αβ(q) from the two components of the displacement field u_α(q) stored on the reciprocal space grid (cf.Fig. 2) and the wavevector q_α according to eqn (B12). As noted in Section IIIB, in NRC there are only two non-vanishing strain fields, namely the longitudinal and transverse strain fields ε_L(q) and ε_T(q) linearly related to the corresponding displacement fields u_L(q) and u_T(q), cf.eqn (16) and (17).

7 Correlation function fields

Using the strain tensor fields ε_αβ(q) in reciprocal space computed according to eqn (C9) for each c and t we obtain the correlation functions c_αβγδ(q) = 〈ε_αβ(q)ε_γδ(−q)〉 averaged over all c and t. For the reported correlation functions in a coordinate system turned by an angle α we first compute the new components and of displacement field and wavevector. (Alternatively, one may also rotate ε_γδ(q).) For the ICFs c_L(q) and c_T(q) obtained using NRC we first get the longitudinal and transverse displacement fields u_L(q) and u_T(q) in NRC and from those using eqn (16) and (17) the longitudinal and transverse strains ε_L(q) and ε_T(q). c_L(q) = 〈ε_L(q)ε_L(−q)〉 and c_T(q) = 〈ε_T(q)ε_T(−q)〉 are computed by averaging over all c and t and all wavevectors q with magnitude |q| within a chosen bin around q. The correlation functions c_αβγδ(r) in real space (either for unrotated or α-rotated coordinate systems) are finally obtained by inverse FFT.

Appendix D: From c_L(q) and c_T(q) to c_αβγδ(r)

As shown in Appendix A4, a forth-order tensor field describing an isotropic achiral system in two dimensions is given by eqn (A17) in terms of four invariants i_n(q). In turn these invariants are expressed in terms of the alternative set of invariants c_L(q), c_T(q), c_N(q) and c_T(q). Due to eqn (18) we have c_N(q) = c_⊥(q) = 0 for strain correlations. The correlation function c_αβγδ(q) in reciprocal space are thus given by the invariantsand

i₄(q) =c_L(q) − 4c_T(q). (D2)

More specifically, this implieswith c = cos(θ) = q̂₁ and s = sin(θ) = q̂₂ being coefficients depending only on the wavevector angle θ. Alternatively, the six relations eqn (D3) may also be obtained using that the components ε_αβ(q) in the original coordinate frame can be expressed asin terms of the longitudinal and transverse strains ε_L(q) and ε_T(q) and the fact that ε_L(q) and ε_T(q) fluctuate independently. For the correlation functions in rotated coordinate systems one simply replaces θ by x = θ − α. Note also that and do in general not vanish for all x in standard (unrotated or rotated) coordinates. The values for NRC are obtained by setting x = 0, i.e. s = 0 and c = 1. We expand the angle-dependent coefficients before c_L(q) and c_T(q) in terms of the planar harmonic functions cos(pθ) and sin(pθ) with p being integers using standard trigonometric relations.⁵⁶ This implies for example
where we have omitted the arguments q on the l.h.s. and q on the r.h.s. The prefactor of cos(4θ) and sin(4θ) is always proportional to i₄(q). Having in mind the equipartition relation eqn (2) and using the creep compliances J₁ and J₂ introduced in eqn (19) one sees thatin the low- q limit. We thus get in reciprocal space
where the dots mark constant terms. These terms are irrelevant for the inverse FT, only leading to contributions at r = 0. We may thus take advantage of the analytical result for the inverse FT eqn (B18). This leads towith θ denoting the polar angle of the field vector r. The correlation functions in rotated coordinate systems generalize the above equations by substituting θ with the angle difference x = θ − α. These results can be rewritten compactly using the general form expected from eqn (A21) for a manifest two-dimensional, isotropic and achiral forth-order tensor field in real space yieldingwhere the invariants ĩ_n(r) in real space are given bySince i₁ = −(2J₁ + J₂)/4, i₂ = (2J₁ + J₂)/8, i₃ = (2J₁ + J₂)/4 and i₄ = −J₁, this is consistent with the more general relation eqn (A22).

The angle dependence for the shear-strain autocorrelation function in real space is investigated in Fig. 7 where we plot using linear coordinates as a function of x for different r and α. (To obtain sufficiently high statistics we need to average over the indicated finite r-bins. This is done by weighting each data entry for a bin with the proper factor 4πr².) The data compare well with the prediction, eqn (D9), confirming thus especially the scaling with angle difference x = θ − α. Naturally, the statistics deteriorates with increasing r due to the faster decay of the correlations as compared to the noise.

Fig. 7 Rescaled shear-strain autocorrelation function in real space as a function of x = θ − α comparing data for different r-intervals and rotation angles α with the prediction (bold solid line).

Acknowledgements

We are indebted to the University of Strasbourg for computational resources.

References

1. A. J. McConnell, Applications of Tensor Analysis, Hassell Street Press, 2021.
2. J. A. Schouten, Tensor Analysis for Physicists, Dover Publications, Oxford, 2015.
3. W. Schultz-Piszachich, Mathematik 11: Tensoralgebra und -analysis. Thun und Frankfurt/Main, Verlag Harri Deutsch, 1977.
4. J. P. Wittmer, A. N. Semenov and J. Baschnagel, Correlations of tensor field components in isotropic systems with an application to stress correlations in elastic bodies, Phys. Rev. E, 2023, 108, 015002.
5. R. J. A. Lambourne, Relativity, Gravitation, and Cosmology, Cambridge University Press, Cambridge, 2010.
6. T. Frankel, The Geometry of Physics: An Introduction, Cambridge University Press, Cambridge, 3rd edn, 2012.
7. P. M. Chaikin and T. C. Lubensky, Principles of condensed matter physics, Cambridge University Press, 1995.
8. L. D. Landau and E. M. Lifshitz, Theory of Elasticity, Pergamon Press, New York, 1959.
9. E. B. Tadmor, R. E. Miller and R. S. Elliot, Continuum Mechanics and Thermodynamics, Cambridge University Press, Cambridge, 2012.
10. C. Klix, F. Ebert, F. Weysser, M. Fuchs, G. Maret and P. Keim, Glass elasticity from particle trajectories, Phys. Rev. Lett., 2012, 109, 178301.
11. C. L. Klix, G. Maret and P. Keim, Discontinuous shear modulus determines the glass transition temperature, Phys. Rev. X, 2015, 5, 041033.
12. M. Maier, A. Zippelius and M. Fuchs, Emergence of long-ranged stress correlations at the liquid to glass transition, Phys. Rev. Lett., 2017, 119, 265701.
13. M. Maier, A. Zippelius and M. Fuchs, Stress auto-correlation tensor in glass-forming isothermal fluids: From viscous to elastic response, J. Chem. Phys., 2018, 149, 084502.
14. F. Vogel, A. Zippelius and M. Fuchs, Emergence of Goldstone excitations in stress correlations of glass-forming colloidal dispersions, Europhys. Lett., 2019, 125, 68003.
15. A. Lemaître, Tensorial analysis of Eshelby stresses in 3d supercooled liquids, J. Chem. Phys., 2015, 143, 164515.
16. A. Lemaître, Stress correlations in glasses, J. Chem. Phys., 2018, 149, 104107.
17. L. Klochko, J. Baschnagel, J. P. Wittmer and A. N. Semenov, Long-range stress correlations in viscoelastic and glass-forming fluids, Soft Matter, 2018, 14, 6835.
18. L. Klochko, J. Baschnagel, J. Wittmer, H. Meyer, O. Benzerara and A. N. Semenov, Theory of length-scale dependent relaxation moduli and stress fluctuations in glass-forming and viscoelastic liquids, J. Chem. Phys., 2022, 156, 164505.
19. D. Steffen, L. Schneider, M. Müller and J. Rottler, Molecular simulations and hydrodynamic theory of nonlocal shear stress correlations in supercooled fluids, J. Chem. Phys., 2022, 157, 064501.
20. A. Kabla and G. Debrégeas, Local Stress Relaxation and Shear Banding in a Dry Foam under Shear, Phys. Rev. Lett., 2003, 90, 258303.
21. J. P. Wittmer, H. Xu, O. Benzerara and J. Baschnagel, Fluctuation-dissipation relation between shear stress relaxation modulus and shear stress autocorrelation function revisited, Mol. Phys., 2015, 113, 2881.
22. K. W. Desmond and E. R. Weeks, Measurement of Stress Redistribution in Flowing Emulsions, Phys. Rev. Lett., 2015, 115, 098302.
23. G. Picard, A. Ajdari, F. Lequeux and L. Bocquet, Elastic consequences of a single plastic event: A step towards the microscopic modeling of the flow of yield stress fluids, Eur. Phys. J. E: Soft Matter Biol. Phys., 2004, 15, 371.
24. L. Bocquet, A. Colin and A. Ajdari, Kinetic theory of plastic flow in soft glassy materials, Phys. Rev. Lett., 2009, 103, 036001.
25. J. Chattoraj and A. Lemaître, Elastic Signature of Flow Events in Supercooled Liquids Under Shear, Phys. Rev. Lett., 2013, 111, 066001.
26. H. Mizuno, S. Mossa and J. L. Barrat, Measuring spatial distribution of the local elastic modulus in glasses, Phys Rev E, 2013, 87, 042306.
27. A. Nicolas, J. Rottler and J. L. Barrat, Spatiotemporal correlations between plastic events in the shear flow of athermal amorphous solids, Eur. Phys. J. E: Soft Matter Biol. Phys., 2014, 37, 50.
28. E. Flenner and G. Szamel, Long-range spatial correlations of particle displacements and the emergence of elasticity, Phys. Rev. Lett., 2015, 114, 025501.
29. B. Illing, S. Fritschi, D. Hajnal, C. Klix, P. Keim and M. Fuchs, Strain pattern in supercooled liquids, Phys. Rev. Lett., 2016, 117, 208002.
30. M. Hassani, E. M. Zirdehi, K. Kok, P. Schall, M. Fuchs and F. Varnik, Long-range strain correlations in 3d quiescent glass forming liquids, Europhys. Lett., 2018, 124, 18003.
31. C. Liu, G. Biroli, D. R. Reichman and G. Szamel, Dynamics of liquids in the large-dimensional limit, Phys Rev E, 2021, 104, 054606.
32. R. N. Chacko, F. P. Landes, G. Biroli, O. D. A. J. Liu and D. R. Reichman, Elastoplasticity Mediates Dynamical Heterogeneity Below the Mode Coupling Temperature, Phys. Rev. Lett., 2021, 127, 048002.
33. M. P. Allen and D. J. Tildesley, Computer Simulation of Liquids, Oxford University Press, Oxford, 2nd edon; 2017.
34. J. P. Hansen and I. R. McDonald, Theory of simple liquids, Academic Press, New York, 3rd edn, 2006.
35. D. Forster, Hydrodynamic Fluctuations, Broken Symmetry, and Correlation Functions, Perseus Books, New York, 1995.
36. J. P. Wittmer, H. Xu, P. Polińska, F. Weysser and J. Baschnagel, Shear modulus of simulated glass-forming model systems: Effects of boundary condition, temperature and sampling time, J. Chem. Phys., 2013, 138, 12A533.
37. L. Klochko, J. Baschnagel, J. P. Wittmer and A. N. Semenov, Relaxation dynamics in supercooled oligomer liquids: From shear-stress fluctuations to shear modulus and structural correlations, J. Chem. Phys., 2019, 151, 054504.
38. G. George, L. Klochko, A. N. Semenov, J. Baschnagel and J. P. Wittmer, Ensemble fluctuations matter for variances of macroscopic variables, Eur. Phys. J. E: Soft Matter Biol. Phys., 2021, 44, 13.
39. G. George, L. Klochko, A. N. Semenov, J. Baschnagel and J. P. Wittmer, Fluctuations of non-ergodic stochastic processes, Eur. Phys. J. E: Soft Matter Biol. Phys., 2021, 44, 54.
40. J. D. Ferry, Viscoelastic properties of polymers, John Wiley & Sons, New York, 1980.
41. M. Rubinstein and R. H. Colby, Polymer Physics, Oxford University Press, Oxford, 2003.
42. M. Doi and S. F. Edwards, The Theory of Polymer Dynamics, Clarendon Press, Oxford, 1986.
43. J. K. G. Donth, The Glass Transition: Relaxation dynamics in liquids and disordered materials, Springer, Berlin-Heidelberg, 2001.
44. A. Argon and H. Kuo, Plastic flow in a disordered bubble raft (an analog of a metallic glass), Mater. Sci. Eng., 1979, 39, 101.
45. D. Rodney, A. Tanguy and D. Vandembroucq, Modeling the mechanics of amorphous solids at different length scale and time scale, Modell. Simul. Mater. Sci. Eng., 2011, 19, 083001.
46. A. Nicolas, E. Ferrero, K. Martens and J. L. Barrat, Deformation and flow of amorphous solids: Insights from elastoplastic models, Rev. Mod. Phys., 2018, 90, 045006.
47. M. L. Falk and J. S. Langer, Dynamics of viscoplastic deformation in amorphous solids, Phys Rev E, 1998, 57, 7192.
48. L. Klochko, J. Baschnagel, J. P. Wittmer and A. N. Semenov, Relaxation moduli of glass-forming systems: temperature effects and fluctuations, Soft Matter, 2021, 17, 7867.
49. A. Ninarello and L. Berthier, Coslovich D. Models and Algorithms for the Next Generation of Glass Transition Studies, Phys. Rev. X, 2017, 7, 021039.
50. J. P. Wittmer, A. N. Semenov and J. Baschnagel Strain fluctuations at finite wavevectors for isotropic colloidal glasses using natural rotated coordinates. in preparation. 2023.
51. J. P. Vest, G. Tarjus and P. Viot, Mode-coupling approach for the slow dynamics of a liquid on a spherical substrate, J. Chem. Phys., 2015, 143, 084505.
52. F. Turci, G. Tarjus and C. P. Royall, From Glass Formation to Icosahedral Ordering by Curving Three-Dimensional Space, Phys. Rev. Lett., 2017, 118, 215501.
53. H. Mizuno and S. Mossa, Mat Phys., 2019, 22, 43604.
54. J. P. Wittmer, A. N. Semenov and J. Baschnagel, Different spatial correlation functions for non-ergodic stochastic processes of macroscopic systems, Eur. Phys. J. E: Soft Matter Biol. Phys., 2022, 45, 65.
55. W. H. Press, S. A. Teukolsky, W. T. Vetterling and B. P. Flannery, Numerical Recipes in FORTRAN: the art of scientific computing, Cambridge University Press, Cambridge, 1992.
56. M. Abramowitz and I. A. Stegun, Handbook of Mathematical Functions, Dover, New York, 1964.
57. J. P. Wittmer, H. Xu, P. Polińska, C. Gillig, J. Helfferich and F. Weysser, et al., Compressibility and pressure correlations in isotropic solids and fluids, Eur. Phys. J. E: Soft Matter Biol. Phys., 2013, 36, 131.
58. J. P. Wittmer, H. Xu and J. Baschnagel, Shear-stress relaxation and ensemble transformation of shear-stress auto-correlation functions revisited, Phys Rev E, 2015, 91, 022107.
59. J. F. Lutsko, Stress and elastic constants in anisotropic solids: Molecular dynamics techniques, J. Appl. Phys., 1988, 64, 1152.
60. J. F. Lutsko, Generalized expressions for the calculation of elastic constants by computer simulation, J. Appl. Phys., 1989, 65, 2991.

---

Footnotes

† See, e.g., the wikipedia entries on quadrupoles and general multipolar expansions as used, say, in electrostatics. For the planar harmonic basis functions cos(pθ) or sin(pθ) a monopole corresponds to p = 0, a dipole to p = 1, a quadrupole to p = 2 and an octupole to p = 4.

‡ Let us impose in a rotated coordinate system at α = −π/4 a symmetric source tensor with a finite “shear” and vanishing diagonal components . Using eqn (A2) this implies S₁₁ = −S₂₂ = s and S₁₂ = 0 in the original coordinate system at α = 0.

§ The omitted constant terms correspond to localized δ(r) contributions to strain correlation functions in real space. For example, such a contribution to c₁₂₁₂(r) is .

¶ The elastic modulus tensor E_αβλδ(q) for isotropic bodies at finite wavevector q is not only characterized by the longitudinal modulus L(q) and the shear modulus G(q) but also by a third modulus M(q) called the “mixed modulus”.¹⁷ Note that , and for isotropic bodies in NRC. It is not possible to determine M(q) solely using strain fluctuations. This requires the additional measurement of stress fields in NRC. M(q) may then be obtained using either the appropriate stress–strain or stress–stress correlation functions in reciprocal space.^17,50

|| The infinite integral eqn (11.4.16) of ref. 56 is expressed in terms of two coefficients μ and ν which in our case take the (real) values μ = 1 and ν = p. It is stated that eqn (11.4.16) holds for , being consistent with the condition p > −2 noted in eqn (B17), but also that , being at first sight in conflict with μ = 1. However, the integral divergence for t → ∞ for μ > 1/2 is fictitious as, e.g., discussed in the Wikipedia entry on “Oscillatory integrals” where it is noted that “Oscillatory integrals make rigorous many arguments that, on a naive level, appear to use divergent integrals”. Note that eqn. (B18) also holds for p = 0 and finite r > 0 which follows from the fact that the FT of a constant is zero everywhere except at the origin. See ref. 4 for an alternative more straight-forward but also more lengthy derivation of eqn. (B18) using the asymptotic behavior of the confluent hypergeometric Kummer function M(a, b, z).

---