<figure><img src="image_1.png" /></figure>

HAL
open science

# Stress correlations and stress memory kernels in viscoelastic fluids

Niklas Grimm, Jörg Baschnagel, Alexander N. Semenov, Annette Zippelius,

Matthias Fuchs

► To cite this version:

Niklas Grimm, Jörg Baschnagel, Alexander N. Semenov, Annette Zippelius, Matthias Fuchs. Stress correlations and stress memory kernels in viscoelastic fluids. Soft Matter, 2025, 21 (21), pp.4256-4274.
10.1039/d5sm00156k . hal-05130155

HAL Id: hal-05130155

https://hal.science/hal-05130155v1

Submitted on 25 Jun 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Cite this: DOI: 00.0000/xxxxxxxxxx

Received Date
Accepted Date

DOI: 00.0000/xxxxxxxxxx

# Stress correlations and stress memory kernels in viscoelastic fluids

N. Grimm,ª J. Baschnagel,ᵇ A. N. Semenov,ᵇ A. Zippelius,ᶜ and Matthias Fuchsª*

We discuss the spatio-temporal correlations of stress fluctuations in viscoelastic fluids, including colloidal dispersions, polymer melts and glass-forming liquids. First, we relate the tensor of stress auto-correlation functions to its corresponding tensor of generalized Onsager transport kernels. The relation is valid at finite wavelengths and frequencies, and serves as basis for continuum mechanics approximations, where it leads to the viscosities in a fluid and to the elastic constants in a solid. Second, we consider the long wavelength limit and re-derive the far-field power-law stress correlations in fluid states. Two theoretical approaches are studied and give equivalent results. One approach is based on the familiar Zwanzig–Mori decomposition of the stress tensor, while the other approach decomposes the stress into a deterministic one due to the history of prior internal flow and a contribution due to “stress noise” that persists even if internal flow is suppressed. The general results are used to predict the distance dependence of the shear stress correlation function both for 2D and 3D systems. An encouraging agreement of these predictions with simulation data on 2D and 3D binary hard-sphere mixtures is observed.

## 1 Introduction

Viscoelastic fluids, particularly supercooled liquids, exhibit a complex interplay between liquid-like flow and solid-like elasticity¹,². When a liquid is cooled towards the glassy state, the arrest of structural rearrangements is marked by the gradual slowdown of stress fluctuations, which reflects the emergent elasticity of the system³. Far away from the glass transition at high temperature or low density, a (simple) liquid behaves as an ideal Newtonian fluid, with shear stress directly proportional to the shear rate, but as the temperature is reduced or the density increased, the fluid becomes more viscous and approaches the (metastable) glassy state. The time-window widens where the stress response resembles that of a Hookean solid with shear stress proportional to shear strain⁴. While this emergence of viscoelasticity on cooling is well known, a fundamental description is still being searched for⁵,⁶.

The role of memory effects in viscoelastic fluids is crucial for understanding this transition. Maxwell was the first to analyze the shear stress that a material develops in linear response to a small applied shear rate⁷. He introduced the “time of relaxation” τ which characterizes the memory in the response and becomes large in “viscous solids”, while being small in “mobile fluids”⁷.

Thereby, Maxwell could explain the crossover from viscous to elastic behavior. This scenario is well supported by rheological studies which report the frequency-dependent shear moduli, the loss part $G''(\omega)$ and the storage part $G'(\omega)^{1,2}$. They are connected to the Fourier-transform of the time-dependent shear relaxation modulus $G(t)$, which decays to zero in fluid states and exhibits a long-lived plateau $G(t \to \infty) \to \mu$ in deeply supercooled (“glassy”) states. This plateau defines the shear elastic constant (Lamé coefficient) $\mu$ of the glassy liquid, while the relaxation time $\tau$ determines its viscosity by $\eta = \mu\tau$. From the view point of modern linear response theory, it is clear that Maxwell⁷ would identify the shear relaxation modulus $G(t)$ with the correlation function of the shear stress $\langle\sigma_{xy}(0)\sigma_{xy}(t)\rangle$. Evans realized that this equality only holds on the macroscopic scale⁸,⁹, viz. for the volume averaged (shear) stress, equivalent to wavevector $q=0$. Inside the viscoelastic fluid (i.e. for $q \neq 0$), transport processes arise from the collective particle motion, implying that spatio-temporal stress fluctuations are coupled to inter alia space and time dependent density, energy, and momentum fluctuations. The fluctuations of these conserved variables become arbitrarily slow for increasing wavelengths¹⁰. Because of symmetry, they cannot couple to the macroscopic shear stress, but locally stress fluctuations couple to hydrodynamic transport processes¹¹,¹². Moreover, mechanical balance has to hold in steady states, which forces some tensor components of the stress to de-correlate even in frozen-in solid states¹³,¹⁴; in solids, strain replaces the velocity-gradient as conjugate variable to stress. It is necessary to account for the coupling of stress to the conserved and Goldstone fields, where

<footnote>ª University of Konstanz, D-78457 Konstanz, Germany</footnote>

<footnote>ᵇ Université de Strasbourg, CNRS, Institut Charles Sadron UPR 22, F-67000 Strasbourg, France</footnote>

<footnote>ᶜ University of Göttingen, D-37077 Göttingen, Germany</footnote>

<footnote>* Corresponding author: Matthias.Fuchs@uni-konstanz.de</footnote>
the latter exist only in solids¹⁵, in order to understand the development of viscoelastic memory in systems where the (local) structure rearranges slowly. This is also required in order to understand the spatial dependence of stress correlations¹⁶⁻¹⁸, and concerns complex fluids of macromolecular constituents as well as simple fluids when structural arrest sets in.

Considering the shear stress, a detailed description of the transition from a Newtonian fluid to a Hookean solid demanded to account for the coupling of stress and transverse momentum¹¹. In the viscoelastic fluid, mechanical memory effects result from the transverse momentum diffusion, leading to the development of long-range, angular dependent shear stress correlations. Considering the stresses around inclusions in force-free isotropic solids, Eshelby had determined their angular and distance behavior in the far field¹⁹. For the description of spatio-temporal stress correlations in viscoelastic fluids, two theoretical approaches were developed¹¹,¹²,¹⁴,¹⁵,²⁰⁻²⁵, whose relationship is currently unclear. The first aim of the present contribution is to show that both approaches give an equivalent framework for the complete tensor of stress correlation functions, bridging from the equilibrium liquid to the metastable glassy state. This will allow us to link the tensor of stress memory kernels, arising within an approach based on the projection operator formalism by Zwanzig and Mori¹¹,¹⁵,²⁰⁻²², to the tensor of generalized relaxation moduli underlying the second approach¹²,¹⁴,²³⁻²⁵. Considering the long-wavelength limit for the stress, we discuss a generalized Maxwell model that allows to describe the emergence of mechanical memory and long-range stress correlations¹¹,²⁰. Finally, we analyze computer simulation data of the (shear) stress correlations in two-dimensional (2D) and three-dimensional (3D) supercooled binary mixtures of hard spherical particles²⁶ to support the theoretical predictions.

Our paper is organized as follows. Section 2 summarizes the definitions of the density, velocity and stress fields as well as the fundamental conservation laws that link them. The section also defines the correlation function of the stress tensor and reviews the symmetries that result from translational invariances, isotropy, achirality, which we assume to hold for viscoelastic fluids. Sections 3 and 4 present and extend the two theoretical approaches to the stress tensor correlations¹¹,¹²,¹⁴,¹⁵,²⁰⁻²⁵ and, as a main result, establish their equivalence. Sections 5 and 6 then focus on the theoretical predictions for the far field stress correlations in real space and examine them from two perspectives, numerically by comparison with recent simulation results²⁶ for 2D and 3D binary mixtures of hard spheres (Sec. 5) and physically via a generalization of the Maxwell model (Sec. 6). The final section 7 discusses the results obtained and presents our conclusions.

## 2 Density, velocity and stress fields: balance laws and correlation functions

To describe the dynamics of amorphous systems at the collective level, basic quantities are the mass density field and the linear momentum density field. For the definition of these fields, Maier, Zippelius and Fuchs (MZF)¹¹,²⁰ and Semenov, Wittmer, Baschnagel (SWB) and coworkers¹²,¹⁴,²³,²⁵ employ different conventions. Although we adopt the SWB convention for the present

paper, we still summarize the different definitions and their conversion because this may serve as a translation catalog to help reading Refs.¹²,¹⁴,²³,²⁵ and¹¹,²⁰; this is done in Appendix A.

### 2.1 Fields encoding liquid dynamics

Let us consider a system of N classical particles in a volume V = L^d where L is the linear size of the system and d is the spatial dimension. The particle i (= 1, ..., N) has the mass m_i and is at position r_i(t) with velocity v_i(t) at time t. We define the mass density field at point r and time t as

$$ \rho(\mathbf{r},t) = \sum_{i=1}^{N} m_i \delta(\mathbf{r} - \mathbf{r}_i(t)). \quad (1) $$

Introducing the Fourier transform of a function f(r) (called variable or field) by

$$ f(\mathbf{q}) = \frac{1}{V} \int d^d \mathbf{r} f(\mathbf{r}) e^{-i\mathbf{q} \cdot \mathbf{r}}, \quad (2) $$

the transformed function f(q) has the same physical dimension as f(r). Here q denotes a (discrete) wavevector compatible with the boundary conditions (BC) of the system (periodic BC are assumed by default). Fourier transformation of Eq. (1) gives the mass density fluctuation

$$ \delta\rho(\mathbf{q},t) = \rho(\mathbf{q},t) - \langle\rho(\mathbf{q})\rangle = \frac{1}{V} \sum_{i=1}^{N} m_i e^{-i\mathbf{q} \cdot \mathbf{r}_i(t)} - \rho_0 \delta_{\mathbf{q},0}, \quad (3) $$

where ⟨...⟩ stands for the canonical ensemble average and ρ₀ is the mean mass density,

$$ \rho_0 = mn_0, \quad n_0 = \frac{N}{V}, \quad m = \frac{1}{N} \sum_{i=1}^{N} m_i, \quad (4) $$

with n₀ and m being respectively the mean number density and the average particle mass.

The rate of change of ρ(q,t) is defined by the momentum density (mass current density or simply current) J(q,t)

$$ \frac{\partial \rho(\mathbf{q},t)}{\partial t} = i \mathcal{L} \rho(\mathbf{q},t) = -i \mathbf{q} \cdot \mathbf{J}(\mathbf{q},t) = -i \mathbf{q} \cdot [\rho_0 \mathbf{v}(\mathbf{q},t)], \quad (5) $$

where L = −i∂_t is the classical Liouville operator which is Hamiltonian, L† = L⁹,²⁷ (here and in the following, linear operators are denoted by script letters), the current is defined by

$$ \mathbf{J}(\mathbf{q},t) = \frac{1}{V} \sum_{i=1}^{N} m_i \mathbf{v}_i(t) e^{-i\mathbf{q} \cdot \mathbf{r}_i(t)} \quad (6) $$

and v(q,t) is the Fourier transform of the velocity field v(r,t),

$$ \mathbf{v}(\mathbf{q},t) = \frac{1}{N} \sum_{i=1}^{N} \left( \frac{m_i}{m} \right) \mathbf{v}_i(t) e^{-i\mathbf{q} \cdot \mathbf{r}_i(t)}. \quad (7) $$

Upon taking the inverse Fourier transform

$$ f(\mathbf{r}) = \sum_{\mathbf{q}} f(\mathbf{q}) e^{i\mathbf{q} \cdot \mathbf{r}}, \quad (8) $$
we get

$$v(\mathbf{r}, t) = \frac{\mathbf{J}(\mathbf{r}, t)}{\rho_0}, \quad \mathbf{J}(\mathbf{r}, t) = \sum_i m_i v_i(t) \delta(\mathbf{r} - \mathbf{r}_i(t)) \quad (9)$$

so that Eq. (5) yields the familiar mass balance equation (mass conservation law) 27–30

$$\frac{\partial \rho(\mathbf{r}, t)}{\partial t} = - \frac{\partial J_{\alpha}(\mathbf{r}, t)}{\partial r_{\alpha}}, \quad (10)$$

where $J_\alpha$ is the $\alpha$th component of the momentum vector and here and in the following summation over repeated Greek indices is implied.

The rate of change of the current introduces the stress tensor $\sigma$ (with components $\sigma_{\alpha\beta}$) through the fundamental momentum equation (momentum balance) 12,25,28

$$\frac{\partial J_{\alpha}(\mathbf{r}, t)}{\partial t} = \frac{\partial \sigma_{\alpha\beta}(\mathbf{r}, t)}{\partial r_{\beta}}, \quad (11)$$

which reads in Fourier space

$$\frac{\partial J_{\alpha}(\mathbf{q}, t)}{\partial t} = i \mathcal{L} J_{\alpha}(\mathbf{q}, t) = i q_{\beta} \sigma_{\alpha\beta}(\mathbf{q}, t). \quad (12)$$

The stress tensor has a kinetic and a potential part. For pairwise interactions between all particles i and j, $U(r_{ij})$ with $r_{ij}$ being the norm of $\mathbf{r}_{ij} = \mathbf{r}_i - \mathbf{r}_j$, the microscopic definition of $\sigma_{\alpha\beta}(\mathbf{q})$ is given by the Irving–Kirkwood 29,31–34 (or Irving–Kirkwood–Noll 28) expression

$$\begin{aligned} \sigma_{\alpha\beta}(\mathbf{q}) = & \frac{1}{V} \sum_{i>j} \frac{r_{ij,\alpha} r_{ij,\beta}}{r_{ij}} \frac{dU(r_{ij})}{dr_{ij}} \frac{1 - e^{-i\mathbf{q} \cdot \mathbf{r}_{ij}}}{i\mathbf{q} \cdot \mathbf{r}_{ij}} e^{-i\mathbf{q} \cdot \mathbf{r}_i} \\ & - \frac{1}{V} \sum_{i=1}^{N} m_i v_{i,\alpha} v_{i,\beta} e^{-i\mathbf{q} \cdot \mathbf{r}_i}. \end{aligned} \quad (13)$$

This expression for the stress tensor is symmetric, $\sigma_{\alpha\beta} = \sigma_{\beta\alpha}$. The ensemble average of Eq. (13) gives 29,30

$$\langle \sigma_{\alpha\beta}(\mathbf{q}) \rangle = -p \delta_{\alpha\beta} \delta_{\mathbf{q},0}, \quad (14)$$

where $p$ is the equilibrium pressure of the system.

## 2.2 Correlations in amorphous systems: Consequences of translational invariance, isotropy and achirality

We consider amorphous systems at equilibrium. Such systems satisfy certain invariances on average: The systems are invariant with respect to translations in time and space and they are isotropic and achiral.

In this section, we address the implications of these invariances for the correlation functions of the stress tensor [cf. Eq. (13)] and density fields [cf. Eq. (1)]. Stress correlations are defined as

$$C_{\alpha\beta\gamma\delta}(\mathbf{r}, t) = \frac{1}{k_B T} \langle \delta\sigma_{\alpha\beta}(\mathbf{r} + \mathbf{r}', t') \delta\sigma_{\gamma\delta}(\mathbf{r}', t + t') \rangle, \quad (15)$$

where $\delta\sigma_{\alpha\beta}(\mathbf{r}, t) = \sigma_{\alpha\beta}(\mathbf{r}, t) - \langle \sigma_{\alpha\beta}(\mathbf{r}, t) \rangle$. Due to translational invariance, $\langle \sigma_{\alpha\beta}(\mathbf{r}, t) \rangle$ is independent of $\mathbf{r}$ and $t$ ($\langle \sigma_{\alpha\beta}(\mathbf{r}, t) \rangle = \langle \sigma_{\alpha\beta} \rangle$) and the correlation functions only depend on the displacement vector $\mathbf{r}$ between $\mathbf{r}'$ and $\mathbf{r} + \mathbf{r}'$ and on the time interval $t$ between $t'$ and $t' + t$. Translational invariance further implies that $\mathbf{r}'$ can be shifted by a constant vector $\mathbf{T}_r$ ($\mathbf{r}' \to \mathbf{r}' + \mathbf{T}_r$) and $t'$ by a constant time $T_t$ ($t' \to t' + T_t$). Choosing $\mathbf{T}_r = -\mathbf{r}$ and $T_t = -t$ one gets

$$C_{\alpha\beta\gamma\delta}(\mathbf{r}, t) = C_{\gamma\delta\alpha\beta}(-\mathbf{r}, -t). \quad (16)$$

Time inversion ($t \to -t$) changes the sign of the velocity of all particles ($v_i \to -v_i$) while their position remains the same ($\mathbf{r}_i \to \mathbf{r}_i$). As the Hamiltonian is an even function of $v_i$, neither the classical equations of motion nor the canonical average are modified. Moreover, $\sigma_{\alpha\beta}$ and $v_\alpha v_\beta$ are also even functions of $v_i$. Therefore, the rhs of Eq. (16) is invariant with respect to time reversal and we have

$$C_{\alpha\beta\gamma\delta}(\mathbf{r}, t) = C_{\gamma\delta\alpha\beta}(-\mathbf{r}, t). \quad (17)$$

