# Correlations of tensor field components in isotropic systems
with an application to stress correlations in elastic bodies

J.P. Wittmer  Institut Charles Sadron, Université de Strasbourg & CNRS, 23 rue
du Loess, 67034 Strasbourg Cedex, France  A.N. Semenov  Institut Charles
Sadron, Université de Strasbourg & CNRS, 23 rue du Loess, 67034 Strasbourg
Cedex, France  J. Baschnagel  Institut Charles Sadron, Université de
Strasbourg & CNRS, 23 rue du Loess, 67034 Strasbourg Cedex, France

###### Abstract

Correlation functions of components of second-order tensor fields in isotropic
systems can be reduced to an isotropic forth-order tensor field characterized
by a few invariant correlation functions (ICFs). It is emphasized that
components of this field depend in general on the coordinates of the field
vector variable and thus on the orientation of the coordinate system. These
angular dependencies are distinct from those of ordinary anisotropic systems.
As a simple example of the procedure to obtain the ICFs we discuss
correlations of time-averaged stresses in isotropic glasses where only one ICF
in reciprocal space becomes a finite constant $e$ for large sampling times and
small wavevectors. It is shown that $e$ is set by the typical size of the
frozen-in stress components normal to the wavevectors, i.e. it is caused by
the symmetry breaking of the stress for each independent configuration. Using
the presented general mathematical formalism for isotropic tensor fields this
finding explains in turn the observed long-range stress correlations in real
space. Under additional but rather general assumptions $e$ is shown to be
given by a thermodynamic quantity, the equilibrium Young modulus $E$. We thus
relate for certain isotropic amorphous bodies the existence of finite Young or
shear moduli to the symmetry breaking of a stress component in reciprocal
space.

##  I Introduction

#### Tensorial foundation of science and engineering.

The fundamental laws of physics and the constitutive relations of engineering
are formulated in terms of tensors and tensor fields (assigning a tensor to
each point of the mathematical space) McConnell (2021); Schultz-Piszachich
(1977); Lambourne (2010); Tadmor et al. (2012) which by construction
guarantees them to hold independently of the specific coordinate system used
for their description. Moreover, also instabilities and failure in materials
science and engineering, e.g., for granular piles and silos Nedderman (1992)
or plastic deformation in soft or glassy materials Argon and Kuo (1979); Falk
and Langer (1998); Picard et al. (2004); Voigtmann (2014); Fielding (2014),
must be described by appropriate tensorial invariants of tensor fields and
this should also be crucial in principle for mesoscopic computer models Rodney
et al. (2011); Nicolas et al. (2018) of localized plastic failure of a broad
range of systems.

#### Isotropic systems.

Isotropic systems, such as generic isotropic elastic bodies Landau and
Lifshitz (1959); Tadmor et al. (2012), simple and complex fluids Hansen and
McDonald (2006); Rubinstein and Colby (2003), amorphous metals and glasses
Donth (2001), polymer foams and networks Rubinstein and Colby (2003) or, as a
matter of fact, our entire universe Lambourne (2010) are described (at least
on some scales) by isotropic tensors and isotropic tensor fields McConnell
(2021); Tadmor et al. (2012); Schultz-Piszachich (1977). Let us assume for
simplicity that the system is not only isotropic but also spatially
homogeneous foo (a), achiral Landau and Lifshitz (1959), stationary in time
and embedded in a $d$-dimensional Euclidean vector space described by an
orthonormal Cartesian tensor basis Tadmor et al. (2012). A point in this
vector space is either called $\bf r$ (real space) or $\bf q$ (reciprocal
space). It is well known that the components of isotropic tensors remain
unchanged under an orthogonal coordinate transformation Tadmor et al. (2012).
For instance, the component $E_{1212}$ of the forth-order elastic modulus
tensor $E_{\alpha\beta\gamma\delta}$ of an isotropic body Landau and Lifshitz
(1959); Tadmor et al. (2012) is always given by the shear modulus $\mu$, i.e.
an invariant material property. (See Sec. II.4 and Sec. V.3 for details.)
Interestingly, this does in general not hold for the components of isotropic
tensor fields Schultz-Piszachich (1977) which may depend (in real space) not
only on the length $r$ of the field vector $\bf r$ but also on the
(normalized) coordinate dependent components $\hat{r}_{\alpha}$ of its
direction $\hat{\mbox{$\bf r$}}=\mbox{$\bf r$}/r$. This implies that
mathematically and physically legitimate isotropic tensor field components may
depend on the orientation of the coordinate system and this, as we shall see,
holds in a related manner both in real and in reciprocal space. As we shall
emphasize these angular dependencies differ from those of ordinary anisotropic
systems with frame-invariant angular-dependent material functions, say for
crystalline solids Landau and Lifshitz (1959).

Figure 1: Autocorrelation function $\bar{c}_{1212}(\mbox{$\bf r$})$ of time-
averaged shear-stress fields $\bar{\sigma}_{12}(\mbox{$\bf r$})$ of a two-
dimensional isotropic elastic body: (a) Unrotated frame with coordinates
$(r_{1},r_{2})$, (b) frame ($r_{1}^{\prime},r_{2}^{\prime})$ rotated by an
angle $\alpha=30^{\circ}$ (rotations marked by “$\prime$”). Albeit the system
is isotropic, the correlation function is strongly angle dependent revealing
an octupolar symmetry and depends, moreover, on the orientation of the
coordinate system. While each pixel corresponds in (a) and (b) to the same
spatial position $\bf r$, the correlation functions differ by the angle
$\alpha$. $\bar{c}_{1212}(\mbox{$\bf r$})$ is negative (blue) along the axes
and positive (red) along the bisection lines of the respective axes. The
theoretical prediction Eq. (1) is indicated by thin lines.

#### Invariant correlation functions.

Tensor fields are probed experimentally or in computer simulations by means of
correlation functions of some of their components. For instance, as shown in
Fig. 1 for an isotropic elastic body discussed in more detail below, one may
investigate the spatial correlations of the shear-stress component of the
stress tensor field Picard et al. (2004); Maier et al. (2017, 2018); Vogel et
al. (2019); Klochko et al. (2018); Lemaître (2015, 2018); Chowdhury et al.
(2016); Klochko et al. (2022); Steffen et al. (2022). For isotropic systems
such correlation functions must be isotropic tensor fields. They may thus
depend on the orientation of the coordinate systems as demonstrated by the
example given in panel (b) of Fig. 1. Another example of current interest are
the correlations of strain tensor field components of a broad range of
isotropic systems which have also been shown to reveal “anisotropic”
correlation functions Picard et al. (2004); Bocquet et al. (2009); Nicolas et
al. (2014); Flenner and Szamel (2015); Desmond and Weeks (2015); Hassani et
al. (2018); Liu et al. (2021); Chacko et al. (2021). Importantly, correlation
functions of tensor components describe the linear response due to a small
perturbation Hansen and McDonald (2006), say an inclusion in an elastic body
Eshelby (1957, 1959). Such a tensorial response field may thus be angle-
dependent for systems and source terms which are both perfectly isotropic.
(See Sec. II.7 for a comparison of tensorial response fields and associated
correlation fields in isotropic systems.) As we shall remind Lemaître (2015)
and emphasize in this study, while the isotropy of the system may not be
manifested by one correlation function of tensor field components it is
nevertheless present in the general mathematical structure of the complete set
of all correlation functions of the investigated tensor field. While this
makes the interpretation of observed tensor field pattern, like the ones given
in Fig. 1, more intricate, the good new is that this complete set of
correlation functions is determined by a small number of “Invariant
Correlation Functions” (ICFs) Lemaître (2015); Maier et al. (2018); Klochko et
al. (2018). We emphasize here that theory and computational studies should
focus on these ICFs and this in a first step in reciprocal space (cf. Appendix
A). The real-space correlations are then obtained by inverse Fourier
transformation (FT).

#### Specific case considered.

We demonstrate the general procedure by an analysis of stress field
correlations in two-dimensional ($d=2$) simulated isotropic elastic bodies
Landau and Lifshitz (1959); Tadmor et al. (2012). While previous work
characterizes correlations of the instantaneous stress field Maier et al.
(2017, 2018); Vogel et al. (2019); Klochko et al. (2018); Chowdhury et al.
(2016); Klochko et al. (2022); Steffen et al. (2022) or of the stress field of
the system’s “inherent states” (local energy minima) Lemaître (2015, 2018), we
rather compute the time-averaged stress fields
$\bar{\sigma}_{\alpha\beta}(\mbox{$\bf q$})$ in reciprocal space ($\bf q$
being the wavevector) for each independent configuration $c$ and analyze their
correlation functions $\bar{c}_{\alpha\beta\gamma\delta}(\mbox{$\bf q$})$.
This is best done using “Natural Rotated Coordinates” (NRC) aligned with $\bf
q$ allowing the precise determination of the ICFs characterizing the isotropic
tensor field. Importantly, only one of these ICFs is shown in the hydrodynamic
limit and for sufficiently large sampling times $\Delta\tau$ to become a
finite constant $e>0$. This phenomenological constant characterizes the
typical size of the frozen stress components normal to each wavevector $\bf
q$, i.e. it measures the (continuous) symmetry breaking of the stress field in
reciprocal space for each independent configuration. This finding and a proper
treatment of tensor field correlations reflecting the material isotropy
directly imply that the stress correlations in real space must be long-ranged.
For instance, the correlation function of the shear-stress field must decay as

$\bar{c}^{\prime}_{1212}(\mbox{$\bf r$})\simeq-\frac{e}{4\pi
r^{2}}\cos[4(\theta-\alpha)]$ (1)

for sufficiently large $r=|\mbox{$\bf r$}|$ with $\theta$ being the angle of
the field vector $\bf r$ in the unrotated physical system, Fig. 1(a), and
$\alpha$ the rotation angle of the coordinate system, Fig. 1(b). (Rotations of
the coordinate system are marked by primes “$\prime$”.) We thus confirm recent
computer simulations on stress correlations in binary Lennard-Jones glasses
Lemaître (2014, 2015, 2017, 2018); Klochko et al. (2022) and more general
theoretical considerations Maier et al. (2017, 2018); Vogel et al. (2019);
Klochko et al. (2018) on supercooled liquids and amorphous elastic bodies.

#### Outline.

We begin by reviewing in Sec. II general features of tensor fields relevant
for isotropic and achiral systems McConnell (2021); Schouten (2015); Schultz-
Piszachich (1977); Tadmor et al. (2012). The main computational points (model
system, data production) are summarized in Sec. III before we turn in Sec. IV
to our central numerical results. Taking advantage of recent theoretical
studies Klochko et al. (2018, 2022); Maier et al. (2017, 2018); Vogel et al.
(2019) we explain in Sec. V why the only phenomenological parameter $e$ needed
to fit our data should be similar to a thermodynamic quantity, the equilibrium
Young modulus $E$. A summary of the presented work and an outlook are given in
Sec. VI. Appendix A summarizes some properties of FTs Abramowitz and Stegun
(1964); Press et al. (1992) while Appendix B addresses the FT of the relevant
correlation functions in $d=2$. More details on our simulation model and on
the computation of the local stress fields are given in Appendix C and
Appendix D.

##  II Review of isotropic tensor fields

###  II.1 Introduction

#### Generalities.

Familiarity with the general ideas and notations of tensor algebra and
analysis, as developed systematically in the standard textbooks McConnell
(2021); Schouten (2015); Schultz-Piszachich (1977); Tadmor et al. (2012), is
taken for granted. We remind that a tensor field assigns a tensor to each
point of the mathematical space, in our case a $d$-dimensional Euclidean
vector space Schultz-Piszachich (1977). An element of this vector space is
denoted by the “spatial position” $\bf r$ in real space or by the “wavevector”
$\bf q$ for the corresponding Fourier transformed reciprocal space. The
relations for tensor fields are formulated below in reciprocal space since
this is more convenient both on theoretical and numerical grounds due to the
assumed spatial translational invariance. The properties of the corresponding
real space tensor field are then obtained by inverse FT.

#### Tensor components and basis.

We assume for simplicity Cartesian coordinates with an orthonormal basis
$\\{\mbox{$\bf e$}_{1},\ldots,\mbox{$\bf e$}_{d}\\}$ McConnell (2021);
Schultz-Piszachich (1977); Tadmor et al. (2012). Greek letters
$\alpha,\beta,\ldots$ are used for the indices of the tensor (field)
components. A twice repeated index $\alpha$ is summed over the values
$1,\ldots,d$, e.g., $\mbox{$\bf q$}=q_{\alpha}\mbox{$\bf e$}_{\alpha}$ with
$q_{\alpha}$ standing for the vector coordinates. This work is chiefly
concerned with tensors

$\mbox{$\bf T$}^{(o)}=T_{\alpha_{1}\ldots\alpha_{o}}\mbox{$\bf
e$}_{\alpha_{1}}\ldots\mbox{$\bf e$}_{\alpha_{o}}$ (2)

of “order” (or “rank”) $o=2$ and $o=4$ and their corresponding tensor fields
with components depending either on $\bf r$ or $\bf q$. As common we refer to
a tensor (field) $\mbox{$\bf T$}^{(o)}$ by indicating its components
$T_{\alpha_{1}\ldots\alpha_{o}}$. The order of a component is given by the
number of suffixes. Note that

$T_{\alpha_{1}\ldots\alpha_{o}}(\mbox{$\bf q$})={\cal
F}[T_{\alpha_{1}\ldots\alpha_{o}}(\mbox{$\bf r$})]$ (3)

for the $d^{o}$ components in real and reciprocal space.

#### Transforms.

We consider linear orthogonal coordinate transformations (marked by “$\ast$”)
$\mbox{$\bf e$}_{\alpha}^{\ast}=c_{\alpha\beta}\mbox{$\bf e$}_{\beta}$ with
matrix coefficients $c_{\alpha\beta}$ given by the direction cosine
$c_{\alpha\beta}\equiv\cos(\mbox{$\bf e$}_{\alpha}^{\ast},\mbox{$\bf
e$}_{\beta})$ Schultz-Piszachich (1977).
$c_{\alpha\beta}=\delta_{\alpha\beta}$ if nothing is changed. For a simple
reflection of, say, the $1$-axis and a rotation in the $12$-plane by an angle
$\theta$ we have, respectively,

reflection $\displaystyle:$ $\displaystyle c_{1\beta}=-\delta_{1\beta}$ (4)
rotation $\displaystyle:$
$\displaystyle\begin{array}[t]{ll}c_{11}=c_{22}&=\cos(\theta)\\\
c_{12}=-c_{21}&=\sin(\theta)\end{array}$ (7)

and $c_{\alpha\beta}=\delta_{\alpha\beta}$ for all other indices. We remind
that Schultz-Piszachich (1977)

$T^{\ast}_{\alpha_{1}\ldots\alpha_{o}}(\mbox{$\bf
q$})=c_{\alpha_{1}\nu_{1}}\ldots
c_{\alpha_{o}\nu_{o}}T_{\nu_{1}\ldots\nu_{o}}(\mbox{$\bf q$})$ (8)

under a general orthogonal transform. For a reflection of the $1$-axis we thus
have, e.g.,

$\displaystyle T_{12}^{\ast}(\mbox{$\bf q$})$ $\displaystyle=$ $\displaystyle-
T_{12}(\mbox{$\bf q$}),T_{11}^{\ast}(\mbox{$\bf q$})=T_{11}(\mbox{$\bf q$}),$
(9) $\displaystyle T^{\ast}_{1222}(\mbox{$\bf q$})$ $\displaystyle=$
$\displaystyle-T_{1222}(\mbox{$\bf q$}),T^{\ast}_{1221}(\mbox{$\bf
q$})=T_{1221}(\mbox{$\bf q$}),$

i.e. quite generally we have sign inversion for an odd number of indices
$\alpha,\beta,\ldots$ equal to the index of the inverted axis. Please note
that the field vector $\mbox{$\bf q$}=q_{\alpha}\mbox{$\bf
e$}_{\alpha}=q_{\alpha}^{\ast}\mbox{$\bf e$}_{\alpha}^{\ast}$ remains
unchanged by these “passive” transforms albeit its coordinates change.

#### Outline.

We come back to this issue in the next subsection, Sec. II.2. The symmetries
of tensors and tensor fields relevant for the present work are summarized in
Sec. II.3. Isotropic tensors are presented in Sec. II.4, general isotropic
tensor fields of order $1\leq o\leq 4$ in Sec. II.5 and, more specifically,
forth-order isotropic tensor fields $T_{\alpha\beta\gamma\delta}(\mbox{$\bf
q$})$ for $d=2$ in Sec. II.6. Finally, Sec. II.7 outlines some general
relations for second-order tensor fields $R_{\alpha\beta}(\mbox{$\bf r$})$
corresponding to the “linear response” caused by a point-like “source term”
$S_{\alpha\beta}(\mbox{$\bf r$})=s_{\alpha\beta}\delta(\mbox{$\bf r$})$ in
real space.

###  II.2 Isotropic tensors and tensor fields

#### Isotropic tensors.

Isotropic systems are described by “isotropic tensors” and “isotropic tensor
fields”. Components of an isotropic tensor remain unchanged by any orthogonal
coordinate transformation Schultz-Piszachich (1977); Tadmor et al. (2012),
i.e.

$T_{\alpha_{1}\ldots\alpha_{o}}^{\ast}=T_{\alpha_{1}\ldots\alpha_{o}}.$ (10)

As noted at the end of Sec. II.1 the sign of tensor components change for a
reflection of one axis if the number of indices equal to the inverted axis is
odd. Consistency with Eq. (10) implies that all tensor components with an odd
number of equal indices must vanish, e.g.,

$T_{12}=T_{1112}=T_{1222}=T_{1234}=T_{1344}=0.$ (11)

#### Isotropic tensor fields.

The corresponding isotropy condition for tensor fields is given by Schultz-
Piszachich (1977)

$T_{\alpha_{1}\ldots\alpha_{o}}^{\ast}(q_{1},\ldots,q_{d})=T_{\alpha_{1}\ldots\alpha_{o}}(q_{1}^{\ast},\ldots,q_{d}^{\ast})$
(12)

which reduces to Eq. (10) for $\mbox{$\bf q$}={\bf 0}$. Please note that the
fields on the left hand side of Eq. (12) are evaluated with the original
coordinates while the fields on the right hand side 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 $\mbox{$\bf
q$}=(q_{1},\ldots,q_{d})$ while the right hand fields are computed at the
“actively transformed” vector $\mbox{$\bf
q$}^{\ast}=(q_{1}^{\ast},\ldots,q_{d}^{\ast})$. It is for this reason that Eq.
(11) does not hold in general for tensor fields, i.e. finite components with
an odd number of equal indices, e.g., $T_{1222}(\mbox{$\bf q$})\neq 0$, are
possible in principle for finite wavevectors.

#### Natural Rotated Coordinates.

Fortunately, there are convenient coordinates, called “Natural Rotated
Coordinates” (NRC), where the nice symmetry Eq. (11) for isotropic tensors can
be also used for tensor fields. Let us assume that the wavevector $\bf q$
points into the direction of one of the axes, say, $\mbox{$\bf
q$}=q\delta_{1\beta}$ with $q=|\mbox{$\bf q$}|$. (This may be achieved by a
first rotation of the coordinate system.) Let us denote by “$\neg$” arbitrary
inversions of axes in this frame. We may thus rewrite Eq. (12) as
$T_{\alpha_{1}\ldots\alpha_{o}}^{\neg}(q_{1})=T_{\alpha_{1}\ldots\alpha_{o}}(q_{1}^{\neg})$
since $q_{2}=q_{2}^{\neg}=\ldots=q_{d}=q_{d}^{\neg}=0$. If we now assume in
addition that $T_{\alpha_{1}\ldots\alpha_{o}}(\mbox{$\bf q$})$ is an even
function of its field variable $\bf q$ this becomes

