# Noether invariance theory for the equilibrium force structure of soft matter

Sophie Hermann  Theoretische Physik II, Physikalisches Institut, Universität
Bayreuth, D-95447 Bayreuth, Germany  Florian Sammüller  Theoretische Physik
II, Physikalisches Institut, Universität Bayreuth, D-95447 Bayreuth, Germany
Matthias Schmidt  Theoretische Physik II, Physikalisches Institut, Universität
Bayreuth, D-95447 Bayreuth, Germany<EMAIL_ADDRESS>
(26 January 2024)

###### Abstract

We give details and derivations for the Noether invariance theory that
characterizes the spatial equilibrium structure of inhomogeneous classical
many-body systems, as recently proposed and investigated for bulk systems [F.
Sammüller et al., Phys. Rev. Lett. 130, 268203 (2023)]. Thereby an intrinsic
thermal symmetry against a local shifting transformation on phase space is
exploited on the basis of the Noether theorem for invariant variations. We
consider the consequences of the shifting that emerge at second order in the
displacement field that parameterizes the transformation. In a natural way the
standard two-body density distribution is generated. Its second spatial
derivative (Hessian) is thereby balanced by two further and different two-body
correlation functions, which respectively introduce coarse-grained force
correlations and force gradients in a systematic and microscopically sharp way
into the framework. Seperate exact self and distinct sum rules hold expressing
this balance. We exemplify the validity of the theory on the basis of computer
simulations for the Lennard-Jones gas, liquid, and crystal, the Weeks-
Chandler-Andersen fluid, monatomic Molinero-Moore water at ambient conditions,
a three-body-interacting colloidal gel former, the Yukawa and soft-sphere
dipolar fluids, and for isotropic and nematic phases of Gay-Berne particles.
We describe explicitly the derivation of the sum rules based on Noether’s
theorem and also give more elementary proofs based on partial phase space
integration following Yvon’s theorem.

##  I Introduction

The spatial structure of soft matter can be very rich and varied. For a fluid
at the most fundamental level the bulk structure is described by the pair (or
radial) distribution function $g(r)$, which is a measure of the probability of
finding two particles separated by a distance $r$ hansen2013 ; frenkel2023book
. The information contained in $g(r)$ is analogously contained in the the
static structure factor $S(k)$, which is directly connected to the Fourier
transform of $g(r)$ and is accessible in scattering experiments. Real-space
measurements of pair correlation functions in colloidal systems give much
insight into the mesoscopic particle structure royall2007gofr ;
thorneywork2014 ; statt2016 .

The behaviour of $g(r)$ can be systematically analyzed and classified by
asymptotic expansion at large interparticle distances, which allows to
identify different types of decay evans1993decay ; evans1994decay ;
dijkstra2000decay ; grodon2004decay ; cats2021decay . The asymptotic analysis
is fundamental to the behaviour of the system at hand, as the results remain
relevant for the decay of the one-body density profile into the bulk far away
from a substrate or solute that created the density inhomogeneity. Furthermore
$g(r)$ plays the prominent role to determine the thermodynamics and equation
of state via spatial integration according to the virial or compressibility
routes to the pressure as are elementary ingredients of liquid integral
equation theory hansen2013 .

Despite these virtues, there is much activity of going beyond $g(r)$ in the
spatial characterization of soft matter. Recent efforts include characterizing
the structure of liquids using four-point correlation functions zhang2020pnas
. Such measures were argued to be relevant for the dynamics in dense colloidal
liquids singh2023pnas . The packing efficiency of binary hard sphere systems
was related to fluid structure at intermediate ranges yuan2021prl . Emergent
structural correlations in dense liquids were characterized up to the six-body
structure factor pihlajamaa2023 . Using machine learning the averaged
structural features centered around nearby particles were used to predict
dynamics in supercooled liquids boattini2011 .

Noether’s theorem of invariant variations noether1918 ; byers1998  was put in
a statistical physics setting in a number of different ways baez2013markov ;
marvian2014quantum ; sasa2016 ; sasa2019 ; revzen1970 ; baez2020bottom ;
bravetti2023 . The recent thermal Noether invariance theory hermann2021noether
; hermann2022topicalReview ; hermann2022variance ; hermann2022quantum ;
tschopp2022forceDFT ; sammueller2023whatIsLiquid ; robitschko2023any is based
on the identification of the symmetry of the statistical many-body systems
against specific shifting and rotation operations on phase space.
Corresponding force and torque correlation functions are generated and exact
statistical mechanical identities (“sum rules”) play the role that in more
conventional uses of the Noether theorem are played by the resulting
conservation laws. The recent hyperforce approach by Robitschko et al.
robitschko2023any carries the thermal Noether concept further in that it
allows to address the behaviour of arbitrary observables in equilibribum.

Figure 1: Illustrative examples of the Noether force correlation functions for
the Lennard-Jones system in the fcc crystal phase (first colunm), the liquid
phase (second column), and the gas phase (third column). Shown are results for
the standard pair correlation function $g(r)$ in all three phases (first row).
The Noether invariance theory introduces two further pair correlation
functions, ${\sf g}_{\nabla\\!f}(r)$ and ${\sf g}_{f\\!f}(r)$, as are defined
in detail in Sec. VI. These tensorial functions possess radial ($\parallel$)
and transversal ($\perp$) components, both of which are plotted. Briefly, the
two-body force-gradient correlator ${\sf g}_{\nabla\\!f}(r)$ (second row) is a
measure of the change of the force on one particle upon changing the position
of a second particle that is located at a distance $r$ from the first
particle. The force-force pair correlator ${\sf g}_{f\\!f}(r)$ (third row)
associates the individual forces in the particle pair with each other. That
the dotted lines (third panel) coincide with the full lines verifies the
3g-rule (1). In the figure the respective vertical scale factor $S$ is given
in the top left corner of each panel and the scaled values for bulk density
$\rho_{b}$ and temperature $T$ are indicated for each statepoint above the
respective column. The vertical gray lines indicate the position of the first
maximum of $g(r)$ as a guide to the eye. Adapted from the Supplementary
Material of Ref. sammueller2023whatIsLiquid .

The classical question “What is liquid?” barker1976 appears in new light when
regarded from the viewpoint of the thermal Noether invariance theory. In
particular the recently formulated force-force and force-gradient pair
correlation functions sammueller2023whatIsLiquid give fresh insight into the
correlation structure. These correlation functions are generated from
canonical “shifting” transformation on phase space goldstein2002 and thereby
considering the second order in the spatial shifting field. The shifting field
is an important formal device for formulating and exploiting the Noether
invariance, yet the resulting identities are free of the shifting field. Only
genuine system properties are correlated with each other and the framework is
entirely mechanical. In particular forces are at the forefront, as these
appear naturally in response to spatial displacement (shifting) of the free
energy. Forces comprise interparticle, external, and diffusive contributions
and they occur both globally as well as locally resolved in arbitrary
spatially inhomogeneous systems.

For bulk systems the second-order Noether invariance reduces to the “3g-rule”
that relates three different types of pair functions with each other (we
reproduce from Sec. VI):

$\displaystyle\nabla\nabla g(r)$ $\displaystyle+{\sf g}_{\nabla\\!f}(r)+{\sf
g}_{f\\!f}(r)=0.$ (1)

Here $\nabla$ indicates the spatial gradient, such that $\nabla\nabla g(r)$ is
the Hessian of the pair correlation function. Furthermore the force-gradient
pair correlator ${\sf g}_{\nabla\\!f}(r)$ and the force-force pair correlator
${\sf g}_{f\\!f}(r)$ are both $3\times 3$-matrices. These correlation
functions also address particle pairs, but they incorporate additional
information over the mere counting in $g(r)$. In particular ${\sf
g}_{f\\!f}(r)$ correlates the interparticle forces that act on the two
particles in the pair with each other. The force-gradient correlation function
${\sf g}_{\nabla\\!f}(r)$ is the mean gradient of the force on one of the
particles with respect to changing the position of the partner. A more
detailed description is given below.

The Noether 3g-rule (1) was verified by Sammüller et al.
sammueller2023whatIsLiquid across a wide range of model Hamiltonians,
including Lennard-Jones, Yukawa, soft-sphere dipolar, Stockmayer, Gay-Berne,
and Weeks-Chandler-Andersen liquids. Also Hamiltonians which incorporate
three-body interaction terms were considered, such as the monatomic Molinero-
Moore water molinero2009 ; coe2022water and a colloidal gel forming model
proposed by Kob and coworkers saw2009 ; saw2011 ; sammueller2023gel . The
identity (1) holds beyond fluids also in liquid crystals and in solids.
Besides certain very generic features, such as a depletion at small
separation, the individual models and phases display differences in shapes and
forms of $g(r)$, as expected. The observed differences in ${\sf g}_{f\\!f}(r)$
and in ${\sf g}_{\nabla\\!f}(r)$ across different models and phases are
significantly more pronounced though and they are indicative of a wide range
of specific features, including the presence of interparticle attraction and
of particle strand formation in colloidal gels jadrich2023preface and in
dipolar fluids klapp2005topRev ; stevens1995 ; teixeira2000topRev ;
tavares1999 ; allen2019topRev .

To illustrate this behaviour we reproduce in Figure 1 results of Ref.
sammueller2023whatIsLiquid for the Lennard-Jones system in all of its three
stable thermodynamic phases: crystal, liquid, and gas. The data is obtained
from equilibribum sampling via adaptive Brownian dynamics sammueller2021 ,
which allows for tight error control of the occurring forces. The behaviour of
$g(r)$ in the gas can be well-understood from the virial expansion where to
lowest order $g(r)=\exp(-\beta\phi(r))$, where $\beta$ denotes inverse
temperature and $\phi(r)$ is the (Lennard-Jones) pair potential.

The force-force and force-gradient pair correlation functions feature
tensorial components that are radial and transversal with respect to the
difference vector between both particle centers. The mean force gradient ${\sf
g}_{\nabla\\!f}(r)$ vanishes beyond the range of the (truncated) Lennard-Jones
potential (second row in Fig. 1) and it features a clear sign of the
interparticle attraction in the transversal component; see the positive peak
of the pink curve. The force-force pair correlation function ${\sf
g}_{f\\!f}(r)$ (third row in Fig. 1) extends beyond the range of truncation
and it again shows very rich correlation structure. In the symmetry-broken
state of the fcc crystal the results are obtained from only resolving
according to scalar distances. While this implies loss of information over the
full inhomogeneous two-body resolution, the results shown in Fig. 1 for the
Lennard-Jones crystal nevertheless indicate intricate and long-ranged force-
force correlation behaviour. In striking contrast, the force-gradient
correlation function ${\sf g}_{\nabla\\!f}(r)$ remains short-ranged and it
strictly vanishes beyond the cutoff for short-ranged or truncated pair
potentials.

With the central mathematical ideas only being sketched in Ref.
sammueller2023whatIsLiquid , we here give a detailed account of the second-
order thermal Noether invariance theory. We present explicit derivations as
well as several new results, such as the density-force correlation sum rule
(35). We also give proofs of the Noether identities from partial integration
on phase space in analogy to the Yvon theorem hansen2013 . This elementary
route is conceptually simple, as it only requires direct manipulation of the
phase space integrals, and it hence provides a valuable consistency check. We
also describe in detail the simulation results for the variety of models
considered.

The paper is organized as follows. In Sec. II we lay out the statistical many-
body physics under investigation and describe the theoretical setup of thermal
Noether invariance against local shifting. The two-body sum rules that follow
from the second-order invariance are presented in Sec. III. A general density-
force correlation sum rule is proven in Sec. IV. The emerging force
correlators are split into separate contributions in Sec. V, which facilitates
to prove the general Noether two-body sum rules from first principles. The
general inhomogeneous two-body framework is specialized to the case of bulk
systems in Sec. VI and we present and discuss our simulation results therein.
Several important yet technical points are shown in five Appendices. This
includes considerations of the kinetic stress autocorrelations presented in
Appendix A. Kinetic energy curvature terms are addressed in Appendix B and
potential energy curvature contributions are presented in Appendix C.
Alternative explicit proofs of the Noether sum rules based on partial
integration are given in Appendix D. Details about measuring correlation
functions in simulations are given in Appendix E. We present our conclusions
in Sec. VII.

##  II Theoretical setup

We consider systems of $N$ particles with position coordinates
${\mathbf{r}}_{1},\ldots,{\mathbf{r}}_{N}\equiv{\mathbf{r}}^{N}$ and linear
momenta ${\bf p}_{1},\ldots,{\bf p}_{N}\equiv{\bf p}^{N}$. The Hamiltonian $H$
is taken to contain kinetic energy, as well as interparticle and external
energy contributions,

$\displaystyle H$ $\displaystyle=\sum_{i}\frac{{\bf
p}_{i}^{2}}{2m}+u({\mathbf{r}}^{N})+\sum_{i}V_{\rm ext}({\mathbf{r}}_{i}),$
(2)

where the summation index $i=1,\ldots,N$ ranges over all particles, $m$
denotes the particle mass, $u({\mathbf{r}}^{N})$ is the interparticle
interaction potential, and the one-body external potential $V_{\rm
ext}({\mathbf{r}})$ is given as a function of a single (generic) position
variable ${\mathbf{r}}$.

The canonical transformation goldstein2002 of the phase space coordinates that
was put forward in Refs. hermann2022quantum ; tschopp2022forceDFT ;
sammueller2023whatIsLiquid maps original onto new coordinates as well as
original onto new momenta according to the following map:

$\displaystyle{\mathbf{r}}_{i}$
$\displaystyle\to{\mathbf{r}}_{i}+{\bm{\epsilon}}({\mathbf{r}}_{i})=\tilde{\mathbf{r}}_{i},$
(3) $\displaystyle{\bf p}_{i}$
$\displaystyle\to[{\mathbbm{1}}+\nabla_{i}{\bm{\epsilon}}({\mathbf{r}}_{i})]^{-1}\cdot{\bf
p}_{i}=\tilde{\bf p}_{i},$ (4)

where the tilde indicates the new variables. The three-dimensional vector
field ${\bm{\epsilon}}({\mathbf{r}})$ parameterizes the transform,
${\mathbbm{1}}$ denotes the $3\times 3$-unit matrix, $\nabla_{i}$ indicates
the derivative with respect to ${\mathbf{r}}_{i}$, and matrix inversion is
indicated by the superscript $-1$. Here we consider systems in three space
dimensions but the theory is general. The transformation (3) and (4) preserves
the phase space volume element as well as thermal averages tschopp2022forceDFT
; hermann2022quantum ; sammueller2023whatIsLiquid ; goldstein2002 and it
constitutes a canonical transformation in the sense of classical mechanics
goldstein2002 . This property can e.g. be seen explicitly by considering a
generator of the form ${\cal G}=\sum_{i=1}^{N}\tilde{\bf
p}_{i}\cdot({\mathbf{r}}_{i}+{\bm{\epsilon}}({\mathbf{r}}_{i}))$, where as
before the tilde indicates the new phase space variables. Applying the generic
transformation relations $\tilde{\mathbf{r}}_{i}=\partial{\cal
G}/\partial\tilde{\bf p}_{i}$ and ${\bf p}_{i}=\partial{\cal
G}/\partial{\mathbf{r}}_{i}$ goldstein2002 and solving for the new variables
then yields the transformation equations (3) and (4).

We Taylor expand the transformation in the shifting field
${\bm{\epsilon}}({\mathbf{r}})$. The coordinate transformation (3) is already
linear in ${\bm{\epsilon}}({\mathbf{r}})$ and it is hence unaffected. The
momentum transformation (4) can be expanded as a Neumann series, which is the
analog of the geometric series for matrices or, more generally, for bounded
operators. Assuming that both the shifting field and its gradient are small,
up to second order the expansion yields:

$\displaystyle[{\mathbbm{1}}+\nabla_{i}{\bm{\epsilon}}({\mathbf{r}}_{i})]^{-1}={\mathbbm{1}}-\nabla_{i}{\bm{\epsilon}}({\mathbf{r}}_{i})+[\nabla_{i}{\bm{\epsilon}}({\mathbf{r}}_{i})]^{2}-\ldots,$
(5)

where the square on the right hand side denotes matrix (self-)multiplication,
i.e.,
$[\nabla_{i}{\bm{\epsilon}}({\mathbf{r}}_{i})]^{2}=[\nabla_{i}{\bm{\epsilon}}({\mathbf{r}}_{i})]\cdot[\nabla_{i}{\bm{\epsilon}}({\mathbf{r}}_{i})]$.
Verifying Eq. (5) to second order in the shifting gradient is straightforward
via matrix multiplying its both sides by
${\mathbbm{1}}+\nabla_{i}{\bm{\epsilon}}({\mathbf{r}}_{i})$, which per
definition of the inverse gives ${\mathbbm{1}}$ on the left hand side. The
right hand side gives the same result up to second order upon carrying out the
matrix multiplications and simplifying.

When expressed in the new variables, the Hamiltonian becomes functionally
dependent on the shifting field, $H\to H[{\bm{\epsilon}}]$. We first consider
the term linear in ${\bm{\epsilon}}({\mathbf{r}})$ hermann2022quantum ;
tschopp2022forceDFT ; sammueller2023whatIsLiquid . One sees that the position-
resolved one-body force density operator $\hat{\bf F}({\mathbf{r}})$ follows
from functional differentiation according to:

$\displaystyle-\frac{\delta
H[{\bm{\epsilon}}]}{\delta{\bm{\epsilon}}({\mathbf{r}})}\Big{|}_{{\bm{\epsilon}}=0}$
$\displaystyle=\hat{\bf F}({\mathbf{r}}),$ (6)

where $\delta/\delta{\bm{\epsilon}}({\mathbf{r}})$ indicates the functional
derivative with respect to the shifting field ${\bm{\epsilon}}({\mathbf{r}})$;
see e.g. Ref. schmidt2022rmp for an overview of functional differentiation
techniques. After the derivative is taken, ${\bm{\epsilon}}({\mathbf{r}})$ is
set to zero, as indicated in the notation on the left hand side of Eq. (6).
The structure of the force density operator $\hat{\bf F}({\mathbf{r}})$
mirrors that of the Hamiltonian (2) in that it contains kinetic,
interparticle, and external contributions:

