Hyperforce Balance via Thermal Noether Invariance of Any Observable

We demonstrate that thermally averaged classical phase space functions, including local and global forces and energies, are associated with exact hyperforce sum rules that follow from translational Noether invariance. Both global and locally resolved identities hold and they relate the mean gradient of a phasespace function to its negative mean product with the total force. Similar to Hirschfelder’s hypervirial theorem [J. Chem. Phys. 33 , 1462 (1960)], hyperforce sum rules apply to arbitrary observables in equilibrium. Exact hierarchies of higher-order sum rules follow iteratively. As applications we investigate via computer simulations the emerging one-body force fluctuation profiles in confined liquids. These local correlators quantify spatially inhomogeneous self organization and they allow to carry out stringent tests for sufficient sampling in complex systems.


I. INTRODUCTION
The task of predicting thermal averages of phase space functions lies at the center of attention in Statistical Mechanics.Prominent examples include correlation functions and order parameters, but also global quantities such as internal and external energies, the entropy, and much more are considered [1,2].Significant progress has been reported for problem-specific order parameters that are tailored to capture intricate correlations effects.Recent examples that address the spatial ordering behaviour of dense liquids include beyond-two-body correlation functions, as advocated by Kob and coworkers [3,4] and by Janssen and her coworkers [5].
In constrast to such freedom of choice, the variables within classical density functional theory [2,6,7] seem to be a priori uniquely determined by the existence of a generating free energy functional and the associated structure of pairs of conjugate variables, which in particular are the external one-body potential energy V ext (r) and the density profile ρ(r).However, there are recent extensions to density functional theory to systematically include the local compressibility [8][9][10], which forms a well-accessible order parameter for local particle number fluctuations.Technically, the local compressibility constitutes either a parametric derivative of the equilibrium density profile with respect to the chemical potential or, analogously, the covariance of the local density and the global particle number.A generalization from such chemical particle number fluctuation to thermal fluctuations has been recently performed [11,12].Working in the grand ensemble, where the particle number fluctuates, is thereby crucial to not impose artificial constraints on the system.
Besides the standard thermodynamical thinking in terms of thermal and chemical equilibrium, there is much recent progress in the force point of view.Highly efficient force sampling techniques allow to obtain reliable results within many-body simulations that outperform more straightforward counting methods [13][14][15][16][17][18][19].Forces are also at the core of power functional theory [20] as a systematic approach to formulate coupled many-body dynamics on the one-body level of dynamical correlation functions.
To be specific, in thermal equilibrium of a spatially inhomogeneous system, the sum of all mean forces necessarily vanishes at each position r.This is expressed by the following exact sum rule: Here F int (r) is the localized force density that acts at position r due to the interparticle interactions with all surrounding particles, ∇ denotes the derivative with respect to r such that −∇V ext (r) is the external force field, k B indicates the Boltzmann constant, and T is absolute temperature.
The sum of the interparticle and external force densities on the left hand side of Eq. (1) balances the thermal diffusive contribution on its right hand side.This is a classical result due to Yvon, Born and Green (YBG) [1], where for particles that mutually interact only via a pair potential, the interparticle force density F int (r) is expressed as an integral over the two-body density multiplied by the pair force [1].Higher-order versions of Eq. ( 1) form a hierarchy.That the first YBG equation (1) has practical consequences for carrying out sampling tasks in simulations is only a quite recent insight.Loosely speaking, force sampling [13][14][15] amounts to obtaining simulation data for the left hand side of Eq. ( 1) and then in a post processing step dividing by k B T and building the inverse operation of the spatial derivative on the right hand side via suitable integration in position.This method yields results for the density profile ρ(r) which feature a significant reduction of statistical noise [13][14][15][16][17][18][19].
Exploiting Noether's Theorem [21,22] in statistical mechanics has been performed in a variety of ways [23][24][25][26][27][28][29][30].Considering the invariance of statistical mechanical functionals leads very naturally to the notion of statistically averaged forces when spatial displacement is imposed on the system; mean torques emerge when invariance against rotations is addressed [31,32].In previous work, we have shown that the statistical Noether concept also applies quantum mechanically [33] and that it gives access to global force fluctutations [34].Generalizing to local invariance [33,35] that is resolved in spatial position facilitated fresh insights into the correlation structure of the liquid state [36].Considering the first order in the displacement field yields the thermal equilibrium force balance relationship according to the YBG equation (1) [33,35].At second order hitherto unkown two-body force-gradient and force-force correlators emerge and these, together with the standard pair correlation function, are constrained by exact Noether identities [36].
This situation of theory development leaves open the question of whether more general observables that serve as important order parameters and quantifiers of spatial structure will also be affected by the statistical Noether invariance, as one could glean from the generality of the thermal invariance concept.Here we demonstrate that any statistical observable Â is intrinsically associated with a corresponding hierarchy of exact identities that emerges from its statistical shifting invariance properties.We validate the corresponding exact local and global sum rules for a range of relevant observables via many-body simulations of a confined Lennard-Jones fluid.The results clarify a very intimate link of global and locally resolved correlators and they suggest a very general statistical mechanical structure.
Our framework can be viewed as a generalization of the YBG equation (1) to systematically include the dependence on a further given observable Â.The equilibrium force balance (1) itself is recovered for the trivial case Â = 1.Such generalization is not uncommon in Statistical Mechanics.The relationship of our theory and the YBG force balance equation is akin Hirschfelder's hypervirial theorem [37] as a generalization of the standard virial theorem [1] to also invoke an additional dependence on a given phase space function Â.Our theory can hence be viewed as a hyperforce balance relationship, and we derive global and local variants below, see Eqs. ( 4) and (10).We further show that the local version simplifies further to Eq. (11) and more explicitly to Eq. (12) in case of Â being independent of momenta; see Fig. 6 for an overview of results.As a specific example, our methodology not only allows to sample density gradients, as is possible in force sampling schemes [13][14][15][16][17][18][19], but also to sample force density gradients.The general method complements existing counting and force-sampling techniques and it gives much inspiration for rigorous statistical mechanical theories based on exact identities.As we lay out, the degree of numerical accuracy to which the Noether sum rules are satisfied can serve as an estimator for sufficient equilibration of slowly converging systems.
The paper is organized as follows.Starting from the description of the underlying Statistical Mechanics, we present the general Noether invariance theory in Sec.II for cases of global and locally resolved shifting invariance.The theory is applied to specific observables, tested in a confined fluid, and its implications for sampling methodology are discussed in Sec.III.We give conclusions and an outlook in Sec.IV.

A. Statistical Mechanics
We consider general thermal many-body systems of particles with identical mass m, coordinates r 1 , . . ., r N ≡ r N , and momenta p 1 , . . ., p N ≡ p N , where N denotes the number of particles.The Hamiltonian is of the standard form H = i p 2 i /(2m) + u(r N ) + i V ext (r i ), where the sums i = 1, . . ., N run over all particles, u(r N ) is the interparticle interaction potential, and V ext (r) is an external potential that depends on position r.Thermal equilibrium is characterized by a statistical equilibrium ensemble with grand canonical probability distribution Ψ eq = Ξ −1 e −β(H−µN ) , where β = 1/(k B T ) and µ indicates the chemical potential.The normalization factor of Ψ eq is the partition sum Ξ = Tr e −β(H−µN ) , where the classical trace is defined as Tr The thermal equilibrium average A of a given phase space function Â is then obtained as A = ⟨ Â⟩ ≡ Tr Ψ eq Â, where we have suppressed the dependence of Â on the phase space variables r N and p N in the notation, i.e., in full notation we have Â(r N , p N ) as well as potentially further parametric dependence such as on a generic position variable r.

B. Global shifting invariance
To develop the Noether invariance theory, we first consider a global coordinate displacement r i → r i + ϵ 0 ≡ ri , where the shifting vector ϵ 0 = const is independent of position and acts on all particles in the same way [31].The tilde indicates the new coordinates.Expressing H as well as the corresponding distribution function Ψ eq in the new coordinates makes averages become formally dependent on the shifting parameter, i.e., A(ϵ 0 ).However, the coordinate change can also be viewed as a mere re-parameterization of the phase space integral which induces no change to its value such that A(ϵ 0 ) = A, with the right hand side denoting ⟨ Â⟩ in the original representation.Partially differentiating both sides of the equation yields: where the second equality arises trivially from ∂A/∂ϵ 0 = 0, as there is no dependence on ϵ 0 .Carrying out the derivative is straightforward upon using the Boltzmann form of the equilibrium distribution function and noting that the partition sum is independent of ϵ 0 .Explicitly we have ∂Ψ eq (ϵ 0 )/∂ϵ 0 = −βΨ eq (ϵ 0 )∂H(ϵ 0 )/∂ϵ 0 .Using the product rule of differentiation and evaluating at vanishing displacement ϵ 0 = 0 then leads from Eq. (2) to Observing that here ∂/∂ϵ 0 | ϵ0=0 ≡ i ∇ i allows to make the derivatives more explicit with ∇ i indicating the partial derivative with respect to r i .In the first term in Eq. ( 3) we have ext , which defines the global external force operator F0 ext .Here the global interparticle force due to the mutual interactions between all particles in the system vanishes, F0 int = − i ∇ i u(r N ) ≡ 0, as is due to Newton's third law or, analogously, to the translational invariance of u(r N ) against global displacement [31].
Re-ordering the two terms in Eq. ( 3) gives the following global hyperforce identity that holds for any given observable Â(r N , p N ): Here Â = Â(r N , p N ) can feature additional parametric dependence, such as on a generic position argument r.
The sum rule (4) relates the correlation of Â with the external force operator (left hand side) to the mean negative global coordinate derivative of Â (right hand side); here we use the term correlation to imply the average of the product of two observables.As announced in the introduction, Eq. ( 4) is similar to Hirschfelder's hypervirial theorem [37] in the inclusion of the phase space function Â(r N , p N ) and our hyperforce terminology parallels his use of the term hypervirial.While Noether invariance enabled us to obtain the global hyperforce identity Eq. ( 4) constructively, one can verify its validity a posteriori by integration by parts in phase space on the right hand side.The derivatives ∇ i then act on the probability distribution and exploiting again the Boltzmann form leads to Â i ∇ i Ψ eq = −βΨ eq Â i ∇ i H = Ψ eq β F0 ext Â, which gives the left hand side of Eq. ( 4) upon building the trace.Alternatively one can start with the Yvon theorem [1], β⟨ Â∇ i H⟩ = ⟨∇ i Â⟩, which itself also follows from partial phase space integration [1].Summing the Yvon theorem over all particles i and noting that i ∇ i u(r N ) ≡ 0 gives Eq. ( 4).A similar derivation can also be based on the hypervirial theorem [37].