Since we also assume the system to be achiral, it must be invariant with respect to space inversion ($\mathbf{r} \to -\mathbf{r}$), and Eq. (17) gives

$$C_{\alpha\beta\gamma\delta}(\mathbf{r}, t) = C_{\gamma\delta\alpha\beta}(\mathbf{r}, t). \quad (18)$$

We see that translational invariance, time-reversal symmetry and achirality lead to important index symmetries: The stress correlations $C_{\alpha\beta\gamma\delta}$ obey the major symmetry ($\alpha\beta \leftrightarrow \gamma\delta$), in addition to the minor symmetries, ($\alpha\beta \leftrightarrow \beta\alpha$) and ($\gamma\delta \leftrightarrow \delta\gamma$), because the stress tensor is symmetric (cf. Ref. 35 and Eq. (13) or Ref. 28 for a general discussion).

Similarly to [Eq. (15)], the density correlation function is defined as:

$$C^{\rho\rho}(\mathbf{r}, t) = \frac{1}{m\rho_0} \langle \delta\rho(\mathbf{r} + \mathbf{r}', t') \delta\rho(\mathbf{r}', t + t') \rangle, \quad (19)$$

where $\rho_0$ is the mean mass density defined in Eq. (4). $C^{\rho\rho}(\mathbf{r}, t)$ shows similar symmetries with respect to time and space inversions as the stress correlation function.

Note that in the above equations $\delta\rho(\mathbf{r}, t)$ and $\delta\sigma(\mathbf{r}, t)$ are locally defined fields. In the equilibrated amorphous systems we consider, static structural correlations (no time shift, $t=0$) are assumed to be local, in the sense that they have finite microscopic correlation length. This means that for any finite $t$ the correlation functions are local as well, as any signal spreads with a finite velocity, so that their Fourier transforms, defined according to Eq. (2), would lead to a vanishing correlation in the thermodynamic limit due to the factor $1/V$. This is physical reason for introducing the factor $V$ in the $q$-dependent correlation functions derived from real-space functions, Eqs. (15) and (19), to render the $q$-dependent functions finite. Accordingly, upon taking the Fourier transform of Eq. (19) we get

$$C^{\rho\rho}(\mathbf{q}, t) = \frac{V}{\rho_0 m} \langle \delta\rho(\mathbf{q}, t') \delta\rho^*(\mathbf{q}, t + t') \rangle = S(\mathbf{q}, t), \quad (20)$$

where (for $q \neq 0$)

$$S(\mathbf{q}, t) = \frac{1}{N} \sum_{i,j} \frac{m_i m_j}{m^2} \langle e^{-i\mathbf{q} \cdot (\mathbf{r}_i(t') - \mathbf{r}_j(t'))} \rangle \quad (21)$$

is the dynamical structure factor of the system. In a similar way,
we get from Eq. (15)

$$C_{\alpha\beta\gamma\delta}(\mathbf{q},t) = \frac{V}{k_B T} \left\langle \delta\sigma_{\alpha\beta}(\mathbf{q},t') \delta\sigma_{\gamma\delta}^*(\mathbf{q},t+t') \right\rangle, \quad (22)$$

where $\delta\sigma_{\alpha\beta}(\mathbf{q},t) = \sigma_{\alpha\beta}(\mathbf{q},t)$ because $\langle\sigma_{\alpha\beta}(\mathbf{q})\rangle = 0$ for $q \neq 0$. Moreover, as $\rho(\mathbf{r},t)$ and $\sigma_{\alpha\beta}(\mathbf{r},t)$ are real, Eqs. (17) and (18) imply that the wavevector dependent correlation functions are real as well,

$$S(\mathbf{q},t) = S^*(\mathbf{q},t) = S(-\mathbf{q},t), \quad (23)$$

$$C_{\alpha\beta\gamma\delta}(\mathbf{q},t) = C_{\alpha\beta\gamma\delta}^*(\mathbf{q},t) = C_{\alpha\beta\gamma\delta}(-\mathbf{q},t).$$

For isotropic systems the scalar density correlation function is independent of the direction of $\mathbf{q}$: $C^{\rho\rho}(\mathbf{q},t) = C^{\rho\rho}(q,t) = S(q,t)$. The tensorial stress correlation function, on the other hand, depends on the direction of $\mathbf{q}$. However, index symmetry of $C_{\alpha\beta\gamma\delta}(\mathbf{q},t)$ and isotropy of the system constrain this dependence; it occurs in two ways$^{20,23,25,36-38}$: first, by four (in $d=2$) or five (for $d \ge 3$) dynamic correlation functions that only depend on the norm $q$, and second, by factors involving the components $q_\alpha$ of $\mathbf{q}$. The latter factors reflect the coordinate frame chosen to represent $\mathbf{q}$ and contain no physical information. Since the coordinate frame is arbitrary, one can take advantage of this freedom to extract the (four or five) basic dynamic correlation functions—which contain the physical information about the spatio-temporal stress relaxation—by a suitable choice of the coordinate system in simulations$^{14,23,26,36,39}$.

### 3 Stress noise, deterministic stress, and suppression of flow: Relation between stress tensor correlations and tensor of generalized relaxation moduli

The approach by SWB and coworkers is based on the concept of stress noise and the constraint to suppress momentum density fluctuations$^{12,25}$. We summarize and extend this approach in the remainder of this section.

Following Refs.$^{12,25}$ the total stress fluctuations $\sigma_{\alpha\beta}(\mathbf{q},t)$ is written as a sum of two parts:

$$\sigma_{\alpha\beta}(\mathbf{q},t) = \sigma_{\alpha\beta}^{\mathrm{n}}(\mathbf{q},t) + \sigma_{\alpha\beta}^{\mathrm{d}}(\mathbf{q},t), \quad (24)$$

where $\sigma^d$ is defined as the flow-induced “deterministic” stress that is determined by the strain deformation history (cf. Eq. (27) below), while $\sigma^n$ is the “stress noise” collecting all contributions to stress other than those caused by deformation history. The stress noise is omnipresent and may be considered as a genuine stress fluctuation inherent in a system kept at zero strain.

The deformation history results from the local flow, i.e. the velocity field $v(\mathbf{r},t)$ that gives rise to a local strain rate defined as$^{25}$

$$\dot{\gamma}_{\alpha\beta}(\mathbf{r},t) = \frac{\partial}{\partial r_{\beta}} v_{\alpha}(\mathbf{r},t) = \sum_{\mathbf{q}} \dot{\gamma}_{\alpha\beta}(\mathbf{q},t) e^{i\mathbf{q}\cdot\mathbf{r}}, \quad (25)$$

where

$$\dot{\gamma}_{\alpha\beta}(\mathbf{q},t) \equiv \mathrm{iv}_{\alpha}(\mathbf{q},t)q_{\beta} = \frac{i}{\rho_0} J_{\alpha}(\mathbf{q},t)q_{\beta}. \quad (26)$$

In the time interval $dt$ the strain rate causes a weak deformation associated with the infinitesimal step strain $d\gamma_{\alpha\beta}(\mathbf{q},t) =$

$\dot{\gamma}_{\alpha\beta}(\mathbf{q},t)dt$. For weak deformations linear response theory is expected to hold. Inspired by the Boltzmann superposition principle$^{40}$ the flow-induced deterministic stress at time $t$ is thus written as a linear superposition of the stresses resulting from each individual step strain applied at time $t' \le t$:

$$\sigma_{\alpha\beta}^{\mathrm{d}}(\mathbf{q},t) = \int_{-\infty}^{t} E_{\alpha\beta\gamma\delta}(\mathbf{q},t-t') \mathrm{iv}_{\gamma}(\mathbf{q},t') q_{\delta} \, dt'. \quad (27)$$

The integral kernel introduces the tensor $E_{\alpha\beta\gamma\delta}(\mathbf{q},t)$ that embodies the space-resolved viscoelasticity of the equilibrium system: It describes how much stress remains at time $t-t'$ after a deformation at time $t'$. As viscoelasticity at the macroscopic ($q \to 0$) level is characterized by relaxation moduli (shear relaxation modulus, bulk compression modulus$^{41}$), the E-tensor is referred to as tensor of generalized (q-dependent) relaxation moduli$^{25}$.

Let us now consider a canonical ensemble of systems, which is subjected to the following deformation history. The ensemble is at thermal equilibrium at $t = -\infty$. For $-\infty < t < 0$ an external force is applied to all particles in each system, the effect of which is to suppress the flow at a given wavevector $\mathbf{q}$, i.e.

$$J(\mathbf{q},t) = 0 \quad (\text{for } -\infty < t < 0). \quad (28)$$

At $t=0$ the external force is switched off so that the flow occurs again for $t > 0$. Inserting this deformation history into Eq. (27) gives$^{25}$