$\displaystyle\hat{\bf F}({\mathbf{r}})$
$\displaystyle=\nabla\cdot\hat{\bm{\tau}}({\mathbf{r}})+\hat{\bf F}_{\rm
int}({\mathbf{r}})-\hat{\rho}({\mathbf{r}})\nabla V_{\rm ext}({\mathbf{r}}),$
(7)

where $\nabla$ denotes the derivative with respect to ${\mathbf{r}}$. The one-
body operators for the kinetic stress $\hat{\bm{\tau}}({\mathbf{r}})$, the
interparticle force density $\hat{\bf F}_{\rm int}({\mathbf{r}})$
schmidt2022rmp , and the microscopic density distribution
$\hat{\rho}({\mathbf{r}})$ evans1979 ; hansen2013  are given respectively by

$\displaystyle\hat{\bm{\tau}}({\mathbf{r}})$
$\displaystyle=-\sum_{i}\frac{{\bf p}_{i}{\bf
p}_{i}}{m}\delta({\mathbf{r}}-{\mathbf{r}}_{i}),$ (8) $\displaystyle\hat{\bf
F}_{\rm int}({\mathbf{r}})$
$\displaystyle=-\sum_{i}\delta({\mathbf{r}}-{\mathbf{r}}_{i})\nabla_{i}u({\mathbf{r}}^{N}),$
(9) $\displaystyle\hat{\rho}({\mathbf{r}})$
$\displaystyle=\sum_{i}\delta({\mathbf{r}}-{\mathbf{r}}_{i}),$ (10)

where $\delta(\cdot)$ indicates the (three-dimensional) Dirac distribution and
${\bf p}_{i}{\bf p}_{i}$ is the momentum of particle $i$ dyadically multiplied
with itself.

The following analysis within Statistical Mechanics holds equally well with
fixed number of particles, i.e., in the canonical ensemble. To be explicit, we
here formulate our theory in the grand ensemble at temperature $T$ and
chemical potential $\mu$. The grand potential $\Omega$ and the grand partition
sum $\Xi$ are given respectively by

$\displaystyle\Omega$ $\displaystyle=-k_{B}T\ln\Xi,$ (11) $\displaystyle\Xi$
$\displaystyle={\rm Tr}\,{\rm e}^{-\beta(H-\mu N)},$ (12)

where $k_{B}$ is the Boltzmann constant, $\beta=1/(k_{B}T)$ denotes inverse
temperature, and the classical “trace” operation in the grand ensemble is
given by ${\rm Tr}\,\cdot=\sum_{N=0}^{\infty}(N!h^{3N})^{-1}\int
d{\mathbf{r}}_{1}\ldots d{\mathbf{r}}_{N}\int d{\bf p}_{1}\ldots d{\bf
p}_{N}\cdot$, where $h$ indicates the Planck constant. The corresponding grand
probability distribution (Gibbs measure) is $\Psi={\rm e}^{-\beta(H-\mu
N)}/\Xi$. Thermal averages are defined via $\langle\cdot\rangle={\rm
Tr}\,\Psi\cdot$, such that the density profile is the average of the one-body
density operator and written as
$\rho({\mathbf{r}})=\langle\hat{\rho}({\mathbf{r}})\rangle$. Correspondingly
the average kinetic stress is
${\bm{\tau}}({\mathbf{r}})=\langle\hat{\bm{\tau}}({\mathbf{r}})\rangle$ and
the mean interparticle force density is ${\bf F}_{\rm
int}({\mathbf{r}})=\langle\hat{\bf F}_{\rm int}({\mathbf{r}})\rangle$.

At face value the grand partition sum (12) carries an apparent functional
dependence on the shifting field ${\bm{\epsilon}}({\mathbf{r}})$, via the
occurrence of ${\bm{\epsilon}}({\mathbf{r}})$ in the transformed Hamiltonian
$H[{\bm{\epsilon}}]$ hermann2022quantum ; tschopp2022forceDFT ;
sammueller2023whatIsLiquid . We hence also have a functional dependence of the
grand partition sum, $\Xi[{\bm{\epsilon}}]$, and consequently of the grand
potential, $\Omega[{\bm{\epsilon}}]$, as it is generated via the standard
relationship (11). We functionally Taylor expand the grand potential in the
shifting field to quadratic order, which gives

$\displaystyle\Omega[{\bm{\epsilon}}]$ $\displaystyle=\Omega+\int
d{\mathbf{r}}\frac{\delta\Omega[{\bm{\epsilon}}]}{\delta{\bm{\epsilon}}({\mathbf{r}})}\Big{|}_{{\bm{\epsilon}}=0}\cdot{\bm{\epsilon}}({\mathbf{r}})$
(13) $\displaystyle\qquad+\frac{1}{2}\int
d{\mathbf{r}}d{\mathbf{r}}^{\prime}\frac{\delta^{2}\Omega[{\bm{\epsilon}}]}{\delta{\bm{\epsilon}}({\mathbf{r}})\delta{\bm{\epsilon}}({\mathbf{r}}^{\prime})}\Big{|}_{\rm{\bm{\epsilon}}=0}:{\bm{\epsilon}}({\mathbf{r}}^{\prime}){\bm{\epsilon}}({\mathbf{r}})+\ldots.$

The colon denotes a double tensor contraction and $\Omega$, with no argument,
indicates the original grand potential (11) with no transformation applied,
$\Omega=\Omega[0]$, where the functional argument $0$ indicates the special
case of vanishing shifting field, ${\bm{\epsilon}}({\mathbf{r}})=0$.

Despite the apparent dependence on ${\bm{\epsilon}}({\mathbf{r}})$, Noether
invariance hermann2021noether ; hermann2022topicalReview  ascertains that the
grand potential remains unchanged upon applying the transformation, and hence

$\displaystyle\Omega[{\bm{\epsilon}}]=\Omega,$ (14)

which holds true for any (permissible) form of
${\bm{\epsilon}}({\mathbf{r}})$. As a consequence, all terms in the functional
expansion (13) need to vanish identically. At first order in the displacement
field, we hence have

$\displaystyle\frac{\delta\Omega[{\bm{\epsilon}}]}{\delta{\bm{\epsilon}}({\mathbf{r}})}\Big{|}_{{\bm{\epsilon}}=0}$
$\displaystyle=0.$ (15)

Explicitly carrying out the functional derivative on the left hand of Eq. (15)
gives tschopp2022forceDFT ; hermann2022quantum ; sammueller2023whatIsLiquid :

$\displaystyle\frac{\delta\Omega[{\bm{\epsilon}}]}{\delta{\bm{\epsilon}}({\mathbf{r}})}\Big{|}_{{\bm{\epsilon}}=0}$
$\displaystyle=\Big{\langle}\frac{\delta
H[{\bm{\epsilon}}]}{\delta{\bm{\epsilon}}({\mathbf{r}})}\Big{|}_{{\bm{\epsilon}}=0}\Big{\rangle}=-\langle\hat{\bf
F}({\mathbf{r}})\rangle,$ (16)

where we recall the total force density operator $\hat{\bf F}({\mathbf{r}})$
in the form (7). The first equality in Eq. (16) follows from applying the
definition (11) of the grand potential. The second equality is the thermal
average of the total force density operator identity (6).

The expression on the right hand side of Eq. (16) constitutes the negative
averaged one-body force density $-{\bf F}({\mathbf{r}})=-\langle\hat{\bf
F}({\mathbf{r}})\rangle$. Due to the invariance that is expressed via Eq. (15)
we can conclude that the mean one-body force density vanishes in equilibribum,

$\displaystyle{\bf F}({\mathbf{r}})$ $\displaystyle=0.$ (17)

The explicit form (7) of $\hat{\bf F}({\mathbf{r}})$ allows to decompose the
left hand side of Eq. (17) into a sum of ideal gas, interparticle, and
external terms according to:

$\displaystyle-k_{B}T\nabla\rho({\mathbf{r}})+{\bf F}_{\rm
int}({\mathbf{r}})-\rho({\mathbf{r}})\nabla V_{\rm ext}({\mathbf{r}})$
$\displaystyle=0.$ (18)

The force density balance (18) is the analog of the first member of the
classical one-body Yvon-Born-Green equation of liquid state theory hansen2013
, in which ${\bf F}_{\rm int}({\mathbf{r}})$ is further expressed in integral
form.

The present functional route tschopp2022forceDFT ; hermann2022quantum ;
sammueller2023whatIsLiquid identifies the equilibribum force density balance
(18) as the first-order Noether invariance sum rule of the grand potential
with respect to the local shifting (3) and (4) of phase space. This
transformation is canonical and hence the associated symmetry forms a deeply
rooted property of the statistical physics generated by the standard
Hamiltonian (2).

As an aside, the recent force-based density functional theory
tschopp2022forceDFT ; sammueller2022forceDFT takes this concept as the
fundamental building block for making predictions based on an excess (over
ideal gas) free energy density functional. Thereby the interparticle force
density is re-written according to the Yvon-Born-Green hierarchy yvon1935 ;
Born1946 ; hansen2013  as an integral over the pair force times the two-body
density, which in turn is obtained from numerical solution of the
inhomogeneous two-body Ornstein-Zernike equation. Force-density functional
theory has been exemplified for inhomogeneous hard sphere systems in planar
geometry by inputting the two-body direct correlation functional from
fundamental measure theory rosenfeld1989 ; roth2010 , see Refs.
tschopp2022forceDFT ; sammueller2022forceDFT .

##  III Local second-order invariance

The previous Sec. II summarizes the theory developed in Refs.
tschopp2022forceDFT ; hermann2022quantum , which we require for considering
the second order. We recall that for global shifting with spatially uniform
displacement, ${\bm{\epsilon}}({\mathbf{r}})={\bm{\epsilon}}_{0}=\rm const$,
addressing the second order in ${\bm{\epsilon}}_{0}$ yields a sum rule which
relates the global force variance with the mean global potential energy
curvature hermann2022variance .

In this global case the interparticle force contributions vanish and only the
external force field contributes, according to the following global second-
order sum rule:

$\displaystyle\int
d{\mathbf{r}}d{\mathbf{r}}^{\prime}H_{2}({\mathbf{r}},{\mathbf{r}}^{\prime})\nabla
V_{\rm ext}({\mathbf{r}})$ $\displaystyle\nabla^{\prime}V_{\rm
ext}({\mathbf{r}}^{\prime})=$ $\displaystyle k_{B}T\int
d{\mathbf{r}}\rho({\mathbf{r}})\nabla\nabla V_{\rm ext}({\mathbf{r}}),$ (19)

where $\nabla^{\prime}$ denotes the derivative with respect to
${\mathbf{r}}^{\prime}$ and $H_{2}({\mathbf{r}},{\mathbf{r}}^{\prime})={\rm
cov}(\hat{\rho}({\mathbf{r}}),\hat{\rho}({\mathbf{r}}^{\prime}))$ is the
standard two-body correlation function of density fluctuations around the
local density profile hansen2013 ; schmidt2022rmp . Thereby the covariance of
two observables, as represented by phase space functions $\hat{A}$ and
$\hat{B}$, is defined in the standard way as ${\rm
cov}(\hat{A},\hat{B})=\langle\hat{A}\hat{B}\rangle-\langle\hat{A}\rangle\langle\hat{B}\rangle$.
The left hand side of Eq. (19) can upon carrying out the position integrals be
alternatively written as $\langle\hat{\bf F}_{\rm ext}^{0}\hat{\bf F}_{\rm
ext}^{0}\rangle$, where the global external force density operator is defined
as $\hat{\bf F}_{\rm ext}^{0}=-\sum_{i}\nabla_{i}V_{\rm
ext}({\mathbf{r}}_{i})$ and the covariance reduces to the mean of the product,
because of the vanishing mean external force in equilibribum, $\langle\hat{\bf
F}_{\rm ext}^{0}\rangle=0$ hermann2021noether .

Following the strategy of Ref. sammueller2023whatIsLiquid we here describe in
detail the consequences of locally resolved shifting at second order in the
displacement field. We recall that the functional expansion (13) of the grand
potential holds irrespective of the precise form of
${\bm{\epsilon}}({\mathbf{r}})$. Addressing the second order term, this
implies that its prefactor in Eq. (13) needs to vanish identically:

$\displaystyle\frac{\delta^{2}\Omega[{\bm{\epsilon}}]}{\delta{\bm{\epsilon}}({\mathbf{r}})\delta{\bm{\epsilon}}({\mathbf{r}}^{\prime})}\Big{|}_{{\bm{\epsilon}}=0}$
$\displaystyle=0.$ (20)

Making the functional derivative that appears on the left hand side more
explicit yields

$\displaystyle\frac{\delta^{2}\Omega[{\bm{\epsilon}}]}{\delta{\bm{\epsilon}}({\mathbf{r}})\delta{\bm{\epsilon}}({\mathbf{r}}^{\prime})}$
$\displaystyle=-\beta\,{\rm cov}\Big{(}\frac{\delta
H[{\bm{\epsilon}}]}{\delta{\bm{\epsilon}}({\mathbf{r}})},\frac{\delta
H[{\bm{\epsilon}}]}{\delta{\bm{\epsilon}}({\mathbf{r}}^{\prime})}\Big{)}+\Big{\langle}\frac{\delta^{2}H[{\bm{\epsilon}}]}{\delta{\bm{\epsilon}}({\mathbf{r}})\delta{\bm{\epsilon}}({\mathbf{r}}^{\prime})}\Big{\rangle}.$
(21)

The functional derivative of the Hamiltonian, $\delta
H[{\bm{\epsilon}}]/\delta{\bm{\epsilon}}({\mathbf{r}})$, can be expressed as
the negative force density operator via Eq. (6). Then inserting Eq. (21) into
Eq. (20), and grouping terms together allows to formulate the following
locally resolved two-body Noether sum rule:

$\displaystyle\beta\langle\hat{\bf F}({\mathbf{r}})\hat{\bf
F}({\mathbf{r}}^{\prime})\rangle$
$\displaystyle=\Big{\langle}\frac{\delta^{2}H[{\bm{\epsilon}}]}{\delta{\bm{\epsilon}}({\mathbf{r}})\delta{\bm{\epsilon}}({\mathbf{r}}^{\prime})}\Big{|}_{{\bm{\epsilon}}=0}\Big{\rangle}.$
(22)

The correlator on the left hand side of Eq. (22) is analogous to ${\rm
cov}(\hat{\bf F}({\mathbf{r}}),\hat{\bf
F}({\mathbf{r}}^{\prime}))=\langle\hat{\bf F}({\mathbf{r}})\hat{\bf
F}({\mathbf{r}}^{\prime})\rangle$, because the local mean force vanishes in
equilibrium, $\langle\hat{\bf F}({\mathbf{r}})\rangle=0$ hermann2022quantum ;
tschopp2022forceDFT ; we recall the individual force contributions in the
explicit form (18) of the force balance relationship.

Equation (22) constitutes an exact sum rule that relates the dyadic force-
force correlations (left hand side) at two different positions with the mean
curvature of the Hamiltonian with respect to local shifting (right hand side).
Hence the left hand side of this sum rule constitutes an immediately
meaningful standalone physical object, with clearcut statistical mechanical
interpretation. We demonstrate below that this is true also for the more
formal looking right hand side of Eq. (22). That the exact equality (22) holds
generally, at arbitrary positions ${\mathbf{r}}$ and ${\mathbf{r}}^{\prime}$,
could have certainly not been guessed easily a priori (we turn to partial
integration methods below in Appendix D).

We re-write Eq. (22) via multiplying by $\beta$ and then treating the kinetic
contributions separately (see below). For brevity of notation it is useful to
combine the interparticle and external forces into the following dedicated
potential force density operator:

$\displaystyle\hat{\bf F}_{U}({\mathbf{r}})=\hat{\bf F}_{\rm
int}({\mathbf{r}})-\hat{\rho}({\mathbf{r}})\nabla V_{\rm ext}({\mathbf{r}}),$
(23)

where the subscript $U$ refers to the total potential energy and we recall the
definition (9) of the interparticle force operator $\hat{\bf F}_{\rm
int}({\mathbf{r}})$. The total potential energy operator is given as the phase
space function $H_{U}=u({\mathbf{r}}^{N})+\sum_{i}V_{\rm
ext}({\mathbf{r}}_{i})$. Using the coordinate transformation (3) the potential
force density is obtained via functional differentiation according to $-\delta
H_{U}/\delta{\bm{\epsilon}}({\mathbf{r}})|_{{\bm{\epsilon}}=0}=\hat{\bf
F}_{U}({\mathbf{r}})$, i.e. the first functional derivative, cf. Eq. (6), of
the potential energy contribution to the Hamiltonian.

We identify distinct contributions (subscript “dist”) that arise on the left
hand side of Eq. (22). These terms are generated from pairs of particles with
unequal indices such that double sums reduce to
$\sum_{ij(\neq)}\equiv\sum_{i=1}^{N}\sum_{j=1,j\neq i}^{N}$. We can then
rewrite Eq. (22) as the following pair of exact distinct and self (see below)
identities:

$\displaystyle\Big{\langle}\beta\hat{\bf F}_{U}({\mathbf{r}})\beta\hat{\bf
F}_{U}({\mathbf{r}}^{\prime})\Big{\rangle}_{\\!\rm
dist}=\nabla\nabla^{\prime}\rho_{2}({\mathbf{r}},{\mathbf{r}}^{\prime})$
$\displaystyle\quad\qquad\qquad\quad+\Big{\langle}{\sum_{ij(\neq)}^{\phantom{space}}}\delta({\mathbf{r}}-{\mathbf{r}}_{i})\delta({\mathbf{r}}^{\prime}-{\mathbf{r}}_{j})\nabla_{i}\nabla_{j}\beta
u({\mathbf{r}}^{N})\Big{\rangle},$ (24)
$\displaystyle\Big{\langle}\beta\hat{\bf F}_{U}({\mathbf{r}})\beta\hat{\bf
F}_{U}({\mathbf{r}})\Big{\rangle}_{\\!\rm
self}=\nabla\nabla\rho({\mathbf{r}})+\rho({\mathbf{r}})\nabla\nabla\beta
V_{\rm ext}({\mathbf{r}})$
$\displaystyle\qquad\quad\qquad\quad\qquad\qquad+\Big{\langle}\sum_{i}\delta({\mathbf{r}}-{\mathbf{r}}_{i})\nabla_{i}\nabla_{i}\beta
u({\mathbf{r}}^{N})\Big{\rangle}.$ (25)