$T_{\alpha_{1}\ldots\alpha_{o}}^{\neg}(q)=T_{\alpha_{1}\ldots\alpha_{o}}(q),$
(13)

i.e. both sides only depend on the same scalar parameter. We may thus use for
each $q$ the same reasoning as for tensor components.

#### Product theorem for isotropic tensor fields.

Let us state a useful theorem for a general tensor field $\mbox{$\bf
C$}(\mbox{$\bf q$})=\mbox{$\bf A$}(\mbox{$\bf q$})\otimes\mbox{$\bf
B$}(\mbox{$\bf q$})$ with $\mbox{$\bf A$}(\mbox{$\bf q$})$ and $\mbox{$\bf
B$}(\mbox{$\bf q$})$ being two isotropic tensor fields and $\otimes$ standing
either for an outer product, e.g. $C_{\alpha\beta\gamma\delta}(\mbox{$\bf
q$})=A_{\alpha\beta}(\mbox{$\bf q$})B_{\gamma\delta}(\mbox{$\bf q$})$, or an
inner product, e.g. $C_{\alpha\beta\gamma\delta}(\mbox{$\bf
q$})=A_{\alpha\beta\gamma\nu}(\mbox{$\bf q$})B_{\nu\delta}(\mbox{$\bf q$})$.
Hence,

$\displaystyle\mbox{$\bf C$}^{\ast}(\mbox{$\bf q$})$ $\displaystyle=$
$\displaystyle\left(\mbox{$\bf A$}(\mbox{$\bf q$})\otimes\mbox{$\bf
B$}(\mbox{$\bf q$})\right)^{\ast}=\mbox{$\bf A$}^{\ast}(\mbox{$\bf
q$})\otimes\mbox{$\bf B$}^{\ast}(\mbox{$\bf q$})$ (14) $\displaystyle=$
$\displaystyle\mbox{$\bf A$}(\mbox{$\bf q$}^{\ast})\otimes\mbox{$\bf
B$}(\mbox{$\bf q$}^{\ast})=\mbox{$\bf C$}(\mbox{$\bf q$}^{\ast})$

using in the second step a general property of tensor (field) products, due to
Eq. (8), and in the third step Eq. (12) for the fields $\mbox{$\bf
A$}(\mbox{$\bf q$})$ and $\mbox{$\bf B$}(\mbox{$\bf q$})$ where $\mbox{$\bf
q$}^{\ast}$ stands for the “actively” transformed field position. We have thus
demonstrated that $\mbox{$\bf C$}(\mbox{$\bf q$})$ is also an isotropic tensor
field. This theorem allows to construct isotropic tensor fields from known
isotropic tensor fields $\mbox{$\bf A$}(\mbox{$\bf q$})$ and $\mbox{$\bf
B$}(\mbox{$\bf q$})$.

#### Multilinear forms.

Alternatively, isotropic tensor fields may be constructed using multilinear
forms Schultz-Piszachich (1977); Tadmor et al. (2012)

$L(\mbox{$\bf v$}^{1},\ldots,\mbox{$\bf v$}^{o};\mbox{$\bf
q$})=T_{\alpha_{1}\ldots\alpha_{o}}(\mbox{$\bf q$})\ v^{1}_{\alpha_{1}}\ldots
v^{o}_{\alpha_{o}}$ (15)

of order $o$ where $L$ stands for a linear and (first order) homogeneous
functional of a $d$-dimensional vector space Schultz-Piszachich (1977); Tadmor
et al. (2012) and $\mbox{$\bf v$}^{1},\ldots,\mbox{$\bf v$}^{o}$ are $o$
arbitrary vectors of this vector space (with superscripts exceptionally used
here for the numbering of these vectors). For tensor fields this functional
also depends on the tensor field vector $\bf q$. The goal is then to construct
generic isotropic tensor fields associated with multilinear forms. This is
done in a first step by means of additive terms of all possible scalars formed
with the vectors $\mbox{$\bf v$}^{1},\ldots,\mbox{$\bf v$}^{o}$ and $\bf q$,
e.g., inner products $\mbox{$\bf v$}^{1}\cdot\mbox{$\bf v$}^{2}$, $\mbox{$\bf
v$}^{1}\cdot\mbox{$\bf q$}$ or $\mbox{$\bf q$}\cdot\mbox{$\bf q$}$ or triple
products such as $[\mbox{$\bf v$}^{1}\mbox{$\bf v$}^{2}\mbox{$\bf q$}]$. In a
second step all contributions are eliminated which are yet incompatible with
Eq. (12) and other imposed symmetries. We shall illustrate this below in Sec.
II.5.

#### Kronecker and Levi-Civita tensors.

We note for later convenience that the Kronecker symbol $\delta_{\alpha\beta}$
is an invariant tensor, $\delta_{\alpha\beta}^{\ast}=\delta_{\alpha\beta}$,
for any orthogonal transform Schultz-Piszachich (1977). As a consequence, any
tensor field, only containing additive terms such as
$i(q)\delta_{\alpha\beta}\delta_{\gamma\delta}$ with $i(q)$ being an invariant
scalar, is an isotropic tensor field. The same applies for tensor fields with
terms containing one or several factors ${\hat{q}}_{\alpha}$ since in
agreement with Eq. (12) this implies, e.g,

$\left(i(q){\hat{q}}_{\alpha}{\hat{q}}_{\beta}\delta_{\gamma\delta}\right)^{\ast}=i(q){\hat{q}}_{\alpha}^{\ast}{\hat{q}}_{\beta}^{\ast}\delta_{\gamma\delta}$
(16)

where Eq. (8) was used. The situation is more intricate for terms containing
the Levi-Civita (“permutation”) tensor $\varepsilon_{\alpha\beta\gamma}$
McConnell (2021) which is only invariant for the rotation subgroup but in
general not for reflections Schultz-Piszachich (1977).

###  II.3 Assumed symmetries

All second-order tensors in this work are symmetric,
$T_{\alpha\beta}=T_{\beta\alpha}$, and the same applies for the corresponding
tensor fields in either $\bf r$\- or $\bf q$-space. This is, e.g., the case
for the stress fields $\sigma_{\alpha\beta}(\mbox{$\bf q$})$. We assume for
all forth-order tensor fields that

$\displaystyle T_{\alpha\beta\gamma\delta}(\mbox{$\bf q$})$ $\displaystyle=$
$\displaystyle T_{\beta\alpha\gamma\delta}(\mbox{$\bf
q$})=T_{\alpha\beta\delta\gamma}(\mbox{$\bf q$})$ (17) $\displaystyle
T_{\alpha\beta\gamma\delta}(\mbox{$\bf q$})$ $\displaystyle=$ $\displaystyle
T_{\gamma\delta\alpha\beta}(\mbox{$\bf q$})\mbox{ and }$ (18) $\displaystyle
T_{\alpha\beta\gamma\delta}(\mbox{$\bf q$})$ $\displaystyle=$ $\displaystyle
T_{\alpha\beta\gamma\delta}(-\mbox{$\bf q$}).$ (19)

Let us remind that forth-order tensor fields are often constructed by taking
outer products Tadmor et al. (2012) of second-order tensor fields. We
consider, e.g., correlation functions

$T_{\alpha\beta\gamma\delta}(\mbox{$\bf
q$})=\left<\hat{T}_{\alpha\beta}(\mbox{$\bf
q$})\hat{T}_{\gamma\delta}(-\mbox{$\bf q$})\right>$ (20)

with $\hat{T}_{\alpha\beta}(\mbox{$\bf q$})$ being an instantaneous (not
ensemble-averaged) second-order tensor field. Eq. (17) then follows from the
symmetry of the second-order tensor fields. Evenness, Eq. (19), is a necessary
condition for achiral systems. It implies that
$T_{\alpha\beta\gamma\delta}(\mbox{$\bf q$})$ is real if
$T_{\alpha\beta\gamma\delta}(\mbox{$\bf r$})$ is real and, moreover, Eq. (18)
for correlation functions since $\langle\hat{T}_{\alpha\beta}(\mbox{$\bf
q$})\hat{T}_{\gamma\delta}(-\mbox{$\bf
q$})\rangle=\langle\hat{T}_{\gamma\delta}(\mbox{$\bf
q$})\hat{T}_{\alpha\beta}(-\mbox{$\bf q$})\rangle$. As already emphasized, it
is assumed that all our systems are isotropic. This implies that Eqs. (10-12)
must hold for the ensemble-averaged tensor fields. Since our systems are also
achiral, Eq. (13) applies and tensor field components with an odd number of
equal indices must vanish if one axis points into the direction of the
wavevector. We consider in the following subsections isotropic tensors and
tensor fields respecting the above symmetries.

###  II.4 Isotropic tensors

Isotropic tensors of different order are discussed, e.g., in Sec. 2.5.6 of
Ref. Tadmor et al. (2012). While all tensors of odd order must vanish, we have

$\displaystyle T_{\alpha\beta}$ $\displaystyle=$ $\displaystyle
k_{1}\delta_{\alpha\beta},$ (21) $\displaystyle T_{\alpha\beta\gamma\delta}$
$\displaystyle=$ $\displaystyle
i_{1}\delta_{\alpha\beta}\delta_{\gamma\delta}+i_{2}\left(\delta_{\alpha\gamma}\delta_{\beta\delta}+\delta_{\alpha\delta}\delta_{\beta\gamma}\right)$
(22)

where $k_{1}$, $i_{1}$ and $i_{2}$ are invariant scalars. Please note that all
symmetries stated above hold, especially also Eq. (11). An example for an
isotropic second-order tensor is the isotropic stress tensor
$\sigma_{\alpha\beta}=-P\delta_{\alpha\beta}$ (with $P$ being the average
normal pressure) which is assumed in the present work. Note that the symmetry
Eq. (17) was used for the second relation, Eq. (22). Importantly, this implies
that only two coefficients are needed for a forth-order isotropic tensor. As
further discussed in Sec. V.3, the elastic modulus tensor
$E_{\alpha\beta\gamma\delta}$ is thus completely described by two elastic
moduli, say $\lambda$ and $\mu$, and the stress relaxation modulus tensor
$E_{\alpha\beta\gamma\delta}(\tau)$ by two relaxation functions, say the
“mixed relaxation function” $M(\tau)$ and the “shear-stress relaxation
function” $G(\tau)$ foo (b).

###  II.5 Tensor fields for isotropic achiral systems

We begin by summarizing the relevant isotropic tensor fields for $1\leq o\leq
4$ compatible with the assumed symmetries (cf. Sec. II.3). With $l_{n}(q)$,
$k_{n}(q)$, $j_{n}(q)$ and $i_{n}(q)$ being invariant scalar functions of $q$
we have

$\displaystyle T_{\alpha}(\mbox{$\bf q$})$ $\displaystyle=$ $\displaystyle
l_{1}(q)\ q_{\alpha}$ (23) $\displaystyle T_{\alpha\beta}(\mbox{$\bf q$})$
$\displaystyle=$ $\displaystyle k_{1}(q)\ \delta_{\alpha\beta}+k_{2}(q)\
q_{\alpha}q_{\beta}$ (24) $\displaystyle T_{\alpha\beta\gamma}(\mbox{$\bf
q$})$ $\displaystyle=$ $\displaystyle j_{1}(q)\
q_{\alpha}\delta_{\beta\gamma}+j_{2}(q)\ q_{\beta}\delta_{\alpha\gamma}$ (25)
$\displaystyle+$ $\displaystyle j_{3}(q)\
q_{\gamma}\delta_{\alpha\beta}+j_{4}(q)\ q_{\alpha}q_{\beta}q_{\gamma}$
$\displaystyle T_{\alpha\beta\gamma\delta}(\mbox{$\bf q$})$ $\displaystyle=$
$\displaystyle i_{1}(q)\ \delta_{\alpha\beta}\delta_{\gamma\delta}$
$\displaystyle+$ $\displaystyle i_{2}(q)\
\left(\delta_{\alpha\gamma}\delta_{\beta\delta}+\delta_{\alpha\delta}\delta_{\beta\gamma}\right)$
$\displaystyle+$ $\displaystyle i_{3}(q)\
\left(q_{\alpha}q_{\beta}\delta_{\gamma\delta}+q_{\gamma}q_{\delta}\delta_{\alpha\beta}\right)$
$\displaystyle+$ $\displaystyle i_{4}(q)\
q_{\alpha}q_{\beta}q_{\gamma}q_{\delta}$ $\displaystyle+$ $\displaystyle
i_{5}(q)\left(q_{\alpha}q_{\gamma}\delta_{\beta\delta}+q_{\alpha}q_{\delta}\delta_{\beta\gamma}+\right.$
$\displaystyle\hskip
25.6073pt\left.q_{\beta}q_{\gamma}\delta_{\alpha\delta}+q_{\beta}q_{\delta}\delta_{\alpha\gamma}\right)$

with $q_{\alpha}\equiv\mbox{$\bf q$}\cdot\mbox{$\bf e$}_{\alpha}$. Let us
first check that the stated relations are reasonable. All relations reduce
(continuously) for $\mbox{$\bf q$}\to{\bf 0}$ to the isotropic tensors stated
in Sec. II.4; all are, according to the discussion in the last paragraph of
Sec. II.2, isotropic tensor fields consistent with Eq. (12) and all symmetries
stated in Sec. II.3 for the second- and forth-order tensor fields are
satisfied. All tensor fields of even (odd) order are even (odd) with respect
to $\bf q$. Hence, tensor fields of odd order vanish for $\mbox{$\bf
q$}\to{\bf 0}$ consistently with Sec. II.4. Note that the terms due to the
invariants $k_{1}(q)$, $i_{1}(q)$ and $i_{2}(q)$ are strictly isotropic and,
hence, independent of the coordinate system. All other terms depend on the
components $q_{\alpha}$ and thus on the coordinate system.

Following Refs. Schultz-Piszachich (1977); Forster (1995) let us first show
that Eq. (24) holds. According to Eq. (15) one may represent a general second-
order tensor field by a bilinear form $L(\mbox{$\bf u$},\mbox{$\bf
v$};\mbox{$\bf q$})$. Invariant with respect to orthogonal transformations are
the scalars $\mbox{$\bf q$}\cdot\mbox{$\bf q$}=q^{2}$, $\mbox{$\bf
q$}\cdot\mbox{$\bf u$}$, $\mbox{$\bf q$}\cdot\mbox{$\bf v$}$ and,
additionally, the triple product $[\mbox{$\bf u$}\mbox{$\bf v$}\mbox{$\bf
q$}]$ for three-dimensional systems. We obtain thus the general bilinear form

$\displaystyle L(\mbox{$\bf u$},\mbox{$\bf v$};\mbox{$\bf q$})$
$\displaystyle=$ $\displaystyle k_{1}\mbox{$\bf u$}\cdot\mbox{$\bf
v$}+k_{2}(\mbox{$\bf q$}\cdot\mbox{$\bf u$})(\mbox{$\bf q$}\cdot\mbox{$\bf
v$})+k_{3}[\mbox{$\bf u$}\mbox{$\bf v$}\mbox{$\bf q$}]$ $\displaystyle=$
$\displaystyle T_{\alpha\beta}(\mbox{$\bf q$})u_{\alpha}v_{\beta}\ \mbox{ with
}$ $\displaystyle T_{\alpha\beta}(\mbox{$\bf q$})$ $\displaystyle=$
$\displaystyle
k_{1}\delta_{\alpha\beta}+k_{2}q_{\alpha}q_{\beta}+k_{3}\varepsilon_{\alpha\beta\gamma}q_{\gamma}$
(27)

with $k_{1}$, $k_{2}$ and $k_{3}$ being scalar coefficients. While the first
two terms of $T_{\alpha\beta}(\mbox{$\bf q$})$ are fine with respect to Eq.
(12), the last term must be eliminated since $\varepsilon_{\alpha\beta\gamma}$
changes sign for a reflection at one axis Schultz-Piszachich (1977).

The indicated relations for the other fields are obtained in a similar manner
Schultz-Piszachich (1977). Note that using the product theorem, Eq. (14), one
may obtain the isotropic tensor fields of third and forth order as sums of
products of lower-order isotropic tensor fields. For instance, let
$A_{\alpha\beta}(\mbox{$\bf q$})$ and $B_{\alpha\beta}(\mbox{$\bf q$})$ be two
second-order isotropic tensor fields according to Eq. (24). It is readily seen
that $A_{\alpha\beta}(\mbox{$\bf q$})B_{\gamma\delta}(\mbox{$\bf
q$})+A_{\gamma\delta}(\mbox{$\bf q$})B_{\alpha\beta}(\mbox{$\bf q$})$
immediately implies the first four terms of Eq. (II.5), i.e. $i_{5}(q)=0$ if
$T_{\alpha\beta\gamma\delta}(\mbox{$\bf q$})$ is only constructed from two
isotropic second-order tensors. For the indicated more general isotropic
forth-order tensors Eq. (II.5) with $i_{5}(q)\neq 0$ we have included
contributions due to products $T_{\alpha}(\mbox{$\bf
q$})T_{\beta\gamma\delta}(\mbox{$\bf q$})$ of isotropic first- and third-order
tensor fields.

Let us check that terms containing the Levi-Civita tensor
$\varepsilon_{\alpha\beta\gamma}$ cannot contribute additional terms to a
forth-order tensor field $T_{\alpha\beta\gamma\delta}(\mbox{$\bf q$})$ obeying
the assumed symmetries. The forth-order multilinear form $L(\mbox{$\bf
u$},\mbox{$\bf v$},\mbox{$\bf s$},\mbox{$\bf t$};\mbox{$\bf
q$})=T_{\alpha\beta\gamma\delta}(\mbox{$\bf
q$})u_{\alpha}v_{\beta}s_{\gamma}t_{\delta}$ may indeed apriori contain in
$d=3$ terms of products of invariants such as

$\displaystyle(\mbox{$\bf u$}\cdot\mbox{$\bf v$})\ \left[\mbox{$\bf
s$}\mbox{$\bf t$}\mbox{$\bf q$}\right]$ $\displaystyle=$
$\displaystyle\delta_{\alpha\beta}\varepsilon_{\gamma\delta\nu}q_{\nu}\ \ \ \
\ \ u_{\alpha}v_{\beta}s_{\gamma}t_{\delta},$ (28)
$\displaystyle\left[\mbox{$\bf u$}\mbox{$\bf v$}\mbox{$\bf q$}\right]\
\left[\mbox{$\bf s$}\mbox{$\bf t$}\mbox{$\bf q$}\right]$ $\displaystyle=$
$\displaystyle\varepsilon_{\alpha\beta\nu}\varepsilon_{\gamma\delta\mu}q_{\nu}q_{\mu}\
\ u_{\alpha}v_{\beta}s_{\gamma}t_{\delta}.$ (29)

Terms of the first type are disallowed for the same reason as argued for the
second-order isotropic tensor, Eq. (27). Note also that such terms would not
be compatible with Eq. (19). Terms of the the second type, Eq. (29), can be
expressed using Sarrus’ law as

$\displaystyle\varepsilon_{\alpha\beta\nu}\varepsilon_{\gamma\delta\mu}q_{\nu}q_{\mu}$
$\displaystyle=$
$\displaystyle\delta_{\alpha\gamma}\delta_{\beta\delta}q^{2}+\delta_{\alpha\delta}q_{\beta}q_{\gamma}+\delta_{\beta\gamma}q_{\alpha}q_{\delta}$
$\displaystyle-$
$\displaystyle\delta_{\alpha\delta}\delta_{\beta\gamma}q^{2}-\delta_{\alpha\gamma}q_{\beta}q_{\delta}-\delta_{\beta\delta}q_{\alpha}q_{\gamma}.$