$$\sigma_{\alpha\beta}^{\mathrm{d}}(\mathbf{q},t) = \sigma_{\alpha\beta}^{\mathrm{d}}(\mathbf{q}) + \int_{0}^{t} E_{\alpha\beta\gamma\delta}(\mathbf{q},t-t') \mathrm{iv}_{\gamma}(\mathbf{q},t) q_{\delta} \, dt'. \quad (29)$$

In addition to changing the lower integral limit ($-\infty \to 0$), the no-flow (or no-current) condition (28) adds the initial value $\sigma_{\alpha\beta}^d(\mathbf{q})$. This is a result of mass conservation. Why?

Due to mass conservation [cf. Eq. (5)], Eq. (28) implies that the mass density fluctuations is arrested at the value it has at $t = -\infty$: $\delta\rho(\mathbf{q},t = -\infty) = \delta\rho(\mathbf{q})$. Since the system is equilibrated at $t = -\infty$, $\delta\rho(\mathbf{q})$ is a canonical density fluctuation. This arrested density fluctuation may be thought of as being created by a small longitudinal strain applied at $t \to -\infty$, just before the flow is stopped. Mass conservation allows to find this strain: With $u_\alpha(\mathbf{q}) = \int v_\alpha(\mathbf{q},t)dt$ being the total displacement at wavevector $\mathbf{q}$, Eq. (5) implies

$$\varepsilon(\mathbf{q}) \equiv i u_{\alpha}(\mathbf{q}) q_{\alpha} = - \frac{\delta\rho(\mathbf{q})}{\rho_0}, \quad (30)$$

where $\varepsilon(\mathbf{q})$ denotes the trace of the strain tensor $\gamma_{\alpha\beta}(\mathbf{q}) = i u_{\alpha}(\mathbf{q})q_{\beta}$. The deformation $\varepsilon(\mathbf{q})$ at $t \to -\infty$ creates a time-dependent stress response of the system, which relaxes towards an equilibrium value, $\sigma_{\alpha\beta}^e(\mathbf{q})$, for times much longer than the finite relaxation time of the fluid system considered. The stress relaxes to $\sigma_{\alpha\beta}^e(\mathbf{q})$ (for a given $\varepsilon(\mathbf{q})$) because the system is conditionally ergodic: just one new invariant ($\delta\rho(\mathbf{q})$) is added to the common conservation laws, and we use the classical hypothesis that all allowed states, compatible with conservation laws, are accessible. We then identify $\sigma_{\alpha\beta}^e(\mathbf{q})$ with the time-averaged stress response $\overline{\sigma_{\alpha\beta}(\mathbf{q},t)}$ (time-averaged over times much longer than the finite relaxation time in the interval $-\infty < t < 0$) to the small
strain $\varepsilon(\mathbf{q})$,

$$ \sigma_{\alpha\beta}^{e}(\mathbf{q}) = \overline{\sigma_{\alpha\beta}(\mathbf{q}, t)} = \kappa_{\alpha\beta}(\mathbf{q})\varepsilon(\mathbf{q}), \quad (31) $$

where (cf. Eqs. (92) and (93) of Ref. 25)

$$ \kappa_{\alpha\beta}(\mathbf{q}) = (L_e(q) - M_e(q)) \hat{q}_\alpha \hat{q}_\beta + M_e(q) \delta_{\alpha\beta}. \quad (32) $$

Here $L_e(q)$ is the equilibrium longitudinal modulus that is related to the static structure factor by the generalized compressibility equation²⁹

$$ L_e(q) = k_B T n_0 / S(q), \quad (33) $$

and $M_e(q)$ is the transverse (mixed) modulus generalizing the Lamé coefficient $\lambda$. Due to the random character of the stress noise, its time-averaged value vanishes and we have $\overline{\sigma_{\alpha\beta}(\mathbf{q}, t)} = \overline{\sigma_{\alpha\beta}^d(\mathbf{q}, t)}$ from Eq. (24). This time-averaged deterministic stress remains constant as long as the flow is suppressed and therefore sets the initial value $\sigma^d(\mathbf{q}, 0) = \sigma^d(\mathbf{q})$ at time $t = 0$ in Eq. (29) when the no-flow constraint is released. Combining these arguments with Eq. (31), we get

$$ \sigma_{\alpha\beta}^{e}(\mathbf{q}) = \sigma_{\alpha\beta}^{d}(\mathbf{q}) = \kappa_{\alpha\beta}(\mathbf{q})\varepsilon(\mathbf{q}) = -\kappa_{\alpha\beta}\delta\rho(\mathbf{q})/\rho_0, \quad (34) $$

where Eq. (30) was used in the last step.

The stress response $\sigma_{\alpha\beta}^e(\mathbf{q})$ does not only determine the initial value, $\sigma_{\alpha\beta}^d(\mathbf{q}, t = 0) = \sigma_{\alpha\beta}^d(\mathbf{q})$, of the deterministic stress, but also serves to derive the equilibrium value $E_{\alpha\beta\gamma\delta}^e(\mathbf{q})$ of the **E**-tensor. This equilibrium value controls the stress response $\sigma_{\alpha\beta}^e(\mathbf{q})$ at long times after a weak impulsive (fast) deformation $\gamma_{\alpha\beta}(\mathbf{q})$ is applied to the system in the infinite past, while no further deformation at the given $\mathbf{q}$ is allowed for all times $-\infty < t < 0$. Then, we can write

$$ \sigma_{\alpha\beta}^{e}(\mathbf{q}) = E_{\alpha\beta\gamma\delta}^{e}(\mathbf{q})\gamma_{\gamma\delta}(\mathbf{q}). \quad (35) $$

Since Eq. (32) implies that

$$ \kappa_{\gamma\delta}q_{\delta} = L_e q_{\gamma}, \quad (36) $$

Eq. (35) with the choice

$$ E_{\alpha\beta\gamma\delta}^{e}(\mathbf{q}) = \kappa_{\alpha\beta}(\mathbf{q})\kappa_{\gamma\delta}(\mathbf{q})/L_e(q) \quad (37) $$

yields the same stress response $\sigma_{\alpha\beta}^e(\mathbf{q})$ as Eq. (34). The $E^e$-tensor is the tensor of static, equilibrium, elastic moduli for a system with no equilibrium shear modulus, as the fluid system we consider. This tensor will later be an important quantity that appears in the stress correlations.

It is possible to relate $E^e$ to the coupling between stress and density at equilibrium. Equation (34) indicates that the equilibrium density fluctuation $\delta\rho(\mathbf{q})$ determines the long-time stress response $\sigma_{\alpha\beta}^e(\mathbf{q})$. This stress response must coincide with a typical stress fluctuation, $\delta\sigma_{\alpha\beta}(\mathbf{q})$, inherent in the equilibrium system. Let us therefore define the equilibrium cross correlation between stress and density as

$$ C_{\alpha\beta}^{\sigma\rho}(\mathbf{q}) = -\frac{V}{k_B T \rho_0} \langle \delta\sigma_{\alpha\beta}(\mathbf{q})\delta\rho^*(\mathbf{q}) \rangle. \quad (38) $$

The prefactor ensures that $C_{\alpha\beta}^{\sigma\rho}(\mathbf{q})$ is dimensionless and the minus sign implies that it corresponds to a pressure-density coupling.

Inserting Eq. (34) we get

$$ C_{\alpha\beta}^{\sigma\rho}(\mathbf{q}) = \frac{S(q)}{k_B T n_0} \kappa_{\alpha\beta}, \quad (39) $$

where $S(q)$ is the static structure factor [cf. Eq. (20)]

$$ S(q) = S(q, t = 0) = \frac{V}{\rho_0 m} \langle \delta\rho(\mathbf{q})\delta\rho^*(\mathbf{q}) \rangle. \quad (40) $$

Since $M_e(q) \to L_e(q)$ in the hydrodynamic limit $q \to 0$²⁵, Eqs. (32) and (33) imply that

$$ C_{\alpha\beta}^{\sigma\rho}(\mathbf{q}) \to \delta_{\alpha\beta} \quad \text{for } q \to 0. \quad (41) $$

Recalling Eq. (33) we see that $E^e$ from Eq. (37) can also be written as

$$ E_{\alpha\beta\gamma\delta}^{e}(\mathbf{q}) = \frac{k_B T n_0}{S(q)} C_{\alpha\beta}^{\sigma\rho}(\mathbf{q}) C_{\gamma\delta}^{\sigma\rho}(\mathbf{q}). \quad (42) $$

This expression will prove useful for the comparison with the Zwanzig–Mori approach in Sec. 4.

The tensor $E^e$ also enters the expression for correlation function of the noise stress defined in analogy to Eq. (22) as

$$ C_{\alpha\beta\gamma\delta}^{n}(\mathbf{q}, t) = \frac{V}{k_B T} \langle \delta\sigma_{\alpha\beta}^{n}(\mathbf{q}, t') \delta\sigma_{\gamma\delta}^{n}(-\mathbf{q}, t+t') \rangle. \quad (43) $$

By arguing that the correlation properties of the stress noise are the same, irrespective of whether the no-flow constraint is applied ($t < 0$) or not ($t > 0$), the analysis of Ref. ²⁵ let us to posit (cf. Sect. 6.2 of²⁵)

$$ C_{\alpha\beta\gamma\delta}^{n}(\mathbf{q}, t) = E_{\alpha\beta\gamma\delta}(\mathbf{q}, |t|) - E_{\alpha\beta\gamma\delta}^{e}(\mathbf{q}). \quad (44) $$

Being a noise correlation function, $C_{\alpha\beta\gamma\delta}^n(\mathbf{q}, t)$ must vanish in the limit $t \to \infty$, which implies that $E_{\alpha\beta\gamma\delta}(\mathbf{q}, |t| \to \infty) = E_{\alpha\beta\gamma\delta}^e(\mathbf{q})$. Equation (44) goes beyond the linear-response-like introduction of the **E**-tensor by Eq. (27). In Eq. (27), the summation involved in the product $E_{\alpha\beta\gamma\delta}q_{\delta}$ implies that not the full **E**-tensor, but only its $(\alpha\beta\gamma)$ components can be obtained from such a linear-response-like relation. (We will encounter a similar situation in Sec. 4.3.) Equation (44), on the other hand, enables us to define all components of the **E**-tensor and moreover shows that **E** obeys the same minor $(\alpha\beta \leftrightarrow \beta\alpha)$ and major $(\alpha\beta \leftrightarrow \gamma\delta)$ symmetries as the stress tensor correlations.

Using Eq. (44) and that $\langle \sigma_{\alpha\beta}^d(\mathbf{q}, 0)\sigma_{\gamma\delta}^n(-\mathbf{q}, t) \rangle = 0$ for $t > 0$ due to the stochastic nature of the stress noise, Eqs. (29), (12), (34), (33) and (37) allow us to establish the following relation between the (full) stress correlation tensor (C-tensor) and the **E**-tensor (cf. Eq. (115) in Ref. ²⁵)

$$ \check{C}_{\alpha\beta\gamma\delta}(\mathbf{q}, s) = \check{E}_{\alpha\beta\gamma\delta}(\mathbf{q}, s) - \frac{q_\theta q_\eta}{\rho_0 s^2} \check{E}_{\alpha\beta\zeta\eta}(\mathbf{q}, s) \check{C}_{\theta\zeta\gamma\delta}(\mathbf{q}, s), \quad (45) $$

where

$$ \check{f}(s) = s \int_{0}^{\infty} dt f(t) e^{-st} \quad (46) $$

is the Laplace–Carson transform ⁴² (or “s-transform”¹⁴,²³,²⁵).

Equations (44) and (45) are the central relations of this section, which will be compared with the results from the Zwanzig–Mori
approach in Sec. 4.

## 4 Projection operator derivation of the viscosity kernel

For evaluating time-dependent correlations of arbitrary fields, such as the stress field, slow transport processes in the system need to be understood first. They are caused by conservation laws, since long-wavelength fluctuations of conserved variables require nonlocal balancing processes and exhibit time-scales which consequently grow with wavelength$^{43}$. For liquid systems the relevant conserved fields are the density and momentum field (energy conservation shall be neglected for simplicity). The smooth fluctuations of the conserved fields exhibit slow diffusional or propagative processes; in the case of momentum and density, they are transverse momentum diffusion and attenuated compression (viz. sound) waves$^{30}$. The quasi-static coupling of the stress to density and momentum causes thermal stress fluctuations to participate in these slow fluctuations of the conserved fields. The remaining stress fluctuations arise from local dissipation processes and should relax quickly in liquid systems. In viscoelastic fluids, however, also local structural processes become slow, and the separation of slow conserved dynamics from slow structural dynamics becomes essential for a correct description of arbitrary fluctuations. In the projection operator formalism this is anticipated as the local dynamics is encoded in memory kernels. For fluid systems (in the normal liquid state), the memory kernels can be replaced by transport coefficients, viz. the viscosities in the case of the stress. In viscoelastic systems, however, the time and (possibly also) space dependence of the memory kernels cannot be neglected. In the following, we address this problem for the stress within the projection operator formalism by Zwanzig and Mori$^{27,29,44,45}$ which is introduced next.

### 4.1 Zwanzig-Mori projection operator in fluids

Let $B(t)$ be a dynamic variable and $\delta B(t) = B(t) - \langle B \rangle$ its fluctuations about the equilibrium average $\langle B \rangle$. The variable $B$ shall depend on $t$ only through the microstate whose time evolution is generated by the Hamiltonian dynamics of the classical system. With the Liouville operator $\mathcal{L}$ the time dependence of $B$ can be written as $B(t) = e^{i\mathcal{L}t} B(0)$. The correlation function of two variables $A$ and $B$ takes the form

$$C^{AB}(t) = \langle \delta A^*(0) \delta B(t) \rangle = \langle \delta A^*(0) \mathcal{R}(t) \delta B(0) \rangle, \quad (47)$$

where the time-evolution operator

$$\mathcal{R}(t) = \exp(i\mathcal{L}t) \quad (48)$$

generates the classical dynamics of the system. The basic idea of the Zwanzig-Mori (ZM) projection operator formalism$^{27,29,44,45}$ is to decompose the dynamics into two contributions, one from the conserved variables and a “rest” that is free of the slow dynamics of the conservation laws. This decomposition is achieved by introducing a projection operator $\mathcal{P}$ on the conserved variables and defining the “rest” via the complementary operator $\mathcal{Q} = 1 - \mathcal{P}$

where $1$ is the identity operator $^*$.

Using the definitions summarized in Sec. 2.1, we write the projection operator $\mathcal{P}$ as (cf. Eq. (12) in Ref.$^{20}$)

$$\mathcal{P} = \frac{V}{m\rho_0 S(q)} \delta\rho(\mathbf{q}) \langle \delta\rho^*(\mathbf{q}) + \frac{mN}{k_B T} v_\gamma(\mathbf{q}) \rangle \langle v_\gamma^*(\mathbf{q}) \rangle, \quad (49)$$

where $\delta\rho(\mathbf{q}) = \delta\rho(\mathbf{q}, t=0)$ and $v_\gamma(\mathbf{q}) = v_\gamma(\mathbf{q}, t=0)$ are respectively the initial fluctuations of the collective density [cf. Eq. (3)] and of component $\gamma$ of the collective velocity [cf. Eq. (7)] $^\dagger$. The normalization constants are given by the static structure factor $S(q)$ of Eq. (40) and by

$$\langle v_{\alpha}^{*}(\mathbf{q}) v_{\beta}(\mathbf{q}) \rangle = \frac{k_B T}{N m} \delta_{\alpha\beta}, \quad (50)$$

with $k_B T/m$ being the square of the thermal velocity of a particle (for a system where all particles have the same mass $m$). Since $\langle \delta\rho^*(\mathbf{q}) v_\gamma(\mathbf{q}) \rangle = \langle v_\gamma^*(\mathbf{q}) \delta\rho(\mathbf{q}) \rangle = 0$ due to different parity of $\delta\rho$ and $v_\gamma$ upon time reversal, the normalization constants ensure the general relations

$$\mathcal{P}^2 = \mathcal{P}, \quad \mathcal{Q}^2 = \mathcal{Q}, \quad \mathcal{P}\mathcal{Q} = \mathcal{Q}\mathcal{P} = 0, \quad (51)$$

in addition to $\mathcal{P} = \mathcal{P}^\dagger$ and $\mathcal{Q} = \mathcal{Q}^\dagger$.

With the projectors $\mathcal{P}$ and $\mathcal{Q}$ the ZM decomposition starts by an exact rewriting of the time-evolution operator $\mathcal{R}(t)$. There are various ways of carrying out this rewriting$^{9,27,46}$. While prior work by MZF on stress correlations$^{20,22}$ employed the ZM decomposition proposed by Götze and Latz$^{27,46}$, we use here an alternative based on the Dyson decomposition$^9$. Both decompositions introduce the reduced time-evolution operator $\mathcal{R}'(t)$ defined as

$$\mathcal{R}'(t) = \mathcal{Q} \exp(i\mathcal{Q}\mathcal{L}\mathcal{Q}t) \mathcal{Q} = \mathcal{Q} \exp(i\mathcal{L}\mathcal{Q}t), \quad (52)$$

where the last equality follows from expanding the exponential operator as a power series and using $\mathcal{Q}^2 = \mathcal{Q}$. The reduced time-evolution operator generates that part of the dynamics of the system, which is not coupled to the conserved quantities. In order to ease the comparison with the prior MZF$^{20,22}$ calculations, the differences in the employed conventions are summarized in the Appendix A.

### 4.2 Zwanzig-Mori decomposition

In the ZM approach, the time evolutions generated by $\mathcal{R}(t)$ and $\mathcal{R}'(t)$ are linked by the fundamental Dyson relation$^9$:

$$\mathcal{R}(t) = \mathcal{R}(t)\mathcal{P} + \mathcal{R}'(t) + \int_0^t dt' \mathcal{R}(t')\mathcal{P}i\mathcal{L}\mathcal{R}'(t-t'). \quad (53)$$

<footnote>* Note that operators like $\mathcal{L}$, $\mathcal{P}$ or $\mathcal{Q}$ are linear transformations applied to the variables $(A, B, \dots)$ which are considered as functions in the phase space of the system, i.e. $A = A(\Gamma)$ where $\Gamma$ is the microstate representing the positions and momenta of all particles in the system. The variables can also depend on additional parameters, like $\mathbf{q}$ and $t$.</footnote>

<footnote>† The projection operator $\mathcal{P}$ depends on the wavevector $\mathbf{q}$ considered as a parameter. Let us denote the first term on the rhs of Eq. (49) by $\mathcal{P}_1$ and $\mathcal{P}_{1A}$ by $B$. Then, $B = (V/m\rho_0 S(q)) \delta\rho(\mathbf{q}, \Gamma) \int \delta\rho^*(\mathbf{q}, \Gamma') A(\Gamma')P(\Gamma')d\Gamma'$, where $\delta\rho(\mathbf{q}, \Gamma)$ is the density fluctuation in the microstate $\Gamma$ and $P(\Gamma)$ is the equilibrium (Gibbs) probability distribution of the ensemble of microstates in phase space.</footnote>
The full time evolution contains the conserved hydrodynamic modes which instantaneously couple into dynamics via the pro- jector P (first term on the rhs) and additionally also in a retarded way via the memory integral (last term on the rhs). Those two terms give the quasi-static coupling of the dynamics to the con- served fields and make the difference between R ( t ) and R $^{′}$( t ).

The Dyson relation can be derived by computing the time derivative of R ( t ) Q to obtain the differential equation

$$
\partial_t \mathcal{R}(t) Q = \mathcal{R}(t) \mathcal{P} i \mathcal{L} Q + \mathcal{R}(t) Q i \mathcal{L} Q,
$$

whose formal solution,

$$
\mathcal{R}(t) Q = Q e^{i \mathcal{L} Q t} + \int_0^t dt' \mathcal{R}(t') \mathcal{P} i \mathcal{L} Q e^{i \mathcal{L} Q(t-t')},
$$

agrees with Eq. (53) when using Eq. (52) and Q = 1 − P .

4.2.1 Decomposition of the deviatoric stress

Let us consider the part of the stress tensor that does not couple statically to the conserved variables, viz. Q σ . Its long wavelength limit will later be identified with the deviatoric stress of contin- uum mechanics (cf. Eq. (109) below). Thus, we refer to

$$
\boldsymbol{\sigma}^{\text{dev}} \equiv \mathcal{Q}\boldsymbol{\sigma} \qquad (54)
$$

as the deviatoric stress. Applying Eq. (53) to Q σ and inserting the momentum equation (12) we have

$$
\begin{equation}
\begin{split}
\sigma_{\alpha\beta}^{\text{dev}}(\mathbf{q}, t) \equiv \mathcal{R}(t) Q \sigma_{\alpha\beta}(\mathbf{q}) &= \mathcal{R}'(t) \sigma_{\alpha\beta}(\mathbf{q}) + \\
&\quad \int_0^t dt' i v_{\zeta}(\mathbf{q}, t') q_{\eta} M_{\zeta\eta\alpha\beta}(\mathbf{q}, t-t'),
\end{split}
\tag{55}
\end{equation}
$$

where we introduced the memory kernel M defined as

$$
M_{\zeta\eta\alpha\beta}(\mathbf{q}, t) \equiv \frac{V}{k_B T} \langle \boldsymbol{\sigma}_{\zeta\eta}^*(\mathbf{q}) \mathcal{R}'(t) \boldsymbol{\sigma}_{\alpha\beta}(\mathbf{q}) \rangle, \quad (56)
$$

The time derivative of this equation would lead to a familiar gen- eralized Langevin equation $^{9,29}$, where the fluctuation of a vari- able contains a noise term and a non-Markovian contribution given by the memory kernel that is the correlation function of the noise. In the following, however, the precursor (55) to the Langevin equation will prove more useful. From Eq. (55) we see that the memory kernel M determines the dependence of the devi- atoric stress Q σ on the flow history, viz. on v ( q , t ). The M -tensor obeys the same symmetries as the stress correlation tensor (cf. Sec. 2.2), since the reduced dynamics is also time reversible $^{20,22}$.

As shown in the next section, M is related to the tensor of gen- eralized relaxation moduli, E . The Fourier transform of M ( q , t ) with respect to time gives the tensor of viscosities in the hydro- dynamic limit. Thus, M will also be referred to as generalized viscosity tensor (or viscosity kernel) in the following.

The relation of M to stress correlations can further be worked out by first multiplying Eq. (55) with σ ∗ from the left and then

taking the thermal average. This gives

$$
\begin{align}
& \frac{V}{k_B T} \langle \sigma_{\gamma\delta}^*(\mathbf{q}) \mathcal{R}(t) Q \sigma_{\alpha\beta}(\mathbf{q}) \rangle = M_{\gamma\delta\alpha\beta}(\mathbf{q}, t) \tag{57} \\
& \quad + \frac{V}{k_B T} \int_0^t dt' \langle \sigma_{\gamma\delta}^*(\mathbf{q}) \mathcal{R}(t') v_{\zeta}(\mathbf{q}) \rangle i q_{\eta} M_{\zeta\eta\alpha\beta}(\mathbf{q}, t-t') \nonumber
\end{align}
$$

The cross-correlations current and stress can be related to the stress correlation function via momentum conservation. Tech- nically, this enters through integration of the time derivative of R ( t ): R ( t ) = 1 + ∫ t 0 dt $^{′}$R ( t $^{′}$) i L . As stress and momentum have dif- ferent time parity, their stationary (equal time) overlap vanishes, and the result reads

$$
\frac{V}{k_B T} \langle \sigma_{\gamma\delta}^*(\mathbf{q}) \mathcal{R}(t) v_{\zeta}(\mathbf{q}) \rangle = \frac{1}{\rho_0} \int_0^t dt' C_{\gamma\delta\theta\zeta}(\mathbf{q}, t') i q_{\theta} \quad (58)
$$

with C ( q , t ) being the stress correlation tensor [Eq. (22)]. Defin- ing the Laplace transform of a function f ( t ) as

$$
\tilde{f}(s) = \int_{0}^{\infty} dt f(t) e^{-st}, \qquad (59)
$$

Eq. (57) leads to the central formula of the ZM approach

$$
\begin{equation}
\begin{split}
\frac{V}{k_B T} \langle \sigma_{\gamma\delta}^*(\mathbf{q}) \tilde{\mathcal{R}}(s) Q \sigma_{\alpha\beta}(\mathbf{q}) \rangle ={}& \tilde{M}_{\gamma\delta\alpha\beta}(\mathbf{q}, s) \\
& - \frac{q_{\theta}q_{\eta}}{s\rho_0} \tilde{C}_{\gamma\delta\theta\zeta}(\mathbf{q}, s) \tilde{M}_{\zeta\eta\alpha\beta}(\mathbf{q}, s).
\end{split}
\tag{60}
\end{equation}
$$

Let us note in passing that overdamping the dynamics by hand, viz. setting m∂$_{t}$ v → ζ$_{0}$ v with friction coefficient ζ$_{0}$ , amounts to replacing ms → ζ$_{0}$ in Eq. (60) and leads to the equivalent result for Brownian dynamics $^{21}$.

As stated at the beginning of this section, Q σ will be identified as the deviatoric part of the stress [cf. Eq. (54)], which needs to be understood in hydrodynamic and elasticity theories. To move on to these coarse-grained descriptions valid at large length and time scales, Eq. (60) is the main starting point derived within the ZM approach.

4.2.2 Stress-density coupling and total stress

Equation (60) is still incomplete because not the full stress, but the deviatoric stress Q σ appears on the rhs. How about the com- plement, the P -projected stress,

$$
\sigma^P \equiv P\sigma, \qquad (61)
$$

and its time evolution? Using Eq. (49) and considering that 〈 v ∗ $_{γ}$σ$_{αβ}$ 〉 = 0 due to different parity under time reversal, σ $^{P}$( q , t ) is given by

$$
\begin{align}
\sigma_{\alpha\beta}^{\mathrm{P}}(\mathbf{q}, t) &\equiv \mathcal{R}(t) \mathcal{P} \sigma_{\alpha\beta}(\mathbf{q}) = \frac{V}{m\rho_0 S(q)} \langle \delta\rho^*(\mathbf{q}) \sigma_{\alpha\beta}(\mathbf{q}) \rangle \mathcal{R}(t) \delta\rho(\mathbf{q}) \nonumber \\
&= -\frac{k_B T}{mS(q)} C_{\alpha\beta}^{\sigma\rho}(\mathbf{q}) \delta\rho(\mathbf{q}, t), \tag{62}
\end{align}
$$

where Eq. (38) was used to obtain the last line. As it is de- sirable to find a closed equation for the stress fluctuations, two time derivatives are applied next. Using the conservation laws,
Eqs. (10) and (12), this leads to

$$ \partial_t^2 \mathcal{R}(t) \mathcal{P} \sigma_{\alpha\beta}(\mathbf{q}) = -\frac{k_B T q_\theta q_\eta}{mS(q)} C_{\alpha\beta}^{\sigma\rho}(\mathbf{q}) \sigma_{\theta\eta}(\mathbf{q}, t). \quad (63) $$

Both Eqs. (60) and (63) could be combined via Laplace transformation. Yet, there is a simpler way to determine the stress correlation function. On multiplying Eq. (62) by $\sigma_{\gamma\delta}^*(\mathbf{q})$ from the left followed by thermal averaging, one arrives at

$$ \frac{V}{k_B T} \langle \sigma_{\gamma\delta}^*(\mathbf{q}) \mathcal{R}(t) \mathcal{P} \sigma_{\alpha\beta}(\mathbf{q}) \rangle = \frac{k_B T n_0}{S(q)} C_{\alpha\beta}^{\sigma\rho}(\mathbf{q}) C_{\gamma\delta}^{\sigma\rho}(\mathbf{q}, t), \quad (64) $$

where Eq. (38) was used to express the rhs in terms of the cross-correlator $C_{\gamma\delta}^{\sigma\rho}(\mathbf{q}, t)$ between density and stress,

$$ C_{\gamma\delta}^{\sigma\rho}(\mathbf{q}, t) = -\frac{V}{k_B T \rho_0} \langle \delta\sigma_{\gamma\delta}(\mathbf{q}) \delta\rho^*(\mathbf{q}, t) \rangle, \quad (65) $$

which has the initial value $C_{\gamma\delta}^{\sigma\rho}(\mathbf{q}, t=0) = C_{\gamma\delta}^{\sigma\rho}(\mathbf{q})$. Next, doing a Laplace transform of Eq. (64) and combining the result with Eq. (60) we find

$$ \begin{aligned} \tilde{C}_{\gamma\delta\alpha\beta}(\mathbf{q}, s) = & \frac{k_B T n_0}{S(q)} C_{\alpha\beta}^{\sigma\rho}(\mathbf{q}) \tilde{C}_{\gamma\delta}^{\sigma\rho}(\mathbf{q}, s) + \tilde{M}_{\gamma\delta\alpha\beta}(\mathbf{q}, s) \\ & - \frac{q_\theta q_\eta}{s \rho_0} \tilde{C}_{\gamma\delta\theta\zeta}(\mathbf{q}, s) \tilde{M}_{\zeta\eta\alpha\beta}(\mathbf{q}, s). \end{aligned} \quad (66) $$

Equation (66) still contains the cross-correlator between stress and density as an unknown. To calculate this correlator we first combine mass and momentum conservation [Eqs. (10), (12)] to get

$$ \partial_t^2 \delta\rho(\mathbf{q}, t) = q_\theta q_\zeta \sigma_{\theta\zeta}(\mathbf{q}, t), $$

then insert the Laplace transform of this equation into Eq. (65), and finally take into account that the equilibrium average of the correlation between $\partial_t \delta\rho(\mathbf{q}, t)|_{t=0}$ ($=-i\epsilon_E J_E(\mathbf{q}, t=0)$ and $\sigma_{\gamma\delta}(\mathbf{q}, 0)$ vanishes due to different parity under time reversal. This leads to

$$ \tilde{C}_{\gamma\delta}^{\sigma\rho}(\mathbf{q}, s) = \frac{1}{s} C_{\gamma\delta}^{\sigma\rho}(\mathbf{q}) - \frac{q_\theta q_\zeta}{\rho_0 s^2} \tilde{C}_{\gamma\delta\theta\zeta}(\mathbf{q}, s), \quad (67) $$

where the major index symmetry of $C_{\gamma\delta\theta\zeta}$ was also used. Substituting Eq. (67) into Eq. (66), recalling Eq. (42),

$$ E_{\alpha\beta\gamma\delta}^{e}(\mathbf{q}) = \frac{k_B T n_0}{S(q)} C_{\alpha\beta}^{\sigma\rho}(\mathbf{q}) C_{\gamma\delta}^{\sigma\rho}(\mathbf{q}), $$

and using that Eqs. (36), (33) and (39) imply

$$ q_{\eta} C_{\zeta\eta}^{\sigma\rho}(\mathbf{q}) = q_{\zeta}, \quad (68) $$

we finally get for the full stress correlator

$$ \begin{aligned} \tilde{C}_{\gamma\delta\alpha\beta}(\mathbf{q}, s) = & \tilde{M}_{\gamma\delta\alpha\beta}(\mathbf{q}, s) + \frac{1}{s} E_{\gamma\delta\alpha\beta}^{e}(\mathbf{q}) \\ & - \frac{q_\theta q_\eta}{s \rho_0} \tilde{C}_{\gamma\delta\theta\zeta}(\mathbf{q}, s) \left[ \tilde{M}_{\zeta\eta\alpha\beta}(\mathbf{q}, s) + \frac{1}{s} E_{\zeta\eta\alpha\beta}^{e}(\mathbf{q}) \right]. \end{aligned} \quad (69) $$

Because of the minor and major index symmetries of the **C**, **M** and **E**^e tensors, this result from the ZM approach agrees with the

result for the stress correlator obtained in Sec. 3 [cf. Eq. (45)], if we relate the memory kernel **M** with the tensor **E** of generalized relaxation moduli in the following way

$$ E_{\alpha\beta\gamma\delta}(\mathbf{q}, t) = M_{\alpha\beta\gamma\delta}(\mathbf{q}, t) + E_{\alpha\beta\gamma\delta}^{e}(\mathbf{q}) \quad (t > 0), \quad (70) $$

where $E_{\alpha\beta\gamma\delta}^{e}(\mathbf{q})$ can also be written in terms of the projector **P** as

$$ E_{\alpha\beta\gamma\delta}^{e}(\mathbf{q}) = \frac{V}{k_B T} \langle \sigma_{\alpha\beta}^{*}(\mathbf{q}) \mathcal{P} \sigma_{\gamma\delta}(\mathbf{q}) \rangle, \quad (71) $$

due to Eqs. (38) and (49). Comparing Eq. (70) with Eq. (44), we can further conclude that the memory kernel **M** coincides with the correlation function of the stress noise for $t > 0$:

$$ M_{\alpha\beta\gamma\delta}(\mathbf{q}, t) = C_{\alpha\beta\gamma\delta}^{\mathrm{n}}(\mathbf{q}, t). \quad (72) $$

Equations (69), (70) and (72) are one of the central results of our study: They establish the equivalence of the approach developed by SWB and coworkers$^{12,25}$ and the ZM approach developed by MZF$^{11,20,22}$ and further extended here. Since stress is not conserved in a fluid, the correlation of stress fluctuations **C**(q,t) is determined by local stress kernels, either **M** or **E**, which are connected via Eq. (70). Macroscopically at q = 0, correlation function and memory kernel agree, as Maxwell anticipated. Yet at finite wavevectors, the conservation laws dictate that the stress correlation function becomes anisotropic and deviates from the memory kernel, Eq. (69).

### 4.2.3 Comparison of the fluctuating stress decompositions of Sections 3 and 4

The approaches of Secs. 3 and 4 both employ a splitting of the total stress into two parts: stress noise $\sigma^n$ and deterministic stress $\sigma^d$ in Sec. 3 [cf. Eq. (24)],

$$ \sigma_{\alpha\beta}(\mathbf{q}, t) = \sigma_{\alpha\beta}^{n}(\mathbf{q}, t) + \sigma_{\alpha\beta}^{d}(\mathbf{q}, t), \quad (73) $$

and deviatoric stress $\sigma^{\mathrm{dev}} = Q\sigma$ and density-related $P$-projected stress $\sigma^P = P\sigma$ in the present section 4,

$$ \sigma_{\alpha\beta}(\mathbf{q}, t) = \mathcal{R}(t) Q \sigma_{\alpha\beta}(\mathbf{q}) + \mathcal{R}(t) P \sigma_{\alpha\beta}(\mathbf{q}). \quad (74) $$

This raises the question of whether both decompositions can be linked. To address this point we first notice that insertion of Eqs. (55) and (62) into Eq. (74) yields

$$ \begin{aligned} \sigma_{\alpha\beta}(\mathbf{q}, t) = & \mathcal{R}'(t) \sigma_{\alpha\beta}(\mathbf{q}) - \frac{k_B T}{mS(q)} C_{\alpha\beta}^{\sigma\rho}(\mathbf{q}) \delta\rho(\mathbf{q}, t) \\ & + \int_0^t dt' M_{\alpha\beta\gamma\delta}(\mathbf{q}, t-t') i v_{\gamma}(\mathbf{q}, t') q_{\delta}, \end{aligned} \quad (75) $$

where the major index symmetry $(\alpha\beta \leftrightarrow \gamma\delta)$ of **M** was utilized. The term involving $\delta\rho(\mathbf{q}, t)$ can be rewritten by integrating Eq. (5) and using Eq. (68),

$$ \delta\rho(\mathbf{q}, t) = \delta\rho(\mathbf{q}) - \rho_0 \int_0^t dt' i v_{\gamma}(\mathbf{q}, t') q_{\delta} C_{\gamma\delta}^{\sigma\rho}(\mathbf{q}), $$
which allows to write the P-projected stress as

$$
\sigma_{\alpha\beta}^{\mathrm{P}}(\mathbf{q}, t) = -\frac{k_{\mathrm{B}}T}{mS(q)} C_{\alpha\beta}^{\sigma\rho}(\mathbf{q}) \delta\rho(\mathbf{q}, t) - \frac{k_{\mathrm{B}}T}{mS(q)} C_{\alpha\beta}^{\sigma\rho}(\mathbf{q}) \delta\rho(\mathbf{q}) \\
+ E_{\alpha\beta\gamma\delta}^{\mathrm{e}}(\mathbf{q}) \int_{0}^{t} \mathrm{d}t' \mathrm{i}v_{\gamma}(\mathbf{q}, t') q_{\delta}, \qquad (76)
$$

if Eq. (42) was also taken into account. Since Eqs. (39) and (34)
imply that (in the case of no flow at $t < 0$)

$$
\sigma_{\alpha\beta}^{\mathrm{P}}(\mathbf{q}, t=0) = -\frac{k_B T}{mS(q)} C_{\alpha\beta}^{\sigma\rho}(\mathbf{q}) \delta\rho(\mathbf{q}) = \sigma_{\alpha\beta}^{\mathrm{d}}(\mathbf{q}), \quad (77)
$$

we can write Eq. (75) with the help of Eq. (70) as

$$
\begin{align}
\sigma_{\alpha\beta}(\mathbf{q}, t) &= \mathcal{R}'(t)\sigma_{\alpha\beta}(\mathbf{q}) + \nonumber \\
&\quad \sigma_{\alpha\beta}^{\mathrm{d}}(\mathbf{q}) + \int_{0}^{t} \mathrm{d}t' E_{\alpha\beta\gamma\delta}(\mathbf{q}, t-t')\mathrm{i}v_{\gamma}(\mathbf{q}, t')q_{\delta}. \tag{78}
\end{align}
$$

As the second line of Eq. (78) is equal to $\sigma_{\alpha\beta}^d(\mathbf{q}, t)$ due to Eq. (29), this shows that Eq. (73) and Eq. (74) are equivalent, if one identifies for $t \ge 0$

$$
\sigma_{\alpha\beta}^{n}(\mathbf{q}, t) = \mathcal{R}'(t) \sigma_{\alpha\beta}(\mathbf{q}) = \mathcal{R}'(t) Q \sigma_{\alpha\beta}(\mathbf{q}), \quad (79)
$$

$$
\begin{equation}
\begin{split}
\sigma_{\alpha\beta}^{\mathrm{d}}(\mathbf{q}, t) = & -\frac{k_{\mathrm{B}}T}{mS(q)} C_{\alpha\beta}^{\sigma\rho}(\mathbf{q}) \delta\rho(\mathbf{q}, t) \\
& + \int_{0}^{t} \mathrm{d}t' M_{\alpha\beta\gamma\delta}(\mathbf{q}, t-t') \mathrm{i}v_{\gamma}(\mathbf{q}, t') q_{\delta},
\end{split}
\tag{80}
\end{equation}
$$

where Eq. (52) and $\mathcal{Q}^2 = \mathcal{Q}$ were employed in Eq. (79). Note that the first term in the rhs of Eq. (80) is equal to $\boldsymbol{\sigma}^{\mathrm{P}}(\mathbf{q}, t)$.

Equations (79) and (80) are another central result of our study. They provide the correspondence of the stress decompositions of both approaches and also show that in general $\boldsymbol{\sigma}^n(\mathbf{q}, t) \neq \boldsymbol{\sigma}^{\mathrm{dev}}(\mathbf{q}, t)$ and $\boldsymbol{\sigma}^d(\mathbf{q}, t) \neq \boldsymbol{\sigma}^\mathrm{P}(\mathbf{q}, t)$, except for $t=0$ where $\mathcal{R}'(0)=\mathcal{Q}$ [cf. Eq. (52)] and Eq. (77) holds. For $t>0$ Eq. (79) implies that the classical Newtonian time evolution of $\boldsymbol{\sigma}^n(\mathbf{q}, t)$ corresponds to the reduced deviatoric stress with a projected dynamics generated by the reduced time-evolution operator $\mathcal{R}'(t)$. As far as $\boldsymbol{\sigma}^d(\mathbf{q}, t)$ is concerned, Eq. (80) implies that its time dependence is not only determined by density fluctuations—and so by $\boldsymbol{\sigma}^\mathrm{P}(\mathbf{q}, t)$—but also by the flow history, i.e. by the full stress due to Eq. (12). This discussion reveals the nontrivial nature of the equality between the memory kernel and the stress noise correlator [cf. Eq. (72)] at the ensemble average level.

4.3 Memory kernels and response functions

In order to interpret the memory kernel **M** a linear response con-
sideration proves useful. It starts from the fluctuation-dissipation
theorem which links the correlation function to a response func-
tion. The latter can be considered as a functional derivative with
respect to an external field. Martin⁴⁷ or respectively Kadanoff
and Martin⁴⁸ suggested that it is also possible to express the re-
sponse function in terms of functional derivatives with respect to
internal fields. Translating the ZM calculation in this style results

in an interpretation of the ZM memory kernel.

4.3.1 Response to an external field

Let us apply a small perturbative external velocity, $v^{\mathrm{ext}}(\mathbf{r}, t)$, to the system. The applied perturbation affects the dynamics of the system and consequently the Liouville equation. Let $v^{\mathrm{ext}}(\mathbf{q}, t)$ be the Fourier mode of the external velocity. Linear response theory³⁰ gives for any variable $A$ the expectation value of its change in linear order in the external perturbation, $\Delta A = \langle A(\mathbf{q}, t) \rangle^{\mathrm{lr}}$, via a Kubo relation⁴⁴ (for an explicit calculation see Ref. ²⁰)

$$
\langle A(\mathbf{q}, t) \rangle^{\mathrm{lr}} = \frac{V}{k_B T} \int_{-\infty}^{t} \mathrm{d}t' \langle \sigma_{\gamma\delta}^{*}(\mathbf{q}) A(\mathbf{q}) \rangle \mathrm{i} q_{\gamma} v_{\delta}^{\mathrm{ext}}(\mathbf{q}, t'), \quad (81)
$$

where 〈 ... 〉<sup>lr</sup> denotes the average over a time dependent distribu-
tion function in the linear approximation, while the average on
the rhs is performed in the unperturbed equilibrium system. ‡
This result is a manifestation of the fluctuation-dissipation theo-
rem. Translational invariance dictates the occurrence of only one
wavevector. Equation (81) motivates the definition of a response
function via a functional derivative. Choosing the stress as vari-
able A in Eq. (81) implies that the integral kernel becomes the
stress correlation function

$$
\langle \sigma_{\alpha\beta}(\mathbf{q}, t) \rangle^{\mathrm{lr}} = \int_{-\infty}^{+\infty} \mathrm{d}t' C_{\gamma\delta\alpha\beta}(\mathbf{q}, t-t') \Theta(t-t') \mathrm{i} q_{\gamma} v_{\delta}^{\mathrm{ext}}(\mathbf{q}, t'), \quad (82)
$$

where $\Theta(t)$ is the Heaviside function. The stress correlation func-
tion may thus be identified with an *after-effect* function<sup>30</sup>, giving
the stress response at time $t$ to an external velocity fluctuations
at time $t'$. The external perturbation can always be considered
as a superposition of monochromatic plane waves, which factor-
izes in linear response<sup>30</sup>. It is therefore sufficient to consider
a single frequency $\omega$. This motivates a Fourier-transformation
$\mathcal{F}[f(t)](\omega) = \int_{-\infty}^{+\infty} dt \, e^{-i\omega t} f(t)$, leading to

$$
\begin{align}
\frac{\partial \langle \sigma_{\alpha\beta}(\mathbf{q}, \omega) \rangle^{\mathrm{lr}}}{\partial v_{\delta}^{\mathrm{ext}}(\mathbf{q}, \omega)} &= \frac{V}{k_{\mathrm{B}}T} \langle \sigma_{\gamma\delta}^{*}(\mathbf{q}) \tilde{\mathcal{R}}(s = i\omega) \sigma_{\alpha\beta}(\mathbf{q}) \rangle \mathrm{i} q_{\gamma} \nonumber \\
&= \tilde{C}_{\gamma\delta\alpha\beta}(\mathbf{q}, s = i\omega) \mathrm{i} q_{\gamma}, \tag{83}
\end{align}
$$

where the Laplace transform, respectively the half-sided Fourier-
transform, appears. Note the partial, rather than functional,
derivative in the frequency domain.

The external velocity field does not only cause a stress but also
an internal velocity. In linear order, the resulting mean internal
velocity reads $v(\mathbf{r}, t) = (v(\mathbf{r}, t))^{\mathrm{lr}}$. The general Kubo relation (81)
can be used with variable change $A \rightarrow v$, which leads to the linear

<footnote>‡ Note that a possible finite average $\langle A \rangle$ would not matter in the following.</footnote>
response relation:

$$
\begin{align}
\frac{\partial v_{\alpha}(\mathbf{q}, \omega)}{\partial v_{\delta}^{\text{ext}}(\mathbf{q}, \omega)} &= \frac{i q_{\beta} i q_{\gamma} V}{i \omega \rho_{0} k_{\text{B}} T} \langle \sigma_{\gamma \delta}^{*}(\mathbf{q}) \tilde{\mathcal{R}}(s = i \omega) \sigma_{\alpha \beta}(\mathbf{q}) \rangle \nonumber \\
&= \frac{i}{\omega \rho_{0}} q_{\beta} q_{\gamma} \tilde{C}_{\gamma \delta \alpha \beta}(\mathbf{q}, s = i \omega), \tag{84}
\end{align}
$$

where momentum conservation [Eq. (12)] and Eq. (58) were
used. Equation (84) and (81) show that the internal velocity field
is caused by previous external flow acting at later time via the
time-integrated stress correlations.

The two linear response relations, Eqs. (83) and (84), will be
used next to interpret the ZM memory kernel.

4.3.2 Viscosity kernel

Now, the central result of the ZM approach [Eq. (60)] can be in-
terpreted in the sense of Martin and Kadanoff$^{47,48}$. Taking the de-
viatoric stress $\sigma_{\alpha\beta}^{\text{dev}} = Q\sigma_{\alpha\beta}$ [cf. Eq. (54)] as variable A in Eq. (81)
and carrying out a Fourier transformation we can identify the lhs
of Eq. (60) as the derivative of the deviatoric stress with respect
to the external velocity field:

$$
\begin{align}
\frac{\partial \langle Q \sigma_{\alpha\beta}(\mathbf{q}, \omega) \rangle^{\text{lr}}}{\partial v_{\delta}^{\text{ext}}(\mathbf{q}, \omega)} &= \frac{V}{k_B T} \langle \sigma_{\gamma\delta}^{*} \tilde{\mathcal{R}}(s = i\omega) Q \sigma_{\alpha\beta} \rangle i q_{\gamma} \tag{85} \\
&= \tilde{M}_{\eta\zeta\alpha\beta}(\mathbf{q}, s = i\omega) i q_{\eta} \left[ \delta_{\zeta\delta} - \frac{q_{\theta}q_{\gamma}}{i\omega\rho_0} \tilde{C}_{\gamma\delta\theta\zeta}(\mathbf{q}, s = i\omega) \right]. \nonumber
\end{align}
$$

Here the minor index symmetry ($\eta\zeta \leftrightarrow \zeta\eta$) of **M** was used to
obtain the last line. The rhs of Eq. (85) may be interpreted by the
chain rule of differentiation,

$$
\frac{\partial \langle Q \sigma_{\alpha \beta}(\mathbf{q}, \omega) \rangle^{\text{lr}}}{\partial v_{\delta}^{\text{ext}}(\mathbf{q}, \omega)} = \frac{\partial \langle Q \sigma_{\alpha \beta}(\mathbf{q}, \omega) \rangle^{\text{lr}}}{\partial v_{\zeta}^{\text{tot}}(\mathbf{q}, \omega)} \frac{\partial v_{\zeta}^{\text{tot}}(\mathbf{q}, \omega)}{\partial v_{\delta}^{\text{ext}}(\mathbf{q}, \omega)}, \quad (86)
$$

where we take the a priori unspecified velocity $\mathbf{v}^{\text{tot}}$ as the mean
local velocity, $\mathbf{v}^{\text{tot}}(\mathbf{r}, t) = \mathbf{v}^{\text{ext}}(\mathbf{r}, t) + \mathbf{v}(\mathbf{r}, t)$, which contains the ex-
ternal and the (mean) internal velocity, and we anticipate that
both velocity contributions build up stresses locally.

This interpretation is supported by the second coefficient in the
rhs of Eq. (86). Looking at Eq. (84), one sees that the term in the
brackets in Eq. (85) equals

$$
\frac{\partial v_{\zeta}^{\text{tot}}(\mathbf{q}, \omega)}{\partial v_{\delta}^{\text{ext}}(\mathbf{q}, \omega)} = \delta_{\zeta\delta} - \frac{q_{\theta}q_{\gamma}}{i\omega\rho_0} \tilde{C}_{\gamma\delta\theta\zeta}(\mathbf{q}, s = i\omega), \quad (87)
$$

where the minor symmetry of **C** was also used. Equation (87)
takes into account that the external velocity field does not only
produce a stress, but also indirectly causes an internal current.

Considering the first factor in Eq. (86), comparison of
Eqs. (85), (86) and (87) shows that we can equate the memory
kernel with the response of the deviatoric stress to the total local
current:

$$
\tilde{M}_{\gamma\delta\alpha\beta}(\mathbf{q}, s = i\omega) i q_{\gamma} = \frac{\partial \langle Q \sigma_{\alpha\beta}(\mathbf{q}, \omega) \rangle^{\text{lr}}}{\partial v_{\delta}^{\text{tot}}(\mathbf{q}, \omega)}. \quad (88)
$$

It is this memory kernel that gets suggested as frequency and
wavevector dependent generalized viscosity $\eta(\mathbf{q}, \omega)$ in the gen-

eralized Maxwell model<sup>11,20</sup>; see also Sect. 6 below.

$$
\tilde{\eta}_{\gamma\delta\alpha\beta}(\mathbf{q}, \omega) \equiv \tilde{M}_{\gamma\delta\alpha\beta}(\mathbf{q}, s = i\omega), \quad (89)
$$

Thus, the memory kernel is a generalized Green–Kubo transport
coefficient which, in the limit of long wavelengths and small fre-
quencies, converges towards the viscosity. This holds as the pro-
jected dynamics simplifies, viz. $\tilde{\mathcal{R}}'(s) \to Q\tilde{\mathcal{R}}(s)Q$ for $q \to 0$, which
follows from Eq. (53); the expected Green–Kubo relation holds<sup>30</sup>.

The interpretation of the viscosity as response to the gradient
of the total local velocity field is fully in accord with the ideas by
Martin and Kadanoff. The ZM calculation disentangles the local
stress-velocity relation from the non-local one. While the former
is encoded in the viscosity kernel, the latter arises from the trans-
port of the conserved momentum over larger distances. These
considerations will be used in Sec. 6 to derive the generalized
hydrodynamic equations.

5 Shear stress correlations and comparison with
simulation data

**5.1 Theoretical predictions for shear stress correlations**

Let us consider the stress correlation function in real space de-
fined in Eq. (15). Here we aim to obtain its shear component
$C_{xyxy}(\mathbf{r}, t)$ which will be used in Sec. 5.3 for comparison with sim-
ulation data. The spatial Fourier transform of $C_{xyxy}(\mathbf{r}, t)$ is defined
in Eq. (22). Predictions for the stress correlator $C_{\alpha\beta\gamma\delta}(\mathbf{q}, t)$ were
provided in Ref. 25 (see Appendix B there) and in Ref. 20,22. The
general expression for $C(\mathbf{q}, t)$ involves five scalar material func-
tions of $q$ and $t$ (see discussion after Eq. (23)). In the regime of
small $q$ (corresponding to distances $r \gg d_l$ where $d_l$ is the particle
size) only two material functions are involved, namely the shear
relaxation modulus $G(q,t)$ and the longitudinal relaxation modu-
lus $L(q,t)$, which are related to the **E**-tensor by (see Eq. (61) in
Ref. 25)

$$
L(q,t) = E_{1111}(q,t) \quad \text{and} \quad G(q,t) = E_{1212}(q,t), \quad (90)
$$

where the 1-axis is oriented along the wavevector **q** and the 2-
axis is perpendicular to **q**$^{23,25}$. In terms of the Laplace–Carson
transform [Eq. (46)] the resulting expression for the shear stress
correlator reads

$$
\begin{equation}
\begin{split}
\tilde{C}_{xyxy}(\mathbf{q}, s) = \tilde{G}(q, s) - \tilde{G}(q, s)^2 & \frac{q_x^2 + q_y^2}{\rho_0 s^2 + q^2 \tilde{G}(q, s)} \\
& + \frac{4}{q^2} \tilde{G}(q, s)^2 q_x^2 q_y^2 \left[ \frac{1}{\rho_0 s^2 + q^2 \tilde{G}(q, s)} - \frac{1}{\rho_0 s^2 + q^2 \tilde{L}(q, s)} \right].
\end{split}
\tag{91}
\end{equation}
$$

This equation is equivalent to that obtained in Ref. 20 (see Eq. (28)
there), provided that $G(q,t)$ is identified with $G_q^\perp(t)$ and $L(q,t)$
with $\tilde{G}_q^\parallel(t)$ (see Sec. III.D of Ref. 20 for the definition of $G_q^\perp(t)$
and $\tilde{G}_q^\parallel(t)$) and the memory kernel **M** is replaced by Eq. (70).
Equation (91) is also in harmony with earlier expressions for the
shear-stress correlation function obtained for incompressible sys-
tems, $L \gg G^{11}$, and for two-dimensional glass-forming liquids $^{12}$.

Noting that $s \to \infty$ in Eq. (91) corresponds to $t \to 0$ (zero time shift), we find $C_{xyxy}(\mathbf{q}, t \to 0) = G(q, t \to 0) = 0$ leading to extremely
short-range correlations, close to $G(t=0)\delta(r)$, in real space. More generally, the short-range character of the equal-time correlations of local stresses at $t=0$ stems from the analytical $q$-dependence of the relaxation moduli. It reflects the disordered nature of well-equilibrated amorphous or liquid systems. This feature implies short correlation lengths for any local variable at $t=0$ and can be considered as a fundamental property of such systems.

Let us consider a simple model with constant—both time and wavevector independent—relaxation moduli,

$$L(q,t) = L_{\text{pl}} = \text{const}, \quad G(q,t) = G_{\text{pl}} = \text{const}. \qquad (92)$$

This approximation is valid for small $q$ for amorphous systems in the glassy plateau regime or for supercooled liquids near the glass transition where their relaxation time exceeds the observation time$^{20,23}$.

Here we shall focus on two time regimes: $\tau_L \ll t \ll \tau_T$ and $t \gg \tau_T$, where $\tau_L = r/c_L$, $\tau_T = r/c_T$, and $c_L = \sqrt{L_{\text{pl}}/\rho_0}$, $c_T = \sqrt{G_{\text{pl}}/\rho_0}$ are the longitudinal and transverse sound velocities, respectively. Using Eq. (91) the result for the shear stress correlator reads

$$C_{xyxy}(\mathbf{q}, s) \approx G_{\text{pl}} - G_{\text{pl}}^2 \frac{q_x^2 + q_y^2}{\rho_0 s^2} + 4G_{\text{pl}} \hat{q}_x^2 \hat{q}_y^2 \left( \frac{G_{\text{pl}}q^2}{\rho_0 s^2} - \frac{G_{\text{pl}}}{L_{\text{pl}}} \right), \quad \tau_L \ll 1/s \ll \tau_T, \quad (93)$$

$$C_{xyxy}(\mathbf{q}, s) \approx G_{\text{pl}} - G_{\text{pl}} \left( \hat{q}_x^2 + \hat{q}_y^2 \right) + 4G_{\text{pl}} \left( 1 - G_{\text{pl}}/L_{\text{pl}} \right) \hat{q}_x^2 \hat{q}_y^2, \quad 1/s \gg \tau_T. \quad (94)$$

Note that $s$ roughly corresponds to $1/t$ and $q$ to $1/r$. Note also that the rhs of Eq. (94) does not depend on the Laplace parameter $s$, i.e. the correlation function is static in the considered long-time regime.

Doing the inverse Fourier–Laplace transform of Eq. (93) in three-dimensional space (3D) and neglecting the first two terms in its rhs (as their contribution is localized at $r=0$) we get

$$C_{xyxy}(\mathbf{r}, t) \approx \frac{3}{2\pi} \frac{G_{\text{pl}}^2 t^2}{\rho_0 r^5} \left[ 1 - 5(\hat{x}^2 + \hat{y}^2) + 35\hat{x}^2\hat{y}^2 \right] - \frac{G_{\text{pl}}^2}{L_{\text{pl}}} \frac{1}{2\pi r^3} \left[ 1 - 3(\hat{x}^2 + \hat{y}^2) + 15\hat{x}^2\hat{y}^2 \right], \quad \tau_L \ll t \ll \tau_T. \quad (95)$$

Note that the last (second) term in the rhs of Eq. (95) vanishes in the incompressible limit $L_{\text{pl}} \to \infty$, $\tau_L \to 0$.

In the second regime, $t \gg \tau_T$, we get from Eq. (94):

$$C_{xyxy}(\mathbf{r}, t) \approx \frac{G_{\text{pl}}}{4\pi r^3} \left[ (1 - G_{\text{pl}}/L_{\text{pl}}) (2 - 6(\hat{x}^2 + \hat{y}^2) + 30\hat{x}^2\hat{y}^2) + 3(\hat{x}^2 + \hat{y}^2) - 2 \right]. \quad (96)$$

Note that the correlation function does not depend on $t$ in the time regime considered.

For two-dimensional (2D) systems the analogous correlation functions have been obtained in Ref. $^{12}$. The results read:

$$C_{xyxy}(\mathbf{r}, t) \approx -\frac{6}{\pi} \frac{G_{\text{pl}}}{r^4} (c_{\text{T}} t)^2 \cos(4\theta), \quad \tau_L \ll t \ll \tau_T, \quad (97)$$

$$C_{xyxy}(\mathbf{r}, t) \approx -\frac{1}{\pi} \left( G_{\text{pl}} - \frac{G_{\text{pl}}^2}{L_{\text{pl}}} \right) \frac{1}{r^2} \cos(4\theta), \quad t \gg \tau_T. \quad (98)$$

Here $\theta$ is the polar angle ($\hat{x} = \cos\theta$, $\hat{y} = \sin\theta$).

## 5.2 Simulation methods

To test the theoretical results we perform event driven Newtonian simulations of a binary hard sphere mixture in two (2D) and three dimensions (3D)$^{26}$. The ratio of the particle numbers is $N_s/N_l = 1$ and the aspect ratio is $d_l/d_s = 1.4$, where the indices $l$ and $s$ denote large and small particle species, respectively. In three dimensions the simulated density is $\phi = 0.59$ which is at the glass transition density of mode-coupling theory (MCT) found in$^{49,50}$. In two dimensions the simulated packing fraction is $\phi = 0.79$, which is slightly below the MCT transition density of $\phi_c = 0.7948$$^{51}$. The kinetic contribution to the stress is neglected and we define the local stress in the hard sphere system as

$$\boldsymbol{\sigma}(\mathbf{r}_B, t) = \frac{1}{V_B \delta t} \sum_{\tau_c \in \Delta t; \mathbf{r}_c \in V_B} \Delta \mathbf{r}^c \Delta \mathbf{p}^c. \quad (99)$$

Here the sum runs over all collisions that take place at time $\tau_c$ in the time window $\Delta t = [t - \delta t/2, t + \delta t/2]$ and at location $\mathbf{r}_c$ in a cubic volume $V_B$ around $\mathbf{r}_B$. $\Delta \mathbf{r}^c$ denotes the vector connecting the spheres' centers and $\Delta \mathbf{p}^c$ denotes momentum that is transferred by the collision. The stress profiles in the following section show the equilibrium correlation of the stresses according to Eq. (99) and Eq. (15).

## 5.3 Simulation results

Figure 1 shows the shear stress correlation function $C_{xyxy}(\mathbf{r}, t)$ [Eq. (15)], which would describe the response to a unit shear strain at the center [cf. Eq. (130)], for two increasing times from left to right and in the $xy$-plane in two dimensions. The octopolar symmetry of alternating negative and positive lobes, predicted by Eqs. (97) and (98), can be observed$^{23}$ as well as the propagation in space. The white circles contain no data due to the binning procedure. The right panel corresponds to the time, where the red bar is located at the bottom edge of the lower panel of Fig. 2. This also coincides with the red curve in the bottom panel of Fig. 3. The two red circles mark the expected positions of the the longitudinal (outer circle) and transversal (inner circle) wavefronts.

Figure 2 § shows the shear relaxation modulus $G(t)$. The stretched final relaxation that is observed over three decades in time is a hallmark of viscoelasticity and the reason why solid-like

<footnote>§ The same data and fits are part of Ref. $^{26}$. There, the 2D relaxation time was mistakenly given as $\tau_{v_0}/d_l = 77$ which is corrected here.</footnote>
<figure><img src="image_3.png" /><figcaption>Fig. 1 Evolution of the shear stress correlation in the xy-plane for two different times in two dimensions. The left panel shows the correlation at $t v_0/d_l = 1$ (here $v_0$ is thermal velocity) and the right one at $t v_0/d_l = 3$, which coincides with the red radial curve in the lower panel of Fig. 3. The red circles in the right panel indicate the positions of longitudinal and transversal sound waves.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>Fig. 2 Shear relaxation modulus of the investigated systems in 3D (top) and 2D (bottom). The colored bars mark the times of the curves in Fig. 3. The red line shows a stretched exponential fit. For 3D the parameters are: $\beta = 0.43$, $\tau v_0/d_l = 90$ and $G_{pl}/(k_B T n_0) = 43.12$. In 2D the parameters are $\tau v_0/d_l = 103$: $\beta = 0.55$ and $G_{pl}/(k_B T n_0) = 43.91$. The red bars on the vertical axes mark $G_{pl}$. In both panels, the blue crosses indicate the time $t^*$ used in Fig. 4.</figcaption></figure>

stress correlations can be observed in the liquid state for intermediate times. Stretched exponential decay, implying the existence of a distribution of relaxation times, which goes beyond Maxwell's original exponential ansatz, is very familiar in supercooled liquids$^{2,4,52}$; because we address only the asymptotic limits, this phenomenon does not affect our work. The final decay of $G(t)$ defines the state to be liquid. The colored bars at the bottom edges mark the times of the similarly colored curves in Fig. 3. Their times all lie in the window where $G(t)$ is close to its plateau value $G_{pl}$. The red curves in both panels show a fit to a stretched exponential, $G(t) = G_{pl} \exp(-(t/\tau)^\beta)$, which we use as a tool to extract $G_{pl}$. For the quantitative analysis presented in Fig. 3 the fitted plateau value $G_{pl}$ will be used to connect our simulation results to the theory of Sec. 5.1. This fit parameter is also used to determine the transversal sound velocity $c_T = \sqrt{G_{pl}/\rho_0}$. For the times marked by the colored bars the temporal evolution of the stress profile is presented in Fig. 3. The blue crosses mark the time for which the far-field behavior of the shear stress correlations is shown in Fig. 4 for two and three dimensions, respectively.

Figure 3 shows the evolution of $C_{xyxy}^{\text{diag}}(r,t)$ for increasing times along the diagonals in the xy-plane for 2D and 3D. The colors of the curves coincide with the times where the bars of the same color at the lower edges of Fig. 2 are located. The red vertical dashed bars denote the position of the transverse sound wave, $r_1 = c_T t$, and of the longitudinal sound wave, $r_2 = c_L t$, respectively. The longitudinal sound velocity $c_L$ was measured from the dispersion relation of the longitudinal velocity correlator. The values are in 2D: $c_T/c_L = 0.21$ and in 3D: $c_T/c_L = 0.28$. The black dashed lines show the expected power-laws derived in section 5.1 in 3D and 2D. For short distances, $r \ll c_T t$, we plot Eq. (96) (for 3D) and Eq. (98) (for 2D) which give for $r$ along the diagonals and incompressible systems ($L_{pl} \gg G_{pl}$):

$$C_{xyxy}(\mathbf{r}) = \begin{cases} \frac{9G_{\text{pl}}}{8\pi r^3}, & 3\text{D}, \\ \frac{G_{\text{pl}}}{\pi r^2}, & 2\text{D}. \end{cases} \qquad (100)$$

The correlation function is static in this case. The agreement of the incompressible model is not perfect. Among the approximations made is the use of $G_{pl}$ from the fit instead of the specific measured $G(t)$ value. These times, however, are so short that force correlations have not decayed yet$^{26}$. This means that the curves in Fig. 3 at short distances still contain fast processes leading to deviations from the idealized model considered in Sec. 5.1 (and in Sec. 6). As these processes have decayed at the time $t^*$ of the blue crosses in Fig. 2, the agreement in Fig. 4 is more convincing. It is seen that, in contrast to the $r^{-(d+2)}$ power law, the amplitude is not strongly time dependent and (closely) follows $G(t)$, which we take as a proxy for the constant $G_{pl}$ in the given time range.

For distances $c_T t \ll r \ll c_L t$ the theoretical predictions are given in Eq. (95) (for 3D) and Eq. (97) (for 2D), along the diagonals

<footnote>¶ The times of the 3D curves are shifted by $\delta \tau v_0/d_l = -0.1$ compared to the same data points in Ref. 26 due to a different time assignment procedure.</footnote>
leading to

$$C_{xyxy}(\mathbf{r}, t) = \begin{cases} \frac{57(c_{\mathrm{T}}t)^2 G_{\mathrm{pl}}}{8\pi r^5}, & 3\mathrm{D}, \\ \frac{6(c_{\mathrm{T}}t)^2 G_{\mathrm{pl}}}{\pi r^4}, & 2\mathrm{D}. \end{cases} \qquad (101)$$

This is confirmed by simulation as can be seen from the two dashed black lines at larger distances in Fig. 3. In each panel, the amplitudes of the power laws are calculated from Eq. (101) at the times of the red and olive curves, while $G_{\mathrm{pl}}$ was again taken from the stretched-exponential fit in Fig. 2.

The dashed-dotted black curve in the bottom panel of Fig. 3 is calculated via Eq. (42) from Ref. 12. Equation (42) represents the inverse Fourier-Laplace transform of Eq. (91) when using the model specified in Eq. (92). In fitting Eq. (42) to our simulation data we included the fact that in hard sphere simulations the stress can only be computed by averaging over a time window of $\delta t$ [cf. Eq. (99)] with $\delta t = d_l/v_0$ for the simulations in 2D. This was done via

$$\frac{1}{(\delta t)^2} \left\langle \left( \int_0^{\delta t} dt' \sigma_{xy}(\mathbf{r}', t') \right) \left( \int_{3\delta t}^{4\delta t} dt \sigma_{xy}(\mathbf{r}, t) \right) \right\rangle = \frac{k_B T}{(\delta t)^2} \int_0^{\delta t} dt' \int_{3\delta t}^{4\delta t} dt C_{xyxy}(t-t', |\mathbf{r}'-\mathbf{r}|), \quad (102)$$

so the theory is time averaged in the same way for consistency. In using Eq. (42) from Ref. 12, a correction of order $O((c_T/c_L)^4)$ is neglected. The integration limits describe the averaging done for the red 2D-curve.

For large times $t = t^*$, indicated by the blue crosses in Fig. 2, the $r^{-d}$ power-laws span the whole simulation box. At this time all force fluctuations have decayed and one observes similar behavior as in a force-free solid state but with a time-dependent shear modulus, $G_{\mathrm{pl}} \to G(t^*)$. The results of this analysis are shown in Fig. 4, i.e. the simulation data together with Eq. (100) when choosing $G_{\mathrm{pl}} = G(t^*)$ and $t = t^*$.

# 6 Generalized Maxwell model

## 6.1 Viscoelastic fluid case

In this section the equations of generalized linear hydrodynamics are derived from the ZM results of Sec. 4. This provides another route to the far-field power-law stress correlations obtained in Sec. 5, which will also be generalized to arbitrary geometry. The resulting generalized Maxwell model has already been derived for incompressible systems 11,20, where the microscopically founded ZM equations for transverse flow were analyzed in the small wavevector regime. The equations generalize the Navier-Stokes equations to include frequency-dependent, viz. non-Markovian, viscosity kernels. By going beyond the familiar Maxwell model, the generalized Maxwell model achieves to connect the limit of linear local (viscous) stress relaxation to the opposite limit of linear nonlocal (Hookean) elasticity. Here a short-cut shall be presented, which has already been used to find the generalized hydrodynamic equations in (compressible) systems of Brownian particles 21.

<figure><img src="image_5.png" /></figure>

<figure><img src="image_6.png" /><figcaption>Fig. 3 Temporal evolution of the shear stress correlations along the diagonals for three (top) and two dimensions (bottom). The dashed lines for shorter distances denote a power law following $r^{-d}$ with an amplitude given by the shear modulus $G_{\mathrm{pl}}$, i.e. Eq. (100). The two dashed lines at larger distances in each panel denote the power law $r^{-(d+2)}$ given by Eq. (101) with the times of the red and olive bars in Fig. 2 for each panel. The vertical dashed lines denote the positions of shear (slower) and longitudinal (faster) sound waves at the time of the red curve in both panels; the r-values correspond to the circles in the right panel of Fig. 1. The black dashed-dotted curve in the lower panel shows Eq. (42) of Ref. 12, adapted by Eq. (102) to the time of the red curve.</figcaption></figure>

Using a coarse-graining approach 53 we define the mass density field of continuum mechanics as

$$\rho(\mathbf{r}, t) = m \sum_{i=1}^{N} \phi(\mathbf{r} - \mathbf{r}_i(t)) \qquad (103)$$

and the corresponding current as

$$\mathbf{J}(\mathbf{r}, t) = m \sum_{i=1}^{N} \mathbf{v}_i(t) \phi(\mathbf{r} - \mathbf{r}_i(t)), \qquad (104)$$

where we assume the system to consist of monodisperse particles with the same mass $m$. The coarse-graining function $\phi(\mathbf{r})$ generalizes the Dirac delta appearing in the definitions of Sec. 2. The reason for introducing $\phi(\mathbf{r})$ is that continuum mechanics involves the concept of a continuum field which implies a spatial averaging over some finite length scale. The function $\phi(\mathbf{r})$ accounts for this averaging and can be taken e.g. as a Gaussian of width $w$ and normalization $\int d^d\mathbf{r} \phi(\mathbf{r}) = 1$ $^{28}$. While the limit $w \to 0$ leads back
<figure><img src="image_7.png" /><figcaption>Fig. 4 Diagonal shear stress correlations for time $t^*$, where the blue crosses in Fig. 2 are located, in $d=2$ and $d=3$. An $r^{-d}$ power-law well describes the simulation data for all $r$ in the whole simulation box. The dashed lines show the theoretical predictions from Eq. (100) with $G_{pl} = G(t^*)$, where $G(t^*)$ is the value of $G(t)$ at the blue crosses in Fig. 2. Here we find the agreement to be much better than in Fig. 3, since the considered time is truly larger than the correlation time of the forces in the system, which is discussed in more detail in Ref. 26.</figcaption></figure>

to the definitions with the Dirac delta given in Sec. 2, here we are interested in the limit of large $w$ where a generalized hydrodynamic approach becomes valid.

According to the momentum balance equation (11), the change in the coarse-grained momentum field is connected to the change in the average velocity and to the force density field

$$ \partial_t J(r,t) = \rho_0 \partial_t v(r,t) = \nabla \cdot \boldsymbol{\sigma}(r,t). \quad (105) $$

The force density is given by the divergence of the coarse grained stress tensor, which follows from the Irving-Kirkwood expression, Eq. (13), upon inverse Fourier transformation and replacing the Dirac delta by $\phi(r)$:

$$ \begin{aligned} \boldsymbol{\sigma}(\mathbf{r},t) = & -\sum_j m \mathbf{v}_j(t) \mathbf{v}_j(t) \phi(\mathbf{r} - \mathbf{r}_j(t)) - \\ & \frac{1}{2} \sum_{k \neq j} \mathbf{r}_{jk}(t) \mathbf{F}_{jk}(t) \int_0^1 ds \phi(\mathbf{r} - \mathbf{r}_j(t) + s \mathbf{r}_{jk}(t)). \end{aligned} \quad (106) $$

It defines the coarse grained stress tensor up to a divergence free term which can be argued to vanish$^{31}$.

We now take the coarse-graining size $w$ to be so large that the associated smoothing in space and time eliminates fluctuations and the stress tensor of a quiescent system is obtained: $\boldsymbol{\sigma}(\mathbf{r},t) \to -p\mathbf{1}$ where $\mathbf{1}$ is the identity matrix. The average stress equals the negative pressure which is constant and homogeneous. As expected, Eq. (105) then states that a quiescent system is force-free and does not flow. $||$

Leaving the uninteresting quiescent system behind, Onsager's regression hypothesis can be used. It assumes local thermodynamic equilibrium holds and that averages are space and time dependent because of the applied external perturbations$^{9,55}$. They are calculated in the linear response formalism of Sect. 4.3.

$$ \boldsymbol{\sigma}(\mathbf{r},t) = \langle \boldsymbol{\sigma}(\mathbf{q} \to 0, t) \rangle_{(\mathbf{r},t)}^{\mathrm{lr}}. \quad (107) $$

Here, the coarse-graining size $w$ is so large that in the evaluation of the linear response the small wavevector limit is taken and spatial variation only remains because of the smooth variation of the external perturbation. Applying the by now familiar decomposition, $\boldsymbol{\sigma} = \mathcal{P}\boldsymbol{\sigma} + \mathcal{Q}\boldsymbol{\sigma}$, the generalized hydrodynamic approximation for the $\mathcal{P}$-projected stress becomes (with the help of Eqs. (62), (41) and using that $S(q \to 0) = n_0 k_B T \kappa_T$ is determined by the isothermal compressibility $\kappa_T$ for monodisperse systems)

$$ \langle \mathcal{P} \boldsymbol{\sigma}(\mathbf{q} \to 0, t) \rangle_{(\mathbf{r},t)}^{\mathrm{lr}} \to \frac{-\delta \rho(\mathbf{r},t)}{\kappa_T \rho_0} \mathbf{1} = -p(\mathbf{r},t) \mathbf{1} \quad (108) $$

where $p(\mathbf{r},t)$ is the pressure fluctuation of an isothermal system (with no temperature fluctuations on all length scales, see e.g. Eq. (8.3.14) of Ref. 30). Equation (108) identifies the $\mathcal{P}$-projected stress with the mean (or spherical) stress of continuum mechanics. Note that the remainder $\mathcal{Q}\boldsymbol{\sigma}$ can possess a finite trace.

The identification of the $\mathcal{Q}$-projected stress as the deviatoric stress of hydrodynamics connects the coarse-grained description to the one based on correlation-functions:

$$ \delta \boldsymbol{\sigma}(\mathbf{r},t) = \langle \mathcal{Q} \boldsymbol{\sigma}(\mathbf{q} \to 0, t) \rangle_{(\mathbf{r},t)}^{\mathrm{lr}}. \quad (109) $$

We again assume local thermodynamic equilibrium and presume that the deviatoric stress is produced by the local current perturbation according to the linear response calculation in the microscopic approach. Importantly, we follow the idea by Martin and Kadanoff, as worked out in section 4.3, that the proper Onsager transport kernel, viz. the viscosity kernel, should be defined in the stress response to the local field, which differs from the external field. Recall from Eq. (87): $\mathbf{v}^{\mathrm{tot}}(\mathbf{r},t) = \mathbf{v}^{\mathrm{ext}}(\mathbf{r},t) + \mathbf{v}(\mathbf{r},t)$. Keeping this in mind, Eq. (88) then gives the viscosity kernel ($\boldsymbol{\eta}(\mathbf{q} = \mathbf{0},t)$) as the response function of the stress to the gradient of the local current. While hydrodynamics assumes rapidly decaying kernels, we will first consider the general non-Markovian case, where the time-dependent kernel $\boldsymbol{\eta}(t) = \boldsymbol{\eta}(\mathbf{q} = \mathbf{0},t)$ encodes retarded response. Because of isotropy and the limit $\mathbf{q} \to 0$, it can be given by two functions of time$^{10,20}$

$$ \eta_{\alpha\beta\gamma\delta}(t) = \zeta(t) \delta_{\alpha\beta} \delta_{\gamma\delta} + G(t) (\delta_{\alpha\gamma} \delta_{\beta\delta} + \delta_{\alpha\delta} \delta_{\beta\gamma}). \quad (110) $$

Here $G(t)$ is the shear relaxation modulus and $\zeta(t)$ is a second (mixed) modulus related to shear and bulk moduli$^{10}$. **

Putting all preceding results together the generalized hydrodynamic equation for the stress tensor reads

$$ \boldsymbol{\sigma}(\mathbf{r},t) = -p(\mathbf{r},t)\mathbf{1} + \delta\boldsymbol{\sigma}(\mathbf{r},t), \quad (111) $$

<footnote>|| Note that continuum mechanics descriptions often do not consider the equilibrium pressure as part of the stress$^{54}$.</footnote>

<footnote>** The modulus $\zeta(t)$ can be related to the longitudinal modulus $L(t) = \zeta(t) + 2G(t)$ or to the transverse (or mixed) relaxation modulus $M(q \to 0, t) = M(t)^{12,25}$ by $\zeta(t) = M(t) - M_e(q \to 0) = M(t) - \frac{1}{\kappa_T}$, since $M_e(q \to 0) = L_e(q \to 0)^{25}$ and $L_e(q \to 0) = 1/\kappa_T$ for monodisperse systems due to Eq. (33). In a solid holds $M(t) = \lambda^{23}$.</footnote>
and

$$ \delta\boldsymbol{\sigma}(\mathbf{r},t) = \eta(t) * \nabla\mathbf{v}^{\text{tot}}(\mathbf{r},t), \quad (112) $$

with “*” denoting the convolution (with respect to time).

We now consider the same situation as in Sec. 4.3: The application of a perturbative external force creates a velocity field $\mathbf{v}^{\text{tot}}(\mathbf{r},t) = \mathbf{v}^{\text{ext}}(\mathbf{r},t) + \mathbf{v}(\mathbf{r},t)$, comprising an external ($\mathbf{v}^{\text{ext}}$) and an internal ($\mathbf{v}$) part, and we anticipate that both parts generate a stress. Using Eqs. (110) and (112) the response of the deviatoric stress to the external velocity field reads

$$ \delta\boldsymbol{\sigma}^{\text{ext}}(\mathbf{r},t) = \zeta(t) * \text{Tr}(\boldsymbol{\kappa}^{\text{ext}}(\mathbf{r},t))\mathbf{1} + 2G(t) * \boldsymbol{\kappa}^{\text{ext}}(\mathbf{r},t), \quad (113) $$

where Tr(...) means the trace of the matrix and we utilized the index symmetries of $\eta$ to replace the velocity gradient by the (symmetrized) strain-rate (or velocity-gradient) tensor

$$ \boldsymbol{\kappa}^{\text{ext}}(\mathbf{r},t) = \frac{1}{2} \left( \nabla \mathbf{v}^{\text{ext}}(\mathbf{r},t) + (\nabla \mathbf{v}^{\text{ext}}(\mathbf{r},t))^{\text{T}} \right), $$

with the superscript “T” denoting the transpose operation. From Eq. (113) we find the external force density

$$ \nabla \cdot \delta\boldsymbol{\sigma}^{\text{ext}}(\mathbf{r},t) = \zeta(t) * \nabla \text{Tr}(\boldsymbol{\kappa}^{\text{ext}}(\mathbf{r},t)) + 2G(t) * \nabla \cdot \boldsymbol{\kappa}^{\text{ext}}(\mathbf{r},t). \quad (114) $$

Similarly, application of Eqs. (110) and (112) to the internal part of the velocity field leads to the internal viscoelastic force density

$$ \nabla \cdot \delta\boldsymbol{\sigma}^{\text{int}}(\mathbf{r},t) = (\zeta(t) + G(t)) * \nabla(\nabla \cdot \mathbf{v}(\mathbf{r},t)) + G(t) * \nabla^2 \mathbf{v}(\mathbf{r},t). \quad (115) $$

Gathering all forces and inserting into the moment equation (11), we obtain the generalization of the linearized Navier-Stokes equation for the momentum density

$$ \rho_0 \partial_t \mathbf{v}(\mathbf{r},t) = -\nabla p(\mathbf{r},t) + G(t) * \nabla^2 \mathbf{v}(\mathbf{r},t) + (\zeta(t) + G(t)) * \nabla(\nabla \cdot \mathbf{v}(\mathbf{r},t)) \\ + 2G(t) * \nabla \cdot \boldsymbol{\kappa}^{\text{ext}}(\mathbf{r},t) + \zeta(t) * \nabla \text{Tr}(\boldsymbol{\kappa}^{\text{ext}}(\mathbf{r},t)). \quad (116) $$

To complete this equation, mass conservation [Eq. (10)] and Eq. (108) are employed to get the pressure field,

$$ \nabla \cdot \mathbf{v}(\mathbf{r},t) = -\frac{1}{\rho_0} \partial_t \delta\rho(\mathbf{r},t) = -\kappa_T \partial_t p(\mathbf{r},t). \quad (117) $$

The case of incompressible fluids is recovered for vanishing compressibility $\kappa_T \to 0$, i.e. $\nabla \cdot \mathbf{v} = 0$.

Equation (116) defines the generalized Maxwell model. In analogy to the Maxwell model of viscoelasticity$^{30,54}$, Eq. (116) generalizes the familiar hydrodynamic equations of liquids by including memory effects via the kernels $\zeta(t)$ and $G(t)$. The assumption of locality is kept because $\zeta(t)$ and $G(t)$ are independent of position. In generalization of Maxwell's macroscopic consideration, it splits the local current into two contributions, the external and the internal one. This enables the model to recover the nonlocality of elasticity in solid states. The time dependence of

the kernels enables the generalized Maxwell model to encompass both the case of viscous fluids and of amorphous solids. These cases are given by two limits of $\zeta(t)$ and $G(t)$: For viscous fluids we take

$$ \zeta(t) = \eta_2 \delta(t), \quad G(t) = \eta \delta(t), \quad (118) $$

where the constant coefficients are the shear viscosity $\eta$ and the second coefficient of viscosity $\eta_2$, which is related to the bulk viscosity$^{55}$. For amorphous solids, on the other hand, the memory kernels adopt constant plateau values, $\zeta_{\text{pl}}$ and $G_{\text{pl}}$, in a time window $\tau_{\text{mic}} \ll t \ll \tau$, where $\tau_{\text{mic}}$ is relaxation time of short-time processes and $\tau$ is the terminal (alpha) relaxation time of the system. Since $\tau$ for amorphous solids is much larger than the longest observation time of the system, the plateau values dominate the behavior of $\zeta(t)$ and $G(t)$ and determine the elastic moduli of the (temporarily frozen) solid$^{23,56}$. Based on the standard notation for isotropic solids we thus choose to express $\zeta_{\text{pl}}$ and $G_{\text{pl}}$ in terms of the Lamé coefficients $\lambda$ and $\mu$. So, we write (for $\tau_{\text{mic}} \ll t \ll \tau$)

$$ \zeta(t) = \zeta_{\text{pl}} \equiv \tilde{\lambda} = \lambda - 1/\kappa_T, \quad G(t) = G_{\text{pl}} \equiv \mu. \quad (119) $$

To find the Green's function of the problem one chooses the external strain-rate tensor as $\boldsymbol{\kappa}^{\text{ext}}(\mathbf{r},t) = \boldsymbol{\kappa}^{\text{ext}}\delta(\mathbf{r}) = [\boldsymbol{\varepsilon}^{\text{ext}}\delta(t)]\delta(\mathbf{r})$, where $\boldsymbol{\varepsilon}^{\text{ext}}$ denotes a constant strain (times a volume) and $\boldsymbol{\kappa}^{\text{ext}}$ is a constant strain rate (times a volume). Upon inserting this choice and Eqs. (118) and (119) into Eq. (113), we have

$$ \delta\boldsymbol{\sigma}^{\text{ext}}(\mathbf{r},t) = \begin{cases} \delta(\mathbf{r}) (\eta_2 \text{Tr}(\boldsymbol{\kappa}^{\text{ext}})\mathbf{1} + 2\eta \boldsymbol{\kappa}^{\text{ext}}), & \text{liquid,} \\ \delta(\mathbf{r})\Theta(t) (\tilde{\lambda} \text{Tr}(\boldsymbol{\varepsilon}^{\text{ext}})\mathbf{1} + 2\mu \boldsymbol{\varepsilon}^{\text{ext}}), & \text{solid,} \end{cases} \quad (120) $$

where $\Theta(t)$ is the Heaviside function. In the incompressible limit, the Green's function was already derived from the ZM equations$^{11,20}$. Here we extend this analysis to obtain the two far-field power-laws derived for the shear stress correlations in Sec. 5.1 and observed in the simulations as discussed in Sec. 5.3.

## 6.2 Limit of incompressible elastic solids

For incompressible solids Eqs. (116), (117), (119) and (120) imply

$$ \rho_0 \partial_t^2 \mathbf{u}(\mathbf{r},t) = -\nabla p(\mathbf{r}) + \mu \nabla^2 \mathbf{u}(\mathbf{r},t) + 2\mu (\boldsymbol{\varepsilon} \cdot \nabla) \delta(\mathbf{r}) \Theta(t), \quad (121) $$

$$ \nabla \cdot \mathbf{u}(\mathbf{r},t) = 0, \quad (122) $$

where $\mathbf{u}(\mathbf{r},t)$ denotes the displacement field entering via $\mathbf{u}(\mathbf{r},t) = \int_0^t \mathrm{d}t' \mathbf{v}(\mathbf{r},t')$. Note that we suppressed the superscript “ext” for the external strain $\boldsymbol{\varepsilon}^{\text{ext}}$ for notational simplicity.

Equation (121) is similar to the equations of incompressible unsteady Stokes flow. In fact, if one considers the overdamped limit of colloids, the second-order time derivative on the rhs becomes a first-order derivative and the resulting equation is the same as for incompressible Stokes flow. The equations of Stokes flow are solved in Ref. 57 (see pp. 149). In two dimensions the problem is solved in Ref. 58. In the following sections, only asymptotic solutions are derived, since this shows more clearly how the power
laws discussed in Sec. 5.1 and Sec. 5.3 emerge.

## 6.3 First power law regime

Inserting the incompressibility constraint (122) into Eq. (121) we find for $t > 0$

$$ \nabla^2 p(\mathbf{r}) = 2\mu(\nabla \cdot \boldsymbol{\varepsilon} \cdot \nabla)\delta(\mathbf{r}). \quad (123) $$

This Poisson problem has the solution (for $r > 0$)

$$ p(\mathbf{r}) = \frac{\mu}{2\pi r^3} (\boldsymbol{\varepsilon} : \mathbf{1} - 3\hat{\mathbf{r}} \cdot \boldsymbol{\varepsilon} \cdot \hat{\mathbf{r}}). \quad (124) $$

Consider now distances much larger than the position of the transverse wavefront, i.e. $r \gg c_T t$. Since $c_T = \sqrt{\mu/\rho_0}$, this is analogous to considering $\mu = 0$. Then, Eq. (121) becomes

$$ \lim_{(c_T t/r) \to 0} \rho_0 \partial_t^2 \mathbf{u}(\mathbf{r}, t) = -\nabla p(\mathbf{r}), \quad (125) $$

which can be integrated over time to read

$$ \mathbf{u}(\mathbf{r}, t) = -\frac{t^2}{2\rho_0} \nabla p(\mathbf{r}). \quad (126) $$

With the result for the pressure from Eq. (124) the final equation for the displacement field is given by

$$ \mathbf{u}(\mathbf{r}, t) = \frac{3(c_T t)^2}{4\pi r^4} (\mathrm{Tr}(\boldsymbol{\varepsilon}\hat{\mathbf{r}} - 5\hat{\mathbf{r}}(\hat{\mathbf{r}} \cdot \boldsymbol{\varepsilon} \cdot \hat{\mathbf{r}}) + 2\boldsymbol{\varepsilon} \cdot \hat{\mathbf{r}})). \quad (127) $$

The stress follows from $\boldsymbol{\sigma} = \mu(\nabla\mathbf{u} + (\nabla\mathbf{u})^T)$. For the special case of a shear inclusion$^{59}$ the strain tensor reads (of infinitesimal constant strength $\boldsymbol{\varepsilon}$)

$$ \boldsymbol{\varepsilon} = \frac{\boldsymbol{\varepsilon}}{2}(\hat{\mathbf{x}}\hat{\mathbf{y}} + \hat{\mathbf{y}}\hat{\mathbf{x}}), \quad (128) $$

and the resulting stress for $r \gg c_T t$ becomes

$$ \sigma_{xy}(\mathbf{r}, t) = \frac{3(c_T t)^2 \mu}{2\pi r^5} (1 - 5(\hat{x}^2 + \hat{y}^2) + 35\hat{x}^2\hat{y}^2) \boldsymbol{\varepsilon}. \quad (129) $$

Let us emphasize two aspects of this result:

(i) Our elastic model shows that the power-law decay $r^{-(d+2)}$ and the amplitude $\sim t^2$ of the shear stress have their origin in the adjustment of the pressure. Interestingly, the model gives back the true far-field behavior in an incompressible system at finite time. As discussed in Sec. 5.3, simulations and the analysis of Eq. (116) show that Eq. (129) also holds in a viscoelastic fluid with slowly decaying amplitude$^{26}$.

(ii) The rhs of Eq. (129) is proportional to the shear stress correlation function given in Eq. (95) in the limit of incompressible systems ($L_{pl} \to \infty$), i.e. the stress response, $\sigma_{xy}$, to the shear perturbation $\boldsymbol{\varepsilon}$ is

$$ \sigma_{xy}(\mathbf{r}, t) = C_{xyxy}(\mathbf{r}, t)\boldsymbol{\varepsilon}. \quad (130) $$

## 6.4 Second power law regime

For $t \gg r/c_T$ one can consider the stationary problem, i.e. Eq. (121) where the time derivative vanishes. This resembles the identical problem solved in Ref.$^{59}$, leading to the Eshelby

inclusion solution

$$ \mathbf{u}(\mathbf{r}) = \frac{3}{4\pi r^2} \left( \frac{1}{3} \mathrm{Tr}(\boldsymbol{\varepsilon})\hat{\mathbf{r}} - (\hat{\mathbf{r}} \cdot \boldsymbol{\varepsilon} \cdot \hat{\mathbf{r}})\hat{\mathbf{r}} \right). \quad (131) $$

For a pure shear transformation, Eq. (128), the shear stress components show the characteristic $r^{-d}$ decay and a four-fold spatial symmetry in the $xy$-plane

$$ \sigma_{xy}(\mathbf{r}) = \frac{3\mu}{4\pi r^3} (10\hat{x}^2\hat{y}^2 - (\hat{x}^2 + \hat{y}^2)) \boldsymbol{\varepsilon}. \quad (132) $$

Again we find that the stress response stated in the above equation agrees with the rhs of Eq. (96) (for $L_{pl} \to \infty$) so that Eq. (130) holds in this long-time regime as well.

# 7 Discussion and conclusions

1. In the present paper we provided a proof of the general relation between the tensor of stress correlations, $C(q,t)$, and the tensor of generalized relaxation moduli, $E(q,t)$, using the rigorous Zwanzig-Mori (ZM) projection operator formalism. The major advance (as compared to the previous work with the same formalism$^{21}$) is that the $C-E$ relation [cf. Eqs. (69), (70)] is now established to be exact for any $q$ in classical systems with Newtonian dynamics, thereby extending prior work on stress correlations with the ZM approach$^{20,22}$.

2. Recently, the $C-E$ relation was also obtained by a different approach using the concepts of stress noise, flow-induced deterministic stress and the fluctuation-dissipation theorem (FDT)$^{25}$. Here we demonstrated full agreement between this approach and the ZM approach (see Sec. 4.2.2 and Sec. 4.2.3). In particular, we proved that the memory kernel $M(q,t)$ of the ZM approach coincides with the correlation function of stress noise $C^n(q,t)$ [cf. Eq. (72)] and that the equilibrium elasticity tensor $E^e(q)$ equals the equilibrium (static) correlation function of the $P$-projected stress, $\sigma^P = P\sigma$ [cf. Eq. (71)].

Physically, these results are gratifying. The memory kernel $M(q,t)$ captures the residual time dependence of the system after removal of the ensemble-averaged dynamics of the conserved quantities, collective density and collective momentum. In the long-wavelength limit, it reduces to the non-Markovian viscosity kernel, $M(q \to 0, t) \to \eta(t)$. The stress-noise correlation function $C^n(q,t)$ allows for a similar interpretation: When constraining the dynamics of the system's microstate in phase space such that $J(q,t) = 0$ for all $t < 0$$^{25}$, the time dependence of the stress is fully given by the stress noise whose ensemble-averaged correlation function then reflects the residual stress relaxation void of contributions from collective density or momentum fluctuations. Thus, at the level of the ensemble-averaged correlation functions, $M(q,t)$ and $C^n(q,t)$, the approaches based on the ZM projection operator formalism and on the no-current constraint become equivalent.

3. This equivalence is remarkable, since projected and constrained dynamics are different. In the constrained case, just like with classical Newtonian dynamics, the time dependence of a variable $A(\Gamma) \equiv A(\Gamma, 0)$ is defined by the evolution of the system's microstate:

$$ A(\Gamma, t) = A(\Gamma'(t)), $$
where $\Gamma'(t)$ is the point in phase space where the system arrives at time $t$, if initially (at $t=0$) its microstate was $\Gamma$. Thus, in this case a change in $A$ is solely related to the variations of the phase-space coordinates. By contrast, the projected dynamics involves, in addition to such variations of $\Gamma$, also small increments of $A$ associated with the projector $Q=1-P$ [due to Eq. (52)], which are proportional to the conserved variables at each time step.

4. Both approaches employ a decomposition of the total stress $\boldsymbol{\sigma}$. The stress-noise approach splits the stress in stress-noise and deterministic parts, $\boldsymbol{\sigma} = \boldsymbol{\sigma}^n + \boldsymbol{\sigma}^d$ [cf. Eq. (73)], while the ZM approach decomposes the stress in a deviatoric part, $\boldsymbol{\sigma}^{\text{dev}} = Q\boldsymbol{\sigma}$, and a density related, $P$-projected part, $\boldsymbol{\sigma}^P = P\boldsymbol{\sigma}:\boldsymbol{\sigma} = Q\boldsymbol{\sigma} + P\boldsymbol{\sigma}$. Under the nonequilibrium condition of no flow at wavevector $\boldsymbol{q}$ for $t < 0$, Eqs. (79) and (80) establish the correspondence between both decompositions for $t \ge 0$ and show that in general $\boldsymbol{\sigma}^n(\boldsymbol{q}, t) \neq \boldsymbol{\sigma}^{\text{dev}}(\boldsymbol{q}, t)$ and $\boldsymbol{\sigma}^d(\boldsymbol{q}, t) \neq \boldsymbol{\sigma}^P(\boldsymbol{q}, t)$. Thus, not only the relevant dynamics are different (the reduced dynamics defined by the evolution operator $\mathcal{R}'(t)$ in the ZM approach vs. the constrained dynamics for the stress-noise approach), but also the respective stress variables are different as well. This notion emphasizes the nontrivial (i.e. physically meaningful) nature of the basic C-M relation, Eq. (72).

5. Returning to Eq. (70) defining the tensor of generalized relaxation moduli, $E(\boldsymbol{q}, t)$, we observe that it consists of a viscoelastic term $M(\boldsymbol{q}, t)$ and a purely elastic term $E^e(\boldsymbol{q})$ which reflects the isothermal compressibility for monodisperse systems in the limit $q \to 0$ (cf. Sec. 3 and Sec. 6.1). The Fourier transform of $M(\boldsymbol{q}, t)\Theta(t)$ is $\eta(\boldsymbol{q}, \omega)$, which is called the generalized viscosity (cf. Sec. 4.3) because the two components of $\eta(\boldsymbol{q} \to 0, \omega \to 0)$ are equal to the steady-state shear $\eta$ and bulk viscosity $\eta_V$, respectively. In the hydrodynamic regime the function $\eta_{xyxy}(\boldsymbol{q} \to 0, \omega) = \eta^*(\omega)$ is known as the complex viscosity$^{54}$. The generalized viscosity $\eta(\boldsymbol{q}, \omega)$ is dominated by its real part, $\Re\eta(\boldsymbol{q}, \omega)$, in liquids at low frequencies $\omega$, reflecting the dissipative nature of this regime. By contrast, in glassy systems (like supercooled liquids below the glass transition) the generalized viscosity is not real even at low $\omega$: $\eta(\boldsymbol{q}, \omega)$ is then typically dominated by its imaginary part, pointing to a largely elastic nature of the strain-induced stress in such systems. Importantly, in Sec. 4.3.2, the viscosity kernel was identified in the response of the stress to the local velocity gradient. It consists of the external and the internal velocity fields, because the external perturbation generates stress but also momentum currents. It is also noteworthy that the Laplace–Carson transform [cf. Eq. (46)] of the generalized relaxation moduli, $E(\boldsymbol{q}, s = i\omega)$, is analogous to the complex dynamical shear modulus $G^*(\omega) = G'(\omega) + iG''(\omega)$, which is one of the main rheological characteristics of a material$^{54}$.

6. In Sec. 5 we illuminate two main regimes for the real-space decay of stress correlations in an amorphous system (near or below the glass transition) for $c_L t \gg r$: $C_{xyxy}(\boldsymbol{r}, t) \propto r^{-d}$ at $r \ll c_T t$ and $C_{xyxy}(\boldsymbol{r}, t) \propto r^{-(d+2)}$ at $r \gg c_T t$ (where $d$ is the space dimension; $c_L$ and $c_T$ are longitudinal and transverse sound velocities, respectively). It is also shown (cf. Sec. 6) that the shear stress response of an incompressible system to a weak (of strength $\epsilon$) local perturbation near the origin (for example, due to a small particle inclusion) follows the same power laws. Moreover, we

demonstrated that this shear stress response is proportional to $C_{xyxy}(\boldsymbol{r}, t)$, see Eq. (130).

Although further work is needed to fully elucidate this point, it appears as if Eq. (130) comes from the FDT (a Kubo equation like Eq. (81)). Another subtle question is whether such a relation would also hold for compressible systems with time-dependent relaxation moduli. Moreover, it is an interesting open problem how to find the perturbation strength $\epsilon$ in the case of a locally strong deformation of the system involving nonlinear effects. We plan to return to these issues elsewhere.

## A Definitions employed by MZF

MZF consider a system of classical particles where all particles have the same mass, $m_i = m^{11,20}$. Therefore, mass density and number density fields can be treated as equivalent and so do momentum and velocity fields.

The density, velocity and stress fields are defined differently by MZF. Their relations to the SWB definitions are given by

$$ \rho^{\text{MZF}}(\boldsymbol{r}, t) = \rho^{\text{SWB}}(\boldsymbol{r}, t)/(m\sqrt{N}), $$

$$ \boldsymbol{v}^{\text{MZF}}(\boldsymbol{r}, t) = \boldsymbol{v}^{\text{SWB}}(\boldsymbol{r}, t)\sqrt{N}/V, $$

$$ \sigma_{\alpha\beta}^{\text{MZF}}(\boldsymbol{r}, t) = -\sigma_{\alpha\beta}^{\text{SWB}}(\boldsymbol{r}, t)/\sqrt{N}. $$

With the convention for the Fourier and inverse Fourier transforms employed by MZF$^{11,20}$,

$$ f(\boldsymbol{q}) = \int d^d \boldsymbol{r} f(\boldsymbol{r}) e^{i\boldsymbol{q}\cdot\boldsymbol{r}}, \qquad (A.1) $$

$$ f(\boldsymbol{r}) = \frac{1}{V} \sum_{\boldsymbol{q}} f(\boldsymbol{q}) e^{-i\boldsymbol{q}\cdot\boldsymbol{r}}, \qquad (A.2) $$

the relations in wavevector space read

$$ \rho^{\text{MZF}}(\boldsymbol{q}, t) = V \rho^{\text{SWB}}(-\boldsymbol{q}, t)/(m\sqrt{N}), \qquad (A.3) $$

$$ \boldsymbol{v}^{\text{MZF}}(\boldsymbol{q}, t) = \boldsymbol{v}^{\text{SWB}}(-\boldsymbol{q}, t)\sqrt{N}, \qquad (A.4) $$

$$ \sigma_{\alpha\beta}^{\text{MZF}}(\boldsymbol{q}, t) = -V\sigma_{\alpha\beta}^{\text{SWB}}(-\boldsymbol{q}, t)/\sqrt{N}. \qquad (A.4) $$

The rate of change of the velocity field introduces the stress tensor $\sigma_{\alpha\beta}$ through the conservation law which reads in Fourier space$^{27,29,30}$

$$ \frac{\partial v_{\alpha}^{\text{MZF}}(\boldsymbol{q}, t)}{\partial t} = i \mathcal{L} v_{\alpha}(\boldsymbol{q}, t) = i q_{\beta} \sigma_{\alpha\beta}^{\text{MZF}}(\boldsymbol{q}, t)/m. \qquad (A.5) $$

The ensemble average of this (symmetric) stress tensor gives

$$ \langle \sigma_{\alpha\beta}^{\text{MZF}}(\boldsymbol{q}) \rangle = \frac{V}{\sqrt{N}} p \delta_{\alpha\beta} \delta_{\boldsymbol{q},0}, \qquad (A.6) $$

where $p$ is the equilibrium pressure of the system.

Note that SWB$^{12,23,25}$ introduce the stress tensor through the fundamental momentum equation (11) which corresponds to a force balance. The resulting stress is the negative of a pressure [cf. Eq. (14)]. This is the usual definition taken in continuum mechanics$^{28}$. MZF$^{11,20}$, on the other hand, introduce the stress
tensor via the momentum conservation law which has the form of a continuity equation. The resulting stress has the same sign as the pressure [cf. Eq. (A.6)]. Identifying the stress tensor with the pressure tensor is the definition often employed in liquid state theory<sup>29,30</sup>.

This concludes our comparative overview of the different def-
initions. In the present work, we adopt the SWB definitions
and thus drop the superscript SWB. There is one difference to
Refs.<sup>12,23,25</sup>, though, as all correlation functions are defined as
⟨δA(t')δB(t+t')⟩ here, while SWB used the equivalent conven-
tion ⟨δA(t+t')δB(t')⟩.

Acknowledgments

We dedicate this work to the memory of Stefan Egelhaaf, a great scientist and a dear colleague. We are deeply grateful for the long, fruitful and always amiable collaboration we could enjoy with him. A.N. and J.B. are particularly grateful to J. Wittmer (Strasbourg) for fruitful collaboration and many helpful discussions on stress correlations. N.G., A. Z. and M.F. thank M. Maier for helpful discussions and earlier collaboration. N.G. acknowledges funding from the Deutsche Forschungsgemeinschaft (DFG) via SFB 1432, project C07. A.N. and J.B. acknowledge funding from the French National Agency for Research (ANR-24-CE30-7668).

Conflicts of interest

There are no conflicts to declare.

Data availability

Data for this article and its evaluation is available at https://doi.org/10.5281/zenodo.15113359.

Notes and references

1 N. Wagner and J. Mewis, *Theory and Applications of Colloidal Suspension Rheology*, Cambridge University Press, 2021.

2 R. G. Larson, *The structure and rheology of complex fluids*, Oxford University Press, New York, 1999.

3 K. Binder and W. Kob, *Glassy Materials and Disordered Solids*, WORLD SCIENTIFIC, Revised edn, 2011.

4 W. Götze, *Complex Dynamics of Glass-Forming Liquids, A Mode-Coupling Theory*, Oxford University Press, 2009.

5 T. Sentjabrskaja, M. Laurati and S. U. Egelhaaf, *Eur. Phys. J.: Spec. Top.*, 2017, **226**, 3023–3037.

6 J. C. Dyre, *J. Phys. Chem. Lett.*, 2024, **15**, 1603–1617.

7 J. C. Maxwell, *Philos. Trans. R. Soc. A*, 1867, **157**, 49.

8 D. Evans, *Phys. Rev. A*, 1981, **23**, 2622.

9 D. J. Evans and G. Morriss, *Statistical mechanics of nonequilibrium liquids*, Cambridge University Press, London, 2008.

10 J.-P. Hansen and I. R. McDonald, *Theory of Simple Liquids*, Academic Press, London, 2nd edn, 1986.

11 M. Maier, A. Zippelius and M. Fuchs, *Phys. Rev. Lett.*, 2017, **119**, 265701.

12 L. Klochko, J. Baschnagel, J. P. Wittmer and A. N. Semenov, *Soft Matter*, 2018, **14**, 6835.

13 A. Lemaître, *J. Chem. Phys.*, 2018, **149**, 104107.

14 L. Klochko, J. Baschnagel, J. P. Wittmer, H. Meyer, O. Benzer-
ara and A. N. Semenov, *J. Chem. Phys.*, 2022, **156**, 164505.

15 F. Vogel, A. Zippelius and M. Fuchs, *EPL*, 2019, **125**, 68003.

16 K. W. Desmond and E. R. Weeks, *Phys. Rev. Lett.*, 2015, **115**, 098302.

17 H. Tong, S. Sengupta and H. Tanaka, *Nat. Comm.*, 2020, **11**, 4863.

18 D. Steffen, L. Schneider, M. Müller and J. Rottler, *J. Chem. Phys.*, 2022, **157**, 064501.

19 J. Eshelby, *Proc. R. Soc. London A*, 1957, **241**, 376.

20 M. Maier, A. Zippelius and M. Fuchs, *J. Chem. Phys.*, 2018, **149**, 084502.

21 F. Vogel and M. Fuchs, *Eur. Phys. J. E*, 2020, **43**, 70.

22 M. Maier, A. Zippelius and M. Fuchs, *J. Chem. Phys.*, 2024, **161**, 199901.

23 J. P. Wittmer, A. N. Semenov and J. Baschnagel, *Phys. Rev. E*, 2023, **108**, 015002.

24 J. P. Wittmer, A. N. Semenov and J. Baschnagel, *Soft Matter*, 2023, **19**, 6140.

25 A. N. Semenov and J. Baschnagel, *Polymers*, 2024, **16**, 2336.

26 N. Grimm, M. von Bischopinck, A. Zumbusch and M. Fuchs, *J. Chem. Phys.*, 2024, **161**, 144118.

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

28 E. B. Tadmor and R. E. Miller, *Modeling Materials: Continuum, Atomistic and Multiscale Techniques*, Cambridge University Press, Cambridge, 2014.

29 U. Balucani and M. Zoppi, *Dynamics of the Liquid State*, Oxford University Press, Oxford, 2003.

30 J. P. Hansen and I. R. McDonald, *Theory of Simple Liquids*, Academic Press, London, 1986.

31 I. Goldhirsch, *Granular Matter*, 2010, **12**, 239–252.

32 J. Farago, H. Meyer, J. Baschnagel and A. N. Semenov, *Phys. Rev. E*, 2012, **85**, 051807.

33 A. N. Semenov, J. Farago and H. Meyer, *J. Chem. Phys.*, 2012, **136**, 244905.

34 K. Shi, E. R. Smith, E. E. Santiso and K. E. Gubbins, *J. Chem. Phys.*, 2023, **158**, 040901.

35 P. C. Martin, O. Parodi and P. S. Pershan, *Phys. Rev. A*, 1972, **6**, 2401.

36 A. Lemaître, *Phys. Rev. E*, 2017, **96**, 052101.

37 A. Lemaître, *J. Chem. Phys.*, 2015, **143**, 164515.

38 J. P. Wittmer and J. Baschnagel, *Phys. Rev. E*, 2024, **110**, 044604.

39 A. Lemaître, *J. Chem. Phys.*, 2018, **149**, 104107.

40 M. Rubinstein and R. H. Colby, *Polymer Physics*, Oxford University Press, Oxford, 2003.

41 L. Klochko, J. Baschnagel, J. P. Wittmer and A. N. Semenov, *Soft Matter*, 2021, **17**, 7867.

42 I. Rubinstein and L. Rubinstein, *Partial differential equations in classical mathematical physics*, Cambridge University Press, Cambridge, UK, 1998.

43 D. Forster, *Hydrodynamic Fluctuations, Broken Symmetry, and Correlation Functions*, CRC Press, Boca Raton, 2018.
44 R. Zwanzig, *Annual Review of Physical Chemistry*, 1965, **16**,

45 H. Mori, *Progr. Theoret. Phys.*, 1965, **33**, 423-455.

46 W. Götze and A. Latz, *J. Phys. Condens. Matter*, 1989, **1**, 4169.

47 P. Martin, *Statistical Mechanics of Equilibrium and Non-Equilibrium*, North-Holland Publ., Amsterdam, 1965, p. 100.

48 L. Kadanoff and P. Martin, *Annals of Physics*, 1963, **24**, 419-469.

49 E. Flenner, M. Zhang and G. Szamel, *Phys. Rev. E*, 2011, **83**, 051501.

50 G. Brambilla, D. El Masri, M. Pierno, L. Berthier, L. Cipelletti, G. Petekidis and A. B. Schofield, *Phys. Rev. Lett.*, 2009, **102**, 085703.

51 F. Weysser and D. Hajnal, *Phys. Rev. E*, 2011, **83**, 041503.

52 J. Song, N. Holten-Andersen and G. H. Mckinley, *Soft Matter*,

2023, **19**, 7885-7906.

53 I. Goldhirsch and C. Goldenberg, *Eur. Phys. Jour. E*, 2002, **9**, 245.

54 R. Bird, R. Armstrong and O. Hassager, *Dynamics of Polymeric Fluids*, Wiley, New York, 1987.

55 P. M. Chaikin and T. C. Lubensky, *Principles of condensed matter physics*, Cambridge University Press, Cambridge, 1997.

56 C. L. Klix, F. Ebert, F. Weysser, M. Fuchs, G. Maret and P. Keim, *Phys. Rev. Lett.*, 2012, **109**, 178301.

57 S. Kim and S. J. Karrila, *Microhydrodynamics*, Dover Publications, 2005.

58 R. B. Guenther and E. A. Thomann, *Journal of Mathematical Fluid Mechanics*, 2007, **9**, 489-505.

59 G. Picard, A. Ajdari, F. Lequeux and L. Bocquet, *Eur. Phys. J. E*, 2004, **15**, 371.