The two-body density on the right hand side of Eq. (24) consists of
contributions only from distinct particle pairs and it has the standard
definition hansen2013 :

$\displaystyle\rho_{2}({\mathbf{r}},{\mathbf{r}}^{\prime})$
$\displaystyle=\langle\hat{\rho}({\mathbf{r}})\hat{\rho}({\mathbf{r}}^{\prime})\rangle_{\rm
dist}$
$\displaystyle=\Big{\langle}\sum_{ij(\neq)}\delta({\mathbf{r}}-{\mathbf{r}}_{i})\delta({\mathbf{r}}^{\prime}-{\mathbf{r}}_{j})\Big{\rangle}.$
(26)

The self average on the left hand side of Eq. (25) only involves a single
delta function in position, i.e., it is defined such that
$\langle\beta\hat{\bf F}_{U}({\mathbf{r}})\beta\hat{\bf
F}_{U}({\mathbf{r}})\rangle_{\rm self}=\langle\sum_{i}[\nabla_{i}\beta
u({\mathbf{r}}^{N})+\nabla_{i}\beta V_{\rm
ext}({\mathbf{r}}_{i})][\nabla_{i}\beta u({\mathbf{r}}^{N})+\nabla_{i}\beta
V_{\rm ext}({\mathbf{r}}_{i})]\delta({\mathbf{r}}-{\mathbf{r}}_{i})\rangle$.

The sum rules (24) and (25) give much concrete relevance to the more abstract
form (22). In particular the curvature terms on the right hand sides are given
as explicit averages, which are accessible in many-body simulations. The
kinetic stress autocorrelations that emerge from Eq. (22) are addressed in
Appendix A and the energy curvature contributions are discussed in Appendix B.
We point the reader at precise points in the following argumentation to these
Appendices.

##  IV Density-force Noether identity

We next consider the correlator of the density operator with the total force
density operator, i.e. $\langle\hat{\rho}({\mathbf{r}})\hat{\bf
F}({\mathbf{r}}^{\prime})\rangle$. Besides being interesting in its own right,
this correlator is relevant for the left hand side of Eq. (22) via the
correlation of external and total force densities, which can be written as
$-\langle\hat{\rho}({\mathbf{r}})\hat{\bf
F}({\mathbf{r}}^{\prime})\rangle\beta\nabla V_{\rm ext}({\mathbf{r}})$.

We demonstrate two distinct routes to formulate a valid sum rule for the
density-force correlator. First we start from the locally resolved equilibrium
force balance relationship (17), $\langle\hat{\bf F}({\mathbf{r}})\rangle=0$.
This identity holds regardless of the form of the external potential $V_{\rm
ext}({\mathbf{r}})$, and we can hence functionally differentiate both sides of
the equation with respect to $V_{\rm ext}({\mathbf{r}})$. Clearly the right
hand side will remain zero. Differentiating the left hand side yields

$\displaystyle\frac{\delta\langle\hat{\bf F}({\mathbf{r}})\rangle}{\delta
V_{\rm ext}({\mathbf{r}}^{\prime})}$ $\displaystyle=-\beta{\rm cov}(\hat{\bf
F}({\mathbf{r}}),\hat{\rho}({\mathbf{r}}^{\prime}))+\Big{\langle}\frac{\delta\hat{\bf
F}({\mathbf{r}})}{\delta V_{\rm ext}({\mathbf{r}}^{\prime})}\Big{\rangle}.$
(27)

That Eq. (27) holds can be seen by writing out explicitly the thermal average
on the left hand side and functionally differentiating all dependencies on the
external potential as we will show. As an aside, Eckert et al.
eckert2023fluctuation have recently pointed out that for an arbitrary phase
space function $\hat{A}$, which is taken to be independent of $V_{\rm
ext}({\mathbf{r}})$, one has $\delta\langle\hat{A}\rangle/\delta V_{\rm
ext}({\mathbf{r}})=-\beta{\rm cov}(\hat{A},\hat{\rho}({\mathbf{r}}))$, which
upon choosing $\hat{A}=\hat{\bf F}({\mathbf{r}})$ generates the first term on
the right hand side of Eq. (27); the second term then stems from taking
account of the explicit dependence of $\hat{\bf F}({\mathbf{r}})$ on $V_{\rm
ext}({\mathbf{r}})$ according to Eq. (7) as we demonstrate in the following.
Considering the functional derivative on the left hand side of Eq. (27), from
the product rule one obtains the following structure:

$\displaystyle{\rm Tr}\,\frac{\delta}{\delta V_{\rm
ext}({\mathbf{r}}^{\prime})}\Xi^{-1}e^{-\beta(H-\mu N)}\hat{\bf
F}({\mathbf{r}})$ $\displaystyle\qquad={\rm Tr}\,\frac{\delta\Xi^{-1}}{\delta
V_{\rm ext}({\mathbf{r}}^{\prime})}e^{-\beta(H-\mu N)}\hat{\bf
F}({\mathbf{r}})$ $\displaystyle\qquad\quad+{\rm Tr}\,\Xi^{-1}e^{-\beta(H-\mu
N)}\Big{[}\frac{\delta(-\beta H)}{\delta V_{\rm
ext}({\mathbf{r}}^{\prime})}\hat{\bf F}({\mathbf{r}})+\frac{\delta\hat{\bf
F}({\mathbf{r}})}{\delta V_{\rm ext}({\mathbf{r}}^{\prime})}\Big{]}.$ (28)

The first term on the right hand side gives
$\beta\rho({\mathbf{r}}^{\prime}){\bf F}({\mathbf{r}})$ upon observing that
$-\Xi^{-2}\delta\Xi/\delta V_{\rm
ext}({\mathbf{r}}^{\prime})=\Xi^{-1}\beta\rho({\mathbf{r}}^{\prime})$, where
the density profile can be taken out of the trace in Eq. (27). The second term
follows straightforwardly as
$-\beta\langle\hat{\rho}({\mathbf{r}}^{\prime}){\bf
F}({\mathbf{r}})\rangle+\langle\delta\hat{\bf F}({\mathbf{r}})/\delta V_{\rm
ext}({\mathbf{r}}^{\prime})\rangle$ when taking account of the fact that
$\delta H/\delta V_{\rm
ext}({\mathbf{r}}^{\prime})=\hat{\rho}({\mathbf{r}}^{\prime})$. Regrouping the
terms then yields Eq. (27).

The covariance on the right hand side of Eq. (27) reduces to the correlation,
i.e., the mean of the product of the two operators, ${\rm cov}(\beta\hat{\bf
F}({\mathbf{r}}),\hat{\rho}({\mathbf{r}}^{\prime}))=\langle\beta\hat{\bf
F}({\mathbf{r}})\hat{\rho}({\mathbf{r}}^{\prime})\rangle$, again because the
average force density vanishes in equilibrium, $\langle\hat{\bf
F}({\mathbf{r}})\rangle={\bf F}({\mathbf{r}})=0$, cf. Eq. (17). The second
term in Eq. (27) can be re-written as follows:

$\displaystyle\Big{\langle}\frac{\delta\hat{\bf F}({\mathbf{r}})}{\delta
V_{\rm ext}({\mathbf{r}}^{\prime})}\Big{\rangle}$
$\displaystyle=-\Big{\langle}\frac{\delta}{\delta V_{\rm
ext}({\mathbf{r}}^{\prime})}\sum_{i}\delta({\mathbf{r}}-{\mathbf{r}}_{i})\nabla_{i}V_{\rm
ext}({\mathbf{r}}_{i})\Big{\rangle}$ (29)
$\displaystyle=-\Big{\langle}\sum_{i}\delta({\mathbf{r}}-{\mathbf{r}}_{i})\nabla_{i}\frac{\delta
V_{\rm ext}({\mathbf{r}}_{i})}{\delta V_{\rm
ext}({\mathbf{r}}^{\prime})}\Big{\rangle}$ (30)
$\displaystyle=\Big{\langle}\sum_{i}\delta({\mathbf{r}}-{\mathbf{r}}_{i})\nabla^{\prime}\delta({\mathbf{r}}_{i}-{\mathbf{r}}^{\prime})\Big{\rangle}$
(31)
$\displaystyle=\nabla^{\prime}\Big{\langle}\sum_{i}\delta({\mathbf{r}}-{\mathbf{r}}_{i})\delta({\mathbf{r}}_{i}-{\mathbf{r}}^{\prime})\Big{\rangle}$
(32) $\displaystyle\equiv\nabla^{\prime}\rho_{2}^{\rm
self}({\mathbf{r}},{\mathbf{r}}^{\prime}).$ (33)

In the last step above we have introduced the self two-body density, defined
as

$\displaystyle\rho_{2}^{\rm self}({\mathbf{r}},{\mathbf{r}}^{\prime})$
$\displaystyle=\Big{\langle}\sum_{i}\delta({\mathbf{r}}-{\mathbf{r}}_{i})\delta({\mathbf{r}}^{\prime}-{\mathbf{r}}_{i})\Big{\rangle}.$
(34)

In typical applications one would re-write the right hand side of Eq. (34) as
$\rho({\mathbf{r}})\delta({\mathbf{r}}-{\mathbf{r}}^{\prime})$, but due to the
present gradient structure some care is required and we thus keep
$\rho_{2}^{\rm self}({\mathbf{r}},{\mathbf{r}}^{\prime})$ in its full form
(34).

Collecting terms and recalling that the right hand side of Eq. (27) vanishes
identically we obtain the following compact force-density correlation sum
rule:

$\displaystyle\Big{\langle}\beta\hat{\bf
F}({\mathbf{r}})\hat{\rho}({\mathbf{r}}^{\prime})\Big{\rangle}$
$\displaystyle=\nabla^{\prime}\rho_{2}^{\rm
self}({\mathbf{r}},{\mathbf{r}}^{\prime}).$ (35)

Equation (35) is a remarkable and highly nontrival result, given the seemingly
complex and a priori possibly highly correlated nature of its left hand side.
However, in reality only the relatively simple local term on the right hand
side remains. Hence for two distinct positions ${\mathbf{r}}$ and
${\mathbf{r}}^{\prime}$ the singular right hand side vanishes and Eq. (35)
reduces trivially to

$\displaystyle\Big{\langle}\beta\hat{\bf
F}({\mathbf{r}})\hat{\rho}({\mathbf{r}}^{\prime})\Big{\rangle}$
$\displaystyle=0,\quad{\mathbf{r}}\neq{\mathbf{r}}^{\prime}.$ (36)

We can also conclude, as the distinct contribution vansishes, that at all
positions ${\mathbf{r}}$ and ${\mathbf{r}}^{\prime}$, even when
${\mathbf{r}}={\mathbf{r}}^{\prime}$, we have

$\displaystyle\Big{\langle}\beta\hat{\bf
F}({\mathbf{r}})\hat{\rho}({\mathbf{r}}^{\prime})\Big{\rangle}_{\rm dist}$
$\displaystyle=0.$ (37)

We recall that the present proof of the density-force sum rule (35) is based
on functionally differentiating ${\bf F}({\mathbf{r}})=0$ with respect to
$V_{\rm ext}({\mathbf{r}}^{\prime})$. An alternative derivation of Eq. (35)
can be based on the Noether invariance of the density profile. This thermal
symmetry is expressed as
$\delta\rho({\mathbf{r}}^{\prime})/\delta{\bm{\epsilon}}({\mathbf{r}})|_{{\bm{\epsilon}}=0}=0$.
As before we denote the transformed positions by $\tilde{\mathbf{r}}_{i}$ and
consider

$\displaystyle\frac{\delta\rho({\mathbf{r}}^{\prime})}{\delta{\bm{\epsilon}}({\mathbf{r}})}$
$\displaystyle=\frac{\delta}{\delta{\bm{\epsilon}}({\mathbf{r}})}\Big{\langle}\sum_{i}\delta({\mathbf{r}}^{\prime}-\tilde{\mathbf{r}}_{i})\Big{\rangle}$
(38) $\displaystyle=\Big{\langle}\beta\hat{\bf
F}({\mathbf{r}})\hat{\rho}({\mathbf{r}}^{\prime})\Big{\rangle}+\Big{\langle}\sum_{i}\frac{\delta}{\delta{\bm{\epsilon}}({\mathbf{r}})}\delta({\mathbf{r}}^{\prime}-{\mathbf{r}}_{i}-{\bm{\epsilon}}({\mathbf{r}}_{i}))\Big{\rangle},$
(39)

where the first term on the right hand side results from functionally
differentiating the probability distribution and it is already evaluated at
${\bm{\epsilon}}({\mathbf{r}})=0$. The second term can be re-written as

$\displaystyle\Big{\langle}\sum_{i}\nabla^{\prime}\delta({\mathbf{r}}^{\prime}-{\mathbf{r}}_{i}-{\bm{\epsilon}}({\mathbf{r}}_{i}))\cdot\frac{\delta[-{\bm{\epsilon}}({\mathbf{r}}_{i})]}{\delta{\bm{\epsilon}}({\mathbf{r}})}\Big{|}_{{\bm{\epsilon}}=0}\Big{\rangle}$
$\displaystyle\qquad\qquad=-\nabla^{\prime}\Big{\langle}\sum_{i}\delta({\mathbf{r}}^{\prime}-{\mathbf{r}}_{i})\delta({\mathbf{r}}_{i}-{\mathbf{r}})\Big{\rangle}$
$\displaystyle\qquad\qquad=-\nabla^{\prime}\rho_{2}^{\rm
self}({\mathbf{r}},{\mathbf{r}}^{\prime}),$ (40)

where we have used the definition (34) of the self two-body density
distribution $\rho_{2}^{\rm self}({\mathbf{r}},{\mathbf{r}}^{\prime})$. The
result (40), when inserted into Eq. (39) and recalling that the left hand side
thereof needs to vanish identically, gives the sum rule (35).

That both derivations give identical results ultimately lies in the fact that
the order of differentation is the only genuine difference, when starting from
the free energy as an overarching object and building a mixed second
derivative, where the order of differentation is interchanged:
$\delta^{2}\Omega[{\bm{\epsilon}}]/[\delta{\bm{\epsilon}}({\mathbf{r}})\delta
V_{\rm ext}({\mathbf{r}}^{\prime})]=\delta^{2}\Omega[{\bm{\epsilon}}]/[\delta
V_{\rm ext}({\mathbf{r}}^{\prime})\delta{\bm{\epsilon}}({\mathbf{r}})]$ where
we recall that $\delta\Omega/\delta V_{\rm
ext}({\mathbf{r}})=\rho({\mathbf{r}})$ and
$\delta\Omega[{\bm{\epsilon}}]/\delta{\bm{\epsilon}}({\mathbf{r}})|_{{\bm{\epsilon}}=0}=-{\bf
F}({\mathbf{r}})$ according to Eq. (16).

Lastly, Eq. (35) can be viewed as a special case of the recent more general
hyperforce correlation theory by Robitschko et al. robitschko2023any . Their
theory applies to general phase space functions $\hat{A}({\mathbf{r}}^{N},{\bf
p}^{N})$ and for configuration-dependent cases, $\hat{A}({\mathbf{r}}^{N})$,
it ascertains the validity of the identity $\langle\beta\hat{\bf
F}({\mathbf{r}})\hat{A}({\mathbf{r}}^{N})\rangle=-\langle\sum_{i}\delta({\mathbf{r}}-{\mathbf{r}}_{i})\nabla_{i}\hat{A}({\mathbf{r}}^{N})\rangle$.
Choosing the density operator as the observable of interest,
$\hat{A}({\mathbf{r}}^{N})=\hat{\rho}({\mathbf{r}}^{\prime})=\sum_{j}\delta({\mathbf{r}}^{\prime}-{\mathbf{r}}_{j})$,
then replacing
$\nabla_{i}\delta({\mathbf{r}}^{\prime}-{\mathbf{r}}_{i})=-\nabla^{\prime}\delta({\mathbf{r}}^{\prime}-{\mathbf{r}}_{i})$,
and identifying $\rho_{2}^{\rm self}({\mathbf{r}},{\mathbf{r}}^{\prime})$
according to Eq. (34) also yields Eq. (35).

In summary, we have obtained the density-force sum rule (35), which is ready
to use in the further investigation of the local second-order Noether
invariance structure. For subsequent use we find it convenient to re-write the
density-force correlation sum rule (35) using the splitting (7) of the total
force density operator $\hat{\bf
F}({\mathbf{r}})=\nabla\cdot\hat{\bm{\tau}}({\mathbf{r}})+\hat{\bf
F}_{U}({\mathbf{r}})$, where we recall the definition (8) of the kinetic
stress density operator $\hat{\bm{\tau}}({\mathbf{r}})$ and the definition
(23) of the potential force density operator $\hat{\bf F}_{U}({\mathbf{r}})$.
After re-arranging Eq. (35) we find

$\displaystyle\Big{\langle}\beta\hat{\bf
F}_{U}({\mathbf{r}})\hat{\rho}({\mathbf{r}}^{\prime})\Big{\rangle}$
$\displaystyle=\nabla^{\prime}\rho_{2}^{\rm
self}({\mathbf{r}},{\mathbf{r}}^{\prime})-\nabla\cdot\Big{\langle}\beta\hat{\bm{\tau}}({\mathbf{r}})\hat{\rho}({\mathbf{r}}^{\prime})\Big{\rangle}$
(41) $\displaystyle=(\nabla^{\prime}+\nabla)\rho_{2}^{\rm
self}({\mathbf{r}},{\mathbf{r}}^{\prime})+\nabla\rho_{2}({\mathbf{r}},{\mathbf{r}}^{\prime}),$
(42)