C. Local shifting invariance
Before presenting explicit applications of Eq. ( 4) to specific forms of Â, we first generalize to the fully position-resolved case.In generalization of the uniform coordinate displacement of Sec.II B, we consider the following local transformation on phase space, as parameterized by a three-dimensional vector field ϵ(r) [33,35]: The gradient ∇ i ϵ(r i ) is a 3×3 matrix, 1 denotes the 3×3 identity matrix and the superscript −1 indicates matrix inversion.Figure 1(a) depicts an illustration of the spatial transformation (5).The momentum transformation (6) has the following Taylor expansion to lowest order in the displacement field: The joint transformation ( 5) and ( 6) is canonical [33,35,38] and it hence preserves the phase space volume element, dr i dp i = dr i dp i , where the tilde indicates the transformed variables [right hand sides of Eqs. ( 5) and (6)].The Hamiltonian also remains unchanged (up to expressing the original via the new variables).Hence the partition sum Ξ is an invariant under the joint transformation ( 5) and ( 6) [33,35,36].Together with invariance of the integration measure, the setup implies that any average A = ⟨ Â⟩ = Tr Ψ eq Â is an invariant.This property holds despite the explicit occurrence of the shifting field ϵ(r) in the integrand, and hence A[ϵ] = A, where the left hand side carries the apparent dependence on the shifting field and the right hand side is the average expressed in the original variables where ϵ(r) is absent.
We use standard notation to express dependence on a function (so-called functional dependence) by bracketed arguments.
From the local Noether invariance we can conclude from functionally differentiating A[ϵ] = A with respect to the shifting field that The right hand side of Eq. ( 7) vanishes trivially due to the average A being independent of ϵ(r) in the original representation and hence δA/δϵ(r) = 0. Carrying out the functional derivative on the left hand side of Eq. ( 7) requires to functionally differentiate the equilibrium distribution, δΨ eq [ϵ]/δϵ(r) = −βΨ eq [ϵ]δH[ϵ]/δϵ(r), as follows from the chain rule.This allows to rewrite Eq. ( 7) upon using the product rule and re-ordering the resulting two terms as where we have evaluated both sides at vanishing shifting field, ϵ(r) = 0. Differentiating the transformed Hamiltonian with respect to the shifting field gives −δH[ϵ]/δϵ(r)| ϵ=0 = F(r) [33,35], where the position-resolved total force operator comprises the following three terms: The right hand side of Eq. ( 9) features the one-body kinematic stress operator τ (r) = − i δ(r − r i )p i p i /m [20], the one-body interparticle force density operator Fint (r) = − i δ(r − r i )∇ i u(r N ) [20], the standard form of the density operator ρ(r) = i δ(r − r i ) [1,2,20], and the external force field −∇V ext (r).(The force density operator F(r) defined in (9) also arises as the time derivative of the one-body current operator [20].In equilibrium the kinematic term in (9) reduces to a diffusive contribution: ∇ • ⟨ τ (r)⟩ = −k B T ∇ρ(r); we refer the Reader to Ref. [20] for details.) Equation ( 8) can then be written upon carrying out the functional derivatives and using Eq. ( 9) (on the left hand side) together with the chain rule (on the right hand side), as the following hyperforce sum rule that holds for a given form of Â = Â(r N , p N ): where Â can again feature additional parametric dependencies, such as on r.
For cases where the observable under consideration is independent of the momenta, i.e., Â ≡ Â(r N ), the second term on the right hand side of Eq. ( 10) vanishes and we obtain the coordinate-based local hyperforce sum rule: Equation ( 11) can be viewed as a generalization of the framework developed by Coles et al. [17], where they consider observables of the specific form Â = i a i δ(r − r i ), where a i is a unique property of particle i only, such as e.g. its charge or, when taking orientational degrees of freedom into account, its polarization [17].
As a consistency check, from integrating the locally resolved sum rules (10) and (11) over position r, i.e., applying dr to both sides of these equations, and observing that drδ(r − r i ) = 1, one retrieves the global Noether identity (4).Here the global external force is the only remaining non-trivial global force contribution, F0 ≡ dr F(r) = F0 ext , as the global interparticle force vanishes due to Newton's third law [31] and there is also no global diffusive effect due to vanishing boundary terms.Hence Eq. ( 4) continues to hold upon replacing F0 ext by the global total force operator F0 .As a further remark, by splitting off the kinetic term in Eq. (11), its left hand side can be re-written as β⟨ F(r) Â(r N )⟩ = β⟨ FU (r) Â(r N )⟩ − ∇⟨ρ(r) Â(r N )⟩, with the potential force density operator being given as the sum of interparticle and external contributions: FU (r) = Fint (r) − ρ(r)∇V ext (r).
In summary, the local force decomposition into ideal, interparticle and external contributions in Eq. (11) allows to obtain the following more explicit form, which holds, as we recall, provided that Â = Â(r N ) is independent of the momenta: For a given explicit form of Â, such as in the concrete examples discussed below, the sum rule (12) connects the three irreducible correlators ⟨ Fint (r) Â⟩, ⟨ρ(r) Â⟩, and ⟨ i δ(r − r i )∇ i Â⟩ in a formally exact and nontrivial way with each other.Setting Â = 1 recovers the YBG equation (which we take to imply thermal averages being taken), as then then last term in Eq. ( 12) vanishes and the remaining terms constitute Eq. (1).Paralleling the naming convention of the hypervirial theorem, which generalizes the standard virial theorem to include a further observable, Eq. ( 12) attains the status of a hyper-YBG equation or hyperforce balance relationship.Concrete applications thereof are shown below in Sec.III B. The correlators on the left hand sides of the sum rules (4), (10), and (11) also constitute covariances.We recall that the covariance of two observables Â and B, as defined via cov( Â, B) = ⟨ Â B⟩−⟨ Â⟩⟨ B⟩ measures the correlation of the fluctuations of the two observables around their respective mean.In the present case the mean force vanishes both globally, ⟨ F0 ext ⟩ = 0, and locally, ⟨ F(r)⟩ = 0. Hence we can formally subtract the vanishing averages and re-express ⟨ F0 ext Â⟩ = cov( F0 ext , Â) as well as ⟨ F(r) Â⟩ = cov( F(r), Â).Besides the conceptual difference between correlation and covariance, in practical sampling schemes it can be beneficial to work with covariances rather than correlations to reduce statistical noise, as we will demonstrate further below.
That Eq. ( 11) holds can again be verified a posteriori by phase space coordinate integration by parts on the right hand side.Due to the product rule two contributions result, one from the Boltzmann factor: i δ(r − r i )∇ i Ψ eq = −βΨ eq i δ(r − r i )∇ i H, and one from the Dirac distribution: Ψ eq i ∇ i δ(r − r i ) = −Ψ eq ∇ρ(r).Together with the factor Â their combination yields the left hand side of Eq. ( 11) upon identifying F(r) via Eq.( 9).The more general Eq. ( 10) follows analogously upon integrating by parts also with respect to the momenta.