This is, however, in conflict with Eq. (17) and Eq. (18). To enforce, e.g.,
the $\alpha\leftrightarrow\beta$-symmetry the multilinear form must also
contain a term $\left[\mbox{$\bf v$}\mbox{$\bf u$}\mbox{$\bf
q$}\right]\left[\mbox{$\bf s$}\mbox{$\bf t$}\mbox{$\bf q$}\right]$ which is
equal to $-\left[\mbox{$\bf u$}\mbox{$\bf v$}\mbox{$\bf
q$}\right]\left[\mbox{$\bf s$}\mbox{$\bf t$}\mbox{$\bf q$}\right]$. All
contributions of the second type needed for symmetry reasons thus exactly
cancel.

Finally, let us note that for physical reasons it is useful to rewrite for
finite wavevectors ($\mbox{$\bf q$}\neq{\bf 0})$ the isotropic forth-order
tensor field Eq. (II.5) in terms of the components
${\hat{q}}_{\alpha}={\hat{\bf q}}\cdot\mbox{$\bf e$}_{\alpha}$ of the
normalized wavevector ${\hat{\bf q}}$. It is thus convenient to bring in
factors of $q$ and to redefine $i_{3}(q)\to i_{3}(q)/q^{2}$, $i_{4}(q)\to
i_{4}(q)/q^{4}$ and $i_{5}(q)\to i_{5}(q)/q^{2}$. We thus rewrite Eq. (II.5)
as

$\displaystyle T_{\alpha\beta\gamma\delta}(\mbox{$\bf q$})$ $\displaystyle=$
$\displaystyle i_{1}(q)\ \delta_{\alpha\beta}\delta_{\gamma\delta}$
$\displaystyle+$ $\displaystyle
i_{2}(q)\left(\delta_{\alpha\gamma}\delta_{\beta\delta}+\delta_{\alpha\delta}\delta_{\beta\gamma}\right)$
$\displaystyle+$ $\displaystyle
i_{3}(q)\left({\hat{q}}_{\alpha}{\hat{q}}_{\beta}\delta_{\gamma\delta}+{\hat{q}}_{\gamma}{\hat{q}}_{\delta}\delta_{\alpha\beta}\right)$
$\displaystyle+$ $\displaystyle i_{4}(q)\
{\hat{q}}_{\alpha}{\hat{q}}_{\beta}{\hat{q}}_{\gamma}{\hat{q}}_{\delta}$
$\displaystyle+$ $\displaystyle
i_{5}(q)\left({\hat{q}}_{\alpha}{\hat{q}}_{\gamma}\delta_{\beta\delta}+{\hat{q}}_{\alpha}{\hat{q}}_{\delta}\delta_{\beta\gamma}+\right.$
$\displaystyle\hskip
25.6073pt\left.{\hat{q}}_{\beta}{\hat{q}}_{\gamma}\delta_{\alpha\delta}+{\hat{q}}_{\beta}{\hat{q}}_{\delta}\delta_{\alpha\gamma}\right).$

Now all $i_{n}(q)$ have the same physical units. As we shall see, the
$i_{n}(q)$ become often constant, $i_{n}(q)\to i_{n}$, or negligibly tiny for
sufficiently small (but finite) $q$.

###  II.6 Isotropic $T_{\alpha\beta\gamma\delta}(\mbox{$\bf q$})$ in two
dimensions

We have stated in Eq. (II.5) the general form of forth-order tensor fields
consistent with the assumed symmetries. As shown here, not all indicated terms
are needed for the two-dimensional systems studied numerically in this work.
To see this let us, following the discussion in Sec. II.2, rotate the
coordinate system such that the $1$-axis points into the direction of $\bf q$,
i.e. $q_{\alpha}^{\prime}=q\delta_{1\alpha}$ with the prime “$\prime$” marking
the rotated frame. Using this coordinate system we define the four functions
Klochko et al. (2018)

$\left.\begin{array}[]{ll}c_{\mathrm{L}}(q)&\equiv
T^{\prime}_{1111}(\mbox{$\bf q$})\\\ c_{\mathrm{N}}(q)&\equiv
T^{\prime}_{2222}(\mbox{$\bf q$})\\\ c_{\perp}(q)&\equiv
T^{\prime}_{1122}(\mbox{$\bf q$})\\\ c_{\mathrm{T}}(q)&\equiv
T^{\prime}_{1212}(\mbox{$\bf q$})\end{array}\right\\}\ \mbox{ for }\
q_{\alpha}^{\prime}=q\delta_{1\alpha}.$ (31)

Since the system is isotropic, these functions depend on the wavelength $q$
but not on the direction $\hat{\mbox{$\bf q$}}$ of the wavector $\bf q$. In
other words, they are invariant under rotation and they do not change either
(being only dependent on $q$) if one of the coordinate axes is inverted.
Importantly, all other components
$T^{\prime}_{\alpha\beta\gamma\delta}(\mbox{$\bf q$})$ are either by Eq. (17)
and Eq. (18) identical to these invariants or must vanish for an odd number of
equal indices due to Eq. (13) as discussed in Sec. II.2. The $d^{o}=16$
components $T^{\prime}_{\alpha\beta\gamma\delta}(\mbox{$\bf q$})$ are thus
completely determined by the four invariants and this for any $\bf q$. The
tensor field $T_{\alpha\beta\gamma\delta}(\mbox{$\bf q$})$ in the original
frame may then be obtained by the inverse rotation of
$T^{\prime}_{\alpha\beta\gamma\delta}(\mbox{$\bf q$})$ to the original
unrotated frame using Eq. (8). Let us define the coefficients $i_{n}(q)$ using

$\displaystyle c_{\mathrm{L}}(q)$ $\displaystyle=$ $\displaystyle
i_{1}(q)+2i_{2}(q)+2i_{3}(q)+i_{4}(q)$ $\displaystyle c_{\mathrm{N}}(q)$
$\displaystyle=$ $\displaystyle i_{1}(q)+2i_{2}(q)$ $\displaystyle
c_{\perp}(q)$ $\displaystyle=$ $\displaystyle i_{1}(q)+i_{3}(q)$
$\displaystyle c_{\mathrm{T}}(q)$ $\displaystyle=$ $\displaystyle i_{2}(q)$
(32)

which is equivalent to the inverse relations

$\displaystyle i_{1}(q)$ $\displaystyle=$ $\displaystyle
c_{\mathrm{N}}(q)-2c_{\mathrm{T}}(q)$ (33) $\displaystyle i_{2}(q)$
$\displaystyle=$ $\displaystyle c_{\mathrm{T}}(q)$ $\displaystyle i_{3}(q)$
$\displaystyle=$ $\displaystyle
c_{\perp}(q)-c_{\mathrm{N}}(q)+2c_{\mathrm{T}}(q)$ $\displaystyle i_{4}(q)$
$\displaystyle=$ $\displaystyle
c_{\mathrm{L}}(q)+c_{\mathrm{N}}(q)-2c_{\perp}(q)-4c_{\mathrm{T}}(q).$

Consistently with Ref. Klochko et al. (2018) it is then seen that

$\displaystyle T_{\alpha\beta\gamma\delta}(\mbox{$\bf q$})$ $\displaystyle=$
$\displaystyle i_{1}(q)\ \delta_{\alpha\beta}\delta_{\gamma\delta}$
$\displaystyle+$ $\displaystyle
i_{2}(q)\left(\delta_{\alpha\gamma}\delta_{\beta\delta}+\delta_{\alpha\delta}\delta_{\beta\gamma}\right)$
$\displaystyle+$ $\displaystyle
i_{3}(q)\left({\hat{q}}_{\alpha}{\hat{q}}_{\beta}\delta_{\gamma\delta}+{\hat{q}}_{\gamma}{\hat{q}}_{\delta}\delta_{\alpha\beta}\right)$
$\displaystyle+$ $\displaystyle i_{4}(q)\
{\hat{q}}_{\alpha}{\hat{q}}_{\beta}{\hat{q}}_{\gamma}{\hat{q}}_{\delta}$

which agrees with Eq. (II.5) if we set $i_{5}(q)\equiv 0$.

###  II.7 Response to point sources

Isotropic tensor fields may also be constructed by taking the (inner or outer)
product of a tensor field and a (constant) tensor. We focus here on the
second-order tensor field

$R_{\alpha\beta}(\mbox{$\bf q$})=\frac{1}{V}\
C_{\alpha\beta\gamma\delta}(\mbox{$\bf q$})\ s_{\gamma\delta}$ (35)

obtained from a forth-order tensor field
$C_{\alpha\beta\gamma\delta}(\mbox{$\bf q$})$ and a second-order tensor
$s_{\alpha\beta}$ and where (as always) summation over repeated indices is
implied. For later convenience we have introduced the volume $V$ of the
system. The results presented below are readily generalized for different
types of products of tensor fields and tensors and for dimensions other than
$d=2$.

If both $C_{\alpha\beta\gamma\delta}(\mbox{$\bf q$})$ and $s_{\alpha\beta}$
are isotropic, the product theorem Eq. (14) discussed in Sec. II.2 implies
that $R_{\alpha\beta}(\mbox{$\bf q$})$ must also be an isotropic tensor field.
Under the additional assumptions stated in Sec. II.3
$R_{\alpha\beta}(\mbox{$\bf q$})$ is then given by Eq. (24) in terms of two
invariants $k_{1}(q)$ and $k_{2}(q)$. These invariants can in turn be
expressed in terms of the invariants of
$C_{\alpha\beta\gamma\delta}(\mbox{$\bf q$})$ and $s_{\alpha\beta}$. It is
important to emphasize that albeit being closely related
$R_{\alpha\beta}(\mbox{$\bf q$})$ and $C_{\alpha\beta\gamma\delta}(\mbox{$\bf
q$})$ have in general different angular dependences. This may be readily seen
by focusing on the specific cases $R_{12}(\mbox{$\bf q$})$ and
$C_{1212}(\mbox{$\bf q$})$. It follows from Eq. (24) that $R_{12}(\mbox{$\bf
q$})\propto{\hat{q}}_{1}{\hat{q}}_{2}\propto\sin(2\theta)$ with
${\hat{q}}_{1}=\cos(\theta)$ and ${\hat{q}}_{2}=\sin(\theta)$ and from Eq.
(II.6) that $C_{1212}(\mbox{$\bf
q$})=i_{2}(q)+i_{4}(q){\hat{q}}_{1}^{2}{\hat{q}}_{2}^{2}$ is given by an
angular-independent scalar plus a term proportional to $\cos(4\theta)$. In
other words, $R_{12}(\mbox{$\bf q$})$ is a quadrupolar field whereas
$C_{1212}(\mbox{$\bf q$})$ is octupolar foo (c). Both fields thus reveal
distinct angular patterns.

For reasons which will become obvious below we shall call
$R_{\alpha\beta}(\mbox{$\bf q$})$ the “response”,
$C_{\alpha\beta\gamma\delta}(\mbox{$\bf q$})$ the “correlation function” or
“propagator” and $s_{\gamma\delta}$ the “source” or “perturbation”. Up to now
we have not used that the field vector $\bf q$ refers to the wavevector
characterizing fields in reciprocal space. Using Eq. (91) it is seen that the
tensor $s_{\alpha\beta}/V$ in reciprocal space corresponds to a “point source”
$S_{\alpha\beta}(\mbox{$\bf r$})=s_{\alpha\beta}\delta(\mbox{$\bf r$})$ in
real space where we have used Dirac’s delta function. The response
$R_{\alpha\beta}(\mbox{$\bf r$})={\cal F}^{-1}[R_{\alpha\beta}(\mbox{$\bf
q$})]$ in real space is then given by

$R_{\alpha\beta}(\mbox{$\bf r$})=C_{\alpha\beta\gamma\delta}(\mbox{$\bf
r$})s_{\gamma\delta}$ (36)

using the correlation function $C_{\alpha\beta\gamma\delta}(\mbox{$\bf
r$})={\cal F}^{-1}[C_{\alpha\beta\gamma\delta}(\mbox{$\bf q$})]$ in real
space. For a more general source term $S_{\alpha\beta}(\mbox{$\bf r$})$ we
have of course a convolution relation

$R_{\alpha\beta}(\mbox{$\bf r$})=\frac{1}{V}\int\mathrm{d}\mbox{$\bf
r$}^{\prime}C_{\alpha\beta\gamma\delta}(\mbox{$\bf r$}-\mbox{$\bf
r$}^{\prime})\ S_{\gamma\delta}(\mbox{$\bf r$}^{\prime})$ (37)

which reduces to Eq. (36) for a point source. Importantly, all statements made
above for the reciprocal space remain valid in real space, i.e. that
$R_{\alpha\beta}(\mbox{$\bf r$})$ is an isotropic tensor field if
$C_{\alpha\beta\gamma\delta}(\mbox{$\bf r$})$ and $S_{\alpha\beta}(\mbox{$\bf
r$})$ are isotropic and, more specifically, that $R_{12}(\mbox{$\bf r$})$ is a
quadrupolar field while $C_{1212}(\mbox{$\bf r$})$ is octupolar.

Importantly, in many physical situations the source is in fact not isotropic
and thus in turn the response field not consistent with Eq. (24). We remind
that, e.g., the mechanical response of amorphous solids under loading proceeds
from local and irreversible rearrangements, resetting disorder locally thus
generating a highly non-trivial mechanical noise (“shear transformation
zones”) Argon and Kuo (1979); Picard et al. (2004); Rodney et al. (2011);
Nicolas et al. (2018); Falk and Langer (1998). According to a popular model of
localized plastic failure two orthogonal twin force dipoles of opposite signs
may be imposed at the origin Picard et al. (2004). This suggests to relax the
isotropy condition for $s_{\alpha\beta}$. Since the source tensor is still
symmetric it may be diagonalized by an appropriate rotation of the coordinate
system where $s_{12}=s_{21}=0$ and $s_{11}$ and $s_{22}$ become the two (in
general not identical) eigenvalues. Hence,

$R_{\alpha\beta}(\mbox{$\bf q$})=s_{11}C_{\alpha\beta 11}(\mbox{$\bf
q$})+s_{22}C_{\alpha\beta 22}(\mbox{$\bf q$})$

with the isotropic correlation tensor field still being given by Eq. (II.6).
Specifically, this implies

$R_{12}(\mbox{$\bf
q$})={\hat{q}}_{1}{\hat{q}}_{2}\left[i_{3}(q)(s_{11}+s_{22})+i_{4}(q)\underline{(s_{11}{\hat{q}}_{1}^{2}+s_{22}{\hat{q}}_{2}^{2})}\right].$

For $s_{11}=s_{22}$ the underlined term becomes a constant and we recover the
isotropic and quadrupolar response field discussed above. Interestingly, for
eigenvalues of opposite sign, $s_{11}=-s_{22}$, we obtain

$R_{12}(\mbox{$\bf q$})=s_{11}i_{4}(q)\
{\hat{q}}_{1}{\hat{q}}_{2}({\hat{q}}_{1}^{2}-{\hat{q}}_{2}^{2})\propto\sin(4\theta).$
(38)

The (non-isotropic) response field $R_{12}(\mbox{$\bf q$})$ thus has in this
case the same multipole pattern as the (isotropic) correlation field
$C_{1212}(\mbox{$\bf q$})$ albeit shifted by an angle $\pi/8$ foo (c). It is
readily seen by inverse FT that the same behavior applies in real space.

In summary, two different types of angular dependence of a response field must
be distinguished. If the angular dependence is consistent with Eq. (24), this
behavior should not be called “anisotropy” since the angle dependence is
basically a trivial consequence of the fact that tensor field components are
measured. If on the other hand $R_{\alpha\beta}(\mbox{$\bf q$})$ is not
consistent with Eq. (24) this suggests that either the correlation tensor
field $C_{\alpha\beta\gamma\delta}(\mbox{$\bf q$})$ and/or the source tensor
$s_{\alpha\beta}$ are not isotropic. In many physical situations this is in
fact expected for the source term while the correlation fields may be assumed
to be isotropic. The physical behavior of response and correlation fields,
albeit closely related, then differ and should thus not be lumped together.

##  III Computational issues

#### Algorithm and systems.

We present below numerical data obtained for amorphous glasses formed by
polydisperse Lennard-Jones (pLJ) particles George et al. (2021a); Klochko et
al. (2022) simulated by means of Monte Carlo (MC) simulations Allen and
Tildesley (2017). See Appendix C for details (Hamiltonian, units,
configuration preparation, data averaging). We focus on systems containing
$n=10000$ particles at a working temperature $T=0.2$ much lower than the glass
transition temperature $T_{\mathrm{g}}\approx 0.26$, i.e. for our largest
sampling time $\Delta\tau=10^{7}$ the systems behave as solid elastic bodies
and all stochastic processes are stationary George et al. (2021a).
Importantly, all $N_{\mathrm{c}}=200$ completely independent configurations
$c$ are quenched and tempered by means of a mix of local and swap MC hopping
moves Ninarello et al. (2017); George et al. (2021a) (being thus effectively
kept adiabatically at thermal equilibrium) while the data production runs are
sampled only using local MC moves. For each $c$ we store several time-series
containing each $N_{\mathrm{t}}=10000$ frames $t$.