where the form (42) arises from performing the momentum integrals over the
standard Maxwell probability distribution in the second term of Eq. (41) and
identifying
$\langle\hat{\rho}({\mathbf{r}})\hat{\rho}({\mathbf{r}}^{\prime})\rangle=\rho_{2}({\mathbf{r}},{\mathbf{r}}^{\prime})+\rho_{2}^{\rm
self}({\mathbf{r}},{\mathbf{r}}^{\prime})$. Building furthermore the gradient
with respect to the position ${\mathbf{r}}^{\prime}$ and transposing [note
that the transpose $(\nabla^{\prime}\nabla)^{\sf T}=\nabla\nabla^{\prime}$]
yields straightforwardly the following tensorial identities:

$\displaystyle\Big{\langle}\beta\hat{\bf
F}_{U}({\mathbf{r}})\nabla^{\prime}\hat{\rho}({\mathbf{r}}^{\prime})\Big{\rangle}$
$\displaystyle=(\nabla^{\prime}\nabla^{\prime}+\nabla\nabla^{\prime})\rho_{2}^{\rm
self}({\mathbf{r}},{\mathbf{r}}^{\prime})$
$\displaystyle\quad+\nabla\nabla^{\prime}\rho_{2}({\mathbf{r}},{\mathbf{r}}^{\prime}),$
(43) $\displaystyle\Big{\langle}\nabla\hat{\rho}({\mathbf{r}})\beta\hat{\bf
F}_{U}({\mathbf{r}}^{\prime})\Big{\rangle}$
$\displaystyle=(\nabla\nabla+\nabla\nabla^{\prime})\rho_{2}^{\rm
self}({\mathbf{r}},{\mathbf{r}}^{\prime})$
$\displaystyle\quad+\nabla\nabla^{\prime}\rho_{2}({\mathbf{r}},{\mathbf{r}}^{\prime}).$
(44)

The version Eq. (44) is thereby obtained from re-transposing (43) and
interchanging the positions arguments ${\mathbf{r}}$ and
${\mathbf{r}}^{\prime}$ upon observing the exchange symmetry $\rho_{2}^{\rm
self}({\mathbf{r}},{\mathbf{r}}^{\prime})=\rho_{2}^{\rm
self}({\mathbf{r}}^{\prime},{\mathbf{r}})$; we recall the definition (34) of
the self two-body density and clearly also
$\rho_{2}({\mathbf{r}},{\mathbf{r}}^{\prime})=\rho_{2}({\mathbf{r}}^{\prime},{\mathbf{r}})$
as is apparent from its definition (26).

The identities (43) and (44) are now in a form ready to be used subsequently
in Sec. V in the proof of the Noether two-body fore correlation sum rules (24)
and (25).

##  V Force correlator splitting

We aim to re-write the full force autocorrelator, as it e.g. appears on the
left hand side of Eq. (22), via the autocorrelator of the potential force
density operator $\hat{\bf F}_{U}({\mathbf{r}})$, see Eq. (23), in order to
simplify the structure of the trivial kinetic contributions. As above, we use
the definition (8) of the kinetic stress operator to express $\hat{\bf
F}({\mathbf{r}})=\nabla\cdot\hat{\bm{\tau}}({\mathbf{r}})+\hat{\bf
F}_{U}({\mathbf{r}})$. We start with the force autocorrelator and obtain by
multiplying out and re-ordering:

$\displaystyle\Big{\langle}\hat{\bf F}({\mathbf{r}})\hat{\bf
F}({\mathbf{r}}^{\prime})\Big{\rangle}$
$\displaystyle=\Big{\langle}\nabla\cdot\hat{\bm{\tau}}({\mathbf{r}})\nabla^{\prime}\cdot\hat{\bm{\tau}}({\mathbf{r}}^{\prime})\Big{\rangle}+\Big{\langle}\hat{\bf
F}_{U}({\mathbf{r}})\hat{\bf F}_{U}^{\prime}({\mathbf{r}})\Big{\rangle}$
$\displaystyle\;+\Big{\langle}\hat{\bf
F}_{U}({\mathbf{r}})\nabla^{\prime}\cdot\hat{\bm{\tau}}({\mathbf{r}}^{\prime})\Big{\rangle}+\Big{\langle}\nabla\cdot\hat{\bm{\tau}}({\mathbf{r}})\hat{\bf
F}_{U}({\mathbf{r}}^{\prime})\Big{\rangle}.$ (45)

We start with the purely kinetic two-body contribution, i.e. the first term in
Eq. (45) and obtain the following result:

$\displaystyle\Big{\langle}\nabla\cdot\beta\hat{\bm{\tau}}({\mathbf{r}})\nabla^{\prime}\cdot\beta\hat{\bm{\tau}}({\mathbf{r}}^{\prime})\Big{\rangle}$
$\displaystyle=(\nabla\nabla^{\prime}+\nabla^{\prime}\nabla+{\mathbbm{1}}\nabla\cdot\nabla^{\prime})\rho_{2}^{\rm
self}({\mathbf{r}},{\mathbf{r}}^{\prime})+\nabla\nabla^{\prime}\rho_{2}({\mathbf{r}},{\mathbf{r}}^{\prime})$
(46)

where we have again used the definition (34) of the self two-body density and
we recall that ${\mathbbm{1}}$ indicates the $3\times 3$-unit matrix. In
Appendix A we lay out the explicit calculation of both the self and distinct
contributions to Eq. (46).

The structure of the two mixed terms in Eq. (45) allows to carry out the
momentum integrals in a simple way. This enables us to simplify the two
kinetic contributions according to

$\displaystyle\Big{\langle}\hat{\bf
F}_{U}({\mathbf{r}})\nabla^{\prime}\cdot\beta\hat{\bm{\tau}}({\mathbf{r}}^{\prime})\Big{\rangle}$
$\displaystyle=-\Big{\langle}\hat{\bf
F}_{U}({\mathbf{r}})\nabla^{\prime}\hat{\rho}({\mathbf{r}}^{\prime})\Big{\rangle},$
(47)
$\displaystyle\Big{\langle}\nabla\cdot\beta\hat{\bm{\tau}}({\mathbf{r}})\hat{\bf
F}_{U}({\mathbf{r}}^{\prime})\Big{\rangle}$
$\displaystyle=-\Big{\langle}\nabla\hat{\rho}({\mathbf{r}})\hat{\bf
F}_{U}({\mathbf{r}}^{\prime})\Big{\rangle}.$ (48)

Collecting all terms and using Eqs. (43) and (44) allows then to re-write Eq.
(45) upon multiplying by $\beta^{2}$ as

$\displaystyle\Big{\langle}\beta\hat{\bf F}({\mathbf{r}})\beta\hat{\bf
F}({\mathbf{r}}^{\prime})\Big{\rangle}=\Big{\langle}\beta\hat{\bf
F}_{U}({\mathbf{r}})\beta\hat{\bf
F}_{U}({\mathbf{r}}^{\prime})\Big{\rangle}-\nabla\nabla^{\prime}\rho_{2}({\mathbf{r}},{\mathbf{r}}^{\prime})$
$\displaystyle\quad+({\mathbbm{1}}\nabla\cdot\nabla^{\prime}+\nabla^{\prime}\nabla-\nabla\nabla^{\prime}-\nabla^{\prime}\nabla^{\prime}-\nabla\nabla)\rho_{2}^{\rm
self}({\mathbf{r}},{\mathbf{r}}^{\prime}).$ (49)

Equation (49) constitutes an explicit form of the left hand side of the
locally resolved two-body sum rule Eq. (22). We recall that the sum rule (22)
balances the force-force correlations on its left hand side with an energy
curvature term on the right hand side.

In order to address this right hand side, we split the Hamiltonian, expressed
in the new coordinates and hence functionally depending on the shifting field
${\bm{\epsilon}}({\mathbf{r}})$, into kinetic and potential energy
contributions, $H[{\bm{\epsilon}}]=H_{\rm
kin}[{\bm{\epsilon}}]+H_{U}[{\bm{\epsilon}}]$. The energy curvature also
consists of kinetic and potential energy contributions. We defer the
calculations for the kinetic energy curvature to Appendix B and for the
potential energy curvature to Appendix C. The results are as follows:

$\displaystyle\Big{\langle}\frac{\delta^{2}H_{\rm
kin}[{\bm{\epsilon}}]}{\delta{\bm{\epsilon}}({\mathbf{r}})\delta{\bm{\epsilon}}({\mathbf{r}}^{\prime})}\Big{|}_{{\bm{\epsilon}}=0}\Big{\rangle}$
$\displaystyle=k_{B}T({\mathbbm{1}}\nabla\cdot\nabla^{\prime}+2\nabla^{\prime}\nabla)\rho_{2}^{\rm
self}({\mathbf{r}},{\mathbf{r}}^{\prime}),$ (50)
$\displaystyle\Big{\langle}\frac{\delta^{2}H_{U}[{\bm{\epsilon}}]}{\delta{\bm{\epsilon}}({\mathbf{r}})\delta{\bm{\epsilon}}({\mathbf{r}}^{\prime})}\Big{|}_{{\bm{\epsilon}}=0}\Big{\rangle}$
$\displaystyle=$
$\displaystyle\Big{\langle}\sum_{ij}\delta({\mathbf{r}}-{\mathbf{r}}_{i})\delta({\mathbf{r}}^{\prime}-{\mathbf{r}}_{j})\nabla_{i}\nabla_{j}u({\mathbf{r}}^{N})\Big{\rangle}$
$\displaystyle+\rho_{2}^{\rm
self}({\mathbf{r}},{\mathbf{r}}^{\prime})\nabla\nabla V_{\rm
ext}({\mathbf{r}}),$ (51)

and we recall that their sum constitutes the right hand side of the
generically expressed Nother sum rule (22).

We can hence express Eq. (22) by rewriting its left hand side via the right
hand side of Eq. (49) and its right hand side via the sum of the right hand
sides of the curvature results (50) and (51). Collecting all terms yields the
following result:

$\displaystyle\Big{\langle}\beta\hat{\bf F}_{U}({\mathbf{r}})\beta\hat{\bf
F}_{U}({\mathbf{r}}^{\prime})\Big{\rangle}-\nabla\nabla^{\prime}\rho_{2}({\mathbf{r}},{\mathbf{r}}^{\prime})$
$\displaystyle\quad+({\mathbbm{1}}\nabla\cdot\nabla^{\prime}+\nabla^{\prime}\nabla-\nabla\nabla^{\prime}-\nabla^{\prime}\nabla^{\prime}-\nabla\nabla)\rho_{2}^{\rm
self}({\mathbf{r}},{\mathbf{r}}^{\prime})=$
$\displaystyle({\mathbbm{1}}\nabla\cdot\nabla^{\prime}+2\nabla^{\prime}\nabla)\rho_{2}^{\rm
self}({\mathbf{r}},{\mathbf{r}}^{\prime})$
$\displaystyle\quad+\Big{\langle}\sum_{ij}\delta({\mathbf{r}}-{\mathbf{r}}_{i})\delta({\mathbf{r}}^{\prime}-{\mathbf{r}}_{j})\nabla_{i}\nabla_{j}\beta
u({\mathbf{r}}^{N})\Big{\rangle}$ $\displaystyle\quad+\rho_{2}^{\rm
self}({\mathbf{r}},{\mathbf{r}}^{\prime})\nabla\nabla\beta V_{\rm
ext}({\mathbf{r}}).$ (52)

Re-arranging and simplifying yields:

$\displaystyle\Big{\langle}\beta\hat{\bf F}_{U}({\mathbf{r}})\beta\hat{\bf
F}_{U}({\mathbf{r}}^{\prime})\Big{\rangle}=$
$\displaystyle\quad\nabla\nabla^{\prime}\rho_{2}({\mathbf{r}},{\mathbf{r}}^{\prime})+\delta({\mathbf{r}}-{\mathbf{r}}^{\prime})[\nabla\nabla\rho({\mathbf{r}})+\rho({\mathbf{r}})\nabla\nabla\beta
V_{\rm ext}({\mathbf{r}})]$
$\displaystyle\quad+\Big{\langle}\sum_{ij}\delta({\mathbf{r}}-{\mathbf{r}}_{i})\delta({\mathbf{r}}^{\prime}-{\mathbf{r}}_{j})\nabla_{i}\nabla_{j}\beta
u({\mathbf{r}}^{N})\Big{\rangle},$ (53)

where we could finally replace the self two-body density as follows. In the
external curvature term we have $\rho_{2}^{\rm
self}({\mathbf{r}},{\mathbf{r}}^{\prime})=\delta({\mathbf{r}}-{\mathbf{r}}^{\prime})\rho({\mathbf{r}})$
and in the ideal contribution we could simplify
$(\nabla+\nabla^{\prime})(\nabla+\nabla^{\prime})\rho_{2}^{\rm
self}({\mathbf{r}},{\mathbf{r}}^{\prime})=\delta({\mathbf{r}}-{\mathbf{r}}^{\prime})\nabla\nabla\rho({\mathbf{r}})$,
which can be verified explicitly by calculating derivatives in
$(\nabla+\nabla^{\prime})(\nabla+\nabla^{\prime})[\delta({\mathbf{r}}-{\mathbf{r}}^{\prime})\rho({\mathbf{r}})]$
using the product rule and exploiting
$\nabla^{\prime}\delta({\mathbf{r}}-{\mathbf{r}}^{\prime})=-\nabla\delta({\mathbf{r}}-{\mathbf{r}}^{\prime})$.
Splitting Eq. (53) into distinct and self contributions, where the latter are
identified as having a common factor
$\delta({\mathbf{r}}-{\mathbf{r}}^{\prime})$, leads respectively to the
distinct sum rule (24) and upon leaving away the common factor
$\delta({\mathbf{r}}-{\mathbf{r}}^{\prime})$ to the self sum rule (25).

This completes our proof of the locally resolved two-body Noether invariance
sum rules (24) and (25). As mentioned above, all derivations continue to hold
in the canonical ensemble. The mechanism for this universality is the
primarily mechanical nature, we recall the shifting transformation (3) and
(4), of the considered thermal invariance, which is oblivious to the presence
of a particle bath.

Our theory is now ready to be applied to real systems and we present explicit
results for a variety of common soft matter models in Sec. VI. These results
demonstrate the validity of the sum rules and they show the prowess of the
force-force and force-gradient correlation functions to systematically
quantify and shed light on the self-structuring mechanisms at play in
equilibrium. Readers who are keen to first follow further formal argumentation
that demonstrates the validity of the Nother sum rules are welcome to go to
Appendix D, where we lay out our partial phase space integration route.

Figure 2: Correlation functions analogous to Fig. 1, with results for the
Lennard-Jones liquid shown again as a reference (first column), the Weeks-
Chandler-Andersen fluid (second column), monatomic water at (approximately)
ambient conditions molinero2009 ; coe2022water (third column), and the three-
body interacting gel saw2009 ; saw2011 ; sammueller2023gel (fourth column).
Shown are results for the standard pair correlation function $g(r)$ for each
system (top row). Further two-body structure is revealed by the radial
($\parallel$) and transversal ($\perp$) components of the force-gradient
correlator ${\sf g}_{\nabla\\!f}(r)$ (middle row) and of the force-force
correlator ${\sf g}_{f\\!f}(r)$ (bottom row). The dotted lines in the middle
row indicate the two components of the force-gradient correlation function as
obtained from the simplifications (64) and (65) for the two pairwise (Lennard-
Jones and Weeks-Chandler-Andersen) Hamiltonians. The dotted lines in the
bottom row indicate the results according to the Noether force sum rules (62)
and (63). Taken from Ref. sammueller2023whatIsLiquid .  Figure 3: Correlation
functions analogous to those shown in Figs. 1 and 2, but for the Yukawa liquid
(first column), the soft-sphere dipolar fluid (second column), the Stockmayer
fluid (third column), and the Gay-Berne model in the isotropic (fourth column)
and nematic phase (fifth column). The results for the anisotropic models are
obtained from canonical Monte Carlo simulations, and they are averaged over
the microscopic orientations; the simulation box volume is $V=(20\sigma)^{3}$
and the long-ranged interactions are cut off at radial distance of $10\sigma$.
Shown are the pair correlation function $g(r)$ (top row), the force-gradient
correlator ${\sf g}_{\nabla\\!f}(r)$ (middle row), and the force-force
correlator ${\sf g}_{f\\!f}(r)$ together with the result from the Noether
identities (62) and (63) (bottom row). Adapted from Ref.
sammueller2023whatIsLiquid (Supplementary Material).

##  VI Noether structure in bulk

We specialize to homogeneous bulk liquid states, where
$\rho({\mathbf{r}})=\rho_{b}=\rm const$ and $V_{\rm ext}({\mathbf{r}})=0$. The
potential forces then only arise from interparticle contributions and hence
$\hat{\bf F}_{U}({\mathbf{r}})=\hat{\bf F}_{\rm int}({\mathbf{r}})$ with the
interparticle force density operator $\hat{\bf F}_{\rm int}({\mathbf{r}})$
being given by Eq. (9). We address the distinct sum rule (24) and use the
standard pair correlation function or “radial distribution function” $g(r)$.
Here $r=|{\mathbf{r}}-{\mathbf{r}}^{\prime}|$ denotes the separation distance
between the two particles. The relationship to the two-body density is simply
via normalization with the squared bulk density $\rho_{b}$:

$\displaystyle
g(|{\mathbf{r}}-{\mathbf{r}}^{\prime}|)=\rho_{2}({\mathbf{r}},{\mathbf{r}}^{\prime})/\rho_{b}^{2}.$
(54)

Furthermore we follow Ref. sammueller2023whatIsLiquid in introducing both the
force-force pair correlation function ${\sf g}_{f\\!f}(r)$ and the force
gradient correlator ${\sf g}_{\nabla\\!f}(r)$. These are rank-two tensorial
correlation functions, as represented by $3\times 3$-matrices in case of a
three-dimensional system, and they are defined via

