Noether invariance theory for the equilibrium force structure of soft matter

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 (Sammüller et al 2023 Phys. Rev. Lett. 130 268203). 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 is thereby balanced by two further and different two-body correlation functions, which respectively introduce thermally averaged force correlations and force gradients in a systematic and microscopically sharp way into the framework. Separate exact self and distinct sum rules express 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 [1,2].The information contained in g(r) is analogously contained in 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 [3][4][5].
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 [6][7][8][9][10].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 [1].
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 [11].Such measures were argued to be relevant for the dynamics in dense colloidal liquids [12].The packing efficiency of binary hard sphere systems was related to fluid structure at intermediate ranges [13].Emergent structural correlations in dense liquids were characterized up to the six-body structure factor [14].Using machine learning the averaged structural features centered around nearby particles were used to predict dynamics in supercooled liquids [15].
Noether's theorem of invariant variations [16,17] was put in a statistical physics setting in a number of different ways [18][19][20][21][22][23][24].The recent thermal Noether invariance theory [25][26][27][28][29][30][31] 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. [31] carries the thermal Noether concept further in that it allows to address the behaviour of arbitrary observables in equilibribum.
The classical question "What is liquid?"[32] 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 [30] give fresh insight into the correlation structure.These correlation functions are generated from canonical "shifting" transformation on phase space [33] 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 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, g ∇f (r) and g ff (r), as are defined in detail in Sec.VI.These tensorial functions possess radial (∥) and transversal (⊥) components, both of which are plotted.Briefly, the two-body force-gradient correlator g ∇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 g ff (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 ρ 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. [30].
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): Here ∇ indicates the spatial gradient, such that ∇∇g(r) is the Hessian of the pair correlation function.Furthermore the force-gradient pair correlator g ∇f (r) and the force-force pair correlator g ff (r) are both 3 × 3-matrices.These correlation functions also address particle pairs, but they incorporate additional information over the mere counting in g(r).In particular g ff (r) correlates the interparticle forces that act on the two particles in the pair with each other.The force-gradient correlation function g ∇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. [30] 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 [34,35] and a colloidal gel forming model proposed by Kob and coworkers [36][37][38].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 g ff (r) and in g ∇f (r) across different models and phases are signifi-cantly 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 [39] and in dipolar fluids [40][41][42][43][44].
To illustrate this behaviour we reproduce in Figure 1 results of Ref. [30] for the Lennard-Jones system in all of its three stable thermodynamic phases: crystal, liquid, and gas.The pair potential is thereby truncated and shifted at a cutoff distance of r c = 2.5σ, where σ denotes the particle size.The data is obtained from equilibribum sampling via adaptive Brownian dynamics [45], 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(−βϕ(r)), where β denotes inverse temperature and ϕ(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.Due to the rotational symmetry of fluid states and the absence of any chirality, these are the only two relevant tensor components to consider.In principle the structure of the crystal is much richer, but we project onto the same two tensor components and hence average over the further spatial and rotational inhomogeneities of the crystal.The mean force gradient g ∇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 g ff (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 g ∇f (r) remains short-ranged and it strictly vanishes beyond the cutoff for short-ranged or truncated pair potentials.The accordance of the dashed and the solid lines in the third row of Fig. 1 confirms numerically that the 3g-rule sum rule (1) is satisfied.
With the central mathematical ideas only being sketched in Ref. [30], 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 [1].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 presented in Sec.III 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 r 1 , . . ., r N ≡ r N and linear momenta p 1 , . . ., p N ≡ p N .The Hamiltonian H is taken to contain kinetic energy, as well as interparticle and external energy contributions, where the summation index i = 1, . . ., N ranges over all particles, m denotes the particle mass, u(r N ) is the interparticle interaction potential, and the one-body external potential V ext (r) is given as a function of a single (generic) position variable r.The canonical transformation [33] of the phase space coordinates that was put forward in Refs.[28][29][30] maps original onto new coordinates as well as original onto new momenta according to the following map: where the tilde indicates the new variables.The threedimensional vector field ϵ(r) parameterizes the transform, 1 denotes the 3 × 3-unit matrix, ∇ i indicates the derivative with respect to 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 [28][29][30]33] and it constitutes a canonical transformation in the sense of classical mechanics [33].This property can e.g.be seen explicitly by considering a generator of the form G = N i=1 pi • (r i + ϵ(r i )), where as before the tilde indicates the new phase space variables.Applying the generic transformation relations ri = ∂G/∂ pi and p i = ∂G/∂r i [33] and solving for the new variables then yields the transformation equations ( 3) and (4).The form of the shifting vector field ϵ(r) must thereby such that the transformation is bijective [28].
We Taylor expand the transformation in the shifting field ϵ(r).The coordinate transformation (3) is already linear in ϵ(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: (5) where the square on the right hand side denotes matrix (self-)multiplication, i.e., [ . Verifying Eq. ( 5) to second order in the shifting gradient is straightforward via matrix multiplying its both sides by 1 + ∇ i ϵ(r i ), which per definition of the inverse gives 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 → H[ϵ].We first consider the term linear in ϵ(r) [28][29][30].One sees that the position-resolved one-body force density operator F(r) follows from functional differentiation according to: where δ/δϵ(r) indicates the functional derivative with respect to the shifting field ϵ(r); see e.g.Ref. [46] for an overview of functional differentiation techniques.After the derivative is taken, ϵ(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 F(r) mirrors that of the Hamiltonian (2) in that it contains kinetic, interparticle, and external contributions: where ∇ denotes the derivative with respect to r.The one-body operators for the kinetic stress τ (r), the interparticle force density Fint (r) [46], and the microscopic density distribution ρ(r) [1,47] are given respectively by where δ(•) indicates the (three-dimensional) Dirac distribution and p i 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 µ.The grand potential Ω and the grand partition sum Ξ are given respectively by where k B is the Boltzmann constant, β = 1/(k B T ) denotes inverse temperature, and the classical "trace" operation in the grand ensemble is given by Tr The corresponding grand probability distribution (Gibbs measure) is Ψ = e −β(H−µN ) /Ξ. Thermal averages are defined via ⟨•⟩ = Tr Ψ•, such that the density profile is the average of the one-body density operator and written as ρ(r) = ⟨ρ(r)⟩.Correspondingly the average kinetic stress is τ (r) = ⟨ τ (r)⟩ and the mean interparticle force density is At face value the grand partition sum ( 12) carries an apparent functional dependence on the shifting field ϵ(r), via the occurrence of ϵ(r) in the transformed Hamiltonian H[ϵ] [28][29][30].We hence also have a functional dependence of the grand partition sum, Ξ[ϵ], and consequently of the grand potential, Ω[ϵ], 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 The colon denotes a double tensor contraction and Ω, with no argument, indicates the original grand potential (11) with no transformation applied, Ω = Ω[0], where the functional argument 0 indicates the special case of vanishing shifting field, ϵ(r) = 0. Despite the apparent dependence on ϵ(r), Noether invariance [25,26] ascertains that the grand potential remains unchanged upon applying the transformation, and hence which holds true for any (permissible) form of ϵ(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 Explicitly carrying out the functional derivative on the left hand of Eq. ( 15) gives [28][29][30]: where we recall the total force density operator F(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 −F(r) = −⟨ F(r)⟩.Due to the invariance that is expressed via Eq.( 15) we can conclude that the mean one-body force density vanishes in equilibribum, The explicit form (7) of F(r) allows to decompose the left hand side of Eq. ( 17) into a sum of ideal gas, interparticle, and external terms according to: 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 [1], in which F int (r) is further expressed in integral form.The present functional route [28][29][30] 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 [29,48] 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 [1,49,50] 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.Forcedensity 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 [51,52], see Refs.[29,48].

III. LOCAL SECOND-ORDER INVARIANCE
The previous Sec.II summarizes the theory developed in Refs.[28,29], which we require for considering the second order.We recall that for global shifting with spatially uniform displacement, ϵ(r) = ϵ 0 = const, addressing the second order in ϵ 0 yields a sum rule which relates the global force variance with the mean global potential energy curvature [27].
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: where ∇ ′ denotes the derivative with respect to r ′ and H 2 (r, r ′ ) = cov(ρ(r), ρ(r ′ )) is the standard two-body correlation function of density fluctuations around the local density profile [1,46].Thereby the covariance of two observables, as represented by phase space functions Â and B, is defined in the standard way as cov( Â, B) = ⟨ Â B⟩ − ⟨ Â⟩⟨ B⟩.Upon carrying out the position integrals the left hand side of Eq. ( 19) can be alternatively written as ⟨ F0 ext F0 ext ⟩, where the global external force density operator is defined as F0 ext = − i ∇ i V ext (r i ) and the covariance reduces to the mean of the product, because of the vanishing mean external force in equilibribum, ⟨ F0 ext ⟩ = 0 [25].Following the strategy of Ref. [30] 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 ϵ(r).Addressing the second order term, this implies that its prefactor in Eq. ( 13) needs to vanish identically: Making the functional derivative that appears on the left hand side more explicit yields The functional derivative of the Hamiltonian, δH[ϵ]/δϵ(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: The correlator on the left hand side of Eq. ( 22) is analogous to cov( F(r), F(r ′ )) = ⟨ F(r) F(r ′ )⟩, because the local mean force vanishes in equilibrium, ⟨ F(r)⟩ = 0 [28,29]; 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 r and r ′ , could have certainly not been guessed easily a priori (we turn to partial integration methods below and present details thereof in Appendix D).
We re-write Eq. ( 22) via multiplying by β and then treating the kinetic contributions separately (see Sec. V).For brevity of notation it is useful to combine the interparticle and external forces into the following dedicated potential force density operator: where the subscript U refers to the total potential energy and we recall the definition (9) of the interparticle force operator Fint (r).The total potential energy operator is given as the phase space function Using the coordinate transformation (3) the potential force density is obtained via functional differentiation according to −δH U /δϵ(r)| ϵ=0 = FU (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 ij(̸ =) ≡ N i=1 N j=1,j̸ =i .As we will show in detail below in Sec.V, we can then obtain from Eq. ( 22) the following pair of exact distinct and self identities: 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 [1]: The self average on the left hand side of Eq. ( 25) only involves a single delta distribution in position, i.e., it is defined such that ⟨β FU (r)β FU (r) of the integral over r i and then drop the common factor δ(r − r ′ ) on both sides of the self equation in order to arrive at the form (25).
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. ⟨ρ(r) F(r ′ )⟩.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 −⟨ρ(r) F(r ′ )⟩β∇V ext (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), ⟨ F(r)⟩ = 0.This identity holds regardless of the form of the external potential V ext (r), and we can hence functionally differentiate both sides of the equation with respect to V ext (r).Clearly the right hand side will remain zero.Differentiating the left hand side yields 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. [53] have recently pointed out that for an arbitrary phase space function Â, which is taken to be independent of V ext (r), one has δ⟨ Â⟩/δV ext (r) = −βcov( Â, ρ(r)), which upon choosing Â = F(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 F(r) on V ext (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: The first term on the right hand side gives βρ(r ′ )F(r) upon observing that −Ξ −2 δΞ/δV ext (r ′ ) = Ξ −1 βρ(r ′ ), where the density profile can be taken out of the trace in Eq. ( 27).The second term follows straightforwardly as −β⟨ρ(r ′ )F(r)⟩ + ⟨δ F(r)/δV ext (r ′ )⟩ when taking account of the fact that δH/δV ext (r ′ ) = ρ(r ′ ).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, cov(β F(r), ρ(r ′ )) = ⟨β F(r)ρ(r ′ )⟩, again because the average force density vanishes in equilibrium, ⟨ F(r)⟩ = F(r) = 0, cf.Eq. ( 17).The second term in Eq. ( 27) can be re-written as follows: In the last step above we have introduced the self twobody density, defined as In typical applications one would re-write the right hand side of Eq. ( 34) as ρ(r)δ(r − r ′ ), but due to the present gradient structure some care is required and we thus keep ρ self 2 (r, r ′ ) 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: 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 r and r ′ the singular right hand side vanishes and Eq. ( 35) reduces trivially to We can also conclude, as the distinct contribution vansishes, that at all positions r and r ′ , even when r = r ′ , we have We recall that the present proof of the density-force sum rule (35) is based on functionally differentiating F(r) = 0 with respect to V ext (r ′ ).An alternative derivation of Eq. ( 35) can be based on the Noether invariance of the density profile.This thermal symmetry is expressed as δρ(r ′ )/δϵ(r)| ϵ=0 = 0.As before we denote the transformed positions by ri and consider where the first term on the right hand side results from functionally differentiating the probability distribution and it is already evaluated at ϵ(r) = 0.The second term can be re-written as where we have used the definition (34) of the self twobody density distribution ρ self 2 (r, r ′ ).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.Starting from the free energy as an overarching object and building a mixed second derivative, where the order of differentation is interchanged yields: where we recall that δΩ/δV ext (r) = ρ(r) and δΩ[ϵ]/δϵ(r)| ϵ=0 = −F(r) according to Eq. ( 16).
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 F(r) = ∇• τ (r)+ FU (r), where we recall the definition (8) of the kinetic stress density operator τ (r) and the definition (23) of the potential force density operator FU (r).After re-arranging Eq. ( 35) we find 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 ⟨ρ(r)ρ(r ′ )⟩ = ρ 2 (r, r ′ ) + ρ self 2 (r, r ′ ).Building furthermore the gradient with respect to the position r ′ and transposing [note that the transpose (∇ ′ ∇) T = ∇∇ ′ ] yields straightforwardly the following tensorial identities: Equation ( 44) is thereby obtained from re-transposing (43) and interchanging the positions arguments r and r ′ upon observing the exchange symmetry ρ self 2 (r, r ′ ) = ρ self 2 (r ′ , r); we recall the definition (34) of the self twobody density and clearly also ρ 2 (r, r ′ ) = ρ 2 (r ′ , 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 FU (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 F(r) = ∇• τ (r)+ FU (r).We start with the force autocor-relator and obtain by multiplying out and re-ordering: We start with the purely kinetic two-body contribution, i.e. the first term in Eq. ( 45) and obtain the following result: where we have again used the definition (34) of the self two-body density and we recall that 1 indicates the 3 × 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 Thereby the momentum averages on the left hand side need to be carried out explicitly, upon using the Maxwell distribution, in order to simplify the divergence of the kinetic stress operator (left hand sides) to yield the gradient of the density operator (right hand sides).Then re-introducing the momentum average on the right hand side, as is implied by the ⟨•⟩ notation, creates no harm due to the independence of the operators on the right hand sides from the momentum degrees of freedom and the correct normalization of the average.Collecting all terms and using Eqs.( 43) and ( 44) allows then to re-write Eq. ( 45) upon multiplying by β 2 as 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 ϵ(r), into kinetic and potential energy contributions, . 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: 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: Re-arranging and simplifying yields: where we could finally replace the self two-body density as follows.In the external curvature term we have ρ self 2 (r, r ′ ) = δ(r − r ′ )ρ(r) and in the ideal contribution we could simplify (∇ + ∇ ′ )(∇ + ∇ ′ )ρ self 2 (r, r ′ ) = δ(r−r ′ )∇∇ρ(r), which can be verified explicitly by calculating derivatives in (∇ + ∇ ′ )(∇ + ∇ ′ )[δ(r − r ′ )ρ(r)] using the product rule and exploiting ∇ ′ δ(r−r ′ ) = −∇δ(r−r ′ ).
One can split Eq. ( 53) into regular and singular contributions, where the latter are identified as having a common factor δ(r−r ′ ).This splitting discriminates between self and distinct contributions on the basis of the standard criterion for pairs of particle indices which are equal i = j (self) and different i ̸ = j (distinct).Hence splitting Eq. ( 53) leads respectively to the distinct sum rule (24) and, upon leaving away the common factor δ(r − r ′ ) from the self terms, to the self sum rule (25).
This completes our proof of the locally resolved twobody 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 selfstructuring 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.

VI. NOETHER STRUCTURE IN BULK
We specialize to homogeneous bulk liquid states, where ρ(r) = ρ b = const and V ext (r) = 0.The potential forces then only arise from interparticle contributions and hence FU (r) = Fint (r) with the interparticle force density operator Fint (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 = |r − r ′ | denotes the separation distance between the two particles.The relationship to the two-body density is simply via normalization with the squared bulk density ρ b : Furthermore we follow Ref. [30] in introducing both the force-force pair correlation function g ff (r) and the force gradient correlator g ∇f (r).These are rank-two tensorial correlation functions, as represented by 3 × 3-matrices in case of a three-dimensional system, and they are defined via The force gradient correlator g ∇f (r) is also the negative mean potential curvature, as is apparent from the right hand side of Eq. ( 56) via ∇ i ∇ j u(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.Recalling the interparticle force on particle i as f i (r N ) = −∇ i u(r N ) and that on particle j as f j (r N ) = −∇ j u(r N ) we can re-write the potential energy curvature as  [34,35] (third column), and the three-body interacting gel [36][37][38] (fourth column); our model parameter choices are inline with the values given in these references.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 (∥) and transversal (⊥) components of the force-gradient correlator g ∇f (r) (middle row) and of the force-force correlator g ff (r) (bottom row).The vertical gray lines indicate the position of the first maximum of g(r) as a guide to the eye.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 agreement with the respective solid lines, as obtained from numerically differentiating the interparticle forces, demonstrates that the force gradient correlation function follows directly from the product of g(r) and the scaled derivatives βϕ ′′ (r) and βϕ ′ (r)/r of the pair potential; see Eqs. (64) and (65).The dotted lines in the bottom row indicate the results according to the Noether force sum rules (62) and (63).The agreement with the respective solid lines demonstrates the validity of these sum rules.Taken from Ref. [30].
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: where in the last term one can alternatively replace T , as laid out above.
In the fully position-resolved case, which is appropriate for the investigation of spatially inhomogeneous sit-r − r ′ and the second and third components being equal and corresponding to two mutually perpendicular directions which are also perpendicular to r − r ′ .
More explicitly, in the considered radial symmetry the second derivative operator in Eq. ( 61) reduces to ∇∇ = e ∥ e ∥ ∂ 2 /∂r 2 + (1 − e ∥ e ∥ )r −1 ∂/∂r, where the direction that connects the two points is e ∥ = (r − r ′ )/|r − r ′ |.This unit vector is complemented by two orthogonal directions e ⊥ , e ′ ⊥ such that all three e ∥ , e ⊥ , e ′ ⊥ are mutually orthogonal to each other.The perpendicular components then generate the g ′ (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 ϕ(r), the force gradient correlator (56) reduces to g ∇f (r) = βg(r)∇∇ϕ(r) such that the two non-trivial components (∥ and ⊥) can be written according to This simplification is due to the reduction of the mixed derivative , for i ̸ = 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.The fact that the separation distance r between the two particles is fixed via the delta distributions in position allows to take the Hessian ∇∇ϕ(r) outside of the average, with the remaining average then generating the bare pair correlation function g(r); we recall its definition via Eqs.( 26) and (54).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: 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(−βϕ(r)), as laid out in Ref. [30].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: Furthermore g ff ⊥ (r) = 0 due to the absence of a third particle at ρ b → 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 [30], 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 r−r ′ over all orientations and translations.We recall the display of simulation results in Fig. 1 and in particular the validity of the sum rules ( 62) and ( 63), indicated by the agreement of the solid and dashed lines in the bottom row of Fig. 1.The parallel component of the force gradient correlator g ∇f ∥ (r) has a strong positive peak at a distance near the Lennard-Jones diameter.The perpendicular component g ∇f ⊥ (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), g ff (r), and g ∇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 [34] of the Stillinger-Weber potential [54], and for the three-body interacting gel [36][37][38] in Fig. 2; the model parameter values are inline with the choices in these references.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 ∇f ⊥ (r) is absent in the Weeks-Chandler-Andersen liquid.This behaviour is very notable due to the common view of interparticle attraction having no significant effect on the microscopic structure of liquids.While our present model comparison explicitly confirms this expectation based solely on the behaviour of the pair correlation function g(r), a much more complete picture emerges from the full force-based pair correlation structure.As it turns out, g ∇f (r) is a suitable means to clearly point to the presence of interparticle attraction.We recall the relationship (65) with the pair potential.
Going to the monatomic Molinero-Moore [34,35] 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 [36][37][38] is, similar to the monatomic water model, a mere reparametrization of the original Stillinger-Weber [54] model.Yet both the forcegradient and the force-force correlation structure (fourth column) changes dramatically when going from the liquid to the gel.The perpendicular mean force-gradient g ∇f ⊥ (r) develops a strong positive signal immediately beyond the first peak of g(r), which is indicated by the vertical gray line as a reference.
Even more strikingly the transverse force-force correlation function g ff ⊥ (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. [38]).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 [44], and the Gay-Berne model [44,55,56] 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.Choices for parameter values are given in the caption of Fig. 3, following Refs.[41,56].In all cases considered the sum rules were found to be satisfied to high numerical accuracy.
The results for the point-Yukawa fluid [57,58] (first column) serve again to illustrate generic liquid-like structuring.No positive signal in g ∇f ⊥ (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 ∇f ⊥ (r) has a strong positive peak directly beyond the first peak of g(r) and the prominent peak of g ff ⊥ (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 ∇f ⊥ (r) has lost its negative peak at distances smaller than the location of the first peak of g(r) (vertical gray line).Instead g ∇f ⊥ (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 [30].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 [25][26][27] and subsequently generalized to local shifting [28,29], the approach yields new insight into the correlated pair structure of soft matter [30].
We have here provided a detailed account of the theory, giving derivations, including the treatment of all selfcontributions 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 [34,35], and for the many-body colloidal gel former [36][37][38].
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 threebody 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 [59][60][61][62] and they share formal similarities with sum rules for interface Hamiltonians [63], see Ref. [64] for recent work, and with the Takahashi-Ward identities [65,66] of quantum field theory.As locally resolved interparticle force measurements have been demonstrated in colloidal systems [67], 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 [68][69][70][71][72][73][74] is another interesting point for future work, as can be relating to force-based density functional theory [29,48].The recent neural functional theory [75][76][77][78] 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 [79][80][81][82][83][84].
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 µ 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 [85][86][87].The Noether approach captures genuinely the mechanical fluctuations, as opposed to local chemical (particle number) [53,[85][86][87][88] and thermal (energy) fluctuations [53,87,88].We lastly point to the recent hyperforce generalization [31] that arises from thermal Noether invariance and that applies to arbitrary observables.
The self part of Eq. (A1) consists only of the cases i = j in the double sum.We can hence simplify Eq. (A1) according to: Switching to index notation and denoting Cartesian components by Greek indices, the αγ-component of this second-rank tensor is: 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 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. (A4) is normalized, dp i e −βp 2 i /(2m) /(2πm/β) 3/2 = 1.In the derivation of Eq. (A4) one requires both the second moment of each (αth) momentum component, dp i (p α i ) 2 e −βp 2 i /(2m) /(2πm/β) 3/2 = m/β, and the fourth moment, Recalling the derivative structure of Eq. (A3) we contract Eq. (A4) with ∇ ν ∇ ′ ξ , which yields Using this result we can simplify the expression (A3) as Hence in summary we obtain the following self contribution: where we recall the definition (34) of the two-body self density: ρ self 2 (r, r ′ ) = ⟨ i δ(r − r i )δ(r ′ − r i )⟩.We next consider the distinct part of Eq. (A1), which consists of the cases i ̸ = j.In index notation its αγcomponent is p α i p γ j p ν i p ξ j m 2 δ(r − r i )δ(r ′ − r j ) .(A9) As before we carry out the momentum average which requires the following momentum integrals of particles i and j: Using this result we can re-write Eq. (A9) upon simplifying according to νξ ∇ ν ∇ ′ ξ δ αν δ γξ = ∇ α ∇ ′ γ to hence We can now express Eq.(D5) 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): (D7) The identity (D6) holds for any r, r ′ and Eq.(D7) 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 (D3) for the distinct correlation identity, we consider the double gradient, but with respect to the same position r taken twice, i.e. the Hessian of the density profile: In this derivation we have exploited twice that ∇δ(r − r i ) = −∇ i δ(r − r i ) and then integrated by parts first according to Eq. (D1) and then according to Eq. (D2).
Identifying FU (r) in Eq. (D11) and re-arranging the different terms gives Eq. (D7).The present phase space integration route to the sum rules (D6) and (D7) [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.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 r i and r j be the positions of two particles, such that r ∥ = r j − r i and e ∥ = r ∥ /|r ∥ |.The calculation of g ff (r) proceeds straightforwardly by incorporating the forces f i and f j of particles i and j as given in Eq. ( 57).The radial and tangential components are respectively obtained by the following projections: where the factor 1/2 in Eq. (E2) accounts for the two equal tangential contributions of g ff (r) in a threedimensional bulk fluid (there is only a single radial component).
For the calculation of g ∇f (r), we employ numerical differentiation to evaluate the force gradients ∇ i 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 g ∇f (r) into its radial and tangential components is performed.Assuming f j ∦ e ∥ , a tangential unit vector e ⊥ = t/|t| ⊥ e ∥ can be constructed by choosing t = f j − (f j • e ∥ )e ∥ .The radial and tangential parts of g ∇f (r) then follow via where D i,e is a directional derivative operator regarding the shifting of particle i along the (generic) direction e.
The evaluation of Eqs.(E3) and (E4) is performed numerically in the simulations, i.e., particle i is shifted by a small amount (∼ 10 −5 σ) to calculate the arising finite difference in the force f j of particle j.
The above procedure applies to isotropic particles.For anisotropic particles a simple alternative to construct a suitable vector e ⊥ can be based on drawing a random unit vector e ran ∦ e ∥ .Then e ⊥ = t/|t| with t = e ran − (e ran • e ∥ )e ∥ , which can be used in Eq. (E4).

0 S 2 S
FIG.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, g ∇f (r) and g ff (r), as are defined in detail in Sec.VI.These tensorial functions possess radial (∥) and transversal (⊥) components, both of which are plotted.Briefly, the two-body force-gradient correlator g ∇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 g ff (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 ρ 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.[30].

0 S 2 S
FIG.2: Correlation functions analogous to Fig.resultsfor the Lennard-Jones liquid shown again as a reference (first column), the Weeks-Chandler-Andersen fluid (second column), monatomic water at (approximately) ambient conditions[34,35] (third column), and the three-body interacting gel[36][37][38] (fourth column); our model parameter choices are inline with the values given in these references.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 (∥) and transversal (⊥) components of the force-gradient correlator g ∇f (r) (middle row) and of the force-force correlator g ff (r) (bottom row).The vertical gray lines indicate the position of the first maximum of g(r) as a guide to the eye.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 agreement with the respective solid lines, as obtained from numerically differentiating the interparticle forces, demonstrates that the force gradient correlation function follows directly from the product of g(r) and the scaled derivatives βϕ ′′ (r) and βϕ ′ (r)/r of the pair potential; see Eqs.(64) and(65).The dotted lines in the bottom row indicate the results according to the Noether force sum rules(62) and(63).The agreement with the respective solid lines demonstrates the validity of these sum rules.Taken from Ref.[30].
Appendix E: Measuring force-force and force-gradient correlations