A. Global shifting applications
We turn to applications and hence consider concrete examples for the general phase space function Â, which has remained so far unspecified in the above generic hyperforce framework.We start with investigating the global invariance (4), which as we demonstrate constitutes a powerful device both if Â is a global object or if it is locally resolved via dependence on a position.We first consider the seemingly trivial case Â = 1, for which of course ⟨1⟩ = 1 due to the correct normalization of Ψ eq .The right hand side of Eq. ( 4) vanishes and we obtain F 0 ext ≡ ⟨ F0 ext ⟩ = 0, i.e. the vanishing of the mean external force in equilibrium [31].This is intuitively expected, as can be seen by contradiction as follows.If the mean external force did not vanish, then the system would start to move on average [32] and hence it would not be in equilibrium.
Adressing the global external force and hence setting Â = F0 ext in Eq. ( 4) leads upon simplifying the right hand side via ext ⟩ = drρ(r)∇∇V ext (r) [34].Here the autocorrelation of the global external force (left hand side) equals up to a factor β the mean external potential energy curvature (right hand side).
By iteratively replacing Â with the composite F0 ext Â in Eq. ( 4), one can systematically generate higher order sum rules, starting with β⟨ , where the first term on the right hand side allows repeated application of Eq. ( 4) and the second term can be written via the external potential curvature.
The result is the following global second order hyperforce sum rule: where Â = Â(r N , p N ).The second term on the right hand of Eq. ( 13) side can alternatively be written as an integral over a correlation function as follows: β dr⟨ Âρ(r)⟩∇∇V ext (r), where ⟨ Âρ(r)⟩ is the correlation of Â and the local density operator.Alternatively to the present route via Eq.( 4), the sum rule (13) can equivalently be derived from second-order invariance of ⟨ Â⟩ against global shifting and hence calculating ∂ 2 A(ϵ 0 )/∂ϵ 0 ∂ϵ 0 = 0. Addressing locally resolved correlation functions on the basis of Eq. ( 4) allows to access a higher degree of spatial resolution.We first consider the case Â = ρ(r), which leads upon re-writing the right hand side of Eq. (4) via to the following identity: where we recall that on the right hand side ρ(r) = ⟨ρ(r)⟩ is the averaged density profile.Hence building the correlation of the density operator with the global external force acts to spatially differentiate the density profile.We recall that the density gradient ∇ρ(r), as it occurs on the right hand side of Eq. ( 14), follows alternatively from the YBG equation (1), which upon multiplication by β attains the form ∇ρ(r) = βF int (r)+βF ext (r), where the external force density is simply given as F ext (r) = −ρ(r)∇V ext (r).
To validate the Noether invariance theory and to investigate its implications for the use in force sampling methods, we turn to many-body simulations and consider the Lennrd-Jones (LJ) fluid as a representative microscopic model.The LJ pair potential ϕ(r) between two particles separated by a distance r has the familiar form ϕ(r) = 4ϵ[(σ/r) 12 − (σ/r) 6 ] with the energy scale ϵ and particle size σ both being constants.
As our above statistical mechanical derivations continue to hold canonically, we sample both via adaptive Brownian dynamics (BD) [41] with fixed number of particles, but also using Monte Carlo simulations in the grand canonical ensemble.Spatial inhomogeneity is induced by confining the system between two planar, parallel LJ walls.Each wall is represented by an external potential contribution V wall (z) that we choose to be identical to the LJ interparticle potential ϕ(r), but instead of the radial distance r evaluated as a function of the distance z perpendicular to the wall, , with L indicating the separation distance between the two walls.The specific choice of 12-6-wall potential is made for convenience only and it differs from the physically motivated 9-3-form (see e.g.Ref. [9]).
The system is periodic in the two directions perpendicular to the z-direction across the slit; a sketch is shown in Fig. 1(b).The wall separation distance is chosen as L = 10σ, with σ denoting the LJ particle size, and the lateral box length is also set to 10σ.The LJ potential is cut and shifted with a cutoff distance of 2.5σ.The reduced temperature is k B T /ϵ = 2 with ϵ denoting the LJ energy scale.We use N = 200 particles.Sampling is started after 10 8 time steps that are used for equilibration.The subsequent sampling runlength is 3×10 8 time steps which corresponds to ∼ 2000τ B , where τ B = γσ 2 /ϵ denotes the Brownian timescale with γ being the friction constant.All results that we show for correlators are obtained from evaluation as covariances.Subtracting the residual contribution from the product of the two mean values helps to remove artifacts that occur due to finite sampling.
The density profile of the confined LJ fluid, resolved as a function of the scaled position z/σ across the planar slit, is shown in Fig. 2(a).The shape of the spatial density variation features structured packing effects that appear adjacent to each wall and that become damped towards the middle of the pore.Turning to the density gradient, we present simulation results for both sides of Eq. ( 14) in Fig. 2(b).Equation ( 14) in the specific planar geometry reduces to β⟨F 0,z ext i δ(z −z i )/L 2 ⟩ = ∂ρ(z)/∂z, where δ(z − z i ) is a one-dimensional Dirac distribution, z i is the component of the vector r i across the pore, L 2 is the lateral system area, and the global external force has only a nonvanishing z-component given by  ext Fint(r)⟩ and β∇Fint(r), see Eq. (15).The former route carries less statisical noise and hence can serve a starting point for a force sampling scheme.(d) Comparison of β 2 ⟨ F0 ext F(r)⟩ and the local external potential curvature density βρ(r)∇∇Vext(r), see Eq. (16).
gives the gradient of the density profile, cf.Eq. ( 14), is surely not only at first glance very counter-intuitive.
We next consider the one-body interparticle force density operator Â = Fint (r), for which Eq. ( 4) yields Quite remarkably, Eq. ( 15) gives access to the gradient of the internal force density (right hand side) via sampling the correlation of the local internal force density with the global external force (left hand side).This relationship could be used in a force sampling scheme [13][14][15][16][17][18][19], where one obtains data for the force correlations and via spatial integration obtains the interparticle force density, which we demonstrate below in Sec.III C. We first present simulation results to illustrate the validity of Eq. ( 15) in Fig. 2(c).Remarkably, the results for β⟨ F0 ext Fint (r)⟩ carry much less statistical noise, as the need for building the numerical derivative ∇F int (r) is avoided.However, this is no panacea, as accurately sampling the correlation of the interparticle force density with the global external force also poses challenges to overall equilibration of the system.
Addressing the total force density operator Â = F(r) requires to complement the above considered interparticle force density Fint (r) with the ideal and external contributions.These two latter terms constitute mere variants of the density operator identity (14).First there is the external force density operator Â = Fext (r) = −ρ(r)∇V ext (r) which yields β⟨ F0 ext Fext (r)⟩ = −[∇V ext (r)]∇ρ(r), as ∇V ext (r) can be taken out of the phase space average on the left hand side of Eq. ( 14).Secondly, the diffusive force density, Â = −k B T ∇ρ(r), leads trivially to the gradient of Eq. ( 14).
Collecting all three terms (ideal, interparticle, and external) allows to obtain for the choice Â = F(r) a mixed global-local Noether identity: We present simulation results that validate the sum rule (16) in Fig. 2(d).We have checked that via spatial integration these results also validate the global variance identity β⟨ [34], which follows from Eq. ( 13) with operator Â = 1.For the presently considered system the global force-force correlation on the left hand side is only marginally (0.4%) smaller than the global mean potential curvature on the right hand side.

B. Local shifting applications
We turn to application of the locally resolved hyperforce sum rules (10) and (11).We recall from Sec. II C that the seemingly trivial case Â = 1 reduces the local hyperforce identity (11) to the locally resolved force balance relationship F(r) = ⟨ F(r)⟩ = 0.The definition of   the total force operator (9) and carrying out the average gives the more explicit form −k B T ∇ρ(r) + F int (r) − ρ(r)∇V ext (r) = 0, i.e., the first member (1) of the Yvon-Born-Green hierarchy [1].Typically this identity is derived from multiplying the equilibrium probability distribution function Ψ eq by the gradient of the interaction potential and integrating over the degrees of freedom of N − 1 particles.We emphasize that arguably the simplest possible application of Eq. ( 11) yields such a central result of liquid state theory with very little effort.From the Noetherian point of view the result was also obtained from applying the locally-resolved transformation ( 5) and ( 6) to the free energy [33,35].At the heart of these Noetherian derivations lies the invariance of the Hamiltonian, of the phase space integration measure, and hence of the partition sum.
We next consider setting Â = F0 ext in Eq. ( 10).This constitutes a valuable consistency check with the above global shifting of F(r) that led to Eq. ( 16) and which we can identically reproduce here.Arguably even more fundamentally the sum rule ( 16) can be obtained by considering mixed local-global shifting invariance at second order, i.e., building the mixed derivative ∂(δΩ/δϵ(r))/∂ϵ 0 = 0, where the grand potential Ω = −k B T ln Ξ is subject to the combined displacement r i → r i + ϵ 0 + ϵ(r i ).Furthermore increasing the spatial resolution and hence selecting Â = F(r) in Eq. ( 10) yields the recent Noether-constrained two-body force-correlation theory, which is discussed in detail in Refs.[36].The theory that is presented therein can hence be viewed as the special case of pure force-dependence within the hyperforce framework.
Reverting back to developing the general theory, here we generalize to higher orders by iteratively replacing Â(r N , p N ) by F(r ′ ) Â(r N , p N ) in Eq. ( 8).This leads to the following second order hyperforce sum rule: where evaluation of the right hand side at ϵ(r) = 0 is suppressed in the notation and δ 2 H[ϵ]/δϵ(r)δϵ(r ′ ) is discussed in Ref. [36].In case of no dependence of Â on momenta, i.e.Â = Â(r N ), the second term on the right hand side of Eq. ( 17) can be made more explicit as We proceed beyond forces by turning to energies with the aim of exploiting their thermal Noether invariance against shifting.consider both the global external potential energy Â = i V ext (r i ) as well as the global interparticle energy Â = u(r N ).Applying Eq. ( 11) yields in these two cases respectively the following sum rules: Simulation results that demonstrate the validity of Eqs.(19) and (20) are shown in Fig. 3 The potential energy identities ( 19) and ( 20) can be supplemented by considering the kinetic energy, Â = i p 2 i /(2m), which leads upon using Eq. ( 10) to the identity β⟨ F(r) i p 2 i /(2m)⟩ = −k B T ∇ρ(r).Treating then the entire Hamiltonian, Â = H, follows from adding up all three energy contributions.The result is the compact identity: β⟨ F(r)H⟩ = F(r) = 0.This possibly unexpected behaviour also holds for the global entropy.We choose the entropy operator Â = Ŝ ≡ −k B ln Ψ eq and obtain from Eq. ( 10) k −1 B ⟨ F(r) Ŝ⟩ = F(r) = 0, i.e., the correlation (as well as the covariance) of the entropy operator with the local force density vanishes.This behaviour is very different from the nontrivial fluctuation profile that is obtained from the covariance of the density operator with the global entropy [11,12].
As a final case we consider the center of mass i r i /N as a purely mechanical entity.We multiply by N , such that Â = i r i and obtain from Eq. ( 11) upon multiplication by −1 the result Integrating Eq. ( 21) over position yields a simple relationship of the correlator of the global external force and the center of mass: ⟨ F0 where the mean number of particles is N = drρ(r) ≡ ⟨N ⟩.Except for an additional sum over all particles, this global relationship is akin to the equipartition theorem.We present simulation results for both sides of the locally resolved Eq. ( 21) in Fig. 3(c).The accurate agreement of the respective profiles confirms that the identity (21) indeed offers a rather unusual route to gain access to the density profile.