$\displaystyle{\sf g}_{f\\!f}(|{\mathbf{r}}-{\mathbf{r}}^{\prime}|)$
$\displaystyle=\frac{\beta^{2}}{\rho_{b}^{2}}\langle\hat{\bf F}_{\rm
int}({\mathbf{r}})\hat{\bf F}_{\rm int}({\mathbf{r}}^{\prime})\rangle_{\rm
dist},$ (55) $\displaystyle{\sf
g}_{\nabla\\!f}(|{\mathbf{r}}-{\mathbf{r}}^{\prime}|)$ $\displaystyle=$
$\displaystyle-\frac{\beta}{\rho_{b}^{2}}\Big{\langle}\sum_{ij(\neq)}$
$\displaystyle\delta({\mathbf{r}}-{\mathbf{r}}_{i})\delta({\mathbf{r}}^{\prime}-{\mathbf{r}}_{j})\nabla_{i}\nabla_{j}u({\mathbf{r}}^{N})\Big{\rangle}.$
(56)

The force gradient correlator ${\sf g}_{\nabla\\!f}(r)$ is also the negative
mean potential curvature, as is apparent from the right hand side of Eq. (56)
via $\nabla_{i}\nabla_{j}u({\mathbf{r}}^{N})$. One can hence refer to the
correlator also as the mean negative Hessian, with respect to differentiation
by particle positions, of the interparticle potential energy.

On the other hand an equally valid interpretation is that of the negative
interparticle force gradient. Defining the interparticle force on particle $i$
as ${\bf f}_{i}({\mathbf{r}}^{N})=-\nabla_{i}u({\mathbf{r}}^{N})$ and that on
particle $j$ as ${\bf f}_{j}({\mathbf{r}}^{N})=-\nabla_{j}u({\mathbf{r}}^{N})$
we can re-write the potential energy curvature as
$\nabla_{i}\nabla_{j}u({\mathbf{r}}^{N})=-\nabla_{i}{\bf
f}_{j}({\mathbf{r}}^{N})=-[\nabla_{j}{\bf f}_{i}({\mathbf{r}}^{N})]^{\sf T}$,
where the possibility of interchange of the particle indices in the two latter
expressions is due to interchange of the order of the two derivatives.

Expressing the forces in this way allows to re-write Eqs. (55) and (56) in the
following more explicit forms:

$\displaystyle{\sf g}_{f\\!f}(|{\mathbf{r}}-{\mathbf{r}}^{\prime}|)=$
$\displaystyle\qquad\Big{\langle}\sum_{ij(\neq)}\delta({\mathbf{r}}-{\mathbf{r}}_{i})\delta({\mathbf{r}}^{\prime}-{\mathbf{r}}_{j})\beta{\bf
f}_{i}({\mathbf{r}}^{N})\beta{\bf
f}_{j}({\mathbf{r}}^{N})\Big{\rangle}\Big{/}\rho_{b}^{2},$ (57)
$\displaystyle{\sf g}_{\nabla\\!f}(|{\mathbf{r}}-{\mathbf{r}}^{\prime}|)=$
$\displaystyle\qquad\Big{\langle}\sum_{ij(\neq)}\delta({\mathbf{r}}-{\mathbf{r}}_{i})\delta({\mathbf{r}}^{\prime}-{\mathbf{r}}_{j})\nabla_{i}\beta{\bf
f}_{j}({\mathbf{r}}^{N})\Big{\rangle}\Big{/}\rho_{b}^{2},$ (58)

where in the last term one can alternatively replace $\nabla_{i}{\bf
f}_{j}({\mathbf{r}}^{N})=[\nabla_{j}{\bf f}_{i}({\mathbf{r}}^{N})]^{\sf T}$,
as laid out above.

In the fully position-resolved case, which is appropriate for the
investigation of spatially inhomogeneous situations, where
$\rho({\mathbf{r}})\neq\rm const$, suitable correlators are the following two-
body force-force and force-gradient density distributions:

$\displaystyle{\sf
G}_{f\\!f}({\mathbf{r}},{\mathbf{r}}^{\prime})=\Big{\langle}\sum_{ij(\neq)}\delta({\mathbf{r}}-{\mathbf{r}}_{i})\delta({\mathbf{r}}^{\prime}-{\mathbf{r}}_{j})\beta{\bf
f}_{i}({\mathbf{r}}^{N})\beta{\bf f}_{j}({\mathbf{r}}^{N})\Big{\rangle},$ (59)
$\displaystyle{\sf
G}_{\nabla\\!f}({\mathbf{r}},{\mathbf{r}}^{\prime})=\Big{\langle}\sum_{ij(\neq)}\delta({\mathbf{r}}-{\mathbf{r}}_{i})\delta({\mathbf{r}}^{\prime}-{\mathbf{r}}_{j})\nabla_{i}\beta{\bf
f}_{j}({\mathbf{r}}^{N})\Big{\rangle}.$ (60)

Inhomogeneous versions of Eqs. (57) and (58) can be obtained via normalization
with the (inhomogeneous) density profile: ${\sf
G}_{f\\!f}({\mathbf{r}},{\mathbf{r}}^{\prime})/[\rho({\mathbf{r}})\rho({\mathbf{r}}^{\prime})]$
and ${\sf
G}_{\nabla\\!f}({\mathbf{r}},{\mathbf{r}}^{\prime})/[\rho({\mathbf{r}})\rho({\mathbf{r}}^{\prime})]$.

As announced in the introduction, using the two-body correlators defined via
Eqs. (54)-(56) the distinct sum (24) can be written for the case of a bulk
fluid in the following appealing compact form sammueller2023whatIsLiquid :

$\displaystyle\nabla\nabla g(r)$ $\displaystyle+{\sf g}_{\nabla\\!f}(r)+{\sf
g}_{f\\!f}(r)=0,$ (61)

which is Eq. (1). Due to the rotational symmetry of a bulk fluid, the only
nontrivial tensor components are parallel ($\parallel$) and transversal
($\perp$) to ${\mathbf{r}}-{\mathbf{r}}^{\prime}$, such that Eq. (61) reduces
to the following two nontrivial tensor components:

$\displaystyle
g^{\prime\prime}(r)+g_{{\nabla\\!f}\parallel}(r)+g_{{f\\!f}\parallel}(r)$
$\displaystyle=0,$ (62)
$\displaystyle\frac{g^{\prime}(r)}{r}+g_{{\nabla\\!f}\perp}(r)+g_{{f\\!f}\perp}(r)$
$\displaystyle=0,$ (63)

with the prime denoting the derivative with respect to $r$. In the chosen
coordinate system the matrices are diagonal, ${\rm
diag}(\parallel,\perp,\perp)$, with the first axis being parallel to
${\mathbf{r}}-{\mathbf{r}}^{\prime}$ and the second and third components being
equal and corresponding to two mutually perpendicular directions which are
also perpendicular to ${\mathbf{r}}-{\mathbf{r}}^{\prime}$.

More explicitly, in the considered radial symmetry the second derivative
operator in Eq. (61) reduces to $\nabla\nabla={\bf e}_{\parallel}{\bf
e}_{\parallel}\partial^{2}/\partial r^{2}+({\mathbbm{1}}-{\bf
e}_{\parallel}{\bf e}_{\parallel})r^{-1}\partial/\partial r$, where the
direction that connects the two points is ${\bf
e}_{\parallel}=({\mathbf{r}}-{\mathbf{r}}^{\prime})/|{\mathbf{r}}-{\mathbf{r}}^{\prime}|$.
This unit vector is complemented by two orthogonal directions ${\bf
e}_{\perp},{\bf e}_{\perp}^{\prime}$ such that all three ${\bf
e}_{\parallel},{\bf e}_{\perp},{\bf e}_{\perp}^{\prime}$ are mutually
orthogonal to each other. The perpendicular components then generate the
$g^{\prime}(r)/r$ term in Eq. (63). These tensor operations can be efficiently
performed in numerical work and we give details about the sampling strategies
used in our simulations in Appendix E.

For simple fluids, with the particles interacting mutually only via a pair
potential $\phi(r)$, the force gradient correlator (56) reduces to ${\sf
g}_{\nabla\\!f}(r)=\beta g(r)\nabla\nabla\phi(r)$ such that the two non-
trivial components ($\parallel$ and $\perp$) can be written according to

$\displaystyle g_{{\nabla\\!f}\parallel}(r)$ $\displaystyle=\beta
g(r)\phi^{\prime\prime}(r),$ (64) $\displaystyle g_{{\nabla\\!f}\perp}(r)$
$\displaystyle=\beta g(r)\frac{\phi^{\prime}(r)}{r}.$ (65)

This simplification is due to the reduction of the mixed derivative
$\nabla_{i}\nabla_{j}u({\mathbf{r}}^{N})=\nabla_{i}\nabla_{j}\sum_{kl(\neq)}\phi(|{\mathbf{r}}_{k}-{\mathbf{r}}_{l}|)/2=\nabla_{i}\nabla_{j}\phi(|{\mathbf{r}}_{i}-{\mathbf{r}}_{j}|)$,
for $i\neq j$, such that strikingly only a single pair potential bond, that
between the particle $i$ and $j$, contributes and the sums over particles $k$
and $l$ have disappeared.

As a consequence of the structure of the curvature correlator expressed in
Eqs. (64) and (65), we can rewrite the sum rules (62) and (63) for the case of
pair-wise interparticle forces in the following form:

$\displaystyle
g^{\prime\prime}(r)+\beta\phi^{\prime\prime}(r)g(r)+g_{{f\\!f}\parallel}(r)=0,$
(66)
$\displaystyle\frac{g^{\prime}(r)}{r}+\frac{\beta\phi^{\prime}(r)}{r}g(r)+g_{{f\\!f}\perp}(r)=0.$
(67)

The validity of the identities (66) and (67) can be analytically verified in
the low-density limit of a simple fluid, where $g(r)=\exp(-\beta\phi(r))$, as
laid out in Ref. sammueller2023whatIsLiquid . We here reproduce the argument.
In the low-density limit the force-force correlations are due to the
antiparallel direct forces between a particle pair:
$g_{{f\\!f}\parallel}(r)=-g(r)[\beta\phi^{\prime}(r)]^{2}$. Furthermore
$g_{{f\\!f}\perp}(r)=0$ due to the absence of a third particle at $\rho_{b}\to
0$ that could mediate a transversal force. Insertion into Eqs. (66) and (67)
verifies this solution.

As an illustration of the validity of these sum rules and for demonstrating
that fresh insight into the force correlation structure can be obtained
sammueller2023whatIsLiquid , we return to the Lennard-Jones system mentioned
in the Introduction. This fundamental model for simple systems comprises gas,
liquid and crystal phases. While our considerations were specific to the
symmetry of fluid phases, the structure of Eqs. (66) and (67) remains intact
upon averaging the distance vector ${\mathbf{r}}-{\mathbf{r}}^{\prime}$ over
all orientations and translations. We recall the display of simulation results
in Fig. 1. The parallel component of the force gradient correlator
$g_{{\nabla\\!f}\parallel}(r)$ has a strong positive peak at a distance near
the Lennard-Jones diameter. The perpendicular component
$g_{{\nabla\\!f}\perp}(r)$ is negative, which indicates anti-correlation as is
mediated by the surrounding particles. The structure of the crystal is
particularly rich, as can be expected. In all phases the sum rules are
satisfied to high numerical precision, which is remarkable as the three
correlation functions $g(r)$, ${\sf g}_{f\\!f}(r)$, and ${\sf
g}_{\nabla\\!f}(r)$ seem to be a priori independent of each other.

We show results for the Weeks-Chandler-Andersen liquid, for monatomic water,
using the parameterization by Molinero and Moore molinero2009 of the
Stillinger-Weber potential stillinger1985 , and for the three-body interacting
gel in Fig. 2. The results for the Lennard-Jones liquid are shown as a
reference. Going from the Lennard-Jones liquid (first column) to the Weeks-
Chandler-Andersen liquid (second column) at identical thermodynamic conditions
is a setup that lets one assess the influence of the interparticle attractive
forces, which are absent in the latter model. As is well-known, the pair
correlation function $g(r)$ is only very mildly affected by interparticle
attraction. In striking contrast the clear positive peak of the Lennard-Jones
liquid in the transversal mean force gradient $g_{{\nabla\\!f}\perp}(r)$ is
absent in the Weeks-Chandler-Andersen liquid. Going to the monatomic Molinero-
Moore molinero2009 ; coe2022water water model (third column) reveals richer
force gradient and force-force correlation structure over a broader range of
distances $r$ as compared to the previous pairwise models. We recall that the
three-body gel saw2009 ; saw2011 ; sammueller2023gel is, similar to the
monatomic water model, a mere reparametrization of the original Stillinger-
Weber stillinger1985 model. Yet both the force-gradient and the force-force
correlation structure (fourth column) changes dramatically when going from the
liquid to the gel. The perpendicular mean force-gradient
$g_{{\nabla\\!f}\perp}(r)$ develops a strong positive signal immediately
beyond the first peak of $g(r)$. Even more strikingly the transverse force-
force correlation function $g_{{f\\!f}\perp}(r)$ has a prominent negative peak
as opposed to the positive peak observed in all considered liquids. Such
negative force-force correlations are indicative of mechanical stability of
particle strands that the gel former develops and which are very apparent in
simulation snapshots (see e.g. Ref. sammueller2023gel ). These features cannot
be discerned from $g(r)$ alone, which despite its quantitatively exaggerated
appearance remains liquid-like, cf. the top right panel in Fig. 2.

In Fig. 3 we show results for the Yukawa liquid, the soft-sphere dipolar
fluid, the Stockmayer fluid allen2019topRev , and the Gay-Berne model gay1981
; brown1998 ; allen2019topRev in the isotropic and nematic phase. Here we use
Monte Carlo simulations to generate the numerical results, as this method is
more straightforward to use for the present models with orientational degrees
of freedom. In all cases considered the sum rules were found to be satisfied
to high numerical accuracy.

The results for the point-Yukawa fluid hamaguchi1997 ; dijkstra1998yukawa
(first column) serve again to illustrate generic liquid-like structuring. No
positive signal in $g_{{\nabla\\!f}\perp}(r)$ occurs as is consistent with the
purely repulsive character of the model and the correlation structure is
spread out over a larger range of distance as opposed to the Lennard-Jones
liquid (recall the first column of Fig. 2), which is consistent with the
longer-ranged decay of the Yukawa pair potential. Both the soft-sphere dipolar
system (second column) and the Stockmayer model (third column) display
prominently the signs of chain formation that we had previously identified for
the colloidal gel former: The perpendicular mean force gradient
$g_{{\nabla\\!f}\perp}(r)$ has a strong positive peak directly beyond the
first peak of $g(r)$ and the prominent peak of $g_{{f\\!f}\perp}(r)$ is
negative rather than positive as we find it to be in the simple fluid phases
and in monatomic (liquid) water.

Addressing the liquid-crystal forming Gay-Berne model reveals further
differences. Even the isotropic phase (fourth column of Fig. 3) displays
pronounced differences to all previous models, in that
$g_{{\nabla\\!f}\perp}(r)$ has lost its negative peak at distances smaller
than the location of the first peak of $g(r)$. Instead
$g_{{\nabla\\!f}\perp}(r)$ displays a prominent entirely positive peak, while
the behaviour of both force-force correlation components is qualitatively
similar to the chain forming models. Curiously, upon increasing the density,
such that the Gay-Berne model spontaneously develops nematic order, pronounced
quantitative increases occur in all observed features of both the mean force-
gradient and the force-force correlator. We re-iterate that we have
deliberately averaged over all orientations for simplicity. We expect that
resolving the dependence on orientation with respect to the nematic director
can reveal much further insight into the spatial structure of liquid crystal
ordering.

##  VII Conclusions

In conclusion we have provided a detailed description of the recent Noether
invariance theory for the two-body structure of liquids and of more general
soft matter systems sammueller2023whatIsLiquid . The theory is based on a
specific spatial displacement invariance of the fundamental degrees of freedom
of the statistical many-body physics, cf. Eqs. (3) and (4). Originally
formulated as a global operation hermann2021noether ; hermann2022topicalReview
; hermann2022variance  and subsequently generalized to local shifting
tschopp2022forceDFT ; hermann2022quantum , the approach yields new insight
into the correlated pair structure of soft matter sammueller2023whatIsLiquid .

We have here provided a detailed account of the theory, giving derivations,
including the treatment of all self-contributions as well as alternative
routes of proof explicitly. We have carefully treated the contributions that
arise from the kinetic energy part of the Hamiltonian and have shown how these
reduce in equilibrium to simple diffusive terms. As illustrations of the
framework, we have discussed computer simulation results for the Lennard-Jones
system in all its three phases, for the Yukawa liquid, a soft-sphere dipolar
fluid, the Stockmayer fluid, as well as the Gay-Berne model in the isotropic
and the nematic phase. We have also shown results for the Weeks-Chandler-
Andersen liquid, for monatomic water molinero2009 ; coe2022water , and for the
many-body colloidal gel former saw2009 ; saw2011 ; sammueller2023gel .

The results indicate that the force-based pair correlation functions reveal
much additional insight into the spatial structure over what is provided by
$g(r)$ alone. Crucially, neither their computational cost in simulations nor
the graphical representations as mere distance-dependent plots pose any
practical difficulties. In our present investigation the computationally most
expensive simulation procedure is to calculate the force gradients of the
three-body interacting Hamiltonians via numerical finite differences. As we
have shown it can be sufficient to only consider the two relevant tensor
components, i.e., the radial and transversal directions relative to the pair
distance vector. Hence in terms of simplicity, much of the appeal of $g(r)$ is
retained, while the basis in the Noether invariance also makes for a well-
grounded statistical mechanical foundation.

Our results add to the body of sum rules in statistical mechanics baus1984 ;
evans1990 ; henderson1992 ; triezenberg1972 and they share formal similarities
with sum rules for interface Hamiltonians mikheev1991 , see Ref. squarcini2022
for recent work, and with the Takahashi-Ward identities takahashi1957 ;
ward1950 of quantum field theory. As locally resolved interparticle force
measurements have been demonstrated in colloidal systems dong2022 , future
experimental use of the Noether force correlation functions looks feasible.
Relating to force-sampling methods that reduce the statistical variance
inherent in sampling results borgis2013 ; rotenberg2020 ;
delasheras2018forceSampling ; purohit2019 ; coles2019 ; coles2021 ;
coles2023revelsMD is another interesting point for future work, as can be
relating to force-based density functional theory tschopp2022forceDFT ;
sammueller2022forceDFT . The recent neural functional theory
delasheras2023perspective ; sammueller2023neural ; sammueller2023whyNeural ;
sammueller2023neuralTutorial is based on the powerful concept of using neural
networks to represent functional relationships which encapsulate the
correlation behaviour of complex systems. Noether sum rules have been shown to
provide valuable consistency checks for these neural functionals and they give
much inspiration for further theoretical developments in the spirit of
physics-informed machine learning karniadakis2021 ; ciarella2022 ; janzen2023
; jung2023piml ; rodrigues2023 ; wu2023review .