Figure 2: Two-dimensional ($d=2$) square lattice with $a_{\mathrm{grid}}$
being the lattice constant and $n_{L}=L/a_{\mathrm{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
$\bf r$ of a microcell is either given by the $r_{1}$\- and
$r_{2}$-coordinates (in the principal box) or by the distance $r=|\mbox{$\bf
r$}|$ from the origin (large circle) and the angle $\theta$.

#### Data sampling and analysis.

In a first step various instantaneous properties are computed for each $t$
which are then “$t$-averaged” over the correlated $t$ and finally
“$c$-averaged” over the independent $c$. We thus characterize, e.g., the
elastic modulus tensor $E_{\alpha\beta\gamma\delta}$ by means of the stress-
fluctuation formalism Lutsko (1988, 1989); Wittmer et al. (2013a, b, 2015a,
2015b) from which a finite Young modulus $E\approx 45$ is obtained (cf. Sec.
V.3). Similarly, we compute in turn

  * •

for each $c$ and $t$ the stress tensor field $\sigma_{\alpha\beta}(\mbox{$\bf
r$},t)|_{c}$ (cf. Appendix D) using a regular square grid as shown in Fig. 2
with a lattice constant $a_{\mathrm{grid}}\approx 0.2$,

  * •

the $t$-averaged fields $\bar{\sigma}_{\alpha\beta}(\mbox{$\bf r$})|_{c}$,

  * •

by Fast-Fourier transform the corresponding stress fields
$\bar{\sigma}_{\alpha\beta}(\mbox{$\bf q$})|_{c}={\cal
F}[\bar{\sigma}_{\alpha\beta}(\mbox{$\bf r$})|_{c}]$ in reciprocal space,

  * •

the correlation functions $\bar{c}_{\alpha\beta\gamma\delta}(\mbox{$\bf
q$})|_{c}$ in reciprocal space for each configuration $c$, cf. Eq. (126),

  * •

the $c$-average $\bar{c}_{\alpha\beta\gamma\delta}(\mbox{$\bf q$})$, cf. Eq.
(127),

  * •

and finally by inverse FT $\bar{c}_{\alpha\beta\gamma\delta}(\mbox{$\bf
r$})={\cal F}^{-1}[\bar{c}_{\alpha\beta\gamma\delta}(\mbox{$\bf q$})]$ the
correlation functions in real space.

To obtain the correlation functions in rotated coordinates we rotate first
$\bar{\sigma}_{\alpha\beta}(\mbox{$\bf
q$})|_{c}\to\bar{\sigma}_{\alpha\beta}^{\prime}(\mbox{$\bf q$})|_{c}$ and
perform then all the subsequent steps as before. We note finally that real and
reciprocal space correlation functions have the same dimension due to our FT
convention (cf. Appendix A).

##  IV Main numerical results

Figure 3: ICFs $\bar{c}_{\mathrm{N}}(q)$, $\bar{c}_{\mathrm{L}}(q)$,
$\bar{c}_{\perp}(q)$ and $\bar{c}_{\mathrm{T}}(q)$ for a large sampling time
$\Delta\tau=10^{7}$ and $n=10000$ particles. Horizontal solid lines indicate
the phenomenological constant $e\approx 45$. Inset:
$\bar{c}_{2222}^{\circ}(q,\theta)$ does not depend on $\theta$ while
$\bar{c}_{1212}(q,\theta)$ is well described by Eq. (45) (dashed line). Main
panel: $\bar{c}_{\mathrm{L}}(q)$, $\bar{c}_{\perp}(q)$ and
$\bar{c}_{\mathrm{T}}(q)$ vanish while $\beta V\bar{c}_{\mathrm{N}}(q)\simeq
e>0$ for $q\ll 1$. A thin solid line emphasizes the peak of
$\bar{c}_{\mathrm{N}}(q)$ at $q\approx 6.5$.

#### Determination of ICFs in NRC.

Before we shall have a closer look at
$\bar{c}_{\alpha\beta\gamma\delta}(\mbox{$\bf q$})$ in standard unrotated or
rotated coordinates let us first characterize the correlations in NRC, i.e.
for each wavevector $\bf q$ the coordinate system is rotated until the
$1$-axis coincides with the $\bf q$-direction. We mark these new tensor field
components by “$\circ$” to distinguish them from standard rotated tensor field
components (marked by primes “$\prime$”) using the same rotation for all $\bf
q$. Note that $q_{\alpha}^{\circ}=q\delta_{1\alpha}$. The ICFs
$\bar{c}_{\alpha\beta\gamma\delta}^{\circ}(\mbox{$\bf
q$})=\langle\bar{\sigma}^{\circ}_{\alpha\beta}(\mbox{$\bf
q$})\bar{\sigma}^{\circ}_{\gamma\delta}(-\mbox{$\bf q$})\rangle$ are thus
obtained using the “Invariant Stress Fields” (ISFs)
$\bar{\sigma}^{\circ}_{\alpha\beta}(\mbox{$\bf q$})|_{c}$ rotated differently
for each $\bf q$. Importantly, for strictly isotropic systems
$\bar{c}_{\alpha\beta\gamma\delta}^{\circ}(\mbox{$\bf q$})$ only depends on
the magnitude $q$ of $\bf q$ but not on its direction ${\hat{\bf q}}$.
Consistently with Eq. (31) and following Ref. Klochko et al. (2018), we define

$\displaystyle\bar{c}_{\mathrm{L}}(q)$ $\displaystyle\equiv$
$\displaystyle\left<\bar{c}^{\circ}_{1111}(\mbox{$\bf q$})\right>_{{\hat{\bf
q}}},$ (39) $\displaystyle\bar{c}_{\mathrm{N}}(q)$ $\displaystyle\equiv$
$\displaystyle\left<\bar{c}^{\circ}_{2222}(\mbox{$\bf q$})\right>_{{\hat{\bf
q}}},$ $\displaystyle\bar{c}_{\perp}(q)$ $\displaystyle\equiv$
$\displaystyle\left<\bar{c}^{\circ}_{1122}(\mbox{$\bf q$})\right>_{{\hat{\bf
q}}}\mbox{ and }$ $\displaystyle\bar{c}_{\mathrm{T}}(q)$ $\displaystyle\equiv$
$\displaystyle\left<\bar{c}^{\circ}_{1212}(\mbox{$\bf q$})\right>_{{\hat{\bf
q}}}$

where we average over all ${\hat{\bf q}}$ with $|\mbox{$\bf q$}|\approx q$
(using a bin width similar to the lattice spacing of the grid in reciprocal
space). The ${\hat{\bf q}}$-averaged four ICFs are shown in the main panel of
Fig. 3 for our largest sampling time $\Delta\tau$. The central observation is
that $\bar{c}_{\mathrm{L}}(q)$, $\bar{c}_{\perp}(q)$ and
$\bar{c}_{\mathrm{T}}(q)$ vanish for sufficiently large $\Delta\tau$ while
$\bar{c}_{\mathrm{N}}(q)$ remains finite. Moreover, as emphasized by the bold
solid line in the main panel

$\beta V\bar{c}_{\mathrm{N}}(q)\simeq e\approx 45\mbox{ for }q\ll 1$ (40)

with $\beta=1/T$ being the inverse temperature, $V$ the system volume and $e$
a phenomenological constant characterizing the plateau in the hydrodynamic
limit. As shown in the inset of Fig. 3 for $\bar{c}^{\circ}_{2222}(q,\theta)$
and $\alpha=0^{\circ}$, we have checked for several tensor field components
$\bar{c}^{\circ}_{\alpha\beta\gamma\delta}(q,\theta)$ that the expected
$\theta$-independence for isotropic systems holds (within statistical
accuracy).

Figure 4: $\Delta\tau$-dependence in reciprocal space using NRC and $q\approx
0.06$. Main panel: $\Delta\tau$-dependence of ICFs demonstrating that
$\bar{c}_{\mathrm{N}}(q,\Delta\tau)$ remains finite while all other ICFs decay
inversely with $\Delta\tau$. Inset: $\bar{\sigma}^{\circ}_{11}(\mbox{$\bf
q$},\Delta\tau)$, $\bar{\sigma}^{\circ}_{12}(\mbox{$\bf q$},\Delta\tau)$ and
$\bar{\sigma}^{\circ}_{22}(\mbox{$\bf q$},\Delta\tau)$ in the complex plane
(horizontal axis for real part, vertical axis for imaginary part) as a
function of $\Delta\tau$ (using a logarithmic scale for the data points) for
one $\bf q$ and one $c$. $\bar{\sigma}^{\circ}_{11}(\mbox{$\bf
q$},\Delta\tau)$ and $\bar{\sigma}^{\circ}_{12}(\mbox{$\bf q$},\Delta\tau)$
vanish for large $\Delta\tau$, i.e. converge to the origin of the complex
plane, while $\bar{\sigma}^{\circ}_{22}(\mbox{$\bf q$},\Delta\tau)$ has a
finite attractor $\bar{\sigma}^{\circ\mathrm{q}}_{22}(\mbox{$\bf q$})$.

#### Sampling time dependence and symmetry breaking.

The dependence of the ICFs on $\Delta\tau$ is summarized in Fig. 4. Let us
first note that the ($t$-averaged) force $\bar{g}_{\alpha}(\mbox{$\bf
r$})=\partial_{\beta}\bar{\sigma}_{\alpha\beta}(\mbox{$\bf r$})$ acting on
each material element becomes in reciprocal space

$\bar{g}^{\circ}_{\alpha}(\mbox{$\bf q$})=iq^{\circ}_{\beta}\
\bar{\sigma}^{\circ}_{\alpha\beta}(\mbox{$\bf q$})$ (41)

using NRC. Since $q^{\circ}_{1}=q$ is finite (for $\mbox{$\bf q$}\neq{\bf
0}$), finite $\bar{\sigma}^{\circ}_{11}(\mbox{$\bf q$})$ and
$\bar{\sigma}^{\circ}_{12}(\mbox{$\bf q$})$ correspond to finite forces
$\bar{g}^{\circ}_{\alpha}(\mbox{$\bf q$})$ which in turn generate fluxes.
Finite $\bar{\sigma}^{\circ}_{11}(\mbox{$\bf q$},\Delta\tau)$ and
$\bar{\sigma}^{\circ}_{12}(\mbox{$\bf q$},\Delta\tau)$ must thus rapidly
vanish, as seen from the corresponding trajectories in the inset of Fig. 4,
and therefore $\bar{c}_{\mathrm{L}}(q,\Delta\tau)$,
$\bar{c}_{\perp}(q,\Delta\tau)$ and $\bar{c}_{\mathrm{T}}(q,\Delta\tau)$ also
vanish for large $\Delta\tau$ (main panel) Klochko et al. (2022). Since
$q^{\circ}_{2}=0$ for all $\bf q$ the normal stress
$\bar{\sigma}^{\circ}_{22}(\mbox{$\bf q$},\Delta\tau)$ may be finite without
violating static mechanical equilibrium. As shown in the inset, the
$\bar{\sigma}^{\circ}_{22}(\mbox{$\bf q$},\Delta\tau)$ thus have in general
finite attractors $\bar{\sigma}^{\circ\mathrm{q}}_{22}(\mbox{$\bf q$})$. As
emphasized by the superscript “q”, these are for realistic $\Delta\tau$
essentially quenched stresses. (Only for $\Delta\tau$ of order of the
$\alpha$-relaxation time $\tau_{\alpha}$ these attractors become weakly time-
dependent diffusively decaying for symmetry reasons towards the origin of the
complex plane. Note that $\Delta\tau\ll\tau_{\alpha}$ for $T=0.2$.) Since the
large-$\Delta\tau$ limit of $\bar{c}_{\mathrm{N}}(q,\Delta\tau)$ is the
typical squared magnitude of the
$\bar{\sigma}^{\circ\mathrm{q}}_{22}(\mbox{$\bf q$})$ in the complex plane,
this implies

$e=\beta V\left<\bar{\sigma}^{\circ\mathrm{q}}_{22}(\mbox{$\bf
q$})\bar{\sigma}^{\circ\mathrm{q}}_{22}(-\mbox{$\bf q$})\right>.$ (42)

This means that $e$ is a finite $\Delta\tau$-independent static property
characterizing the typical size of the (continuous) symmetry breaking
associated with the stress components normal to the wavevectors for each
configuration $c$. Note that $e$ does thus not depend on whether we use, e.g.,
a momentum conserving or an overdamped simulation scheme Allen and Tildesley
(2017). In general, it is a fitting parameter depending on the distribution of
the $\bar{\sigma}^{\circ\mathrm{q}}_{22}(\mbox{$\bf q$})$ caused by the
preparation history. Interestingly, it is observed that $e$ is similar to the
Young modulus $E$. This result is in fact expected from recent studies on
equilibrium viscoelastic fluids (including supercooled liquids and glasses)
Klochko et al. (2018, 2022) showing that the ICFs may be expressed in the
small-$q$ limit in terms of invariant macroscopic relaxation functions
Rubinstein and Colby (2003); Klochko et al. (2018, 2022); Maier et al. (2017,
2018); Vogel et al. (2019). Naturally, this requires additional assumptions
the crucial point being here that the systems must be at thermal equilibrium.
See Sec. V for more details.

#### Reciprocal space correlation functions.

We turn now to a coordinate system rotated as in Fig. 1(b) by the same angle
$\alpha$ for all $\bf q$. As shown in Sec. II.6 and using the form-invariance
of isotropic tensor fields, cf. Eq. (12), the four ICFs determine in $d=2$ the
isotropic forth-order tensor field Klochko et al. (2018)

$\displaystyle\bar{c}^{\prime}_{\alpha\beta\gamma\delta}(\mbox{$\bf q$})$
$\displaystyle=$
$\displaystyle\left[\bar{c}_{\mathrm{N}}-2\bar{c}_{\mathrm{T}}\right]\ \
\delta_{\alpha\beta}\delta_{\gamma\delta}$ $\displaystyle+$
$\displaystyle\bar{c}_{\mathrm{T}}\
\left(\delta_{\alpha\gamma}\delta_{\beta\delta}+\delta_{\alpha\delta}\delta_{\beta\gamma}\right)$
$\displaystyle+$
$\displaystyle\left[\bar{c}_{\perp}-\bar{c}_{\mathrm{N}}+2\bar{c}_{\mathrm{T}}\right]\
\left(\hat{q}^{\prime}_{\alpha}\hat{q}^{\prime}_{\beta}\delta_{\gamma\delta}+\hat{q}^{\prime}_{\gamma}\hat{q}^{\prime}_{\delta}\delta_{\alpha\beta}\right)$
$\displaystyle+$
$\displaystyle\left[\bar{c}_{\mathrm{L}}+\bar{c}_{\mathrm{N}}-2\bar{c}_{\perp}-4\bar{c}_{\mathrm{T}}\right]\
\
\hat{q}^{\prime}_{\alpha}\hat{q}^{\prime}_{\beta}\hat{q}^{\prime}_{\gamma}\hat{q}^{\prime}_{\delta}$

with ${\hat{q}}^{\prime}_{\alpha}$ being the $\alpha$-component of ${\hat{\bf
q}}^{\prime}={\hat{\bf q}}$ in rotated coordinates. The $\alpha$-rotation
changes the terms in the last two lines of Eq. (IV). (Minor generalizations
are needed for $d>2$.) Using the known values of the ICFs of our system Eq.
(IV) reduces to

$\displaystyle\beta V\bar{c}^{\prime}_{\alpha\beta\gamma\delta}(\mbox{$\bf
q$})$ $\displaystyle\simeq$ $\displaystyle e\times$ (44)
$\displaystyle(\delta_{\alpha\beta}\delta_{\gamma\delta}$ $\displaystyle-$
$\displaystyle{\hat{q}}^{\prime}_{\alpha}{\hat{q}}^{\prime}_{\beta}\delta_{\gamma\delta}-{\hat{q}}^{\prime}_{\gamma}{\hat{q}}^{\prime}_{\delta}\delta_{\alpha\beta}+{\hat{q}}^{\prime}_{\alpha}{\hat{q}}^{\prime}_{\beta}{\hat{q}}^{\prime}_{\gamma}{\hat{q}}^{\prime}_{\delta})$

for $q\ll 1$ and large $\Delta\tau$. According to Eq. (44) we thus obtain,
e.g., for the shear-stress autocorrelation function

$8\beta V\bar{c}^{\prime}_{1212}(\mbox{$\bf q$})=e\
(1-\cos[4(\theta-\alpha)])$ (45)

for all $\alpha$. As seen from the inset of Fig. 3 for $\alpha=0^{\circ}$,
this prediction (bold dashed line) agrees perfectly with our data (triangles).

Figure 5: Rescaled shear-stress autocorrelation function
$\beta\bar{c}^{\prime}_{1212}(r,\theta)r^{2}$ as a function of
$x=\theta-\alpha$ comparing data for different $r$ and $\alpha$ with the
prediction (bold solid line).

#### Real space correlation functions.

We return now to the correlations
$\bar{c}^{\prime}_{\alpha\beta\gamma\delta}(\mbox{$\bf r$})={\cal
F}^{-1}[\bar{c}^{\prime}_{\alpha\beta\gamma\delta}(\mbox{$\bf q$})]$ in real
space. As already stated in the Introduction, cf. Eq. (1), inverse Fourier
transformation (cf. Appendix B) implies

$-\beta\bar{c}^{\prime}_{1212}(\mbox{$\bf r$})\simeq\frac{e}{4\pi
r^{2}}\cos[4(\theta-\alpha)]\mbox{ for }r\gg 1$ (46)

with $\theta$ being the angle indicated in Fig. 1(a). The same
large-$\Delta\tau$ limit holds for $-\beta\bar{c}^{\prime}_{1122}(\mbox{$\bf
r$})$ and for $\beta(\bar{c}^{\prime}_{1111}(\mbox{$\bf
r$})+\bar{c}^{\prime}_{2222}(\mbox{$\bf r$}))/2$. Moreover,

$\beta(\bar{c}^{\prime}_{1111}(\mbox{$\bf
r$})-\bar{c}^{\prime}_{2222}(\mbox{$\bf r$}))/2\simeq 2\frac{e}{4\pi
r^{2}}\cos[2(\theta-\alpha)]$ (47)

for $r\gg 1$ and for large $\Delta\tau$. The angle dependence for the shear-
stress autocorrelation function in real space is investigated in Fig. 5 where
we plot using linear coordinates $\beta\bar{c}^{\prime}_{1212}(r,\theta)r^{2}$
as a function of $x=\theta-\alpha$ for different $r$ and $\alpha$. The data
compare well with the prediction, Eq. (46). Naturally, the statistics
deteriorates with increasing $r$.

Figure 6: $P[f(\mbox{$\bf r$}),p](r)$ for rescaled correlation functions
$f(\mbox{$\bf r$})$ and mode number $p$ as indicated ($\alpha=0^{\circ}$). A
double-logarithmic representation is used. Only data for $1<r<L/2$ are given.
The bold solid lines represent $e/4\pi r^{2}$. Main panel: Various projections
for $n=10000$. Inset: Projection for $f(\mbox{$\bf
r$})=-\beta\bar{c}_{1212}(\mbox{$\bf r$})$ and $p=4$ for different particle
numbers $n$.

#### Asymptotic $r$-dependence.

A more precise check of the $r$-dependence is obtained using the
$\theta$-average

$P[f,p](r)\equiv\frac{2}{\pi}\int_{0}^{\pi}\mathrm{d}\theta\
f(r,\theta)\cos(p\theta)$ (48)

for $p=2$ and $p=4$. For convenience the prefactor of the integral is chosen
such that $P[\cos(2\theta),2]=P[\cos(4\theta),4]=1$. (On the discrete grid the
integral is replaced by the sum over all grid points in a distance interval
$[r-\delta r/2,r+\delta r/2]$ which is finally normalized by half the number
of grid points.) All correlation functions $f(\mbox{$\bf r$})$ presented in
Fig. 6 are rescaled to make their projections collapse on the same power law
$e/4\pi r^{2}$ (bold solid lines). The main panel presents different
correlation functions for $n=10000$. The negative signs for $\bar{c}_{1212}$
and $\bar{c}_{1122}$ are implied by Eq. (46). Consistently with Eq. (47) the
projection for $(\bar{c}_{1111}-\bar{c}_{2222})/2$ is additionally rescaled
with a prefactor $\beta/2$. Focusing on $f(\mbox{$\bf
r$})=-\beta\bar{c}_{1212}(\mbox{$\bf r$})$ and $p=4$ we verify in the inset
the system-size independence of these results. Data for a broad range of
particle numbers $n$ are presented. As can be seen, all data nicely collapse
on $e/4\pi r^{2}$, confirming thus the predicted long-range correlations for
asymptotically large simulation boxes. Similar results (not shown) have been
found for the projections of other
$\bar{c}_{\alpha\beta\gamma\delta}(\mbox{$\bf r$})$.

##  V Why and when $e\approx E$ holds

###  V.1 Introduction

The phenomenological parameter $e$ was defined in Sec. IV by the limit Eq.
(40). Note that $e$ has the same dimension energy/volume as the stress or an
elastic modulus. We have verified (cf. Fig. 4) that $e$ indeed becomes
$\Delta\tau$-independent for sufficiently large sampling times $\Delta\tau$
for our pLJ particle glasses. Having fitted the value $e\approx 45$ and using
that the ICFs $\bar{c}_{\mathrm{L}}(q)$, $\bar{c}_{\perp}(q)$ and
$\bar{c}_{\mathrm{T}}(q)$ vanish, all numerical results for sufficiently large
$\Delta\tau$, small $q$ or large $r$ can be explained without any additional
physical insight. We have observed, however, that $e$ is similar to the
(equilibrium static) Young modulus $E$. Albeit not strictly necessary for the
main thrust of this work, a demonstration that both constants should be
similar or even equal — under to be specified assumptions and approximations —
must be an important finding allowing to estimate apriori the angular
dependence of the correlation functions. This can indeed be done following
Refs. Klochko et al. (2018, 2022) and in agreement with a different and
complementary approach developed in Refs. Maier et al. (2017, 2018); Vogel et
al. (2019). Naturally, this requires additional physical input. We remind
first in Sec. V.2 how the correlation functions $\bar{c}(\Delta\tau)$ of time-
averaged fields, the focus of the present work, are related to the correlation
functions $c(\tau)$ of instantaneous fields under the assumption that the
relevant stochastic processes are stationary foo (b). A short recap of linear
viscoelasticity is given in Sec. V.3. Using the general theoretical
predictions for the ICFs of instantaneous stress fields Klochko et al. (2018,
2022), reminded in Sec. V.4, it is shown (cf. Sec. V.5) that $e\approx E$ for
equilibrated viscoelastic fluids with a sufficiently broad elastic plateau.

###  V.2 Instantaneous and time-averaged fields

#### General connection for stationary processes.

We discuss in this work spatial correlation functions $\bar{c}$ (both in real
and reciprocal space) of time-averaged stress fields computed over
$N_{\mathrm{t}}=\Delta\tau/\delta\tau$ instantaneous configurations. The
various correlation functions thus depend in principle on the sampling time
$\Delta\tau$ as we have seen in Fig. 4. As shown elsewhere George et al.
(2021a) assuming a stationary stochastic process (both for global properties
as for fields) the $\Delta\tau$-dependence of $\bar{c}(\Delta\tau)$ can be
traced back via

$\bar{c}(\Delta\tau)=\frac{2}{\Delta\tau^{2}}\int_{0}^{\Delta\tau}\mathrm{d}\tau\
(\Delta\tau-\tau)\ c(\tau)$ (49)

to the time dependence of the corresponding correlation function $c(\tau)$ of
the instantaneous fields. Please note that Eq. (49) is closely related to the
equivalence of the Einstein relation, corresponding to $\bar{c}(\Delta\tau)$,
and the Green-Kubo relation, corresponding to $c(\tau)$, for transport
coefficients Hansen and McDonald (2006); Allen and Tildesley (2017); George et
al. (2021a). We thus study in this work correlations within the Einstein
picture. This has the advantage that the integral Eq. (49) filters irrelevant
high frequencies, i.e. $\bar{c}(\Delta\tau)$ is a natural smoothing function
of the instantaneous field correlation function $c(\tau)$ on which previous
work has focused on Maier et al. (2017, 2018); Vogel et al. (2019); Klochko et
al. (2018); Maloney and Lemaître (2004); Lemaître (2014, 2015, 2017, 2018);
Chowdhury et al. (2016); Klochko et al. (2022); Steffen et al. (2022).

#### Plateau values and asymptotic behavior.

Obviously, Eq. (49) implies that $c(\tau)$ is constant iff
$\bar{c}(\Delta\tau)$ is constant and both constants are equal. Importantly,
this even holds if $c(\tau)$ and $\bar{c}(\Delta\tau)$ are only constant for a
finite but sufficiently large time window George et al. (2021a). Hence, if
$c(\tau)$ becomes constant in the large-time limit $\bar{c}(\Delta\tau)$ as
well becomes constant, i.e.

$\lim_{\tau\to\infty}c(\tau)=\lim_{\Delta\tau\to\infty}\bar{c}(\Delta\tau)\equiv
c_{\infty}$ (50)

with $c_{\infty}$ being the common asymptote foo (d). Thus,
$\bar{c}_{\mathrm{L}}(q,\Delta\tau)\approx\bar{c}_{\perp}(q,\Delta\tau)\approx\bar{c}_{\mathrm{T}}(q,\Delta\tau)\simeq
0$ for large $\Delta\tau$ implies $c_{\mathrm{L}}(q,\tau)\approx
c_{\perp}(q,\tau)\approx c_{\mathrm{T}}(q,\tau)\simeq 0$ for large $\tau$ and
visa versa. Similarily,

$\beta Vc_{\mathrm{N}}(q,\tau)\simeq e\ \ \mbox{\Large$\Leftrightarrow$}\ \
\beta V\bar{c}_{\mathrm{N}}(q,\Delta\tau)\simeq e$ (51)

for the ICFs of the transverse normal stresses in the low-$q$ limit for,
respectively, $\tau\to\infty$ and $\Delta\tau\to\infty$ foo (d).

#### Generalized Maxwell model.

It follows directly from Eq. (49) for $c(\tau)=c_{p}\exp(-\tau/\tau_{p})$ that
George et al. (2021a)

$\displaystyle\bar{c}(\Delta\tau)$ $\displaystyle=$ $\displaystyle
c_{p}D(\Delta\tau/\tau_{p})\mbox{ with }$ $\displaystyle D(x)$
$\displaystyle=$ $\displaystyle 2[\exp(-x)-1+x]/x^{2}$ (52)

being the “Debye function” well known in polymer science Doi and Edwards
(1986); Rubinstein and Colby (2003). For systems with overdamped dynamics,
such as for our MC simulations, one expects the relaxation dynamics to be
described by a linear superposition of exponentially decaying Maxwell modes
Rubinstein and Colby (2003); Ferry (1980). For such a “Generalized Maxwell
model” Eq. (52) generalizes to the superposition George et al. (2021a)

$\bar{c}(\Delta\tau)=c_{\infty}+\sum_{p=1}^{p_{\mathrm{max}}}c_{p}D(\Delta\tau/\tau_{p})$
(53)

with $c_{p}$ being the amplitude and $\tau_{p}$ (with $0<\tau_{p}<\infty$) the
relaxation time of a mode $p$. Note that $c_{\infty}$ corresponds to the modes
with virtually infinite relaxation times. We remind that given a sufficiently
high number of modes $p$ (or a distribution of modes) it is in principle
always possible to fit any reasonable $\bar{c}(\Delta\tau)$ using the standard
numerical techniques Ferry (1980); Provencher (1982); Klochko et al. (2022).
In any case, for $\Delta\tau\gg\tau_{\star}$ with $\tau_{\star}=\tau_{1}$
being the largest (Maxwell) relaxation time Eq. (53) leads to

$\bar{c}(\Delta\tau)\simeq
c_{\infty}+\frac{2}{\Delta\tau}\sum_{p=1}^{p_{\mathrm{max}}}c_{p}\tau_{p},$
(54)

i.e. $\bar{c}(\Delta\tau)-c_{\infty}$ ultimately decays inversely with
$\Delta\tau$. This decay is emphasized by the dashed-dotted line in Fig. 4.

###  V.3 Linear elasticity and viscoelasticity

Figure 7: Comparison of the standard MC algorithm only including local hopping
moves (circles) Wittmer et al. (2013a) with the dramatically more efficient
variant including swap MC hopping moves (squares) Ninarello et al. (2017);
George et al. (2021a) for pLJ systems with $n=10000$ particles focusing on the
shear-stress relaxation function $G(\tau)$ at $T=0.2$. $G(\tau)$ is computed
using the shear-stress autocorrelation function and an additive constant
Wittmer et al. (2015a, b); George et al. (2021a); Klochko et al. (2022).
$G(\tau)\approx\mu\approx 14$ for $\tau\gg 10^{3}$ (solid horizontal line) if
only local moves are used while $G(\tau)\to 0$ (dashed horizontal line) if the
unphysical swap MC moves are added. Including swap MC moves for the quenching
and tempering of the systems thus allows to bring the systems close to
(liquid) equilibrium at $T=0.2$. The final production runs are sampled only
using local moves.

#### Static elastic moduli.

Consistently with Sec. II.4, the elastic modulus tensor
$E_{\alpha\beta\gamma\delta}$ for isotropic systems is completely described by
two invariants, say the two Lamé moduli $\lambda$ and $\mu$. Using Eq. (22)
$E_{\alpha\beta\gamma\delta}$ may thus be written as Landau and Lifshitz
(1959); Tadmor et al. (2012)

$E_{\alpha\beta\gamma\delta}=\lambda\delta_{\alpha\beta}\delta_{\gamma\delta}+\mu\left(\delta_{\alpha\gamma}\delta_{\beta\delta}+\delta_{\alpha\delta}\delta_{\beta\gamma}\right).$
(55)

We have determined $\lambda$ and $\mu$ by means of the stress-fluctuation
formalism described elsewhere Lutsko (1988, 1989); Wittmer et al. (2013a, b,
2015a, 2015b). This shows that for our pLJ particle systems at $T=0.2$ we have
$\lambda\approx 39$ and $\mu\approx 14$. The latter value is indicated by the
bold horizontal line in Fig. 7. (We have verified that the fluctuations of
$\lambda$ and $\mu$ between different independent configurations $c$ are
negligible.) Alternatively, one may describe the elastic response of a body by
means of the creep compliance tensor $J_{\alpha\beta\gamma\delta}$ being the
inverse of the elastic modulus tensor:

$J_{\alpha\beta\gamma\delta}E_{\gamma\delta\alpha^{\prime}\beta^{\prime}}=\frac{1}{2}\left(\delta_{\alpha\alpha^{\prime}}\delta_{\beta\beta^{\prime}}+\delta_{\alpha\beta^{\prime}}\delta_{\alpha^{\prime}\beta}\right).$
(56)

For isotropic bodies Landau and Lifshitz (1959); Tadmor et al. (2012)

$J_{\alpha\beta\gamma\delta}=\frac{1+\nu}{2E}\left(\delta_{\alpha\gamma}\delta_{\beta\delta}+\delta_{\alpha\delta}\delta_{\beta\gamma}\right)-\frac{\nu}{E}\delta_{\alpha\beta}\delta_{\gamma\delta}$
(57)

with $E$ being the Young modulus and $\nu$ Poisson’s ratio. Consistently with
Eq. (56) the two sets of invariants $(\lambda,\mu)$ and $(E,\nu)$ are related
in $d$ dimensions by

$\displaystyle\nu$ $\displaystyle=$
$\displaystyle\frac{\lambda}{2\mu+\lambda(d-1)}$ (58) $\displaystyle E$
$\displaystyle=$ $\displaystyle\lambda+2\mu-(d-1)\lambda\nu.$ (59)

For $d=3$ the usual formulae given in the standard textbooks Landau and
Lifshitz (1959); Tadmor et al. (2012) are recovered while

$\nu=\frac{\lambda}{\lambda+2\mu}\mbox{ and
}E=4\mu\frac{\lambda+\mu}{\lambda+2\mu}\mbox{ for }d=2.$ (60)

Using the known values for $\lambda$ and $\mu$ for our pLJ particle systems at
$T=0.2$ this implies $E\approx 45$ and $\nu\approx 0.6$. Please note that the
values $\lambda_{c}$ and $\mu_{c}$ obtained for each independent configuration
$c$ only differ very weakly for $n>2000$ from the $c$-ensemble averages
$\lambda=\sum_{c}\lambda_{c}/N_{\mathrm{c}}$ and
$\mu=\sum_{c}\mu_{c}/N_{\mathrm{c}}$. We may thus determine $E$ either by
first obtaining for each $c$ a Young modulus $E_{c}$ using $\lambda_{c}$ and
$\mu_{c}$ and from this $E=\sum_{c}E_{c}/N_{\mathrm{c}}$ or by directly
applying Eq. (60) for the $c$-averages $\lambda$ and $\mu$.

#### Stress relaxation tensor for isotropic systems.

The static elastic modulus tensor $E_{\alpha\beta\gamma\delta}$ can be
generalized in the time domain as discussed systematically and in more detail
in the literature Ferry (1980); Doi and Edwards (1986); Rubinstein and Colby
(2003); Balucani and Zoppi (1995); Evans and Morriss (1990). This is important
in general for viscoelastic fluids but also for the dynamical response of
(crystalline or amorphous) solid bodies. The forth-order stress relaxation
tensor $E_{\alpha\beta\gamma\delta}(\tau)$ characterizes the time-dependence
of the stress increment $\delta\sigma_{\alpha\beta}(\tau)$ as a function of an
imposed strain $\epsilon_{\alpha\beta}(\tau)$ foo (b). Generalizing Eq. (55)
for isotropic systems we may write

$\displaystyle E_{\alpha\beta\gamma\delta}(\tau)$ $\displaystyle=$
$\displaystyle M(\tau)\ \delta_{\alpha\beta}\delta_{\gamma\delta}$ (61)
$\displaystyle+$ $\displaystyle G(\tau)\
\left(\delta_{\alpha\gamma}\delta_{\beta\delta}+\delta_{\alpha\delta}\delta_{\beta\gamma}\right)$

with $M(\tau)$ being the “mixed relaxation modulus” and $G(\tau)$ the “shear-
stress relaxation modulus”. The macroscopic linear response tensor is thus
fully determined by (again) two invariant response functions. It is useful to
additionally introduce the “longitudinal relaxation modulus”
$L(\tau)=M(\tau)+2G(\tau)$. Note that $G(\tau)=E_{1212}(\tau)$,
$M(\tau)=E_{1122}(\tau)$ and $L(\tau)=E_{1111}(\tau)=E_{2222}(\tau)$. The
relaxation modulus $E_{\alpha\beta\gamma\delta}(\tau)$ reduces (continuously)
to the static modulus $E_{\alpha\beta\gamma\delta}$ for large times, i.e. the
invariant relaxation functions become

$L(\tau)\to\lambda+2\mu,G(\tau)\to\mu,M(\tau)\to\lambda\mbox{ for
}\tau\to\infty$ (62)

with $\lambda$ and $\mu$ being the known Lamé moduli. That this is the case
can be seen for $G(\tau)$ in Fig. 7 for two different MC variants. As expected
for equilibrium liquids the shear modulus vanishes, $\mu=0$, for the variant
with swap MC hopping moves (squares) used for the equilibration of our
systems. The total ensemble of independent configurations $c$ thus corresponds
to a liquid system while each configuration $c$ being confined in a metabasin
if only local moves are included is an (isotropic) elastic body. We note for
later convenience that, hence,

$L(\tau)-\frac{M(\tau)^{2}}{L(\tau)}\to\lambda+2\mu-\frac{\lambda}{\lambda+2\mu}=4\mu\frac{\lambda+\mu}{\lambda+2\mu},$
(63)

being due to Eq. (60) equal to the Young modulus $E$ in two dimension.

#### Experimental relevant time scales.

Strictly speaking, $\mu=E=0$ holds for all systems including amorphous glasses
or even standard solids, if the mathematical limit “$\tau\to\infty$” is read
“for times much larger than any instrinsic relaxation time of the system”
Rubinstein and Colby (2003). In practice, the largest “terminal” relaxation
time $\tau_{\star}$ (for glasses called “$\alpha$-relaxation time” Ferry
(1980)) for many viscoelastic systems is commonly much larger than any
reasonable typical experimental or computational measurement time
$\tau_{\mathrm{exp}}$. Moreover, for many viscoelastic systems the relaxation
moduli are approximatively constant over many orders of magnitude in a time
window $\tau_{1}\ll\tau\ll\tau_{2}$. This means especially that the system can
support a finite shear stress in this time window. It is common that
$\tau_{1}$ is given by the time $\tau_{\mathrm{A}}$ needed to relax the
“affine strains” applied at $\tau=0$ or by the relaxation time
$\tau_{\mathrm{b}}$ within a meta-basin in glassy materials Wittmer et al.
(2013a, b, 2015a, 2015b); George et al. (2021a); Klochko et al. (2022) and
$\tau_{2}$ by the already mentioned terminal relaxation time $\tau_{\star}$.
(Obviously, much more complicated scenarios exist Ferry (1980); Rubinstein and
Colby (2003); Balucani and Zoppi (1995).) We assume in the following for the
simplicity of the discussion that this intermediate plateau regime becomes
very broad such that $\tau_{2}\approx\tau_{\star}$ exceeds
$\tau_{\mathrm{exp}}$ by many orders of magnitude. This is indeed the case,
e.g., for the shear-stress relaxation modulus $G(\tau)$ shown in Fig. 7 for
local MC moves (circles). The meaning of “$\tau\to\infty$” is thus (here as
elsewhere in this work) “for times much larger than $\tau_{1}$ but yet much
smaller than $\tau_{2}$”. In this sense Eq. (62) and Eq. (63) are still valid
with $\lambda$, $\mu$ or $E$ referring now to the (finite) plateau values. It
is these plateau values which are computed by means of out-of-equilibrium
methods or by means of stress or strain fluctuation formulae Allen and
Tildesley (2017) using times series with sampling times
$\tau_{1}\ll\Delta\tau\ll\tau_{2}$.

###  V.4 Key relations for overdamped motion

Previous studies Klochko et al. (2018, 2022); Maier et al. (2017, 2018); Vogel
et al. (2019) have focused on the correlation functions
$c^{\circ}_{\alpha\beta\gamma\delta}(\mbox{$\bf q$},\tau)$ of the
instantaneous “Invariant Stress Fields” (ISFs)
$\hat{\sigma}^{\circ}_{\alpha\beta}(\mbox{$\bf q$},\tau)$ in NRC. Importantly,
it has been shown Klochko et al. (2018) that the four “Invariant Correlation
Functions” (ICFs)

$\displaystyle c_{\mathrm{L}}(q,\tau)$ $\displaystyle\equiv$
$\displaystyle\left<c_{1111}^{\circ}(\mbox{$\bf q$},\tau)\right>_{{\hat{\bf
q}}},$ $\displaystyle c_{\perp}(q,\tau)$ $\displaystyle\equiv$
$\displaystyle\left<c_{1122}^{\circ}(\mbox{$\bf q$},\tau)\right>_{{\hat{\bf
q}}},$ $\displaystyle c_{\mathrm{T}}(q,\tau)$ $\displaystyle\equiv$
$\displaystyle\left<c_{1212}^{\circ}(\mbox{$\bf q$},\tau)\right>_{{\hat{\bf
q}}}\mbox{ and }$ $\displaystyle c_{\mathrm{N}}(q,\tau)$ $\displaystyle\equiv$
$\displaystyle\left<c_{2222}^{\circ}(\mbox{$\bf q$},\tau)\right>_{{\hat{\bf
q}}},$ (64)

averaged over all $\bf q$ of similar magnitude $q$ (within a given bin width
$\delta q$), can be expressed for small wavevectors in terms of two
independent invariant relaxation functions (“material functions”) Klochko et
al. (2018) characterizing the system. Naturally, these relations are
formulated in Fourier-Laplace space. As in Ref. Klochko et al. (2022) we use
here the modified Laplace transform

$f(s)={\cal L}[f(\tau)]=s\int_{0}^{\infty}\mathrm{d}\tau\ f(\tau)e^{-\tau s},$
(65)

called “$s$-transform” Klochko et al. (2022) or “Laplace-Carson transform”
Rubinstein and Rubinstein (1998), for which due to the prefactor $s$ the
original function $f(\tau)$ and its Laplace transform $f(s)$ have the same
dimension. We note by $L(s)={\cal L}[L(\tau)]=M(s)+2G(s)$, $M(s)={\cal
L}[M(\tau)]$ and $G(s)={\cal L}[G(\tau)]$ the Laplace transforms of the
invariant relaxation moduli introduced in Sec. V.3. In the present work MC
simulations have been used, i.e. not a momentum conserving dynamics as assumed
in Refs. Klochko et al. (2018, 2022) but an overdamped simulation dynamics
where the drift velocity of particles is proportional to the imposed body
force. Due to this the slightly modified relations become

$\displaystyle\beta Vc_{\mathrm{L}}(q,s)$ $\displaystyle=$
$\displaystyle\frac{sL(s)}{s+q^{2}L(s)/\zeta}$ (66) $\displaystyle\beta
Vc_{\perp}(q,s)$ $\displaystyle=$
$\displaystyle\frac{sM(s)}{s+q^{2}L(s)/\zeta}$ (67) $\displaystyle\beta
Vc_{\mathrm{T}}(q,s)$ $\displaystyle=$
$\displaystyle\frac{sG(s)}{s+q^{2}G(s)/\zeta}$ (68) $\displaystyle\beta
Vc_{\mathrm{N}}(q,s)$ $\displaystyle=$ $\displaystyle
L(s)-\frac{q^{2}M(s)^{2}/\zeta}{s+q^{2}L(s)/\zeta}$ (69)

with $\zeta$ being the friction constant of the overdamped dynamics foo (e).
All relations hold for momentum conserving dynamics if $\zeta$ is substituted
by the term $\rho s$ due to inertia with $\rho$ being the mass density. For
momentum conserving dynamics, Eq. (68) is well-established Balucani and Zoppi
(1995); Evans and Morriss (1990); its derivation based on the “Fluctuation-
Dissipation Theorem” (FDT) Doi and Edwards (1986); Forster (1995); Evans and
Morriss (1990); Chaikin and Lubensky (1995) is given in Ref. Semenov et al.
(2012) in the context of viscoelastic hydrodynamic interactions in polymer
liquids. The first relation Eq. (66) was already mentioned in Ref. Ruscher et
al. (2017). We note finally that these findings agree (where a comparison is
possible) with an independent and complementary approach to long-wavelength
stress-correlations in viscoelastic fluids at equilibrium based on the
Zwanzig-Mori projection operator formalism Maier et al. (2017, 2018); Vogel et
al. (2019).

###  V.5 Long-time limit for viscoelastic plateau

#### Asymptotic limit.

As noted in Sec. V.3 we focus on the long-time limit of the four key relations
where in a large time window $\tau_{1}\ll\tau\ll\tau_{2}$ all relaxation
moduli are approximatively constant. Due to the Laplace-Carson transform, Eq.
(65), Eq. (62) implies

$L(s)\to\lambda+2\mu,G(s)\to\mu,M(s)\to\lambda$ (70)

for small $s$ (with $\tau_{1}\ll 1/|s|\ll\tau_{2}$). Using the final-value
theorem of Laplace transforms we get

$\displaystyle c_{\mathrm{L}}(q,\tau)\approx c_{\perp}(q,\tau)\approx
c_{\mathrm{T}}(q,\tau)$ $\displaystyle\simeq$ $\displaystyle 0$ (71)
$\displaystyle\beta Vc_{\mathrm{N}}(q,\tau)\to
L(\tau)-\frac{M(\tau)^{2}}{L(\tau)}$ $\displaystyle\simeq$ $\displaystyle E$
(72)

for small $q$ and large $\tau$ foo (d). (Eq. (63) was used in the last step of
Eq. (72).) Using in addition Eq. (50) or Eq. (51) we have thus demonstrated
that for general viscoelastic fluids at thermal equilibrium the
phenomenological parameter $e$, cf. Eq. (40), is indeed given by the Young
modulus $E$ in two dimensions. As mentioned above (cf. Sec. IV), these results
are qualitatively understood by the facts that finite
$\bar{c}_{\mathrm{L}}(q)$, $\bar{c}_{\perp}(q)$ and $\bar{c}_{\mathrm{T}}(q)$
for large $\Delta\tau$ would violate mechanical equilibrium while a finite
$\bar{c}_{\mathrm{N}}(q)$ is possible since a finite transverse normal ISF
$\bar{\sigma}^{\circ}_{22}(\mbox{$\bf q$})$ cannot induce a force on a volume
element —due to Eq. (41) and $q_{2}^{\circ}=0$ in NRC for all $\bf q$ — and,
hence, no deterministic flux.

Figure 8: Scaling of ICFs $\bar{c}_{\mathrm{L}}(q,\Delta\tau)$,
$\bar{c}_{\mathrm{N}}(q,\Delta\tau)$, $\bar{c}_{\perp}(q,\Delta\tau)$ and
$\bar{c}_{\mathrm{T}}(q,\Delta\tau)$ for large $\Delta\tau$. We present the
rescaled ICFs $y=\beta V\bar{c}_{\mathrm{L}}/(\lambda+2\mu)$, $(\beta
V\bar{c}_{\mathrm{N}}-E)/\lambda\nu$ and $\beta V\bar{c}_{\perp}/\lambda$ vs.
$x=\Delta\tau/\tau_{\mathrm{L}}$ and $y=\beta V\bar{c}_{\mathrm{T}}/\mu$ vs.
$x=\Delta\tau/\tau_{\mathrm{T}}$ assuming $\zeta=900$ for $q=q_{\mathrm{min}}$
and $q\approx 0.18$ (filled symbols). The bold dashed line indicates the Debye
relation, the dashed-dotted line shows the expected asymptotic power-law
decay.

#### Predictions for leading deviations.

Importantly, the limits Eq. (71) and Eq. (72) do not depend (to leading order)
on the simulation dynamics, e.g. of whether we have used a momentum conserving
or overdamped dynamics. As emphasized by the dash-dotted line in Fig. 4
$\bar{c}_{\mathrm{L}}(\Delta\tau)$, $\bar{c}_{\perp}(\Delta\tau)$ and
$\bar{c}_{\mathrm{T}}(\Delta\tau)$ vanish inversely with $\Delta\tau$ as one
expects for uncorrelated fluctuations. We show in the remainder of this
subsection that these leading deviations for finite $\Delta\tau$ are
quantitatively described using Eqs. (66-69). To do this let us introduce the
characteristic times

$\tau_{\mathrm{L}}\equiv\zeta/q^{2}(\lambda+2\mu)\mbox{ and
}\tau_{\mathrm{T}}\equiv\zeta/q^{2}\mu$ (73)

characterizing, respectively, the longitudinal and the transverse overdamped
relaxation. Using Eq. (70) this leads to

$\displaystyle\beta Vc_{\mathrm{L}}(q,s)$ $\displaystyle\simeq$
$\displaystyle(\lambda+2\mu)\ \frac{s}{s+1/\tau_{\mathrm{L}}}$ (74)
$\displaystyle\beta Vc_{\perp}(q,s)$ $\displaystyle\simeq$
$\displaystyle\lambda\ \frac{s}{s+1/\tau_{\mathrm{L}}}$ (75)
$\displaystyle\beta Vc_{\mathrm{T}}(q,s)$ $\displaystyle\simeq$
$\displaystyle\mu\ \frac{s}{s+1/\tau_{\mathrm{T}}}$ (76) $\displaystyle\beta
Vc_{\mathrm{N}}(q,s)-E$ $\displaystyle\simeq$ $\displaystyle\lambda\nu\
\frac{s}{s+1/\tau_{\mathrm{L}}}.$ (77)

The Laplace-Carson transformation ${\cal L}[\exp(-t/\tau)]=s/(s+1/\tau)$ thus
gives

$\displaystyle\beta Vc_{\mathrm{L}}(q,t)$ $\displaystyle\simeq$
$\displaystyle(\lambda+2\mu)\exp(-t/\tau_{\mathrm{L}})$ (78)
$\displaystyle\beta Vc_{\perp}(q,t)$ $\displaystyle\simeq$
$\displaystyle\lambda\exp(-t/\tau_{\mathrm{L}})$ (79) $\displaystyle\beta
Vc_{\mathrm{T}}(q,t)$ $\displaystyle\simeq$
$\displaystyle\mu\exp(-t/\tau_{\mathrm{T}})$ (80) $\displaystyle\beta
Vc_{\mathrm{N}}(q,t)-E$ $\displaystyle\simeq$
$\displaystyle\lambda\nu\exp(-t/\tau_{\mathrm{L}}).$ (81)

for the ICFs of the instantaneous ISFs. Hence, Eq. (52) leads to

$\displaystyle\beta V\bar{c}_{\mathrm{L}}(q,\Delta\tau)$ $\displaystyle\simeq$
$\displaystyle(\lambda+2\mu)\ D(\Delta\tau/\tau_{\mathrm{L}})$ (82)
$\displaystyle\beta V\bar{c}_{\perp}(q,\Delta\tau)$ $\displaystyle\simeq$
$\displaystyle\lambda\ D(\Delta\tau/\tau_{\mathrm{L}})$ (83)
$\displaystyle\beta V\bar{c}_{\mathrm{T}}(q,\Delta\tau)$ $\displaystyle\simeq$
$\displaystyle\mu\ D(\Delta\tau/\tau_{\mathrm{T}})$ (84) $\displaystyle\beta
V\bar{c}_{\mathrm{N}}(q,\Delta\tau)-E$ $\displaystyle\simeq$
$\displaystyle\lambda\nu D(\Delta\tau/\tau_{\mathrm{L}})$ (85)

for the ICFs of the time-averaged ISFs.

#### Numerical test of predicted deviations.

These relations should allow a reasonable fit for large $\Delta\tau$ while
additional modes may be relevant in general. We test in Fig. 8 the suggested
scaling for $q=q_{\mathrm{min}}\approx 0.06$ and $q\approx 0.18$. The scaling
variable $x$ of the horizontal axis is $\Delta\tau/\tau_{\mathrm{L}}$ for
$\bar{c}_{\mathrm{L}}$, $\bar{c}_{\mathrm{N}}$ and $\bar{c}_{\perp}$ and
$\Delta\tau/\tau_{\mathrm{T}}$ for $\bar{c}_{\mathrm{T}}$. For the scaling of
$\bar{c}_{\mathrm{N}}$ we have to subtract $E$ and, unfortunately, the scaling
function $y=(\beta V\bar{c}_{\mathrm{N}}-E)/\lambda\nu$ (circles) is less
accurate than the scaling functions of the other three ICFs. Using the known
elastic moduli of our systems and fitting for only one parameter, the friction
coefficient $\zeta\approx 900$, we obtain for sufficiently large $x$ and small
$q$ a good data collapse on the Debye function (bold dashed line). Due to the
limited system size we have only achieved good data collapse for
$q_{\mathrm{min}}\leq q<0.2$. Larger system sizes are clearly warranted in
future work to verify over more than a decade the expected $q$-dependence of
the relaxation times $\tau_{L,T}\propto 1/q^{2}$. See Ref. Klochko et al.
(2022) for the discussion of other possible caveats leading to small
deviations.

##  VI Conclusion

#### Summary.

After presenting in Sec. II a survey of general useful mathematical relations
for isotropic tensor fields we have focused in Sec. IV on the numerical
description of spatial correlations of stress tensor fields in isotropic pLJ
particle systems as described in Sec. III and Appendix C. The phenomenological
parameter $e$ was related to the Young modulus $E$ in Sec. V. Several general
results deserve to be emphasized.

#### Angular dependences.

As demonstrated by our computational example, correlation functions of tensor
field components of perfectly isotropic systems must generally depend on the
components of the field vector and not only on its magnitude, cf. Eq. (IV).
They thus depend both on the direction of this vector in the physical system
(angle $\theta$) and on the orientation of the coordinate system (angle
$\alpha$) as shown in Fig. 1. Importantly, the angular dependence of all
correlation functions boils down to a dependence on the difference
$\theta-\alpha$ for all $\theta$ and $\alpha$ as shown in panel (b) of Fig. 1
and Fig. 5. Obviously, this simple scaling (without characteristic angles)
cannot hold for true anisotropic systems which have material functions
depending explicitly on the direction of the field vector. Unfortunately, many
recent studies do not clarify, e.g., simply by rotating the coordinate system
or equivalently the orientations of the experimental devices measuring the
local tensor field components, which of these two mathematically and
physically very different types of “anisotropy” is involved.

#### Invariant correlation functions.

Just as the invariance of a tensor under orthogonal transformations is only
revealed by the full set of tensor components, one component of the
correlation tensor field may appear to be inconsistent with the assumed
isotropy which becomes only manifested for the total set of correlation
functions. Fortunately, this set is completely characterized by a few ($4$ in
$d=2$ and $5$ in $d>2$) ICFs, cf. Eq. (II.5). Only if a measured subset of
correlation functions cannot be consistently expressed by ICFs a legitimate
claim of “anisotropy” can be made. As we have seen, these ICFs are best
described both theoretically and numerically (cf. Fig. 3) in reciprocal space
using NRC. The investigated system was deliberately simple since only one ICF,
namely the normal stress ICF $\bar{c}_{\mathrm{N}}(q)$, remains finite for
large $\Delta\tau$ and, moreover, constant for sufficiently small $q$, Eq.
(40). This plateau in reciprocal space explains directly (cf. Appendix B) the
long-ranged $1/r^{2}$-decay of the stress correlations in real space
demonstrated in Fig. 6 foo (f). Only one phenomenological parameter $e$ is
needed for the fitting of our data in the large-$\Delta\tau$ limit, cf. Eq.
(40).

#### Time-averaged stress fields.

We have studied in this work correlation functions
$\bar{c}_{\alpha\beta\gamma\delta}(\mbox{$\bf q$},\Delta\tau)$ of time-
averaged stress fields $\bar{\sigma}_{\alpha\beta}(\mbox{$\bf q$},\Delta\tau)$
and not the correlation functions $c_{\alpha\beta\gamma\delta}(\mbox{$\bf
q$},t)$ of instantaneous fields $\sigma_{\alpha\beta}(\mbox{$\bf q$},t)$. As
reminded in Sec. V.2, for stationary processes
$\bar{c}_{\alpha\beta\gamma\delta}(\mbox{$\bf q$},\Delta\tau)$ and
$c_{\alpha\beta\gamma\delta}(\mbox{$\bf q$},t)$ are closely related, cf. Eq.
(49), and both correlation functions have the same asymptote for large
(sampling) times, cf. Eq. (50), i.e. the same static properties are probed.
From the numerical side $\bar{c}_{\alpha\beta\gamma\delta}(\mbox{$\bf
q$},\Delta\tau)$ basically has the advantage that irrelevant high frequencies
are filtered off (as it is usually the case for Einstein relations compared to
Green-Kubo relations Allen and Tildesley (2017)). Compared with the common
method to determine the system’s inherent quenched stresses by quenching to
the local energy minima Lemaître (2015, 2018) this has the advantage that we
measure the relevant stresses at the given finite temperature and not at a
different thermodynamic state. We believe that this method is therefore not
only simpler (no additional quench required) but thermodynamically and
conceptionally better defined. The technical downside is that time-series of
many frames need to be stored.

#### Symmetry breaking.

The main physical reason for focusing on $t$-averaged stress fields
$\bar{\sigma}_{\alpha\beta}(\mbox{$\bf q$},\Delta\tau)$ is, however, that the
systematic time-averaging naturally allows to focus on the quenched stresses
$\sigma^{\mathrm{q}}_{\alpha\beta}(\mbox{$\bf q$})$ for each independent
configuration $c$ by projecting out the trivial instantaneous thermal stress
fluctuations using that $\bar{\sigma}_{\alpha\beta}(\mbox{$\bf
q$},\Delta\tau)\to\sigma^{\mathrm{q}}_{\alpha\beta}(\mbox{$\bf q$})$ for large
$\Delta\tau$. Importantly, the parameter $e$ is set by Eq. (42) in terms of
the typical size of the quenched normal stresses
$\bar{\sigma}^{\circ\mathrm{q}}_{22}(\mbox{$\bf q$})$ in reciprocal space
using NRC. Due to Eq. (41) these stresses may be finite without violating
static mechanical equilibrium as shown by the finite attractor in Fig. 4.

#### Connection to thermodynamic response properties.

Interestingly, $e$ was found to be numerically similar to the macroscopic
Young modulus $E$ at thermal equilibrium in $d=2$. As discussed in Sec. V,
this finding was theoretically anticipated Klochko et al. (2022) since the
predictions for equilibrium viscoelastic fluids, must hold approximatively for
our glasses due to the swap MC moves used for the system preparation. (The
same identification is obtained by an independent and complementary approach
to long-wavelength stress-correlations in glass-forming liquids based on the
Zwanzig-Mori projection operator formalism Maier et al. (2017, 2018); Vogel et
al. (2019).) Altogether this implies

$E=e\simeq\beta V\left<\bar{\sigma}^{\circ\mathrm{q}}_{22}(\mbox{$\bf
q$})\bar{\sigma}^{\circ\mathrm{q}}_{22}(-\mbox{$\bf q$})\right>$ (86)

for sufficiently small wavevectors. The existence of a finite equilibrium
Young modulus $E$ (and, hence, of a finite shear modulus $\mu$) and the broken
(continuous) symmetry of the stress field in reciprocal space, characterized
by the typical size of the quenched stresses
$\bar{\sigma}^{\circ\mathrm{q}}_{22}(\mbox{$\bf q$})$, are thus intimately
connected.

#### Outlook.

The presented work suggests several natural extensions:

  * •

The given tensor fields relations for isotropic systems naturally generalize
to higher spatial dimensions, to tensor fields of different order and to
crosscorrelation functions of different tensorial fields (say, between stress
and strain fields) for which the major suffix symmetry Eq. (18) may not hold
(which merely introduces one additional ICF).

  * •

We have focused in the present work on Euclidean spaces and Carthesian
coordinates. It is possible to generalize our relations for systems embedded
in non-Euclidean spaces, say for glasses on spheres Vest et al. (2015); Turci
et al. (2017), and more general curvilinear coordinate systems McConnell
(2021); Lambourne (2010). Assuming the system to be locally isotropic the goal
is always to construct from the available correlation functions the local
(tangent) isotropic invariants and to check whether this can be done
consistently with all available data sets.

  * •

The proposed methodology could, e.g., also be used for irreversibly
crosslinked polymer networks Rubinstein and Colby (2003) or active non-
equilibrium systems Grosberg and Joanny (2015); Elismaili et al. (2020) which
are both isotropic and highly viscous or even jammed. $e$ is in such cases
just a fitting parameter characterizing the typical size of the quenched
stresses.

  * •

The presented work is also of relevance for the correlations of strain tensor
fields Nicolas et al. (2014); Liu et al. (2021); Chacko et al. (2021) and for
the characterization of plastic deformations Voigtmann (2014); Fielding
(2014); foo (f). It can be shown that the strain correlations of any isotropic
elastic body are characterized by two ICFs which are, moreover, set by means
of the general equipartition theorem Chaikin and Lubensky (1995) in the large-
wavelength limit by the Lamé coefficients $\lambda$ and $\mu$. This leads
again to an octupolar correlation field pattern as seen in Fig. 1 for the
stress.

  * •

The low-$q$ limit $e(\Delta\tau)$ of $\beta
V\bar{c}_{\mathrm{N}}(q,\Delta\tau)$ becomes only a strictly
$\Delta\tau$-independent constant for permanently quenched invariant stress
fields $\bar{\sigma}^{\circ\mathrm{q}}_{22}(\mbox{$\bf q$})$, e.g., for
permanent polymer networks Rubinstein and Colby (2003). By contrast, the
quenched attractors $\bar{\sigma}^{\circ\mathrm{q}}_{22}(\mbox{$\bf q$})$ of
low-temperature glasses become for $\Delta\tau$ similar to the
$\alpha$-relaxation time slow $\Delta\tau$-dependent fields
$\sigma^{\circ\mathrm{s}}_{22}(\mbox{$\bf q$},\Delta\tau)$. This is, e.g., of
relevance for higher temperatures $T$ of the presented pLJ particle model
where after a transient plateau $e(\Delta\tau)$ ultimately vanishes, i.e. the
long-range stress correlations disappear Klochko et al. (2022). Future work
should focus on the distributions of the
$\bar{\sigma}^{\circ}_{\alpha\beta}(\mbox{$\bf q$},\Delta\tau,T)$ as a direct
means to systematically project out fast relaxation modes independently of the
coordinate system.

#### Acknowledgments.

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

##  Appendix A Fourier transformation

#### Continuous Fourier transform.

Following Refs. Wittmer et al. (2022); Klochko et al. (2022) but at variance
with Ref. Klochko et al. (2018) we define in this work the Fourier transform
(FT) $f(\mbox{$\bf q$})={\cal F}[f(\mbox{$\bf r$})]$ of a (real-valued)
function $f(\mbox{$\bf r$})$ in real space by

$f(\mbox{$\bf q$})=\frac{1}{V}\int\mathrm{d}\mbox{$\bf r$}\ f(\mbox{$\bf
r$})\exp(-i\mbox{$\bf q$}\cdot\mbox{$\bf r$})$ (87)

with $V$ being the volume of the system and $\bf q$ a wavevector commensurate
with the simulation box. The inverse FT is then given by

$f(\mbox{$\bf r$})={\cal F}^{-1}[f(\mbox{$\bf
q$})]=\frac{V}{(2\pi)^{d}}\int\mathrm{d}\mbox{$\bf q$}\ f(\mbox{$\bf
q$})\exp(i\mbox{$\bf q$}\cdot\mbox{$\bf r$}).$ (88)

Note that $f(\mbox{$\bf r$})$ and $f(\mbox{$\bf q$})$ have the same dimension.
For notational simplicity the function names remain unchanged. Let us denote
by $f$ without argument the macroscopic field average. Due to Eq. (87) we thus
have the “sum rule”

$f(\mbox{$\bf q$}={\bf 0})=f\equiv\frac{1}{V}\int\mathrm{d}\mbox{$\bf r$}\
f(\mbox{$\bf r$}).$ (89)

We note for later convenience the FTs

$\displaystyle{\cal F}\left[\frac{\partial}{\partial r_{\alpha}}f(\mbox{$\bf
r$})\right]$ $\displaystyle=$ $\displaystyle iq_{\alpha}f(\mbox{$\bf q$})$
(90) $\displaystyle{\cal F}\left[\delta(\mbox{$\bf r$}-\mbox{$\bf v$})\right]$
$\displaystyle=$ $\displaystyle\frac{1}{V}\exp(-i\mbox{$\bf q$}\cdot\mbox{$\bf
v$})\mbox{ and }$ (91) $\displaystyle{\cal F}\left[\int_{0}^{1}\mathrm{d}s\
\delta(\mbox{$\bf r$}-\mbox{$\bf v$}s)\right]$ $\displaystyle=$
$\displaystyle\frac{1}{V}\frac{1-\exp(-i\mbox{$\bf q$}\cdot\mbox{$\bf
v$})}{i\mbox{$\bf q$}\cdot\mbox{$\bf v$}}$ (92)

with $\delta(\mbox{$\bf r$})$ being Dirac’s delta function. Let us consider
the spatial correlation function

$c(\mbox{$\bf r$})=\frac{1}{V}\int\mathrm{d}\mbox{$\bf
r$}^{\prime}g(\mbox{$\bf r$}+\mbox{$\bf r$}^{\prime})h(\mbox{$\bf
r$}^{\prime})$ (93)

where the fields $g(\mbox{$\bf r$})$ and $h(\mbox{$\bf r$})$ are assumed to be
real. According to the “correlation theorem” Press et al. (1992) this becomes

$c(\mbox{$\bf q$})=g(\mbox{$\bf q$})h^{\star}(\mbox{$\bf q$})=g(\mbox{$\bf
q$})h(-\mbox{$\bf q$})$ (94)

in reciprocal space (with $\star$ marking the complex conjugate). For auto-
correlation functions, i.e. for $g(\mbox{$\bf r$})=h(\mbox{$\bf r$})$, this
simplifies to (“Wiener-Khinchin theorem”)

$c(\mbox{$\bf q$})=g(\mbox{$\bf q$})g^{\star}(\mbox{$\bf q$})=|g(\mbox{$\bf
q$})|^{2},$ (95)

i.e. the Fourier transformed auto-correlation functions are real and $\geq 0$
for all $\bf q$. Moreover, we shall consider correlation functions
$c(\mbox{$\bf r$})$, Eq. (93), being even in real space, $c(\mbox{$\bf
r$})=c(-\mbox{$\bf r$})$, and thus also in reciprocal space, $c(\mbox{$\bf
q$})=c(-\mbox{$\bf q$})=c^{\star}(\mbox{$\bf q$})$, i.e. $c(\mbox{$\bf q$})$
is real.

For many reasons it is convenient to consider correlation functions
$c(\mbox{$\bf r$})$ which vanish for large $r=|\mbox{$\bf r$}|$. This is
achieved here by replacing in Eq. (93) $g(\mbox{$\bf r$})$ by
$\tilde{g}(\mbox{$\bf r$})=g(\mbox{$\bf r$})-g$ and $h(\mbox{$\bf r$})$ by
$\tilde{h}(\mbox{$\bf r$})=h(\mbox{$\bf r$})-h$ with $g$ and $h$ being field
averages according to Eq. (89). We thus probe in this work correlation
functions

$\displaystyle c(\mbox{$\bf r$})$ $\displaystyle=$
$\displaystyle\frac{1}{V}\int\mathrm{d}\mbox{$\bf
r$}^{\prime}\tilde{g}(\mbox{$\bf r$}+\mbox{$\bf
r$}^{\prime})\tilde{h}(\mbox{$\bf r$}^{\prime})$ (96) $\displaystyle=$
$\displaystyle\frac{1}{V}\int\mathrm{d}\mbox{$\bf r$}^{\prime}g(\mbox{$\bf
r$}+\mbox{$\bf r$}^{\prime})h(\mbox{$\bf r$}^{\prime})-gh$

using a real offset $gh$. Eq. (94) thus becomes

$c(\mbox{$\bf q$})=\tilde{g}(\mbox{$\bf q$})\tilde{h}(-\mbox{$\bf
q$})=\left\\{\begin{array}[]{ll}g(\mbox{$\bf q$})h(-\mbox{$\bf q$})&\mbox{ for
}\mbox{$\bf q$}\neq{\bf 0}\\\ 0&\mbox{ for }\mbox{$\bf q$}={\bf
0}\end{array}\right.$ (97)

where we have used that $g=g(\mbox{$\bf q$}={\bf 0})$ and $h=h(\mbox{$\bf
q$}={\bf 0})$.

#### Discrete FT on microcell grid.

Numerically, all fields $f(\mbox{$\bf r$})$ are stored on a regular
equidistant grid as shown in Fig. 2. Periodic boundary conditions are assumed
Allen and Tildesley (2017). The discrete FT and its inverse become

$\displaystyle f(\mbox{$\bf q$})$ $\displaystyle=$
$\displaystyle\frac{1}{n_{V}}\sum_{\mbox{$\bf r$}}f(\mbox{$\bf
r$})\exp(-i\mbox{$\bf q$}\cdot\mbox{$\bf r$})$ (98) $\displaystyle
f(\mbox{$\bf r$})$ $\displaystyle=$ $\displaystyle\sum_{\mbox{$\bf
q$}}f(\mbox{$\bf q$})\exp(i\mbox{$\bf q$}\cdot\mbox{$\bf r$})$ (99)

with $\sum_{\mbox{$\bf r$}}$ and $\sum_{\mbox{$\bf q$}}$ being discrete sums
over $n_{V}=V/a_{\mathrm{grid}}^{d}$ grid points in, respectively, real or
reciprocal space. With $n_{L}$ being the number of grid points in each spatial
direction $\alpha$, i.e. $n_{V}=n_{L}^{d}$ and $L=n_{L}a_{\mathrm{grid}}$, we
have

$\frac{r_{\alpha}}{a_{\mathrm{grid}}}=n_{\alpha}\mbox{ and
}q_{\alpha}a_{\mathrm{grid}}=\frac{2\pi}{n_{L}}n_{\alpha}$ (100)

with integers $n_{\alpha}=0,1,\ldots,n_{L}-1$. To take advantage of the
implemented Fast-Fourier transform (FFT) routines Press et al. (1992) $n_{L}$
is an integer-power of $2$. The real space correlation function $c(\mbox{$\bf
r$})$ on the discrete grid is

$c(\mbox{$\bf r$})=\frac{1}{n_{V}}\sum_{\mbox{$\bf r$}^{\prime}}g(\mbox{$\bf
r$}+\mbox{$\bf r$}^{\prime})h(\mbox{$\bf r$}^{\prime})-gh$ (101)

being an operation of order ${\cal O}(n_{V}^{2})$. Importantly, the
periodicity of all fields must be taken into account. It is obviously much
more efficient to first FFT the discrete fields $g(\mbox{$\bf r$})$ and
$h(\mbox{$\bf r$})$ or $\tilde{g}(\mbox{$\bf r$})$ and $\tilde{h}(\mbox{$\bf
r$})$ and to apply then Eq. (94) Allen and Tildesley (2017); Press et al.
(1992). In this manner $c(\mbox{$\bf q$})$ is automatically periodic in all
spatial directions of the simulation box and the same applies to the
correlation function $c(\mbox{$\bf r$})={\cal F}^{-1}[c(\mbox{$\bf q$})]$ in
real space computed using Eq. (99).

##  Appendix B Fourier transforms in two dimensions

#### Correlation functions in reciprocal space.

As shown in Fig. 3 all four ICFs in $d=2$ become constant in reciprocal space
for sufficiently small wavevectors $q$ and large sampling times $\Delta\tau$,
but only $\beta V\bar{c}_{\mathrm{N}}(q)\simeq e$ remains finite. According to
Eq. (33) of Sec. II.6 we thus have

$\displaystyle\beta V\bar{c}_{\alpha\beta\gamma\delta}(\mbox{$\bf q$})$
$\displaystyle\simeq$ $\displaystyle e\times$ (102)
$\displaystyle(\delta_{\alpha\beta}\delta_{\gamma\delta}$ $\displaystyle-$
$\displaystyle{\hat{q}}_{\alpha}{\hat{q}}_{\beta}\delta_{\gamma\delta}-{\hat{q}}_{\gamma}{\hat{q}}_{\delta}\delta_{\alpha\beta}+{\hat{q}}_{\alpha}{\hat{q}}_{\beta}{\hat{q}}_{\gamma}{\hat{q}}_{\delta})$

for the correlation functions in the “old” (unrotated) reference frame. Hence,
we have, e.g.,

$\displaystyle\beta V\bar{c}_{1212}(\mbox{$\bf q$})$ $\displaystyle=$
$\displaystyle\beta V\bar{c}_{1122}(\mbox{$\bf q$})$ $\displaystyle\simeq$
$\displaystyle e\ {\hat{q}}_{1}^{2}{\hat{q}}_{2}^{2}$ (103)
$\displaystyle\beta V\bar{c}_{1111}(\mbox{$\bf q$})$ $\displaystyle\simeq$
$\displaystyle
e\left(1-2{\hat{q}}_{1}^{2}+{\hat{q}}_{1}^{4}\right)=e{\hat{q}}_{2}^{4}$ (104)
$\displaystyle\beta V\bar{c}_{2222}(\mbox{$\bf q$})$ $\displaystyle\simeq$
$\displaystyle
e\left(1-2{\hat{q}}_{2}^{2}+{\hat{q}}_{2}^{4}\right)=e{\hat{q}}_{1}^{4}$ (105)
$\displaystyle\beta V\bar{c}_{1112}(\mbox{$\bf q$})$ $\displaystyle\simeq$
$\displaystyle-e\ {\hat{q}}_{1}{\hat{q}}_{2}^{3}$ (106) $\displaystyle\beta
V\bar{c}_{2221}(\mbox{$\bf q$})$ $\displaystyle\simeq$ $\displaystyle-e\
{\hat{q}}_{1}^{3}{\hat{q}}_{2}.$ (107)

Let us rewrite these cases in terms of the angle $\theta_{q}$ of the
normalized wavevector ${\hat{\bf q}}=\ (\cos(\theta_{q}),\sin(\theta_{q}))$.
Using standard trigonometric relations it is readily seen that

$\displaystyle-\beta V\bar{c}_{1212}(\mbox{$\bf q$})$ $\displaystyle=$
$\displaystyle-\beta V\bar{c}_{1122}(\mbox{$\bf q$})$ $\displaystyle\simeq$
$\displaystyle\frac{e}{8}\left(\cos(4\theta_{q})-1\right)$ (108)
$\displaystyle\beta V(\bar{c}_{1111}(\mbox{$\bf q$})+\bar{c}_{2222}(\mbox{$\bf
q$}))/2$ $\displaystyle\simeq$
$\displaystyle\frac{e}{8}\left(\cos(4\theta_{q})+3\right)$ (109)
$\displaystyle\beta V(\bar{c}_{1111}(\mbox{$\bf q$})-\bar{c}_{2222}(\mbox{$\bf
q$}))/2$ $\displaystyle\simeq$ $\displaystyle-\frac{e}{2}\cos(2\theta_{q})$
(110) $\displaystyle\beta V(\bar{c}_{1112}(\mbox{$\bf
q$})+\bar{c}_{2221}(\mbox{$\bf q$}))/2$ $\displaystyle\simeq$
$\displaystyle-\frac{e}{4}\sin(2\theta_{q})$ (111) $\displaystyle\beta
V(\bar{c}_{1112}(\mbox{$\bf q$})-\bar{c}_{2221}(\mbox{$\bf q$}))/2$
$\displaystyle\simeq$ $\displaystyle\frac{e}{8}\sin(4\theta_{q})$ (112)

in the above-mentioned limits. More generally, all correlation functions in
two dimensions can be expressed by a linear superposition of the orthogonal
basis functions $\cos(p\theta_{q})$ and $\sin(p\theta_{q})$ with $p=0$, $2$
and $4$.

#### Inverse Fourier transform.

To obtain the correlation functions in real space we thus have to compute the
inverse FT in $d=2$ dimensions for

$f(\mbox{$\bf q$})=\frac{1}{V}\ v(q)\ b_{p}(\theta_{q})$ (113)

with $b_{p}(\theta_{q})$ standing for the basis functions $\cos(p\theta_{q})$
or $\sin(p\theta_{q})$ with $p=2$ and $4$. The constant terms ($p=0$) in Eq.
(108) and Eq. (109) are irrelevant for our considerations leading merely to
$\delta$-contributions at $r=0$. The yet unspecified auxiliary function
$v(q)$, with $v(q)\to e$ for $q\to 0$, will be chosen below in a convenient
manner to take advantage of special mathematical functions Abramowitz and
Stegun (1964). Hence, using Eq. (88) we have

$\displaystyle f(\mbox{$\bf r$})$ $\displaystyle=$
$\displaystyle\frac{1}{4\pi^{2}}\int\mathrm{d}q\ qv(q)\times$ (114)
$\displaystyle\int_{0}^{2\pi}\mathrm{d}\theta_{q}\
b_{p}(\theta_{q})\exp[iqr\cos(\theta_{q}-\theta_{r})]$

with $\theta_{r}$ being the angle of ${\hat{\bf
r}}=(\cos(\theta_{r}),\sin(\theta_{r}))$. We make now the substitution
$\theta=\theta_{q}-\theta_{r}$ and use that Abramowitz and Stegun (1964)

$\displaystyle\cos(p\theta+p\theta_{r})+\cos(-p\theta+p\theta_{r})$
$\displaystyle=$ $\displaystyle 2\cos(p\theta)\cos(p\theta_{r})$
$\displaystyle\sin(p\theta+p\theta_{r})+\sin(-p\theta+p\theta_{r})$
$\displaystyle=$ $\displaystyle 2\cos(p\theta)\sin(p\theta_{r}).$

We remind that following Eq. (9.1.21) of Ref. Abramowitz and Stegun (1964) the
integer Bessel function $J_{p}(z)$ may be written

$J_{p}(z)=\frac{i^{-p}}{\pi}\int_{0}^{\pi}\mathrm{d}\theta\cos(p\theta)\exp[iz\cos(\theta)].$
(115)

This leads to

$f(\mbox{$\bf r$})=\frac{i^{p}}{2\pi}b_{p}(\theta_{r})\ \int\mathrm{d}q\
qv(q)\ J_{p}(rq).$ (116)

We use next Eq. (11.4.28) of Ref. Abramowitz and Stegun (1964)

$\displaystyle\int_{0}^{\infty}e^{-(at)^{2}}t^{\mu-1}J_{\nu}(bt)\mathrm{d}t$
$\displaystyle=$ $\displaystyle\frac{\Gamma[(\nu+\mu)/2]\
(b/2a)^{\nu}}{2a^{\mu}\Gamma(\nu+1)}$ (117) $\displaystyle\times\
M[(\nu+\mu)/2,\nu+1,-(b/2a)^{2}]$

for $\Re(\mu+\nu)>0$ and $\Re(a^{2})>0$ relating the general Bessel function
$J_{\nu}(bt)$ to the confluent hypergeometric Kummer function $M(a,b,z)$.
($\Gamma(x)$ denotes the standard Gamma function Abramowitz and Stegun
(1964).) To take advantage of Eq. (117) we finally set
$v(q)=e\exp[-(aq)^{2}]$. Note that $v(q)\to e$ for $aq\to 0$, i.e. the
auxiliary variable $a$ becomes irrelevant for small $q$. We thus rewrite Eq.
(116) as

$\displaystyle f(\mbox{$\bf r$})$ $\displaystyle=$ $\displaystyle e\
\frac{i^{p}}{2\pi}\frac{\Gamma(p/2+1)(r/2a)^{p}}{2a^{2}\Gamma(p+1)}\ \
b_{p}(\theta_{r})$ (118) $\displaystyle\times$ $\displaystyle
M[p/2+1,p+1,-(r/2a)^{2}]$

in terms of Kummer’s function $M$. Following Eq. (13.1.5) of Ref. Abramowitz
and Stegun (1964)

$M(a,b,z)=\frac{\Gamma(b)}{\Gamma(b-a)}(-z)^{-a}\left(1+{\cal
O}(|z|^{-1}\right)$ (119)

for a real part $\Re(z)<0$. Using this expansion in Eq. (118) finally leads to

$f(\mbox{$\bf r$})\simeq e\ \frac{i^{p}p}{2\pi r^{2}}\ b_{p}(\theta_{r})$
(120)

for sufficiently large $r$ and $p>0$. We note that the auxiliary variable $a$
indeed drops out, that $f(\mbox{$\bf r$})=f(-\mbox{$\bf r$})$ and that
$f(\mbox{$\bf r$})$ is real for even $p$.

#### Correlation functions in real space.

Using Eq. (120) we can finally restate Eqs. (108-112) in real space as

$\displaystyle-\beta\bar{c}_{1122}(\mbox{$\bf
r$})=-\beta\bar{c}_{1212}(\mbox{$\bf r$})$ $\displaystyle=$
$\displaystyle\beta\left(\bar{c}_{1111}(\mbox{$\bf
r$})+\bar{c}_{2222}(\mbox{$\bf r$})\right)/2$ $\displaystyle\simeq$
$\displaystyle\frac{e}{4\pi r^{2}}\cos(4\theta_{r})$ (121)
$\displaystyle\beta\left(\bar{c}_{1111}(\mbox{$\bf
r$})-\bar{c}_{2222}(\mbox{$\bf r$})\right)/2$ $\displaystyle\simeq$
$\displaystyle\frac{e}{2\pi r^{2}}\cos(2\theta_{r})$ (122)
$\displaystyle\beta\left(\bar{c}_{1112}(\mbox{$\bf
r$})+\bar{c}_{2221}(\mbox{$\bf r$})\right)/2$ $\displaystyle\simeq$
$\displaystyle\frac{e}{4\pi r^{2}}\sin(2\theta_{r})$ (123)
$\displaystyle\beta\left(\bar{c}_{1112}(\mbox{$\bf
r$})-\bar{c}_{2221}(\mbox{$\bf r$})\right)/2$ $\displaystyle\simeq$
$\displaystyle\frac{e}{4\pi r^{2}}\sin(4\theta_{r})$ (124)

for $r>0$ and $\Delta\tau\to\infty$. Some of these relations are shown in Fig.
5 and Fig. 6 in the main part of this work.

##  Appendix C Additional computational details

#### 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 Wittmer et al. (2013a); Klochko
et al. (2019); George et al. (2021a, b); Klochko et al. (2021, 2022)

$u(s)=4\epsilon\left(\frac{1}{s^{12}}-\frac{1}{s^{6}}\right)\mbox{ with
}s=\frac{r}{(D_{i}+D_{j})/2}$ (125)

being the reduced distance according to the Lorentz rule Hansen and McDonald
(2006). This potential is truncated and shifted Allen and Tildesley (2017)
with a cutoff $s_{\mathrm{cut}}=2s_{\mathrm{min}}$ given by the minimum
$s_{\mathrm{min}}$ of $u(s)$. Lennard-Jones units Allen and Tildesley (2017)
are used throughout this study, i.e. $\epsilon=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_{\mathrm{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 Allen and Tildesley
(2017). Time is measured in units of MC cycles throughout this work.

#### Operational parameters.

We focus in the present work on configurations with $n=10000$ particles albeit
we have sampled a broad range of particle numbers $n$ between $n=100$ and
$n=160000$ (cf. inset of Fig. 6). $N_{\mathrm{c}}=200$ independent
configurations $c$ have been obtained by quenching configurations equilibrated
at a high temperature $T=0.55$ in the liquid limit. This is done using a
combination of local MC moves Allen and Tildesley (2017) and swap MC moves
exchanging the sizes of pairs of particles Ninarello et al. (2017); George et
al. (2021a). In addition an MC barostat Allen and Tildesley (2017) imposes an
average normal stress $P=2$ Wittmer et al. (2013a); George et al. (2021a).

#### Working temperature.

At the working temperature $T=0.2$ we first thoroughly temper over
$\Delta\tau=10^{7}$ all configurations with switched-on local, swap and
barostat MC moves and then again over $\Delta\tau=10^{7}$ with switched-on
local and swap moves and switched-off barostat moves. The shear-stress
relaxation function $G(\tau)$ Wittmer et al. (2015a, b); George et al.
(2021a); Klochko et al. (2022) for this tempering run is shown in Fig. 7. As
can be seen (squares), $G(\tau)$ rapidly decays to a shear modulus $\mu=0$
(dashed line) as expected for an equilibrium liquid. The final production runs
are carried out at constant volume $V$ only keeping the slow local moves.
Importantly, $T=0.2$ is well below the known glass transition temperature
$T_{\mathrm{g}}\approx 0.26$ assuming only local MC hopping moves Wittmer et
al. (2013a); George et al. (2021a). As shown by the circles in Fig. 7 only
keeping the local moves essentially traps the (equilibrated) configurations in
local metabasins where they behave as elastic bodies with a finite shear
modulus $\mu\approx 14$ (solid line). Due to the barostat used for the
quenching of the configuration the box volume $V=L^{d}$ differs slightly
between different configurations $c$ while $V$ is identical for all
(correlated) configurations $t$ of the time-series of the same independent
configuration $c$. In all cases the number density is of order unity, i.e. the
particle number $n$ and the volume $V$ are similar. This implies that the
ideal pressure $P_{\mathrm{id}}=Tn/V\approx 0.2$ is much smaller than the
imposed total pressure $P=P_{\mathrm{id}}+P_{\mathrm{ex}}=2$. As shown
elsewhere Klochko et al. (2019); George et al. (2021a, b); Klochko et al.
(2021, 2022) our systems are homogeneous and isotropic and crystallization is
irrelevant.

#### Data sampling.

For each of the $N_{\mathrm{c}}=200$ independent configurations $c$ we sample
and store four ensembles of time series containing each $N_{\mathrm{t}}=10000$
instantaneous “frames” $t$. These are obtained using the equidistant time
intervals $\delta\tau=1$, $10$, $100$ and $1000$ foo (b). This implies for
each time series a largest sampling time
$\Delta\tau=N_{\mathrm{t}}\delta\tau$, i.e. $\Delta\tau=10^{7}$ for the
largest time interval. This was chiefly done to check that all correlation
functions $\bar{c}_{\alpha\beta\gamma\delta}(\mbox{$\bf r$})$ become indeed
$\Delta\tau$-independent as shown in Fig. 4 and Fig. 8. The time-averaged
stress tensors $\bar{\sigma}_{\alpha\beta}|_{c}$ and stress tensor fields
$\bar{\sigma}_{\alpha\beta}(\mbox{$\bf r$})|_{c}$ are obtained by averaging
over the corresponding instantaneous $\sigma_{\alpha\beta}(t)|_{c}$ and
$\sigma_{\alpha\beta}(\mbox{$\bf r$},t)|_{c}$ using each frame $t$ of a time-
series. We thus obtain first

$\left.\bar{c}_{\alpha\beta\gamma\delta}(\mbox{$\bf
q$})\right|_{c}=\left\\{\begin{array}[]{ll}\left.\bar{\sigma}_{\alpha\beta}(\mbox{$\bf
q$})\right|_{c}\left.\bar{\sigma}_{\gamma\delta}(-\mbox{$\bf
q$})\right|_{c}&\mbox{ for }\mbox{$\bf q$}\neq{\bf 0}\\\ 0&\mbox{ for
}\mbox{$\bf q$}={\bf 0}\end{array}\right.$ (126)

for each $c$ according to Eq. (97), take then the $c$-average

$\bar{c}_{\alpha\beta\gamma\delta}(\mbox{$\bf
q$})=\frac{1}{N_{\mathrm{c}}}\sum_{c=1}^{N_{\mathrm{c}}}\left.\bar{c}_{\alpha\beta\gamma\delta}(\mbox{$\bf
q$})\right|_{c},$ (127)

and perform finally the inverse FFT to real space

$\bar{c}_{\alpha\beta\gamma\delta}(\mbox{$\bf r$})={\cal
F}^{-1}[\bar{c}_{\alpha\beta\gamma\delta}(\mbox{$\bf q$})].$ (128)

##  Appendix D Construction of stress field

#### Macroscopic stresses.

Carets “$\hat{a}$” mark here instantaneous properties and the argument $t$ is
dropped. Let us first remind the Irving-Kirkwood formula for the total
macroscopic ($\mbox{$\bf q$}={\bf 0}$) stress tensor
$\hat{\sigma}_{\alpha\beta}$ Irving and Kirkwood (1950); Allen and Tildesley
(2017). Note that the total system Hamiltonian is the sum of an ideal
contribution, depending only on the momenta, and an excess contribution,
depending only on the particle positions. This implies Allen and Tildesley
(2017) that the total stress tensor
$\hat{\sigma}_{\alpha\beta}=\hat{\sigma}^{\mathrm{id}}_{\alpha\beta}+\hat{\sigma}^{\mathrm{ex}}_{\alpha\beta}$
is a sum of an ideal stress $\hat{\sigma}^{\mathrm{id}}_{\alpha\beta}$ and an
excess stress $\hat{\sigma}^{\mathrm{ex}}_{\alpha\beta}$. The ideal
contribution is Allen and Tildesley (2017)

$\hat{\sigma}^{\mathrm{id}}_{\alpha\beta}=-\frac{1}{V}\sum_{a=1}^{n}\frac{p^{a}_{\alpha}p^{a}_{\beta}}{m}$
(129)

with $p^{a}_{\alpha}$ being the $\alpha$-component of the momentum of particle
$a$. This contribution is naturally not accessible in an MC simulation. Since
we want anyway ultimately to time-average all instantaneous stresses, Eq.
(129) may be replaced without loss of information by its ensemble average
$\langle\hat{\sigma}^{\mathrm{id}}_{\alpha\beta}\rangle=-P_{\mathrm{id}}\delta_{\alpha\beta}$
with $P_{\mathrm{id}}=Tn/V$ Allen and Tildesley (2017). The excess
contribution is given by the sum
$\hat{\sigma}^{\mathrm{ex}}_{\alpha\beta}=\sum_{l}w^{l}_{\alpha\beta}/V$ over
all interacting pairs $l$ of particles $a<b$. Let us denote by $\mbox{$\bf
r$}^{l}=\mbox{$\bf r$}^{b}-\mbox{$\bf r$}^{a}$ the vector from the position
$\mbox{$\bf r$}^{a}$ of $a$ to the position $\mbox{$\bf r$}^{b}$ of $b$, by
$r^{l}=|\mbox{$\bf r$}^{l}|$ its length and by $\hat{r}^{l}_{\alpha}$ a
component of its unit vector. The contribution $w^{l}_{\alpha\beta}$ of the
interaction $l$ to the excess stress is then Allen and Tildesley (2017)

$w^{l}_{\alpha\beta}\equiv
r_{l}u^{\prime}(r_{l})\hat{r}^{l}_{\alpha}\hat{r}^{l}_{\beta}$ (130)

with $u^{\prime}(r)$ being the first derivative of the pair potential.

#### Stress tensor fields in reciprocal space.

The corresponding stress tensor field $\hat{\sigma}_{\alpha\beta}(\mbox{$\bf
q$})=\hat{\sigma}^{\mathrm{id}}_{\alpha\beta}(\mbox{$\bf
q$})+\hat{\sigma}^{\mathrm{ex}}_{\alpha\beta}(\mbox{$\bf q$})$ in reciprocal
space may be directly obtained from the local momentum equation Landau and
Lifshitz (1959); Hansen and McDonald (2006)

$\frac{\partial}{\partial t}\hat{P}_{\alpha}(\mbox{$\bf
r$})=\frac{\partial}{\partial r_{\beta}}\hat{\sigma}_{\alpha\beta}(\mbox{$\bf
r$})$ (131)

as shown long ago Lutsko (1988); Hansen and McDonald (2006). Note that the
momentum density $\hat{P}_{\alpha}(\mbox{$\bf r$})$ is defined by Hansen and
McDonald (2006)

$\hat{P}_{\alpha}(\mbox{$\bf
r$})=\sum_{a=1}^{n}p^{a}_{\alpha}\delta(\mbox{$\bf r$}-\mbox{$\bf r$}^{a})$
(132)

in terms of the individual particle momenta and positions. Rewritten in
Fourier space Eq. (131) becomes

$\displaystyle iq_{\beta}\hat{\sigma}_{\alpha\beta}(\mbox{$\bf q$})$
$\displaystyle=$ $\displaystyle\frac{\partial}{\partial
t}\hat{P}_{\alpha}(\mbox{$\bf q$})$ $\displaystyle=$ $\displaystyle-
iq_{\beta}\frac{1}{V}\sum_{a=1}^{n}\frac{p^{a}_{\alpha}p^{a}_{\beta}}{m}\exp(-i\mbox{$\bf
q$}\cdot\mbox{$\bf r$}^{a})$ $\displaystyle+$
$\displaystyle\frac{1}{V}\sum_{a=1}^{n}\dot{p}^{a}_{\alpha}\exp(-i\mbox{$\bf
q$}\cdot\mbox{$\bf r$}^{a}).$ (134)

The first term in the previous relation gives the ideal stress contribution
Lutsko (1988); Hansen and McDonald (2006)

$\hat{\sigma}^{\mathrm{id}}_{\alpha\beta}(\mbox{$\bf
q$})=-\frac{1}{V}\sum_{a=1}^{n}\frac{p^{a}_{\alpha}p^{a}_{\beta}}{m}\exp(-i\mbox{$\bf
q$}\cdot\mbox{$\bf r$}^{a}).$ (135)

As before for the macroscopic ideal stress we may integrate out the momenta
$\\{\mbox{$\bf p$}\\}$ and replace Eq. (135) by

$\hat{\sigma}^{\mathrm{id}}_{\alpha\beta}(\mbox{$\bf
q$})=-\frac{T}{V}\delta_{\alpha\beta}\sum_{a=1}^{n}\exp(-i\mbox{$\bf
q$}\cdot\mbox{$\bf r$}^{a})$ (136)

which is used for our MC simulations. Using Newton’s second law for
$\dot{p}^{a}_{\alpha}$ and a total potential energy given by pair interactions
such as Eq. (125) the excess stress in reciprocal space is obtained from the
second term in Eq. (134) as shown in Refs. Lutsko (1988); Hansen and McDonald
(2006). This gives

$\hat{\sigma}^{\mathrm{ex}}_{\alpha\beta}(\mbox{$\bf
q$})=\frac{1}{V}\sum_{l}w^{l}_{\alpha\beta}\frac{e^{-i\mbox{$\bf
q$}\cdot\mbox{$\bf r$}^{b}}-e^{-i\mbox{$\bf q$}\cdot\mbox{$\bf
r$}^{a}}}{-i\mbox{$\bf q$}\cdot(\mbox{$\bf r$}^{b}-\mbox{$\bf r$}^{a})}.$
(137)

Let us define $\mbox{$\bf m$}^{l}=(\mbox{$\bf r$}^{a}+\mbox{$\bf r$}^{b})/2$
and $\mbox{$\bf h$}^{l}=(\mbox{$\bf r$}^{b}-\mbox{$\bf r$}^{a})/2$. As noted
by Lemaître Lemaître (2014) the previous relation may be rewritten more
conveniently as

$\hat{\sigma}^{\mathrm{ex}}_{\alpha\beta}(\mbox{$\bf
q$})=\frac{1}{V}\sum_{l}w^{l}_{\alpha\beta}\frac{\sin(\mbox{$\bf
q$}\cdot\mbox{$\bf h$}^{l})}{\mbox{$\bf q$}\cdot\mbox{$\bf
h$}^{l}}e^{-i\mbox{$\bf q$}\cdot\mbox{$\bf m$}^{l}}.$ (138)

We remark first that $\hat{\sigma}^{\mathrm{id}}_{\alpha\beta}(\mbox{$\bf
q$})$ and $\hat{\sigma}^{\mathrm{ex}}_{\alpha\beta}(\mbox{$\bf q$})$ reduce
properly for $\mbox{$\bf q$}\to{\bf 0}$ to the macroscopic stress
contributions noted above. Note also that the term $\sin(x)/x\to 1$ for
$x=\mbox{$\bf q$}\cdot\mbox{$\bf h$}^{l}\ll 1$, i.e. for wavelengths $2\pi/q$
larger than the typical interaction range of order unity in the present model.
In this limit Eq. (138) further simplifies to

$\hat{\sigma}^{\mathrm{ex}}_{\alpha\beta}(\mbox{$\bf
q$})\approx\frac{1}{V}\sum_{l}w^{l}_{\alpha\beta}e^{-i\mbox{$\bf
q$}\cdot\mbox{$\bf m$}^{l}},$ (139)

i.e. only the mean position $\mbox{$\bf m$}^{l}$ of two particles matters, not
their relative orientation. Since this study anyway focuses on the universal
large-wavelength limit, the approximation Eq. (139) should be sufficient.

#### Computation in reciprocal space.

Using Eq. (136) and Eq. (138) or Eq. (139) one may directly compute
$\hat{\sigma}_{\alpha\beta}(\mbox{$\bf q$})$ for a given configuration using a
discrete grid of linear length $n_{L}$ as shown in Fig. 2. Since we anyway
need the stress fields in reciprocal space to obtain the correlation
functions, Eq. (97), this is the most direct method and a useful exercise to
test less direct definitions and approximations. Unfortunately, this is
numerically not the most efficient procedure since it evolves $n_{L}^{d}\times
n\simeq V^{2}$ operations. It is computationally much faster to first obtain
the stress fields in real space with a number of operations of order $n\propto
V$ and then to FFT transform with a number of operations of order $V\log(V)$.

#### Stress tensor fields in real space.

We thus need to state the corresponding relations in real space. Inverse
Fourier transformation yields for Eq. (135)

$\hat{\sigma}^{\mathrm{id}}_{\alpha\beta}(\mbox{$\bf
r$})=-\sum_{a=1}^{n}\frac{p^{a}_{\alpha}p^{a}_{\beta}}{m}\ \delta(\mbox{$\bf
r$}-\mbox{$\bf r$}^{a})$ (140)

and for the preaveraged ideal stress, Eq. (136),

$\hat{\sigma}^{\mathrm{id}}_{\alpha\beta}(\mbox{$\bf
r$})=-T\delta_{\alpha\beta}\sum_{a=1}^{n}\delta(\mbox{$\bf r$}-\mbox{$\bf
r$}^{a}).$ (141)

One confirms using Eq. (92) that Eq. (137) becomes

$\displaystyle\hat{\sigma}^{\mathrm{ex}}_{\alpha\beta}(\mbox{$\bf r$})$
$\displaystyle=$ $\displaystyle\sum_{l}w_{\alpha\beta}^{l}I_{l}(\mbox{$\bf
r$})\mbox{ with }$ (142) $\displaystyle I_{l}(\mbox{$\bf r$})$
$\displaystyle=$ $\displaystyle\int_{0}^{1}\mathrm{d}s\ \delta(\mbox{$\bf
r$}-(\mbox{$\bf r$}^{a}+\mbox{$\bf r$}^{l}s))$

being a line integral between the two particle positions $\mbox{$\bf r$}^{a}$
and $\mbox{$\bf r$}^{b}$ of the interaction $l$. This relation has an old
history going at least back to the work by Kirkwood and Buff Kirkwood and Buff
(1949) for the layer-resolved slabs of microcells Allen and Tildesley (2017).
Various rediscoveries, reformulations and generalizations (e.g., for multibody
potentials and constraint dynamics) of the inverse Fourier transform Eq. (142)
are discussed elsewhere Tadmor and Miller (2011); Kirkwood and Buff (1949);
Irving and Kirkwood (1950); Noll (1955); Varnik et al. (2000); Goldhirsch and
Goldenberg (2002); Vanegas et al. (2014); Lemaître (2014); Behbahani et al.
(2021). If we use instead the approximation Eq. (139) the line integral is
replaced by

$I_{l}(\mbox{$\bf r$})\approx\delta(\mbox{$\bf r$}-\mbox{$\bf m$}^{l}).$ (143)

We emphasize that this coarse-grained expression is completely sufficient for
most applications. (The $\delta$-functions in Eq. (142) or Eq. (139) are
sometimes “blurred” using more general weighting distributions Lemaître
(2014).) Please note that for all stated expressions the volume averages are
identical to the macroscopic relation, i.e.

$\displaystyle\hat{\sigma}^{\mathrm{id}}_{\alpha\beta}$ $\displaystyle=$
$\displaystyle\frac{1}{V}\int\mathrm{d}\mbox{$\bf r$}\
\hat{\sigma}^{\mathrm{id}}_{\alpha\beta}(\mbox{$\bf
r$})=\hat{\sigma}^{\mathrm{id}}_{\alpha\beta}(\mbox{$\bf q$}={\bf 0})$ (144)
$\displaystyle\hat{\sigma}^{\mathrm{ex}}_{\alpha\beta}$ $\displaystyle=$
$\displaystyle\frac{1}{V}\int\mathrm{d}\mbox{$\bf r$}\
\hat{\sigma}^{\mathrm{ex}}_{\alpha\beta}(\mbox{$\bf
r$})=\hat{\sigma}^{\mathrm{ex}}_{\alpha\beta}(\mbox{$\bf q$}={\bf 0})$ (145)

hold exactly without approximation.

Figure 9: Rescaled ICF $\beta V\bar{c}_{\mathrm{N}}(q)$ comparing different
methods and $a_{\mathrm{grid}}=L/n_{L}$. The method “qA” uses Eq. (138)
directly in reciprocal space (solid line), the method “rA” the line
integration Eq. (142) in real space (open symbols), the method “rB” the
coarse-graining approximation Eq. (143) in real space (filled symbols). The
same result is obtained in all cases for sufficiently small $q$. The position
of the peak of $\beta V\bar{c}_{\mathrm{N}}(q)$ at $q\approx 6.5$ is very
similar to the position of the main peak of the structure factor $S(q)$ as
shown in Ref. Klochko et al. (2021).

#### Discrete grid implementations.

The numerical calculation of the above continuous space relations on a
discrete grid of linear length $n_{L}$ (cf. Fig. 2) is obvious for the ideal
contributions, Eq. (136) and Eq. (141), and also for the excess contribution
Eq. (138) in reciprocal space or the coarse-graining approximations Eq. (139)
or Eq. (143). Slightly less trivial is the implementation of the line integral
Eq. (142). Different variants exist for distributing the information of an
interaction $l$ on the grid as discussed in the literature Allen and Tildesley
(2017); Tadmor and Miller (2011). In our view this issue is not crucial since
the corresponding small wavelengths have no universal physical meaning being
due to an artificial computer model and are, moreover, readily renormalized
away as shown by Eq. (139). An important technical point is only that any
reasonable method must strictly obey Eq. (145), i.e. that the contributions to
all grid points from the line integral must be properly weighted. We use a
simple numerical rendering of Eq. (142) on the grid: for each of $k$