C. Hyperforce sampling and equilibration testing
Besides the unexpected insights into the general correlation structure of equilibrium many-body systems that the thermal Noether invariance delivers, our results are useful for the careful assessment and construction of computer sampling schemes.Force sampling [13][14][15][16][17][18][19] in perhaps its most intuitive form [15] rests on spatial integration of the YBG equation ( 1) such that the density profile is obtained via ρ where ρ 0 = const is an integration constant and ∇ −1 is an inverse ∇ operator.The data input on the right hand side is obtained via sampling F int (r) = ⟨ Fint (r)⟩ and either The averages denote those that are being carried out in the simulation.In the present planar geometry ∇ −1 reduces to carrying out a simple position integral, which we make explicit below.Shown is data from the standard counting method according to Eq. ( 22) (orange lines), from force sampling FU (r) according to Eq. ( 23) (green dash-dotted lines), from global external force correlation sampling of ⟨ F0 ext ρ(r)⟩ according to Eq. ( 24) (blue solid lines), and from center-of-mass correlation sampling of ⟨ F(r) i ri⟩ according to Eq. ( 25) (red symbols).Results from the latter route are only displayed from the longest run [panel (d)] and they still display considerable scatter, whereas the results from the remaining three methods already agree very satisfactorily.
Summarizing, we can compare the results from four different routes: i) counting of particle occurrences in a position-resolved histogram, which constitutes the standard method; ii) force sampling [15] according to Eq. (1), iii) hyperforce sampling according to the global external force correlation in Eq. ( 14) with spatial integration post processing, and iv) center-of-mass-based hyperforce sampling according to Eq. ( 21).These routes are respectively given by the following explicit expressions: Here we choose the integration constant as ρ 0 = ρ(0) = 0 due to the divergent wall potential at z = 0.The averages on the above right hand sides denote the actual simulation data, all vectors have been projected onto the z-direction across the pore, and z i denotes the zcomponent of the particle position r i .In more detail, the operators on the right hand sides of Eqs. ( 22) and ( 23) are explicitly given as ρ(z , where we recall that L 2 is the lateral system size and the z-component of the interparticle force on particle i is f int i,z = −∂u(r N )/∂z i .Furthermore the operators on the right hand sides of Eqs. ( 24) and ( 25 Results for the density profile from the four routes ( 22)- (25) are shown in Fig. 4. The simulation parameters are identical as before [Fig.2].We display the four different statistical estimators for the density profile, as obtained after increasing runlength of (a) 10 5 , (b) 10 6 , (c) 10 7 , and (d) 3 × 10 8 simulation steps.We recall that as demonstrated above both in the numerical examples as well as in the formal statistical mechanical derivations, the results from all routes are formally identical.In practice, pronounced differences can be observed and these are due to the simulation averages being mere approximations for the true statistical mechanical equilibrium.
For example the routes ( 23) and ( 24) yield less statistical noise due to the spatial integration, but they however can instead accumulate systematic deviations.The expexted differences between the four methods are also consistently demonstrated by the fact that the results from the different routes mutually agree better for increasing runlengths.Nevertheless, in particular the routes (24) and ( 25) that involve global quantities are very sensitive to choice of runlength and they can hence serve as indicators of the overall quality of the sampling routine, even when quantities beyond the density profile are the very aim of the simulation.Having such tools for quality control can be particularly useful when investigating capillary and wetting phenomena [42][43][44][45][46] where surface phase transitions pose significant challenges for reliable prediction.
The Noetherian hyperforce framework allows us to easily go beyond the density profile and we wish to address the interparticle force density as a target, rather than the mere source that it played in contributing to Eq. (23) above for the force sampling.As a demonstration we use and contrast different estimators for the interparticle force density profile F int (r).The traditional counting method of filling a position-resolved histogram forms the baseline and Eq. ( 15) provides an alternative.These methods are respectively given by Figure 5 presents both canonical averages obtained via adaptive BD [41] as well as grand canonical Monte Carlo data.The chemical potential is chosen as µ/ϵ = 1 and the resulting average number of particles is ⟨N ⟩ = 136.5.In the corresponding adaptive BD simulation runs we have set N = 136, which remains sharply fixed in the course of time.The agreement between both sets of results confirms the expectation of independence of the sum rule validity on the choice of ensemble.This is based on the fact that the theoretical derivations continue to hold with fixed N , as we have also explicitly verified.Hence the mechanical effects that the Noether invariance against spatial displacement captures are oblivious to the presence of global particle number fluctuations.We recall that the later are precisely captured and quantified by the local compressibility [8][9][10].
The comparison of lower and higher quality statistical data, as obtained from sampling every step (top row) or only every 1000th step (bottom row) demonstrates that the force correlation method is a sensitive measure of the degree of sampling quality.

IV. CONCLUSIONS
In conclusion we have developed a statistical mechanical hyperforce framework in generalization of the YBG equilibrium force balance relationship (1).Our theory is based on previously developed global [31,32,34] and local [33,35] shifting transformations on phase space.These variable transformations leave the thermal physics invariant despite an apparent dependence on the transformation parameter.The parameter is a threedimensional globally constant vector in case of global symmetry, which applies to the entirety of the system, and a position-dependent three-dimensional vector field for the locally resolved case.Treating the corresponding phase space transformations according to Noether's invariant variational calculus [21] allows to systematically generate exact identities.
Here we have generalized this Noetherian concept to the equilibrium average of an arbitrary given phase space function Â.The resulting Noether identities couple in a specific manner the forces, which the underlying Hamiltonian generates, to the observable Â and its gradient with respect to the phase space variables.In the positionresolved case we obtain localized correlation functions, with the Dirac distribution generating microscopically sharp, but statistically coarse-grained and hence wellaccessible correlators.In detail, we have presented the global hyperforce sum rule (4) that applies to any given phase space function Â.The local versions comprise Eqs.(10) and Eq. ( 11), where the latter version applies to momentum-independent observables.Decomposing the total force operator into its ideal, interparticle, and external contributions leads to Eq. ( 12), which generalizes the equilibrium force density balance (1).An overview of these general identities is shown in Fig. 6(a).
For a variety of relevant concrete choices of the form of Â, we have demonstrated explicitly their validity via carrying out many-body simulations.This includes forces, energies, and entirely mechanical quantities such as the center of mass; see Fig. 6(b) for a summary of specific examples.We have shown that the sampling quality and equilibration properties depend significantly on the type of underlying sum rule.We argue that this behaviour forms a valuable asset for systematic assessment of simulation quality.
Our hyperforce identities complement the virial [1], hy-mutual interactions between all particles in the system vanishes, F0 int = P i r i u(r N ) ⌘ 0, as is due to Newton's third law or, analogously, to the translational invariance of u(r N ) against global displacement [30].
Re-ordering the two terms in Eq. ( 3) gives the following global hyperforce identity that holds for any given observable Â(r N , p N ): Here Â = Â(r N , p N ) can feature additional parametric dependence, such as on a generic position argument r.
The sum rule ( 4) is striking, as it relates the correlation of Â with the external force operator (left hand side) to the mean negative global coordinate derivative of Â (right hand side).
While Noether invariance enabled us to obtain the global hyperforce identity Eq. ( 4) constructively, one can verify its validity a posteriori by integration by parts in phase space on the right hand side.The derivatives r i then act on the probability distribution and exploiting again the Boltzmann form leads to Â P i r i eq = eq Â P i r i H = eq F0 ext Â, which gives the left hand side of Eq. ( 4) upon building the trace.Alternatively one can start with the Yvon theorem [1], h Âr i Hi = hr i Âi, which itself also follows from partial phase space integration [1].Summing the Yvon theorem over all particles i and noting that P i r i u(r N ) ⌘ 0 gives Eq. ( 4).A similar derivation can also be based on the hypervirial theorem [36].( The right hand side of Eq. ( 9) features the one-bo kinematic stress operator ⌧ (r) = P i (r r i )p i p i / [19], the one-body interparticle force density operat Fint (r) = P i (r r i )r i u(r N ) [19], the standard for of the density operator ⇢(r) = P i (r r i ) [1,2,19], an the external force field rV ext (r).Equation ( 8) can th be written upon carrying out the functional derivativ and using Eq. ( 9) together with the chain rule, as t following hyperforce sum rule that holds for a given for of Â = Â(r N , p N ): where Â can again feature additional parametric depe dencies, such as on r.
For cases where the observable under consideration independent of the momenta, i.e., Â ⌘ Â(r N ), the secon term on the right hand side of Eq. ( 10) vanishes and w obtain the coordinate-based local hyperforce sum rule Balance against Coles et al., Ref. [18].As a consistency check, from integrating the loca resolved sum rules ( 10) and (11) over position r, i. applying R dr to both sides of these equations, and o serving that R dr (r r i ) = 1, one retrieves the glob Noether identity (4).Here the global external force the only remaining non-trivial global force contributio F0 ⌘ R dr F(r) = F0 ext , as the global interparticle for vanishes due to Newton's third law [30] and there also no global di↵usive e↵ect due to vanishing bounda terms.Hence Eq. ( 4) continues to hold upon replaci F0 ext by the global total force operator F0 .As a further remark, by splitting o↵ the kinetic ter in Eq. ( 11), its left hand side can be re-written h F(r) Â(r N )i = h FU (r) Â(r N )i rh⇢(r) Â(r N )i, wi the potential force density operator being given as t sum of interparticle and external contributions: FU (r) Fint (r) ⇢(r)rV ext (r).
In summary, the local force decomposition into ide interparticle and external contributions in Eq. ( 11) allow to obtain the following more explicit form, which hold as we recall, provided that Â = Â(r N ) is independent the momenta: vanishes due to Newton's third law [30] and there is also no global di↵usive e↵ect due to vanishing boundary terms.Hence Eq. ( 4) continues to hold upon replacing F0 ext by the global total force operator F0 .
As a further remark, by splitting o↵ the kinetic term in Eq. ( 11), its left hand side can be re-written as h F(r) Â(r N )i = h FU (r) Â(r N )i rh⇢(r) Â(r N )i, with the potential force density operator being given as the sum of interparticle and external contributions: FU (r) = Fint (r) ⇢(r)rV ext (r).
In summary, the local force decomposition into ideal, interparticle and external contributions in Eq. ( 11) allows to obtain the following more explicit form, which holds, as we recall, provided that Â = Â(r N ) is independent of the momenta: We turn examples fo has remain hyperforce global inva tutes a pow locally reso sition.We for which o ization of and we obt the mean e tuitively ex general global local the system would start to move on average [31] and hence it would not be in equilibrium.
Adressing the global external force and hence setting Â = F0 ext in Eq. ( 4) leads upon simplifying the right hand side via h [33].Here the autocorrelation of the global external force (left hand side) equals up to a factor the mean external potential energy curvature (right hand side).
By iteratively replacing Â with the composite F0 ext Â in Eq. ( 4), one can systematically generate higher order sum rules, starting with h F0 , where the first term on the right hand side allows repeated application of Eq. ( 4) and the second term can be written via the external potential curvature.The result is the following global second order hyperforce sum rule: where Â = Â(r N , p N ).The second term on the right hand of Eq. ( 13) side can alternatively be written as an integral over a correlation function as follows: R drh Â⇢(r)irrV ext (r), where h Â⇢(r)i is the correlation of Â and the local density operator.Alternatively to the present route via Eq.( 4), the sum rule (13) can equivalently be derived from second-order invariance of h Âi against global shifting and hence calculating @ 2 A(✏ 0 )/@✏ 0 @✏ 0 = 0.
Addressing locally resolved correlation functions on the basis of Eq. ( 4) allows to access a higher degree of spatial resolution.We first consider the case Â = ⇢(r), which leads upon re-writing the right hand side of Eq. ( 4) via P i r i ⇢(r) = P i r i (r r i ) = P i r (r r i ) = r⇢(r) to the following remarkable identity: Hence building the correlation of the density operator with the global external force acts to spatially di↵erentiate the density profile.We recall that the density gradient r⇢(r), as it occurs on the right hand side of Eq. ( 14), follows alternatively from the YBG equation ( 1), which upon multiplication by attains the form r⇢(r) = F int (r) + F ext (r), where the external force density is simply given as F ext (r) = ⇢(r)rV ext (r).
To validate the Noether invariance theory and to in vestigate its implications for the use in force sampling methods, we turn to many-body simulations and con sider the Lennard-Jones (LJ) fluid as a representative microscopic model.As our above statistical mechanica derivations continue to hold canonically, we sample both via adaptive Brownian dynamics (BD) [40] with fixed number of particles, but also using Monte Carlo simula tions in the grand canonical ensemble.Spatial inhomo geneity is induced by confining the system between two planar, parallel LJ walls.Each wall is represented by an external potential contribution that is identical to the LJ interparticle potential (r), but instead of the radia distance r evaluated as a function of the distance z per pendicular to the wall.Hence the joint potential of both walls is V ext (z) = (z) + (L z), with L indicating the separation distance between the two walls.
The system is periodic in the two directions perpendic ular to the z-direction across the slit; a sketch is shown in Fig. 1(b).The wall separation distance is chosen as L = 10 , with denoting the LJ particle size, and the lateral box length is also set to 10 .The LJ poten tial is cut and shifted with a cuto↵ distance of 2.5 The reduced temperature is k B T/✏ = 2 with ✏ denot ing the LJ energy scale.We use N = 200 particles Sampling is started after 10 8 time steps that are used for equilibration.The subsequent sampling runlength is 3⇥10 8 time steps which corresponds to ⇠ 2000⌧ B , where 2 /✏ denotes the Brownian timescale with be ing the friction constant.All results that we show for correlators are obtained from evaluation as covariances Subtracting the residual contribution from the product o the two mean values helps to remove artifacts that occur due to finite sampling.
The density profile of the confined LJ fluid, resolved as a function of the scaled position z/ across the pla nar slit, is shown in Fig. 2(a).The shape of the spatia density variation features structured packing e↵ects that appear adjacent to each wall and that become damped towards the middle of the pore.Turning to the density gradient, we present simulation results for both sides o Eq. ( 14) in Fig. 2(b).The comparison of these a pri ori very di↵erent data sets indicates excellent agreement That the correlation of the density operator with the global external force operator indeed gives the gradient of the density profile, cf.Eq. ( 14), is surely a priori very counter-intuitive.
We next consider the one-body interparticle force den

Hyperforces (a)
i of the density operator ⇢(r) = P i (r r i ) [1,2,19], and the external force field rV ext (r).Equation ( 8) can then be written upon carrying out the functional derivatives and using Eq. ( 9) together with the chain rule, as the following hyperforce sum rule that holds for a given form of Â = Â(r N , p N ): where Â can again feature additional parametric dependencies, such as on r.
For cases where the observable under consideration is independent of the momenta, i.e., Â ⌘ Â(r N ), the second term on the right hand side of Eq. ( 10) vanishes and we obtain the coordinate-based local hyperforce sum rule: Balance against Coles et al., Ref. [18].As a consistency check, from integrating the locally resolved sum rules ( 10) and ( 11) over position r, i.e., applying R dr to both sides of these equations, and observing that R dr (r r i ) = 1, one retrieves the global Noether identity (4).Here the global external force is the only remaining non-trivial global force contribution, F0 ⌘ R dr F(r) = F0 ext , as the global interparticle force vanishes due to Newton's third law [30] and there is also no global di↵usive e↵ect due to vanishing boundary terms.Hence Eq. ( 4) continues to hold upon replacing F0 ext by the global total force operator F0 .As a further remark, by splitting o↵ the kinetic term in Eq. ( 11), its left hand side can be re-written as h F(r) Â(r N )i = h FU (r) Â(r N )i rh⇢(r) Â(r N )i, with the potential force density operator being given as the sum of interparticle and external contributions: FU (r) = Fint (r) ⇢(r)rV ext (r).
In summary, the local force decomposition into ideal, interparticle and external contributions in Eq. ( 11) allows to obtain the following more explicit form, which holds, as we recall, provided that Â = Â(r N ) is independent of the momenta:  [33].Here the autocorrelation of the global external force (left hand side) equals up to a factor the mean external potential energy curvature (right hand side).
By iteratively replacing Â with the composite F0 ext Â in Eq. ( 4), one can systematically generate higher order sum rules, starting with h , where the first term on the right hand side allows repeated application of Eq. ( 4) and the second term can be written via the external potential curvature.The result is the following global second order hyperforce sum rule: where Â = Â(r N , p N ).The second term on the right hand of Eq. ( 13) side can alternatively be written as an integral over a correlation function as follows: R drh Â⇢(r)irrV ext (r), where h Â⇢(r)i is the correlation of Â and the local density operator.Alternatively to the present route via Eq.( 4), the sum rule (13) can equivalently be derived from second-order invariance of h Âi against global shifting and hence calculating @ 2 A(✏ 0 )/@✏ 0 @✏ 0 = 0.
Addressing locally resolved correlation functions on the basis of Eq. ( 4) allows to access a higher degree of spatial resolution.We first consider the case Â = ⇢(r), which leads upon re-writing the right hand side of Eq. ( 4) via Hence building the correlation of the density operator with the global external force acts to spatially di↵erentiate the density profile.We recall that the density gradient r⇢(r), as it occurs on the right hand side of Eq. ( 14), follows alternatively from the YBG equation ( 1), which upon multiplication by attains the form r⇢(r) = F int (r) + F ext (r), where the external force density is simply given as F ext (r) = ⇢(r)rV ext (r).
It is interesting to note that Eq. ( 14), when written in covariance form as cov( F0 ext , ⇢(r)) = r⇢(r), mirrors closely the structure of the thermodynamic identity cov(N, ⇢(r)) = @⇢(r)/@µ ⌘ µ (r), with the local compressibility µ (r) [8][9][10].Here rather than the spatial gradient, the thermodynamic parametric derivative with respect to chemical potential occurs.Equation ( 14) ever not required in th form (14). Furthermore, pairs with unequal indice h F0 ext ⇢(r)i dist = F int (r), very di↵erent physical ob To validate the Noeth vestigate its implications methods, we turn to ma sider the Lennard-Jones microscopic model.As o derivations continue to h via adaptive Brownian d number of particles, but tions in the grand canon geneity is induced by con planar, parallel LJ walls an external potential con LJ interparticle potential distance r evaluated as a pendicular to the wall.H walls is V ext (z) = (z) + separation distance betw The system is periodic ular to the z-direction ac in Fig. 1 We next consider the o i i i i [19], the one-body interparticle force density operator Fint (r) = P i (r r i )r i u(r N ) [19], the standard form of the density operator ⇢(r) = P i (r r i ) [1,2,19], and the external force field rV ext (r).Equation ( 8) can then be written upon carrying out the functional derivatives and using Eq. ( 9) together with the chain rule, as the following hyperforce sum rule that holds for a given form of Â = Â(r N , p N ): where Â can again feature additional parametric dependencies, such as on r.
For cases where the observable under consideration is independent of the momenta, i.e., Â ⌘ Â(r N ), the second term on the right hand side of Eq. ( 10) vanishes and we obtain the coordinate-based local hyperforce sum rule: Balance against Coles et al., Ref. [18].
As a consistency check, from integrating the locally resolved sum rules ( 10) and (11) over position r, i.e., applying R dr to both sides of these equations, and observing that R dr (r r i ) = 1, one retrieves the global Noether identity (4).Here the global external force is the only remaining non-trivial global force contribution, F0 ⌘ R dr F(r) = F0 ext , as the global interparticle force vanishes due to Newton's third law [30] and there is also no global di↵usive e↵ect due to vanishing boundary terms.Hence Eq. ( 4) continues to hold upon replacing F0 ext by the global total force operator F0 .As a further remark, by splitting o↵ the kinetic term in Eq. ( 11), its left hand side can be re-written as h F(r) Â(r N )i = h FU (r) Â(r N )i rh⇢(r) Â(r N )i, with the potential force density operator being given as the sum of interparticle and external contributions: FU (r) = Fint (r) ⇢(r)rV ext (r).
In summary, the local force decomposition into ideal, interparticle and external contributions in Eq. ( 11) allows to obtain the following more explicit form, which holds, as we recall, provided that Â = Â(r N ) is independent of the momenta:  10) and (11) over position r, i.e., applying R dr to both sides of these equations, and observing that R dr (r r i ) = 1, one retrieves the global Noether identity (4).Here the global external force is the only remaining non-trivial global force contribution, F0 ⌘ R dr F(r) = F0 ext , as the global interparticle force vanishes due to Newton's third law [30] and there is also no global di↵usive e↵ect due to vanishing boundary terms.Hence Eq. ( 4) continues to hold upon replacing F0 ext by the global total force operator F0 .As a further remark, by splitting o↵ the kinetic term in Eq. ( 11), its left hand side can be re-written as h F(r) Â(r N )i = h FU (r) Â(r N )i rh⇢(r) Â(r N )i, with the potential force density operator being given as the sum of interparticle and external contributions: FU (r) = Fint (r) ⇢(r)rV ext (r).
In summary, the local force decomposition into ideal, interparticle and external contributions in Eq. ( 11) allows to obtain the following more explicit form, which holds, as we recall, provided that Â = Â(r N ) is independent of the momenta: Eq. ( 4), one can systematically generate higher order sum rules, starting with h , where the first term on the right hand side allows repeated application of Eq. ( 4) and the second term can be written via the external potential curvature.The result is the following global second order hyperforce sum rule: where Â = Â(r N , p N ).The second term on the right hand of Eq. ( 13) side can alternatively be written as an integral over a correlation function as follows: R drh Â⇢(r)irrV ext (r), where h Â⇢(r)i is the correlation of Â and the local density operator.Alternatively to the present route via Eq.( 4), the sum rule (13) can equivalently be derived from second-order invariance of h Âi against global shifting and hence calculating @ 2 A(✏ 0 )/@✏ 0 @✏ 0 = 0.
Addressing locally resolved correlation functions on the basis of Eq. ( 4) allows to access a higher degree of spatial resolution.We first consider the case Â = ⇢(r), which leads upon re-writing the right hand side of Eq. ( 4) via Hence building the correlation of the density operator with the global external force acts to spatially di↵erentiate the density profile.We recall that the density gradient r⇢(r), as it occurs on the right hand side of Eq. ( 14), follows alternatively from the YBG equation ( 1), which upon multiplication by attains the form r⇢(r) = F int (r) + F ext (r), where the external force density is simply given as F ext (r) = ⇢(r)rV ext (r).
It is interesting to note that Eq. ( 14), when written in covariance form as cov( F0 ext , ⇢(r)) = r⇢(r), mirrors closely the structure of the thermodynamic identity cov(N, ⇢(r)) = @⇢(r)/@µ ⌘ µ (r), with the local compressibility µ (r) [8][9][10].Here rather than the spatial gradient, the thermodynamic parametric derivative with respect to chemical potential occurs.Equation ( 14) methods, we turn to many-body sider the Lennard-Jones (LJ) flui microscopic model.As our above derivations continue to hold canon via adaptive Brownian dynamics number of particles, but also using tions in the grand canonical ensem geneity is induced by confining th planar, parallel LJ walls.Each w an external potential contribution LJ interparticle potential (r), bu distance r evaluated as a function pendicular to the wall.Hence the walls is V ext (z) = (z) + (L z), separation distance between the tw The system is periodic in the two ular to the z-direction across the s in Fig. 1(b).The wall separation L = 10 , with denoting the LJ lateral box length is also set to tial is cut and shifted with a cu The reduced temperature is k B T/ ing the LJ energy scale.We use Sampling is started after 10 8 tim for equilibration.The subsequent 3⇥10 8 time steps which correspond ⌧ B = 2 /✏ denotes the Brownian ing the friction constant.All resu correlators are obtained from eval Subtracting the residual contributi the two mean values helps to remo due to finite sampling.
The density profile of the confin as a function of the scaled positio nar slit, is shown in Fig. 2(a).Th density variation features structure appear adjacent to each wall and towards the middle of the pore.T gradient, we present simulation re Eq. ( 14) in Fig. 2

(b). The comp ori very di↵erent data sets indicate
That the correlation of the densi global external force operator inde of the density profile, cf.Eq. ( 14), counter-intuitive.
We next consider the one-body i  14), (15), and ( 16).(a) The density profile ⇢(r) of the confined system is shown as a reference.(b) Comparison of the correlator h F0 ext ⇢(r)i and the density gradient r⇢(r), see Eq. ( 14); the zoomed inset demonstrates the respective noise levels and it also shows the scaled force density sum FU (r) = Fint(r) + Fext(r) which equals r⇢(r) due to the local force balance.(c) Comparison of 2 h F0 ext Fint(r)i and rFint(r), see Eq. ( 15).The former route carries less statisical noise and hence can serve a starting point for a force sampling scheme.(d) Comparison of 2 h F0 ext F(r)i and the local external potential curvature density ⇢(r)rrV (r), see sity operator Â = Fint (r), for which Eq. ( 4) yields Quite remarkably, Eq. ( 15) gives access to the gradient of the internal force density (right hand side) via sampling the correlation of the local internal force density with the global external force (left hand side).This relationship could be used in a force sampling scheme [13][14][15][16][17][18], where one obtains data for the force correlations and via spatial integration obtains the interparticle force density, which we demonstrate below in Sec.III C. We first present simulation results to illustrate the validity of Eq. ( 15) in Fig. 2(c).Remarkably, the results for h F0 ext Fint (r)i carry much less statistical noise, as the need for building the numerical derivative rF int (r) is avoided.However, this is no panacea, as accurately sampling the correlation of the interparticle force density with the global external force also poses challenges to overall equilibration of the system.
Addressing the total force density operator Â = F(r) requires to complement the above considered interparticle force density Fint (r) with the ideal and external contributions.These two latter terms constitute mere variants of the density operator identity (14).First there is the external force density operator Â = Fext (r) = ⇢(r)rV ext (r) which yields h F0 ext Fext (r)i = [rV ext (r)]r⇢(r), as rV ext (r) can be taken out of the phase space average on the left hand side of Eq. ( 14).Secondly, the di↵usive force density, Â = k B T r⇢(r), leads trivially to the gradient of Eq. ( 14).
Collecting all three terms (ideal, interparticle, and external) allows to obtain for the choice Â = F(r) a mixed global-local Noether identity: We present simulation results that validate the sum rule (16) in Fig. 2(d).We have checked that via spatial integration these results also validate the global variance identity h F0 ext F0 ext i = R dr⇢(r)rrV ext (r) [33].For the presently considered system the global force-force correlation on the left hand side is only marginally (0.4%) smaller than the global mean potential curvature on the right hand side.

B. Local shifting applications
We turn to application of the locally resolved hyperforce sum rules (10) and (11).We recall from Sec. II C that the seemingly trivial case Â = 1 reduces the local   ⇢(r)rV ext (r) = 0, i.e., the first member (1) of the Yvon-Born-Green hierarchy [1].Typically this identity is derived from multiplying the equilibrium probability distribution function eq by the gradient of the interaction potential and integrating over the degrees of freedom of N 1 particles.We find it remarkable that arguably the simplest possible application of Eq. ( 11) yields such a central result of liquid state theory with very little e↵ort.
the spatial resolution and hence selecting Â = F(r) in Eq. ( 10) yields the recent Noether-constrained two-body force-correlation theory, which is discussed in detail in Refs.[35].The theory that is presented therein can hence be viewed as the special case of pure force-dependence within the hyperforce framework.
Here we generalize to higher orders by iteratively replacing Â(r N , p N ) by F(r 0 ) Â(r N , p N ) in Eq. ( 8).This leads to the following second order hyperforce sum rule: where evaluation of the right hand side at ✏(r) = 0 is suppressed in the notation and 2 H[✏]/ ✏(r) ✏(r 0 ) is discussed in Ref. [35].In case of no dependence of Â on momenta, i.e.Â = Â(r N ), the second term on the right hand side of Eq. ( 17) can be made more explicit as We proceed beyond forces by turning to energies with the aim of exploiting their thermal Noether invariance against shifting.We consider both the global external potential energy Â = P i V ext (r i ) as well as the global interparticle energy Â = u(r N ).Applying Eq. ( 11) yields in these two cases respectively the following sum rules: Simulation results that demonstrate the validity of Eqs.(19) and (20) are shown in Fig. 3    ⇢(r)rV ext (r) = 0, i.e., the first member (1) of the Yvon-Born-Green hierarchy [1].Typically this identity is derived from multiplying the equilibrium probability distribution function eq by the gradient of the interaction potential and integrating over the degrees of freedom of N 1 particles.We find it remarkable that arguably the simplest possible application of Eq. ( 11) yields such a central result of liquid state theory with very little e↵ort.ment r i !r i + ✏ 0 + ✏(r i ).Furthermore increasing the spatial resolution and hence selecting Â = F(r) in Eq. ( 10) yields the recent Noether-constrained two-body force-correlation theory, which is discussed in detail in Refs.[35].The theory that is presented therein can hence be viewed as the special case of pure force-dependence within the hyperforce framework.
Here we generalize to higher orders by iteratively replacing Â(r N , p N ) by F(r 0 ) Â(r N , p N ) in Eq. ( 8).This leads to the following second order hyperforce sum rule: where evaluation of the right hand side at ✏(r) = 0 is suppressed in the notation and 2 H[✏]/ ✏(r) ✏(r 0 ) is discussed in Ref. [35].In case of no dependence of Â on momenta, i.e.Â = Â(r N ), the second term on the right hand side of Eq. ( 17) can be made more explicit as We proceed beyond forces by turning to energies with the aim of exploiting their thermal Noether invariance against shifting.We consider both the global external potential energy Â = P i V ext (r i ) as well as the global interparticle energy Â = u(r N ).Applying Eq. ( 11) yields in these two cases respectively the following sum rules: D Simulation results that demonstrate the validity of Eqs. ( 19) and ( 20) are shown in Fig. 3 23) (green dash-dotted h F0 ext ⇢(r)i according to Eq. ( 24) (blue solid lines), and from center Eq. ( 25) (red symbols).Results from the latter route are only displa considerable scatter, whereas the results from the remaining three m The potential energy identities ( 19) and ( 20) can be supplemented by considering the kinetic energy, Â = P i p 2 i /(2m), which leads upon using Eq. ( 10) to the identity h F(r) Treating then the entire Hamiltonian, Â = H, follows from adding up all three energy contributions.The result is the compact identity: h F(r)Hi = F(r) = 0.This rather striking behaviour also holds for the global entropy.We choose the entropy operator Â = Ŝ ⌘ k B ln eq and obtain from Eq. ( 10) k 1 B h F(r) Ŝi = F(r) = 0, i.e., the correlation (as well as the covariance) of the entropy operator with the local force density vanishes.This behaviour is very di↵erent from the nontrivial fluctuation profile that is obtained from the covariance of the density operator with the global entropy [11,12].
As a final case we consider the center of mass P i r i /N as a purely mechanical entity.We multiply by N , such that Â = P i r i and obtain from Eq. ( 11) upon multiplication by 1 the result ulation results to illustrate the validity of Eq. ( 15) in Fig. 2(c).Remarkably, the results for h F0 ext Fint (r)i carry much less statistical noise, as the need for building the numerical derivative rF int (r) is avoided.However, this is no panacea, as accurately sampling the correlation of the interparticle force density with the global external force also poses challenges to overall equilibration of the system.
Addressing the total force density operator Â = F(r) requires to complement the above considered interparticle force density Fint (r) with the ideal and external contributions.These two latter terms constitute mere variants of the density operator identity (14).First there is the external force density operator Â = Fext (r) = ⇢(r)rV ext (r) which yields h F0 ext Fext (r)i = [rV ext (r)]r⇢(r), as rV ext (r) can be taken out of the phase space average on the left hand side of Eq. ( 14).Secondly, the di↵usive force density, Â = k B T r⇢(r), leads trivially to the gradient of Eq. ( 14).
Collecting all three terms (ideal, interparticle, and external) allows to obtain for the choice Â = F(r) a mixed global-local Noether identity: We present simulation results that validate the sum rule (16) in Fig. 2(d).We have checked that via spatial integration these results also validate the global variance identity h F0 ext F0 ext i = R dr⇢(r)rrV ext (r) [33].For the presently considered system the global force-force correlation on the left hand side is only marginally (0.4%) smaller than the global mean potential curvature on the right hand side.

B. Local shifting applications
We turn to application of the locally resolved hyperforce sum rules (10) and (11).We recall from Sec. II C that the seemingly trivial case Â = 1 reduces the local hyperforce identity (11)  ext at first and second order with an observable Â(r N , p N ).The local sum rules (10), (11), and (12) contain the localized force density operator F(r) and its interparticle contribution Fint(r) and they apply for full phase space-dependence Â(r N , p N ) (upper bubble) or coordinate-only dependence Â(r N ) (lower bubble).In the latter case alternative forms involve the total local force density F(r) or the splitting (1) into its three constituent contributions.(b) Sum rules that result from specific observable choices.Global examples set Â(r N , p N ) as the density operator ρ(r), the interparticle force density operator Fint(r), and the total local force density operator F(r), as corresponds to the sum rules ( 14), (15), and ( 16), respectively.The local examples include the total external potential energy i Vext(ri), the interparticle energy u(r N ), and the sum of all particle coordinates i ri, corresponding to the sum rules ( 19), (20), and (21), respectively.pervirial [37], equipartition [1] and Yvon [1] theorems.Despite certain formal similarities, we emphasize that the underlying phase space invariance is more fundamental than derivations based on ad hoc partial integration.Fur-thermore, the considered invariance operations naturally lead to correlations with either global or locally resolved forces, which are both simple to interpret and straightforward to acquire in simulations.
We have shown how the global hyperforce identity (4) can alternatively be obtained from the Yvon theorem [1].Hence, as anticipated in the discussion by Rotenberg in Ref. [13], the Yvon theorem can indeed be a relevant tool for force sampling.However, the general localized Noether sum rule (10) reaches beyond the Yvon theorem in terms of the momentum effects that are included.We have shown that the derivation of the momentumindependent sum rule (11) based on the Yvon theorem requires to apply the ad hoc localized choice Âδ(r − r i ) and summing over i.As a second step, treating the ideal contribution explicitly allows to identify the sequence of emerging terms as the one-body force operator F(r) at any position.Conversely, the Yvon theorem can be derived as a limit case from Noether invariance upon shifting only one given particle i according to r i → r i + ϵ 0 , and keeping unchanged all other particles coordinates r j with j ̸ = i.
Besides the theoretical connections that the Noether hyperforce sum rules establish, they can serve to carry out tests in theoretical and simulation approaches, with possible fruitful connections to the mapped-averaging force sampling framework [16].The hyperforce sum rules can also provide a starting point, together with the existing body of equilibrium sum rules [42][43][44][45][46], for the construction of new inhomogeneous liquid state approximations.We have exemplified the use of sum rules in providing gauges for equilibration quality of simulation data and we are confident in their future beneficial use in machine-learning approaches such as the recent neural functional theory [47,48].
In the context of the use of machine-learning in Statistical Mechanics [47][48][49][50][51][52][53][54] sum rules were shown to provide tests for the successful construction of neural functionals both in [47,48] and out of equilibrium [49].These sum rules amount to specific force properties, such as the vanishing of the global interparticle force [49] and the interrelation of different orders of direct correlation functions [47,48].The present much more general hyperforce framework can form much inspiration for such approaches as well as for the recent force-based density functional theory [35,59], which was compared [59] to standard fundamental measure theory [60].
In our theoretical derivations we have relied on the grand canonical ensemble, with fixed chemical potential µ and fluctuating number of particles N .Carrying out formal manipulations in this way is often more straightforward than working with fixed N , as is appropriate for a canonical treatment.(Temperature is constant in both ensembles.)A prominent example is to obtain the density profile as a functional derivative ρ(r) = δΩ/δV ext (r) where crucially µ is kept fixed, rather than N , upon building the functional derivative.This prototypical example demonstrates the elegance of working grand canonically, and one could expect that a similar situation applies for the thermal Noether invariance.This, however, is not the case.Rather the phase space shifting transformation, whether global by a constant ϵ 0 or locally resolved in position via a three-dimensional vector field ϵ(r), is an entirely mechanical operation that applies equally well canonically.The shifting invariance gives a powerful new route to correlation functions and and their sum rules, alternative to the traditional method of integrating over degrees of freedom, as pioneered by Yvon [61] and Born and Green [62].
A detailed account of global shifting in the canonical ensemble is provided in Ref. [32].The resulting Noether force identities are analogous in form to the results from a grand canonical treatment [31], with the sole (and trivial) difference of the definition of the respective ensemble averages.Here we find that the analogous situation holds for the hyperforce identities.As they originate from phase space transformations only, they are insensitive to the ensemble differences between the canonical and the grand ensemble.This theoretical fact is corroborated by our computer simulation results, where we have explicitly compared grand canonical Monte Carlo data and canonical results, with the latter obtained via sampling under adaptive BD time evolution [41].
We have used overdamped Brownian time evolution as a means to sample in thermal equilibrium.We find the adaptive Brownian dynamics time stepping algorithm of Ref. [41] to be a convenient choice for our present purposes.The principle validity of the hyperforce sum rules is nevertheless independent thereof and we expect careful use of either the simpler Euler-Maruyama method [41,63] or indeed Molecular Dynamics [63] to yield identical results.In our BD simulations, we have sampled all correlation functions at equal time as is the appropriate limit in equilibrium time evolution to recover static thermal ensemble averages.As an aside, Ref. [31] exploits Noether invariance in nonequilibrium dynamics.
As a final note, we return to classical density functional theory and its prowess in the description and incorporation of the fundamental force correlators that emerge from the hyperforce concept.We recall that classical density functional theory is based on a formally exact variational principle which amounts to solving the following Euler-Lagrange equation: Here c 1 (r, [ρ]) is the one-body direct correlation function of inhomogeneous liquid state theory.This is expressed as a density functional via c 1 (r, [ρ]) = −βδF exc [ρ]/δρ(r), where F exc [ρ] is the intrinsic excess Helmholtz free energy functional, which contains the interparticle interactions, and δ/δρ(r) denotes the functional derivative with respect to the density profile.Solving Eq. ( 28) for given T, µ, and V ext (r) yields results for the equilibrium density profile ρ(r), which is hence the central variable of density functional theory.
A multitude of connections with the current invariance theory emerge naturally.From the density profile, using the hyperforce identities one can obtain results for ⟨ F(r) i r i ⟩ via Eq.( 21), for ⟨ F0 ext F(r)⟩ via Eq.( 16), for ⟨ F(r) i V ext (r i )⟩ via Eq.( 19), and upon building the gradient of the density profile for ⟨ F0 ext ρ(r)⟩ via Eq.( 14).We can make further progress by noting that within density functional theory the interparticle force density is given by F int (r) = k B T ρ(r)∇c 1 (r), where we have suppressed the functional dependence on the density profile in the notation.That this relationship holds can be seen from building the gradient of Eq. ( 28) whereby −∇µ vanishes as the chemical potential is constant, multiplying by ρ(r), and comparing term-wise with the force density balance relationship (1).
Having obtained F int (r) in this (density functional) way gives access to ⟨ F(r)u(r N )⟩ via Eq.( 20).Building the gradient of the interparticle force density via the product rule yields ∇F int (r) = k B T [∇ρ(r)]∇c 1 (r) + k B T ρ(r)∇∇c 1 (r).Again in principle the right hand side is straightforward to evaluate in a typical numerical density functional study as data for both ρ(r) and c 1 (r) is accessible.As a result the correlator ⟨ F0 ext Fint (r)⟩ is available via the hyperforce sum rule (15).
Hence standard results that are obtained within the density functional framework allow to access a wealth of nontrivial force correlation structure.This additional information is not redundant.We compare with Evans and his coworkers' local compressibility χ µ (r) = ∂ρ(r)/∂µ, where similar to the present force setup, one obtains χ µ (r) from processing the density profile.In practice then analyzing χ µ (r) can shed significantly more light on the physics than what is apparent from the density profile alone, as has been demonstrated in a range of insightful studies on drying at substrates and the important phenomenon of hydrophobicity [8][9][10].

FIG. 1 :
FIG.1: Illustrations of the relevant geometries.(a) The shifting field (orange arrows) displaces the coordinates ri of all particles i by a vector ϵ(ri); the specific particle i is highlighted in red and all particles are identical.A corresponding change in momenta, see Eq. (6), compensates the spatial distortion such that the differential phase space volume element (integration measure) remains unchanged.(b) Planar geometry of the confined Lennard-Jones fluid between two smooth parallel soft attractive Lennard-Jones walls; σ is the particle size, z measures the distance across the planar pore, and L is the distance between the two walls.
The comparison of the a priori very different data sets shown in Fig. 2(b) indicates excellent agreement.That the correlation of the density operator with the global external force operator indeed FIG.2: Illustration of the sum rules(14),(15), and(16).The simulation results obtained from adaptive BD sampling of the LJ fluid confined between two parallel planar LJ walls.The profiles are shown as a function of scaled distance z/σ across the planar slit.(a) The density profile ρ(r) of the confined system is shown as a reference.(b) Comparison of the correlator β⟨ F0 ext ρ(r)⟩ and the density gradient ∇ρ(r), see Eq. (14); the zoomed inset demonstrates the respective noise levels and it also shows the scaled force density sum βFU (r) = βFint(r) + βFext(r) which equals ∇ρ(r) due to the local force balance.(c) Comparison of β 2 ⟨ F0ext Fint(r)⟩ and β∇Fint(r), see Eq.(15).The former route carries less statisical noise and hence can serve a starting point for a force sampling scheme.(d) Comparison of β 2 ⟨ F0 ext F(r)⟩ and the local external potential curvature density βρ(r)∇∇Vext(r), see Eq. (16).
(a) and (b), respectively.Here we have increased the overall density by reducing the lateral box size to 5σ and we sample N = 128 particles over time periods of 2000τ [data shown in Figs.3(a) and (b)] and of 8000τ [data shown in Figs.3(c)].

7 3 ×FIG. 4 :
FIG.4: Comparison of standard counting against force sampling.We show results from four different routes towards the density profile ρ(r) as obtained after (a) 10 5 , (b) 10 6 , (c) 10 7 , (d) 3 × 10 8 simulation steps.Shown is data from the standard counting method according to Eq. (22) (orange lines), from force sampling FU (r) according to Eq. (23) (green dash-dotted lines), from global external force correlation sampling of ⟨ F0 ext ρ(r)⟩ according to Eq. (24) (blue solid lines), and from center-of-mass correlation sampling of ⟨ F(r) i ri⟩ according to Eq. (25) (red symbols).Results from the latter route are only displayed from the longest run [panel (d)] and they still display considerable scatter, whereas the results from the remaining three methods already agree very satisfactorily.

FIG. 5 :
FIG.5: Comparison of different statistical estimators for the interparticle one-body force density profile according to Eq. (15).The results are obtained from the standard counting histogram method, Fint(r) = ⟨ Fint(r)⟩ (orangle lines), and from hyperforce sampling and spatial integration of ∇Fint(r) = β⟨ F0ext Fint(r)⟩ (blue lines) with data for the right hand side forming the basis.These methods are explicitly spelled out in Eqs.(26) and (27).Results are shown from grand canonical Monte Carlo simulations [panels (a) and (c)] and from adaptive BD simulations [(b) and (d)].The results stem from sampling 10 5 strongly correlated microstates at every simulation step [(a) and (b)] compared to better decorrelated configurations obtained from also 10 5 configurations, with the samples being taken only every 1000th simulation step [(c) and (d)].The results from sampling are symmetrized with respect to mirroring at the center of the pore, i.e., via building the arithmetic mean [Fint(z) − Fint(L − z)]/2, as is common practice in force sampling schemes.
(b).The wall s L = 10 , with denotin lateral box length is als tial is cut and shifted w The reduced temperatur ing the LJ energy scale Sampling is started afte for equilibration.The su 3⇥10 8 time steps which co ⌧ B = 2 /✏ denotes the ing the friction constant correlators are obtained Subtracting the residual c the two mean values help due to finite sampling.The density profile of as a function of the scal nar slit, is shown in Fig. density variation features appear adjacent to each towards the middle of th gradient, we present simu Eq. (14) in Fig. 2(b).T ori very di↵erent data set That the correlation of global external force ope of the density profile, cf.counter-intuitive.

FIG. 2 :
FIG.2: Simulation results obtained from adaptive BD sampling of the LJ fluid confined between two parallel planar LJ walls.The profiles are shown as a function of scaled distance z/ across the planar slit and they illustrate the sum rules (14),(15), and (16).(a) The density profile ⇢(r) of the confined system is shown as a reference.(b) Comparison of the correlator h F0 ext ⇢(r)i and the density gradient r⇢(r), see Eq. (14); the zoomed inset demonstrates the respective noise levels and it also shows the scaled force density sum FU (r) = Fint(r) + Fext(r) which equals r⇢(r) due to the local force balance.(c) Comparison of 2 h F0 ext Fint(r)i and rFint(r), see Eq. (15).The former route carries less statisical noise and hence can serve a starting point for a force sampling scheme.(d) Comparison of 2 h F0 ext F(r)i and the local external potential curvature density ⇢(r)rrV (r), see

FIG. 3 :
FIG. 3: Demonstration of Noether sum rules (19), (20), and (21), respectively based on the global external and interparticle energies, and on the center of mass.Shown are results from adaptive BD simulations for the LJ fluid between parallel LJ walls as a function of the scaled distance z/ .(a) Comparison of 2 h F(r) P i Vext(ri)i and Fext(r) = ⇢(r)rVext(r), see Eq.(19).(b) Comparison of 2 h F(r)u(r N )i and Fint(r), see Eq.(20).(c) Comparison of h F(r) P i rii and ⇢(r), see Eq. (21).
FIG. 3: Demonstration of Noether sum rules (19), (20), and (21), respectively based on the global external and interparticle energies, and on the center of mass.Shown are results from adaptive BD simulations for the LJ fluid between parallel LJ walls as a function of the scaled distance z/ .(a) Comparison of 2 h F(r) P i Vext(ri)i and Fext(r) = ⇢(r)rVext(r), see Eq.(19).(b) Comparison of 2 h F(r)u(r N )i and Fint(r), see Eq.(20).(c) Comparison of h F(r) P i rii and ⇢(r), see Eq. (21).
FIG. 3: Demonstration of Noether sum rules (19), (20), and (21), respectively based on the global external and interparticle energies, and on the center of mass.Shown are results from adaptive BD simulations for the LJ fluid between parallel LJ walls as a function of the scaled distance z/ .(a) Comparison of 2 h F(r) P i Vext(ri)i and Fext(r) = ⇢(r)rVext(r), see Eq.(19).(b) Comparison of 2 h F(r)u(r N )i and Fint(r), see Eq.(20).(c) Comparison of h F(r) P i rii and ⇢(r), see Eq. (21).
(a) and (b), respectively.Here we have increased the overall density by reducing the lateral box size to 5 and we sample N = 128 particles over time periods of 2000⌧ [data shown in Figs.3(a) and (b)] and of 8000⌧ [data shown in Figs.3(c)].

FIG. 3 :
FIG. 3: Demonstration of Noether sum rules (19), (20), and (21), respectively based on the global external and interparticle energies, and on the center of mass.Shown are results from adaptive BD simulations for the LJ fluid between parallel LJ walls as a function of the scaled distance z/ .(a) Comparison of 2 h F(r) P i Vext(ri)i and Fext(r) = ⇢(r)rVext(r), see Eq.(19).(b) Comparison of 2 h F(r)u(r N )i and Fint(r), see Eq.(20).(c) Comparison of h F(r) P i rii and ⇢(r), see Eq. (21).
FIG. 3: Demonstration of Noether sum rules (19), (20), and (21), respectively based on the global external and interparticle energies, and on the center of mass.Shown are results from adaptive BD simulations for the LJ fluid between parallel LJ walls as a function of the scaled distance z/ .(a) Comparison of 2 h F(r) P i Vext(ri)i and Fext(r) = ⇢(r)rVext(r), see Eq.(19).(b) Comparison of 2 h F(r)u(r N )i and Fint(r), see Eq.(20).(c) Comparison of h F(r) P i rii and ⇢(r), see Eq. (21).
FIG. 3: Demonstration of Noether sum rules (19), (20), and (21), respectively based on the global external and interparticle energies, and on the center of mass.Shown are results from adaptive BD simulations for the LJ fluid between parallel LJ walls as a function of the scaled distance z/ .(a) Comparison of 2 h F(r) P i Vext(ri)i and Fext(r) = ⇢(r)rVext(r), see Eq.(19).(b) Comparison of 2 h F(r)u(r N )i and Fint(r), see Eq.(20).(c) Comparison of h F(r) P i rii and ⇢(r), see Eq. (21).

FIG. 4 :
FIG.4: Illustration of results from four di↵erent routes towards the 10 7 , (d) 3 ⇥ 10 8 simulation steps.Shown is data from the standa from force sampling FU (r) according to Eq. (23) (green dash-dotted h F0 ext ⇢(r)i according to Eq. (24) (blue solid lines), and from center Eq. (25) (red symbols).Results from the latter route are only displa considerable scatter, whereas the results from the remaining three m F0 ext F(r)i ⇢(r)rrV ext (r) : Simulation results obtained from adaptive BD samf the LJ fluid confined between two parallel planar LJ The profiles are shown as a function of scaled dis-/ across the planar slit and they illustrate the sum 14), (15), and (16).(a) The density profile ⇢(r) of fined system is shown as a reference.(b) Comparthe correlator h F0 ext ⇢(r)i and the density gradient see Eq. (14); the zoomed inset demonstrates the ree noise levels and it also shows the scaled force density U (r) = F int (r) + F ext (r) which equals r⇢(r) due local force balance.(c) Comparison of 2 h F0 ext Fint (r)i F int (r), see Eq. (15).The former route carries less al noise and hence can serve a starting point for a force g scheme.(d) Comparison of 2 h F0 ext F(r)i and the loernal potential curvature density ⇢(r)rrV ext (r), see ).

FIG. 6 :
FIG.6: Overview of theoretical results.(a) Hyperforce sum rules that relate a given observable with the forces that act in the system.The global sum rules (4) and (13) couple the global external force operator F0ext at first and second order with an observable Â(r N , p N ).The local sum rules(10),(11), and (12) contain the localized force density operator F(r) and its interparticle contribution Fint(r) and they ap-