For the reason of being fully explicit, we have formulated the theory in the
grand ensemble, as is certainly common for the present type of considerations.
However, the entirety of the argumentation is applicable to a canonical
treatment, where the particle number is fixed. A simple way to see the
equivalence is to recognize that we have never relied on or indeed exploited
the grand ensemble structure; no derivatives with respect to $\mu$ are taken
or similar operations being carried out. Hence the Noether correlation theory
does not interfere with approaches that specifically target on particle number
effects, such as the local compressibility introduced by Evans and coworkers
evans2015jpcm ; evans2019pnas ; coe2022prl . The Noether approach captures
genuinely mechanical fluctuations, as opposed to local chemical (particle
number) evans2015jpcm ; evans2019pnas ; coe2022prl ; eckert2020 ;
eckert2023fluctuation and thermal (energy) fluctuations coe2022prl ;
eckert2020 ; eckert2023fluctuation . We lastly point to the recent hyperforce
generalization robitschko2023any that arises from thermal Noether invariance
and that applies to arbitrary observables.

###### Acknowledgements.

We thank Silas Robitschko and Daniel de las Heras for useful discussions.

##  Appendix A Kinetic stress autocorrelation

The Noether identity (22) contains a singular self contribution which arises
from the autocorrelation of the kinetic stress (8) on the left hand side and
from the kinetic energy on the right hand side. We treat this singular
contribution to the sum rule in the following by first considering the kinetic
stress autocorrelator, as it is contained in the form of a force
autocorrelator on the left hand side of Eq. (22) and correspondingly occuring
(multiplied by $\beta$) on the left hand side of Eq. (46), i.e.,

$\displaystyle\beta^{2}\Big{\langle}\nabla\cdot\sum_{i}\frac{{\bf p}_{i}{\bf
p}_{i}}{m}\delta({\mathbf{r}}-{\mathbf{r}}_{i})\nabla^{\prime}\cdot\sum_{j}\frac{{\bf
p}_{j}{\bf
p}_{j}}{m}\delta({\mathbf{r}}^{\prime}-{\mathbf{r}}_{j})\Big{\rangle}.$ (68)

The self part of Eq. (68) consists only of the cases $i=j$ in the double sum.
We can hence simplify Eq. (68) according to:

$\displaystyle\nabla\nabla^{\prime}:\beta^{2}\Big{\langle}\sum_{i}\frac{{\bf
p}_{i}{\bf p}_{i}{\bf p}_{i}{\bf
p}_{i}}{m^{2}}\delta({\mathbf{r}}-{\mathbf{r}}_{i})\delta({\mathbf{r}}^{\prime}-{\mathbf{r}}_{i})\Big{\rangle}.$
(69)

Switching to index notation and denoting Cartesian components by Greek
indices, the $\alpha\gamma$-component of this second-rank tensor is:

$\displaystyle\sum_{\nu\xi}\nabla_{\nu}\nabla^{\prime}_{\xi}\beta^{2}\Big{\langle}\sum_{i}\frac{p_{i}^{\alpha}p_{i}^{\gamma}p_{i}^{\nu}p_{i}^{\xi}}{m^{2}}\delta({\mathbf{r}}-{\mathbf{r}}_{i})\delta({\mathbf{r}}^{\prime}-{\mathbf{r}}_{i})\Big{\rangle},$
(70)

where the sums over the Greek indices run over all Cartesian components.
Carrying out the momentum average over the momentum tetradic in the integrand
requires the following integral over the Maxwell distribution

$\displaystyle\frac{\beta^{2}}{m^{2}}\int d{\bf p}_{i}$
$\displaystyle\frac{{\rm e}^{-\beta{\bf p}_{i}^{2}/(2m)}}{(2\pi
m/\beta)^{3/2}}p_{i}^{\alpha}p_{i}^{\gamma}p_{i}^{\nu}p_{i}^{\xi}$
$\displaystyle=\delta_{\alpha\nu}\delta_{\gamma\xi}+\delta_{\alpha\xi}\delta_{\gamma\nu}+\delta_{\alpha\gamma}\delta_{\nu\xi},$
(71)

where the combination of the six Kronecker delta symbols on the right hand
side is the isotropic tensor of rank four. The Maxwellian in Eq. (71) is
normalized, $\int d{\bf p}_{i}{\rm e}^{-\beta{\bf p}_{i}^{2}/(2m)}/(2\pi
m/\beta)^{3/2}=1$. In the derivation of Eq. (71) one requires both the second
moment of each ($\alpha$th) momentum component, $\int d{\bf
p}_{i}(p_{i}^{\alpha})^{2}{\rm e}^{-\beta{\bf p}_{i}^{2}/(2m)}/(2\pi
m/\beta)^{3/2}=m/\beta$, and the fourth moment, $\int d{\bf
p}_{i}(p_{i}^{\alpha})^{4}{\rm e}^{-\beta{\bf p}_{i}^{2}/(2m)}/(2\pi
m/\beta)^{3/2}=3m^{2}/\beta^{2}$.

Recalling the derivative structure of Eq. (70) we contract Eq. (71) with
$\nabla_{\nu}\nabla_{\xi}^{\prime}$, which yields

$\displaystyle\sum_{\nu\xi}\nabla_{\nu}\nabla^{\prime}_{\xi}$
$\displaystyle(\delta_{\alpha\nu}\delta_{\gamma\xi}+\delta_{\alpha\xi}\delta_{\gamma\nu}+\delta_{\alpha\gamma}\delta_{\nu\xi})$
$\displaystyle=(\nabla_{\alpha}\nabla^{\prime}_{\gamma}+\nabla^{\prime}_{\alpha}\nabla_{\gamma}+\delta_{\alpha\gamma}\nabla\cdot\nabla^{\prime})$
(72)
$\displaystyle=(\nabla\nabla^{\prime}+\nabla^{\prime}\nabla+{\mathbbm{1}}\nabla\cdot\nabla^{\prime})_{\alpha\gamma}.$
(73)

Using this result we can simplify the expression (70) as

$\displaystyle(\nabla\nabla^{\prime}+\nabla^{\prime}\nabla+{\mathbbm{1}}\nabla\cdot\nabla^{\prime})\Big{\langle}\sum_{i}\delta({\mathbf{r}}-{\mathbf{r}}_{i})\delta({\mathbf{r}}^{\prime}-{\mathbf{r}}_{i})\Big{\rangle}.$
(74)

Hence in summary we obtain the following self contribution:

$\displaystyle\Big{\langle}\nabla\cdot\beta\hat{\bm{\tau}}({\mathbf{r}})$
$\displaystyle\nabla^{\prime}\cdot\beta\hat{\bm{\tau}}({\mathbf{r}}^{\prime})\Big{\rangle}_{\rm
self}$
$\displaystyle=(\nabla\nabla^{\prime}+\nabla^{\prime}\nabla+{\mathbbm{1}}\nabla\cdot\nabla^{\prime})\rho_{2}^{\rm
self}({\mathbf{r}},{\mathbf{r}}^{\prime}),$ (75)

where we recall the definition (34) of the two-body self density:
$\rho_{2}^{\rm
self}({\mathbf{r}},{\mathbf{r}}^{\prime})=\langle\sum_{i}\delta({\mathbf{r}}-{\mathbf{r}}_{i})\delta({\mathbf{r}}^{\prime}-{\mathbf{r}}_{i})\rangle$.

We next consider the distinct part of Eq. (68), which consists of the cases
$i\neq j$. In index notation its $\alpha\gamma$-component is

$\displaystyle\sum_{\nu\xi}\nabla_{\nu}\nabla^{\prime}_{\xi}\beta^{2}\Big{\langle}\sum_{ij(\neq)}\frac{p^{\alpha}_{i}p^{\gamma}_{j}p^{\nu}_{i}p^{\xi}_{j}}{m^{2}}\delta({\mathbf{r}}-{\mathbf{r}}_{i})\delta({\mathbf{r}}^{\prime}-{\mathbf{r}}_{j})\Big{\rangle}.$
(76)

As before we carry out the momentum average which requires the following
momentum integrals of particles $i$ and $j$:

$\displaystyle\frac{\beta^{2}}{m^{2}}\int d{\bf p}_{i}\frac{e^{-\beta{\bf
p}_{i}^{2}/(2m)}}{(2\pi m/\beta)^{3/2}}p^{\alpha}_{i}p^{\nu}_{i}\int d{\bf
p}_{j}\frac{e^{-\beta{\bf p}_{j}^{2}/(2m)}}{(2\pi
m/\beta)^{3/2}}p^{\gamma}_{j}p^{\xi}_{j}$
$\displaystyle\quad=\delta_{\alpha\nu}\delta_{\gamma\xi}.$ (77)

Using this result we can re-write Eq. (76) upon simplifying according to
$\sum_{\nu\xi}\nabla_{\nu}\nabla^{\prime}_{\xi}\delta_{\alpha\nu}\delta_{\gamma\xi}=\nabla_{\alpha}\nabla^{\prime}_{\gamma}$
to hence determine the distinct contribution of the kinetic force
autocorrelator as:

$\displaystyle\Big{\langle}\nabla\cdot\beta\hat{\bm{\tau}}({\mathbf{r}})\nabla^{\prime}\cdot\beta\hat{\bm{\tau}}({\mathbf{r}}^{\prime})\Big{\rangle}_{\text{dist}}$
$\displaystyle\qquad=\nabla\nabla^{\prime}\Big{\langle}\sum_{ij(\neq)}\delta({\mathbf{r}}-{\mathbf{r}}_{i})\delta({\mathbf{r}}^{\prime}-{\mathbf{r}}_{j})\Big{\rangle}$
(78)
$\displaystyle\qquad=\nabla\nabla^{\prime}\rho_{2}({\mathbf{r}},{\mathbf{r}}^{\prime}),$
(79)

where we have used the definition (26) of the distinct two-body density.

Hence adding the self contribution (75) and the distinct results (79) we
obtain the desired total kinetic force autocorrelator, multiplied by
$\beta^{2}$, as the following expression:

$\displaystyle\Big{\langle}\nabla\cdot\beta\hat{\bm{\tau}}({\mathbf{r}})\nabla^{\prime}\cdot\beta\hat{\bm{\tau}}({\mathbf{r}}^{\prime})\Big{\rangle}$
$\displaystyle=(\nabla\nabla^{\prime}+\nabla^{\prime}\nabla+{\mathbbm{1}}\nabla\cdot\nabla^{\prime})\rho_{2}^{\rm
self}({\mathbf{r}},{\mathbf{r}}^{\prime})+\nabla\nabla^{\prime}\rho_{2}({\mathbf{r}},{\mathbf{r}}^{\prime}),$
(80)

which we reproduce as Eq. (46) in Sec. V of the main text.

##  Appendix B Kinetic energy curvature

To second order the momentum transformation (4) can be written as

$\displaystyle{\bf
p}_{i}\to\\{{\mathbbm{1}}-\nabla_{i}{\bm{\epsilon}}({\mathbf{r}}_{i})+[\nabla_{i}{\bm{\epsilon}}({\mathbf{r}}_{i})]^{2}\\}\cdot{\bf
p}_{i}.$ (81)

As a consequence, the squared momentum, as is relevant for the kinetic energy
$H_{\rm kin}=\sum_{i}{\bf p}_{i}^{2}/(2m)$, becomes

$\displaystyle{\bf p}_{i}\cdot{\bf p}_{i}\to$
$\displaystyle[{\mathbbm{1}}-2\nabla_{i}{\bm{\epsilon}}_{i}+(\nabla_{i}{\bm{\epsilon}}_{i})^{\sf
T}\cdot\nabla_{i}{\bm{\epsilon}}_{i},+2(\nabla_{i}{\bm{\epsilon}}_{i})^{2}]:{\bf
p}_{i}{\bf p}_{i},$ (82)

where the colon indicates a double tensor contraction of two matrices ${\sf
A}$ and ${\sf B}$, i.e. ${\sf A}:{\sf B}=\sum_{\alpha\gamma}{\sf
A}_{\alpha\gamma}{\sf B}_{\gamma\alpha}$, with $\alpha,\gamma$ denoting the
Cartesian components; we have used the shorthand
${\bm{\epsilon}}_{i}={\bm{\epsilon}}({\mathbf{r}}_{i})$ and the superscript
$\sf T$ indicates transposition.

We can use the following identities for the linear and quadratic contribution
in the expansion (82). For the linear contribution, we have

$\displaystyle[\nabla_{i}{\bm{\epsilon}}({\mathbf{r}}_{i})]^{\sf T}$
$\displaystyle=-\int
d{\mathbf{r}}{\bm{\epsilon}}({\mathbf{r}})\nabla\delta({\mathbf{r}}-{\mathbf{r}}_{i}),$
(83)
$\displaystyle\Big{(}\frac{\delta}{\delta{\bm{\epsilon}}({\mathbf{r}})}\nabla_{i}{\bm{\epsilon}}({\mathbf{r}}_{i})\Big{)}_{\alpha\gamma\nu}$
$\displaystyle=-\delta_{\alpha\nu}\nabla_{\gamma}\delta({\mathbf{r}}-{\mathbf{r}}_{i}),$
(84)

where as above the Greek indices denote the Cartesian components. We first
consider the following quadratic contribution:

$\displaystyle[\nabla_{i}{\bm{\epsilon}}({\mathbf{r}}_{i})]^{\sf
T}\cdot\nabla_{i}{\bm{\epsilon}}({\mathbf{r}}_{i})$ $\displaystyle=$
$\displaystyle\int
d{\mathbf{r}}d{\mathbf{r}}^{\prime}{\bm{\epsilon}}({\mathbf{r}}){\bm{\epsilon}}({\mathbf{r}}^{\prime})$
$\displaystyle\nabla\cdot\nabla^{\prime}\delta({\mathbf{r}}-{\mathbf{r}}_{i})\delta({\mathbf{r}}^{\prime}-{\mathbf{r}}_{i}),$
(85)

of which the second functional derivative is:

$\displaystyle\Big{(}\frac{\delta^{2}}{\delta{\bm{\epsilon}}({\mathbf{r}})\delta{\bm{\epsilon}}({\mathbf{r}}^{\prime})}\big{\\{}[\nabla_{i}{\bm{\epsilon}}({\mathbf{r}}_{i})]^{\sf
T}\cdot\nabla_{i}{\bm{\epsilon}}({\mathbf{r}}_{i})\big{\\}}\Big{)}_{\alpha\gamma\nu\xi}$
$\displaystyle=$ $\displaystyle
2\nabla\cdot\nabla^{\prime}\delta({\mathbf{r}}-{\mathbf{r}}_{i})$
$\displaystyle\delta({\mathbf{r}}^{\prime}-{\mathbf{r}}_{i})\delta_{\alpha\nu}\delta_{\gamma\xi},$
(86)

where as before $\nabla^{\prime}$ indicates the derivative with respect to
${\mathbf{r}}^{\prime}$. The further quadratic term is:

$\displaystyle[\nabla_{i}{\bm{\epsilon}}({\mathbf{r}}_{i})]^{2}=\int
d{\mathbf{r}}d{\mathbf{r}}^{\prime}{\bm{\epsilon}}({\mathbf{r}}){\bm{\epsilon}}({\mathbf{r}}^{\prime})\cdot\nabla\nabla^{\prime}\delta({\mathbf{r}}-{\mathbf{r}}_{i})\delta({\mathbf{r}}^{\prime}-{\mathbf{r}}_{i}).$
(87)

and its second functional derivative with respect to the shifting field is:

$\displaystyle\Big{(}\frac{\delta^{2}}{\delta{\bm{\epsilon}}({\mathbf{r}})\delta{\bm{\epsilon}}({\mathbf{r}}^{\prime})}[\nabla_{i}{\bm{\epsilon}}({\mathbf{r}}_{i})]^{2}\Big{)}_{\alpha\gamma\nu\xi}=$
$\displaystyle\quad\big{(}\delta_{\alpha\xi}\nabla_{\gamma}\nabla_{\nu}^{\prime}+\delta_{\gamma\xi}\nabla_{\nu}\nabla_{\alpha}^{\prime}\big{)}\delta({\mathbf{r}}-{\mathbf{r}}_{i})\delta({\mathbf{r}}^{\prime}-{\mathbf{r}}_{i}).$
(88)

We next build the second functional derivative of the transformed kinetic
energy:

$\displaystyle\frac{\delta^{2}H_{\rm
kin}[{\bm{\epsilon}}]}{\delta{\bm{\epsilon}}({\mathbf{r}})\delta{\bm{\epsilon}}({\mathbf{r}}^{\prime})}$
$\displaystyle=\frac{\delta^{2}}{\delta{\bm{\epsilon}}({\mathbf{r}})\delta{\bm{\epsilon}}({\mathbf{r}}^{\prime})}\sum_{i}[(\nabla_{i}{\bm{\epsilon}}_{i})^{\sf
T}\cdot\nabla_{i}{\bm{\epsilon}}_{i}+2(\nabla_{i}{\bm{\epsilon}}_{i})^{2})]:\frac{{\bf
p}_{i}{\bf p}_{i}}{2m}$ (89) $\displaystyle=\sum_{i}\Big{(}\frac{{\bf
p}_{i}{\bf
p}_{i}}{m}\nabla\cdot\nabla^{\prime}+\nabla^{\prime}\nabla\cdot\frac{{\bf
p}_{i}{\bf p}_{i}}{m}$ $\displaystyle\qquad\qquad+\frac{{\bf p}_{i}{\bf
p}_{i}}{m}\cdot\nabla^{\prime}\nabla\Big{)}\delta({\mathbf{r}}-{\mathbf{r}}_{i})\delta({\mathbf{r}}^{\prime}-{\mathbf{r}}_{i}),$
(90)

where we have used Eqs. (86) and (88) in the second step and the factors of 2
have cancelled. The double contraction of a fourth- and a second-rank tensor
is thereby defined as $({\sf A}:{\sf B})_{\alpha\gamma}=\sum_{\nu\xi}{\sf
A}_{\alpha\gamma\nu\xi}{\sf B}_{\xi\nu}$. Building the thermal equilibrium
average of Eq. (90) and calculating thereby the momentum integrals explicitly
yields upon simplifying the following final compact result:

$\displaystyle\Big{\langle}$ $\displaystyle\frac{\delta^{2}H_{\rm
kin}[{\bm{\epsilon}}]}{\delta{\bm{\epsilon}}({\mathbf{r}})\delta{\bm{\epsilon}}({\mathbf{r}}^{\prime})}\Big{|}_{{\bm{\epsilon}}=0}\Big{\rangle}$
$\displaystyle=k_{B}T({\mathbbm{1}}\nabla\cdot\nabla^{\prime}+2\nabla^{\prime}\nabla)\rho_{2}^{\rm
self}({\mathbf{r}},{\mathbf{r}}^{\prime}),$ (91)

which is reproduced as Eq. (50) in the main text.

##  Appendix C Potential energy curvature

We consider the potential energy contribution to the Hamiltonian
$H_{U}=u({\mathbf{r}}^{N})+\sum_{i}V_{\rm ext}({\mathbf{r}}_{i})$, such that
the total Hamiltonian (2) is $H=H_{\rm kin}+H_{U}$. Expressed in the
transformed variables we have

$\displaystyle H_{U}[{\bm{\epsilon}}]$
$\displaystyle=u\big{(}{\mathbf{r}}_{1}+{\bm{\epsilon}}({\mathbf{r}}_{1}),\ldots,{\mathbf{r}}_{N}+{\bm{\epsilon}}({\mathbf{r}}_{N})\big{)}$
$\displaystyle\quad+\sum_{i}V_{\rm
ext}({\mathbf{r}}_{i}+{\bm{\epsilon}}({\mathbf{r}}_{i})).$ (92)

The first functional derivative of Eq. (92) is

$\displaystyle\frac{\delta
H_{U}[{\bm{\epsilon}}]}{\delta{\bm{\epsilon}}({\mathbf{r}}^{\prime})}$
$\displaystyle=\sum_{j}\delta({\mathbf{r}}^{\prime}-{\mathbf{r}}_{j})\nabla_{j}u\big{(}{\mathbf{r}}_{1}+{\bm{\epsilon}}({\mathbf{r}}_{1}),\ldots,{\mathbf{r}}_{N}+{\bm{\epsilon}}({\mathbf{r}}_{N})\big{)}$
$\displaystyle\quad+\sum_{i}\delta({\mathbf{r}}^{\prime}-{\mathbf{r}}_{i})\nabla_{i}V_{\rm
ext}({\mathbf{r}}_{i}+{\bm{\epsilon}}({\mathbf{r}}_{i})),$ (93)

where the Dirac distribution is generated from the functional derivative
$\delta[{\mathbf{r}}_{j}+{\bm{\epsilon}}({\mathbf{r}}_{j})]/\delta{\bm{\epsilon}}({\mathbf{r}}^{\prime})=\delta({\mathbf{r}}^{\prime}-{\mathbf{r}}_{j}){\mathbbm{1}}$
and the unit matrix then disapears upon
${\mathbbm{1}}\cdot\nabla_{j}=\nabla_{j}$. The functional curvature with
respect to the shifting field is then obtained as the functional Hessian
(second derivative) of Eq. (92) or analogously as the first functional
derivative of Eq. (93):

$\displaystyle\frac{\delta^{2}H_{U}[{\bm{\epsilon}}]}{\delta{\bm{\epsilon}}({\mathbf{r}})\delta{\bm{\epsilon}}({\mathbf{r}}^{\prime})}\Big{|}_{{\bm{\epsilon}}=0}$
$\displaystyle=\sum_{ij}\delta({\mathbf{r}}-{\mathbf{r}}_{i})\delta({\mathbf{r}}^{\prime}-{\mathbf{r}}_{j})\nabla_{i}\nabla_{j}u({\mathbf{r}}^{N})$
$\displaystyle\quad+\sum_{i}\delta({\mathbf{r}}-{\mathbf{r}}_{i})\delta({\mathbf{r}}^{\prime}-{\mathbf{r}}_{i})\nabla_{i}\nabla_{i}V_{\rm
ext}({\mathbf{r}}_{i}),$ (94)

which follows similarly as Eq. (93). Building then the equilibrium average of
Eq. (94) yields

$\displaystyle\Big{\langle}\frac{\delta^{2}H_{U}[{\bm{\epsilon}}]}{\delta{\bm{\epsilon}}({\mathbf{r}})\delta{\bm{\epsilon}}({\mathbf{r}}^{\prime})}$
$\displaystyle\Big{|}_{{\bm{\epsilon}}=0}\Big{\rangle}=$
$\displaystyle\Big{\langle}\sum_{ij}\delta({\mathbf{r}}-{\mathbf{r}}_{i})\delta({\mathbf{r}}^{\prime}-{\mathbf{r}}_{j})\nabla_{i}\nabla_{j}u({\mathbf{r}}^{N})\Big{\rangle}$
$\displaystyle+\rho_{2}^{\rm
self}({\mathbf{r}},{\mathbf{r}}^{\prime})\nabla\nabla V_{\rm
ext}({\mathbf{r}}),$ (95)

where for the external potential term we have used that here
$\nabla_{i}\nabla_{i}=\nabla\nabla$. Together with the corresponding result
(91) for the kinetic energy, Eq. (95) forms the functional Hessian of the full
Hamiltonian (2) with respect to the shifting field. We reproduce the result
(95) as Eq. (51) in the main text.

##  Appendix D Sum rules via partial integration

As a means of independent verification, we describe an alternative method to
obtain the second-order sum rules, based on more elementary, but very explicit
calculations that do not rely on Noether invariance. On the downside the
thermal symmetry remains hidden.

The argumentation simply rests on suitably chosen integration by parts in
phase space. For a given observable $\hat{A}({\mathbf{r}}^{N},{\bf p}^{N})$
the underlying mechanism is based on two identities. One of them is the Yvon
theorem hansen2013 ; rotenberg2020 :

$\displaystyle\Big{\langle}\nabla_{i}\hat{A}\Big{\rangle}=\Big{\langle}\hat{A}\nabla_{i}\beta
H\Big{\rangle},$ (96)

which can be shown simply by writing out the expression on the left hand side
as $\langle\nabla_{i}\hat{A}\rangle={\rm Tr}\,\Xi^{-1}e^{-\beta(H-\mu
N)}\nabla_{i}\hat{A}=-{\rm Tr}\,\hat{A}\Xi^{-1}\nabla_{i}e^{-\beta(H-\mu
N)}={\rm Tr}\,\hat{A}\Xi^{-1}e^{-\beta(H-\mu N)}\nabla_{i}\beta
H=\langle\hat{A}\nabla_{i}\beta H\rangle$, which is the right hand side of Eq.
(96). In the second step we have assumed that any boundary terms from the
partial integration vanish, because the external potential contains e.g.
contributions from impenetrable container walls.

The second identity involves an additional phase space function $\hat{B}$. We
have

$\displaystyle\Big{\langle}\hat{A}\nabla_{i}\hat{B}\Big{\rangle}=\Big{\langle}\hat{A}\hat{B}\nabla_{i}\beta
H\Big{\rangle}-\Big{\langle}\hat{B}\nabla_{i}\hat{A}\Big{\rangle}.$ (97)

The proof of Eq. (97) rests on analogous argumentation as above, with
additionally taking into account the product rule of differentiation that
generates the second term on the right hand side of Eq. (97). Explicitly, we
start with the left hand side $\langle\hat{A}\nabla_{i}\hat{B}\rangle={\rm
Tr}\,\Xi^{-1}e^{-\beta(H-\mu N)}\hat{A}\nabla_{i}\hat{B}$ and integrate by
parts, which yields $-{\rm Tr}\,\hat{B}\Xi^{-1}\nabla_{i}e^{-\beta(H-\mu
N)}\hat{A}$. Carrying out the derivative gives two terms: ${\rm
Tr}\,\hat{B}\Xi^{-1}e^{-\beta(H-\mu N)}(\nabla_{i}\beta H)\hat{A}-{\rm
Tr}\,\hat{B}\Xi^{-1}e^{-\beta(H-\mu
N)}\nabla_{i}\hat{A}=\langle\hat{A}\hat{B}\nabla_{i}\beta
H\rangle-\langle\hat{B}\nabla_{i}\hat{A}\rangle$, where we have resorted back
to the compact thermal average notation. Upon ordering of the factors in the
first term, we have obtained the right hand side of Eq. (97), as desired.

We apply these formal results to the many-body physics under consideration. We
first consider the distinct case. Identifying the correct starting point is
thereby crucial and we take this to be the double gradient of the (distinct)
two-body density (26):

$\displaystyle\nabla\nabla^{\prime}\rho_{2}({\mathbf{r}},{\mathbf{r}}^{\prime})$
$\displaystyle=\nabla\nabla^{\prime}\Big{\langle}{{\sum_{ij(\neq)}}}\delta_{i}\delta_{j}^{\prime}\Big{\rangle}.$
(98)

As before $\nabla^{\prime}$ indicates the derivative with respect to
${\mathbf{r}}^{\prime}$, the double sum is denoted using the compact notation
$\sum_{ij(\neq)}=\sum_{i=1}^{N}\sum_{j=1,j\neq i}^{N}$ and the delta functions
are abbreviated as $\delta_{i}=\delta({\mathbf{r}}-{\mathbf{r}}_{i})$ and
$\delta_{j}^{\prime}=\delta({\mathbf{r}}^{\prime}-{\mathbf{r}}_{j})$. Due to
the identity
$\nabla\delta({\mathbf{r}}-{\mathbf{r}}_{i})=-\nabla_{i}\delta({\mathbf{r}}-{\mathbf{r}}_{i})$
we can re-write the right hand side of Eq. (98) as

$\displaystyle\Big{\langle}\sum_{ij(\neq)}\nabla_{i}\delta_{i}\nabla_{j}\delta_{j}^{\prime}\Big{\rangle}$
$\displaystyle=\Big{\langle}\sum_{ij(\neq)}(\nabla_{i}\beta
H)\delta_{i}\nabla_{j}\delta_{j}^{\prime}\Big{\rangle}$ (99)
$\displaystyle=\Big{\langle}\sum_{ij(\neq)}(\nabla_{i}\beta H)(\nabla_{j}\beta
H)\delta_{i}\delta_{j}^{\prime}\Big{\rangle}$
$\displaystyle\quad-\Big{\langle}{\sum_{ij(\neq)}}(\nabla_{i}\nabla_{j}\beta
H)\delta_{i}\delta_{j}^{\prime}\Big{\rangle},$ (100)

where we have integrated by parts in the first step with respect to
${\mathbf{r}}_{i}$ according to Eq. (96) and in the second step with respect
to ${\mathbf{r}}_{j}$ according to Eq. (97).

We regroup the double sum inside of the first average in Eq. (100) by
factorizing $\sum_{ij(\neq)}(\nabla_{i}\beta H)(\nabla_{j}\beta
H)\delta_{i}\delta_{j}^{\prime}=\sum_{i}\delta_{i}(\nabla_{i}\beta
H)\sum_{j\neq i}\delta_{j}^{\prime}\nabla_{j}\beta H$. Taking account of the
structure of the Hamiltonian (2) then allows to identify the expression
$-\sum_{i}\delta_{i}\nabla_{i}H=\hat{\bf F}_{U}({\mathbf{r}})$ as the local
potential force density operator that arises from both interparticle and
external effects, see Eq. (23) for the definition of $\hat{\bf
F}_{U}({\mathbf{r}})$.

We can now express Eq. (100) succinctly and hence obtain the distinct part
(24) of the curvature sum rule, which for completeness we supplement by re-
writing the self part (25) (to be proven below):

$\displaystyle\Big{\langle}\beta\hat{\bf F}_{U}({\mathbf{r}})\beta\hat{\bf
F}_{U}({\mathbf{r}}^{\prime})\Big{\rangle}_{\\!\rm dist}=$
$\displaystyle\qquad\qquad\nabla\nabla^{\prime}\rho_{2}^{\rm
dist}({\mathbf{r}},{\mathbf{r}}^{\prime})+\Big{\langle}{\sum_{ij(\neq)}}(\nabla_{i}\nabla_{j}\beta
u)\delta_{i}\delta_{j}^{\prime}\Big{\rangle},$ (101)
$\displaystyle\Big{\langle}\beta\hat{\bf F}_{U}({\mathbf{r}})\beta\hat{\bf
F}_{U}({\mathbf{r}})\Big{\rangle}_{\\!\rm self}=$
$\displaystyle\qquad\quad\nabla\nabla\rho({\mathbf{r}})+\rho({\mathbf{r}})\nabla\nabla\beta
V_{\rm ext}({\mathbf{r}})+\Big{\langle}\sum_{i}(\nabla_{i}\nabla_{i}\beta
u)\delta_{i}\Big{\rangle}.$ (102)

The identity (101) holds for any ${\mathbf{r}},{\mathbf{r}}^{\prime}$ and Eq.
(102) is the corresponding self sum rule, as derived in the following by a
proof based on partial phase space integration, which is similar to the
distinct case above.

Analogously to the starting point (98) for the distinct correlation identity,
we consider the double gradient, but with respect to the same position
${\mathbf{r}}$ taken twice, i.e. the Hessian of the density profile:

$\displaystyle\nabla\nabla\rho({\mathbf{r}})$
$\displaystyle=\nabla\nabla\Big{\langle}\sum_{i}\delta_{i}\Big{\rangle}$ (103)
$\displaystyle=\Big{\langle}\sum_{i}\nabla_{i}\nabla_{i}\delta_{i}\Big{\rangle}$
(104) $\displaystyle=\Big{\langle}\sum_{i}(\nabla_{i}\beta
H)\nabla_{i}\delta_{i}\Big{\rangle}$ (105)
$\displaystyle=\Big{\langle}\sum_{i}(\nabla_{i}\beta H)(\nabla_{i}\beta
H)\delta_{i}\Big{\rangle}$
$\displaystyle\quad-\Big{\langle}\sum_{i}(\nabla_{i}\nabla_{i}\beta
H)\delta_{i}\Big{\rangle}.$ (106)

In this derivation we have exploited twice that
$\nabla\delta({\mathbf{r}}-{\mathbf{r}}_{i})=-\nabla_{i}\delta({\mathbf{r}}-{\mathbf{r}}_{i})$
and then integrated by parts first according to Eq. (96) and then according to
Eq. (97). Identifying $\hat{\bf F}_{U}({\mathbf{r}})$ in Eq. (106) and re-
arranging the different terms gives Eq. (102).

The present phase space integration route to the sum rules (101) and (102) [or
analogously Equations (24) and (25)] is free of any functional calculus as
required in the derivations in the main text. However, the required chain of
these individual steps is not easy to guess and apparently these, as we argue,
very fundamental results have not been written down in the existing
literature. Hence, while the Noether route comes at the expense of having to
engage with some functional calculus, the very significant benefit is the
simplicity of the starting point, which here is the second-order invariance
(20) of the grand potential against spatially inhomogeneous displacement.

##  Appendix E Measuring force-force and force-gradient correlations

We here provide technical details for the sampling of the force-force and
force-gradient correlations as given in Eqs. (57) and (58). Similar to the
measurement of the standard pair correlation function $g(r)$, a loop over
particle pairs is performed in order to calculate an average of a bulk two-
body quantity. However, one not only counts the relative number of particles
separated by a certain distance $r$, but instead also involves the forces and
force gradients that act on the particles in the calculation of the average.

Let ${\mathbf{r}}_{i}$ and ${\mathbf{r}}_{j}$ be the positions of two
particles, such that
${\mathbf{r}}_{\parallel}={\mathbf{r}}_{j}-{\mathbf{r}}_{i}$ and ${\bf
e}_{\parallel}={\mathbf{r}}_{\parallel}/|{\mathbf{r}}_{\parallel}|$. The
calculation of ${\sf g}_{f\\!f}(r)$ proceeds straightforwardly by
incorporating the forces ${\bf f}_{i}$ and ${\bf f}_{j}$ of particles $i$ and
$j$ as given in Eq. (57). The radial and tangential components are
respectively obtained by the following projections:

$\displaystyle({\bf f}_{i}{\bf f}_{j})_{\parallel}$ $\displaystyle=({\bf
f}_{i}\cdot{\bf e}_{\parallel})({\bf f}_{j}\cdot{\bf e}_{\parallel}),$ (107)
$\displaystyle({\bf f}_{i}{\bf f}_{j})_{\perp}$ $\displaystyle=\frac{[{\bf
f}_{i}-({\bf f}_{i}\cdot{\bf e}_{\parallel}){\bf e}_{\parallel}]\cdot[{\bf
f}_{j}-({\bf f}_{j}\cdot{\bf e}_{\parallel}){\bf e}_{\parallel}]}{2},$ (108)

where the factor 1/2 in Eq. (108) accounts for the two equal tangential
contributions of ${\sf g}_{f\\!f}(r)$ in a three-dimensional bulk fluid (there
is only a single radial component).

For the calculation of ${\sf g}_{\nabla\\!f}(r)$, we employ numerical
differentiation to evaluate the force gradients $\nabla_{i}{\bf f}_{j}$ that
appear in Eq. (58), which is simpler to implement than analytic Hessians of
the interaction potential, in particular for the more complex models such as
the Stillinger-Weber or Gay-Berne potentials.

As before, a splitting of the force-gradient correlation function ${\sf
g}_{\nabla\\!f}(r)$ into its radial and tangential components is performed.
Assuming ${\bf f}_{j}\nparallel{\bf e}_{\parallel}$, a tangential unit vector
${\bf e}_{\perp}={\bf t}/|{\bf t}|\perp{\bf e}_{\parallel}$ can be constructed
by choosing ${\bf t}={\bf f}_{j}-({\bf f}_{j}\cdot{\bf e}_{\parallel}){\bf
e}_{\parallel}$. The radial and tangential parts of ${\sf g}_{\nabla\\!f}(r)$
then follow via

$\displaystyle(\nabla_{i}{\bf f}_{j})_{\parallel}$ $\displaystyle=D_{i,{\bf
e}_{\parallel}}({\bf f}_{j}\cdot{\bf e}_{\parallel}),$ (109)
$\displaystyle(\nabla_{i}{\bf f}_{j})_{\perp}$ $\displaystyle=D_{i,{\bf
e}_{\perp}}({\bf f}_{j}\cdot{\bf e}_{\perp}),$ (110)

where $D_{i,{\bf e}}$ is a directional derivative operator regarding the
shifting of particle $i$ along the (generic) direction ${\bf e}$. The
evaluation of Eqs. (109) and (110) is performed numerically in the
simulations, i.e., particle $i$ is shifted by a small amount ($\sim
10^{-5}\sigma$) to calculate the arising finite difference in the force ${\bf
f}_{j}$ of particle $j$.

The above procedure applies to isotropic particles. For anisotropic particles
a simple alternative to construct a suitable vector ${\bf e}_{\perp}$ can be
based on drawing a random unit vector ${\bf e}_{\rm ran}\nparallel{\bf
e}_{\parallel}$. Then ${\bf e}_{\perp}={\bf t}/|{\bf t}|$ with ${\bf t}={\bf
e}_{\rm ran}-({\bf e}_{\rm ran}\cdot{\bf e}_{\parallel}){\bf e}_{\parallel}$,
which can be used in Eq. (110).

## References

  * (1) J. P. Hansen and I. R. McDonald, Theory of Simple Liquids, 4th ed. (Academic Press, London, 2013).
  * (2) D. Frenkel and B. Smit, Understanding Molecular Simulation: From Algorithms to Applications, 3rd ed. (Academic Press, 2023).
  * (3) C. P. Royall, A. A. Louis, and H. Tanaka, Measuring colloidal interactions with confocal microscopy, J. Chem. Phys. 127, 044507 (2007).
  * (4) A. L. Thorneywork, R. Roth, D. G. A. L. Aarts, and R. P. A. Dullens, Communication: Radial distribution functions in a two-dimensional binary colloidal hard sphere system, J. Chem. Phys. 140, 161106 (2014).
  * (5) A. Statt, R. Pinchaipat, F. Turci, R. Evans, and C. P. Royall, Direct observation in 3d of structural crossover in binary hard sphere mixtures, J. Chem. Phys. 144, 144506 (2016).
  * (6) R. Evans, J. R. Henderson, D. C. Hoyle, A. O. Parry, and Z. A. Sabeur, Asymptotic decay of liquid structure: oscillatory liquid-vapour density profiles and the Fisher-Widom line, Mol. Phys. 80, 755 (1993).
  * (7) R. Evans, R. J. F. Leote de Carvalho, J. R. Henderson, and D. C. Hoyle, Asymptotic decay of correlations in liquids and their mixtures, J. Chem. Phys. 100, 591 (1994).
  * (8) M. Dijkstra and R. Evans, A simulation study of the decay of the pair correlation function in simple fluids, J. Chem. Phys. 112, 1449 (2000).
  * (9) C. Grodon, M. Dijkstra, R. Evans, and R. Roth, Decay of correlation functions in hard-sphere mixtures: Structural crossover, J. Chem. Phys. 121, 7869 (2004).
  * (10) P. Cats, R. Evans, A. Härtel, and R. van Roij, Primitive model electrolytes in the near and far field: Decay lengths from DFT and simulations, J. Chem. Phys. 154, 124504 (2021).
  * (11) Z. Zhang and W. Kob, Revealing the three-dimensional structure of liquids using four-point correlation functions, Proc. Nat. Acad. Sci. 117, 14032 (2020).
  * (12) N. Singh, Z. Zhang, A. K. Sood, W. Kob, and R. Ganapathy, Intermediate-range order governs dynamics in dense colloidal liquids, Proc. Nat. Acad. Sci. 120, e2300923120 (2023).
  * (13) H. Yuan, Z. Zhang, W. Kob, and Y. Wang, Connecting packing efficiency of binary hard sphere systems to their intermediate range structure, Phys. Rev. Lett. 127, 278001 (2021).
  * (14) I. Pihlajamaa, C. C. L. Laudicina, C. Luo, L. M. C. Janssen, Emergent structural correlations in dense liquids, PNAS Nexus 2, pgad184 (2023).
  * (15) E. Boattini, F. Smallenburg, and L. Filion, Averaging local structure to predict the dynamic propensity in supercooled liquids, Phys. Rev. Lett. 127, 088007 (2021).
  * (16) E. Noether, Invariante Variationsprobleme, Nachr. d. König. Gesellsch. d. Wiss. zu Göttingen, Math.-Phys. Klasse, 235, 183 (1918). English translation by M. A. Tavel: Invariant variation problems. Transp. Theo. Stat. Phys. 1, 186 (1971); for a version in modern typesetting see: Frank Y. Wang, arXiv:physics/0503066v3 (2018).
  * (17) N. Byers, E. Noether’s discovery of the deep connection between symmetries and conservation laws, arXiv:physics/9807044 (1998).
  * (18) J. C. Baez and B. Fong, A Noether theorem for Markov processes, J. Math. Phys. 54, 013301 (2013).
  * (19) I. Marvian and R. W. Spekkens, Extending Noether’s theorem by quantifying the asymmetry of quantum states, Nat. Commun. 5, 3821 (2014).
  * (20) S. Sasa and Y. Yokokura, Thermodynamic entropy as a Noether invariant, Phys. Rev. Lett. 116, 140601 (2016).
  * (21) S. Sasa, S. Sugiura, and Y. Yokokura, Thermodynamical path integral and emergent symmetry, Phys. Rev. E 99, 022109 (2019).
  * (22) M. Revzen, Functional integrals in statistical physics, Am. J. Phys. 38, 611 (1970).
  * (23) J. C. Baez, Getting to the bottom of Noether’s theorem, arXiv:2006.14741.
  * (24) A. Bravetti, M. A. Garcia-Ariza, and D. Tapias, Thermodynamic entropy as a Noether invariant from contact geometry, Entropy 25, 1082 (2023).
  * (25) S. Hermann and M. Schmidt, Noether’s theorem in statistical mechanics, Commun. Phys. 4, 176 (2021).
  * (26) S. Hermann and M. Schmidt, Why Noether’s theorem applies to statistical mechanics, J. Phys.: Condens. Matter 34, 213001 (2022) (Topical Review).
  * (27) S. Hermann and M. Schmidt, Variance of fluctuations from Noether invariance, Commun. Phys. 5, 276 (2022).
  * (28) S. Hermann and M. Schmidt, Force balance in thermal quantum many-body systems from Noether’s theorem, J. Phys. A: Math. Theor. 55, 464003 (2022). (Special Issue: Claritons and the Asymptotics of ideas: the Physics of Michael Berry).
  * (29) S. M. Tschopp, F. Sammüller, S. Hermann, M. Schmidt, and J. M. Brader, Force density functional theory in- and out-of-equilibrium, Phys. Rev. E 106, 014115 (2022).
  * (30) F. Sammüller, S. Hermann, D. de las Heras, and M. Schmidt, Noether-constrained correlations in equilibrium liquids, Phys. Rev. Lett. 130, 268203 (2023).
  * (31) S. Robitschko, F. Sammüller, M. Schmidt, and S. Hermann, Hyperforce balance from thermal Noether invariance of any observable, Commun. Phys. (resubmitted) arxiv:2308.12098.
  * (32) J. A. Barker and D. Henderson, What is “liquid”? Understanding the states of matter, Rev. Mod. Phys. 48, 587 (1976).
  * (33) H. Goldstein, C. Poole, and J. Safko, Classical Mechanics (Addison-Wesley, New York, 2002).
  * (34) V. Molinero and E. B. Moore, Water modeled as an intermediate element between carbon and silicon, J. Phys. Chem. B 113, 4008 (2009).
  * (35) M. K. Coe, R. Evans, and N. B. Wilding, The coexistence curve and surface tension of a monatomic water model, J. Chem. Phys. 156, 154505 (2022).
  * (36) S. Saw, N. L. Ellegaard, W. Kob, and S. Sastry, Structural relaxation of a gel modeled by three body interactions, Phys. Rev. Lett. 103, 248305 (2009).
  * (37) S. Saw, N. L. Ellegaard, W. Kob, and S. Sastry, Computer simulation study of the phase behavior and structural relaxation in a gel-former modeled by three-body interactions, J. Chem. Phys. 134, 164506 (2011).
  * (38) F. Sammüller, D. de las Heras, and M. Schmidt, Inhomogeneous steady shear dynamics of a three-body colloidal gel former, J. Chem. Phys. 158, 054908 (2023). (Special Topic on Colloidal Gels).
  * (39) R. B. Jadrich, D. J. Milliron, T. M. Truskett, Colloidal gels, J. Chem. Phys. 159, 090401 (2023).
  * (40) S. H. L. Klapp, Topical Review: Dipolar fluids under external perturbations, J. Phys.: Condens. Matter 17, R525 (2005).
  * (41) M. J. Stevens and G. S. Grest, Structure of soft-sphere dipolar fluids, Phys. Rev. E 51, 5962 (1995).
  * (42) P. I. C. Teixeira, J. M. Tavares, and M. M. Telo da Gama, Review Article: The effect of dipolar forces on the structure and thermodynamics of classical fluids, J. Phys.: Condens. Matter 12, R411 (2000).
  * (43) J. M. Tavares, J. J. Weis, and M. M. Telo da Gama, Strongly dipolar fluids at low densities compared to living polymers, Phys. Rev. E 59, 4388 (1999).
  * (44) M. P. Allen, Topical Review: Molecular simulation of liquid crystals Mol. Phys. 117, 2391 (2019).
  * (45) F. Sammüller and M. Schmidt, Adaptive Brownian dynamics, J. Chem. Phys. 155, 134107 (2021).
  * (46) M. Schmidt, Power functional theory for many-body dynamics, Rev. Mod. Phys. 94, 015007 (2022).
  * (47) R. Evans, The nature of the liquid-vapour interface and other topics in the statistical mechanics of non-uniform, classical fluids, Adv. Phys. 28, 143 (1979).
  * (48) F. Sammüller, S. Hermann, and M. Schmidt, Comparative study of force-based classical density functional theory, Phys. Rev. E 107, 034109 (2023).
  * (49) J. Yvon, La théorie statistique des fluides et l’équation d’état (in French), Actualités Scientifiques et Industrielles, (Hermann & Cie., Paris, 1935).
  * (50) M. Born and H. S. Green, A general kinetic theory of liquids I. The molecular distribution functions, Proc. R. Soc. London, Ser. A 188, 10 (1946).
  * (51) Y. Rosenfeld, Free-energy model for the inhomogeneous hard-sphere fluid mixture and density-functional theory of freezing, Phys. Rev. Lett. 63, 980 (1989).
  * (52) R. Roth, Fundamental measure theory for hard-sphere mixtures: a review, J. Phys.: Condens. Matter 22, 063102 (2010).
  * (53) T. Eckert, N. C. X. Stuhlmüller, F. Sammüller, and M. Schmidt, Local measures of fluctuations in inhomogeneous liquids: Statistical mechanics and illustrative applications, J. Phys.: Condens. Matter 35, 425102 (2023).
  * (54) F. H. Stillinger and T. A. Weber, Computer simulation of local order in condensed phases of silicon, Phys. Rev. B 31, 5262 (1985).
  * (55) J. G. Gay and B. J. Berne, Modification of the overlap potential to mimic a linear site-site potential J. Chem. Phys. 74, 3316 (1981).
  * (56) J. T. Brown, M. P. Allen, E. Martín del Río, and E. de Miguel, Effects of elongation on the phase behavior of the Gay-Berne fluid, Phys. Rev. E 57, 6685 (1998).
  * (57) S. Hamaguchi, R. T. Farouki, and D. H. E. Dubin, Triple point of Yukawa systems, Phys. Rev. E 56, 4671 (1997).
  * (58) M. Dijkstra and R. van Roij, Vapour-liquid coexistence for purely repulsive point-Yukawa fluids J. Phys.: Condens. Matter 10, 1219 (1998).
  * (59) M. Baus, Broken symmetry and invariance properties of classical fluids, Mol. Phys. 51, 211 (1984).
  * (60) R. Evans and A. O. Parry, Liquids at interfaces: what can a theorist contribute? J. Phys.: Condens. Matter 2, SA15 (1990).
  * (61) J. R. Henderson, Statistical mechanical sum rules, Chap. 2 in Fundamentals of Inhomogeneous Fluids, edited by D. Henderson (Dekker, New York, 1992).
  * (62) D. G. Triezenberg and R. Zwanzig, Fluctuation theory of surface tension, Phys. Rev. Lett. 28, 1183 (1972).
  * (63) L. V. Mikhheev and J. D. Weeks, Sum rules for interface Hamiltonians, Physica A 177, 495 (1991).
  * (64) A. Squarcini, J. M. Romero-Enrique, and A. O. Parry, Casimir contribution to the interfacial Hamiltonian for 3D wetting, Phys. Rev. Lett. 128, 195701 (2022).
  * (65) Y. Takahashi, On the generalized Ward identity, Il Nuovo Cimento 6, 371 (1957).
  * (66) J. C. Ward, An identity in quantum electrodynamics, Phys. Rev. 78, 182 (1950).
  * (67) J. Dong, F. Turci, R. L. Jack, Ma. A. Faers, and C. P. Royall, Direct imaging of contacts and forces in colloidal gels, J. Chem. Phys. 156, 214907 (2022).
  * (68) D. Borgis, R. Assaraf, B. Rotenberg, and R. Vuilleumier, Computation of pair distribution functions and three-dimensional densities with a reduced variance principle, Mol. Phys. 111, 3486 (2013).
  * (69) B. Rotenberg, Use the force! Reduced variance estimators for densities, radial distribution functions, and local mobilities in molecular simulations, J. Chem. Phys. 153, 150902 (2020).
  * (70) D. de las Heras and M. Schmidt, Better than counting: Density profiles from force sampling, Phys. Rev. Lett. 120, 218001 (2018).
  * (71) A. Purohit, A. J. Schultz, and D. A. Kofke, Force-sampling methods for density distributions as instances of mapped averaging, Mol. Phys. 117, 2822 (2019).
  * (72) S. W. Coles, D. Borgis, R. Vuilleumier, and B. Rotenberg, Computing three-dimensional densities from force densities improves statistical efficiency, J. Chem. Phys. 151, 064124 (2019).
  * (73) S. W. Coles, E. Mangaud, D. Frenkel, and B. Rotenberg, Reduced variance analysis of molecular dynamics simulations by linear combination of estimators, J. Chem. Phys. 154, 191101 (2021).
  * (74) S. W. Coles, B. J. Morgan, and B. Rotenberg, RevelsMD: Reduced variance estimators of the local structure in Molecular Dynamics, arXiv:2310.06149.
  * (75) D. de las Heras, T. Zimmermann, F. Sammüller, S. Hermann, and M. Schmidt, Perspective: How to overcome dynamical density functional theory, J. Phys.: Condens. Matter 35, 271501 (2023). (Invited Perspective).
  * (76) F. Sammüller, S. Hermann, D. de las Heras, and M. Schmidt, Neural functional theory for inhomogeneous fluids: Fundamentals and applications, Proc. Nat. Acad. Sci. 120, e2312484120 (2023).
  * (77) F. Sammüller, S. Hermann, and M. Schmidt, Why neural functionals suit statistical mechanics, arXiv:2312.04681.
  * (78) F. Sammüller, Neural functional theory for inhomogeneous fluids – Tutorial; for online access see the URL: https://github.com/sfalmo/NeuralDFT-Tutorial .
  * (79) G. E. Karniadakis, I. G. Kevrekidis, L. Lu, P. Perdikaris, S. Wang, and L. Yang, Physics-informed machine learning, Nat. Rev. Phys. 3, 422 (2021).
  * (80) S. Ciarella, M. Chiappini, E. Boattini, M. Dijkstra, and L. M. C. Janssen, Dynamics of supercooled liquids from static averaged quantities using machine learning, Mach. Learn.: Sci. Technol. 4, 025010 (2023).
  * (81) G. Janzen, C. Smit, S. Visbeek, V. E. Debets, C. Luo, C. Storm, S. Ciarella, and L. M. C. Janssen, Classifying the age of a glass based on structural properties: A machine learning approach, arXiv:2303.00636.
  * (82) G. Jung, G. Biroli, and L. Berthier, Predicting dynamic heterogeneity in glass-forming liquids by physics-inspired machine learning, Phys. Rev. Lett. 130, 238202 (2023).
  * (83) F. A. Rodrigues, Machine learning in physics: a short guide, EPL 144, 22001 (2023).
  * (84) J. Wu and M. Gu, Perfecting liquid-state theories with machine intelligence, J. Phys. Chem. Lett. 14, 10545 (2023).
  * (85) R. Evans and M. C. Stewart, The local compressibility of liquids near non-adsorbing substrates: a useful measure of solvophobicity and hydrophobicity?, J. Phys.: Condens. Matter 27, 194111 (2015).
  * (86) R. Evans, M. C. Stewart, and N. B. Wilding, A unified description of hydrophilic and superhydrophobic surfaces in terms of the wetting and drying transitions of liquids, Proc. Nat. Acad. Sci. 116, 23901 (2019).
  * (87) M. K. Coe, R. Evans, and N. B. Wilding, Density depletion and enhanced fluctuations in water near hydrophobic solutes: identifying the underlying physics, Phys. Rev. Lett. 128, 045501 (2022).
  * (88) T. Eckert, N. C. X. Stuhlmüller, F. Sammüller, and M. Schmidt, Fluctuation profiles in inhomogeneous fluids, Phys. Rev. Lett. 125, 268004 (2020).