⬇ Download Markdown  ·  ⬇ Download LaTeX  ·  ⬇ original  ·  view original  ·  CSV: table 1

Unified Gauge–Geometry Symmetry for Equilibrium Statistical Mechanics

arXiv:2604.00354v1 [cond-mat.stat-mech]

1

Apr 2026

Department of Physics, Hanoi National University of Education,

136 Xuan Thuy, Cau Giay, Hanoi, Vietnam and Institute of Natural Sciences, Hanoi National University of Education, 136 Xuan Thuy, Cau Giay, Hanoi, Vietnam We present a symmetry-based framework for equilibrium statistical mechanics that formulates a single Lie group combining conventional spacetime symmetries with a recently identified phase-space gauge-shifting invariance [Müller et al., Phys. Rev. Lett. 133, 217101 (2024)]. Using Noether’s theorem, we obtain a set of general Ward identities together with previously unexplored cross-relations arising from the noncommutation of different symmetry generators. The approach extends standard many-body symmetries–such as translations, rotations, Galilean boosts, dilations, and particle exchange–by incorporating an internal gauge-shift symmetry within a unified group structure. The resulting Lie algebra suggests a hierarchy of exact identities that encompass established sum rules and indicate possible cross-coupling relations between distinct response and correlation functions. We also identify a Wigner–Eckart–Ward reduction that simplifies tensor–hyperforce correlators to two scalar radial spectra in isotropic fluids, and we outline an equivariant gauge-constrained-DFT formulation whose Euler–Lagrange equations are constructed to satisfy the corresponding Ward and cross-Ward constraints. This framework provides a consistent organizational basis for phenomena in liquids, mixtures, and interfaces, and may offer a symmetry-based perspective connecting structure, mechanics, and dynamics in many-body systems.

I. INTRODUCTION

Gauge invariance plays a foundational role across physics, linking symmetries to conservation laws via Noether’s theorem [1, 2]. In quantum field theory, local gauge invariance imposes fundamental constraints represented by Yang–Mills charge conservation [3, 4]. Gauge concepts are usually developed in the context of fields and fundamental forces [3, 5], yet recent work shows analogous structures in many-body statistical mechanics via a phase-space ‘shifting’ symmetry that leaves equilibrium averages unchanged. Schmidt and co-workers applied Noether’s calculus to classical statistical functionals, deriving exact balance identities (sum rules) for interparticle forces and introducing higher-order “hyperforce” fields that generalize force correlations to arbitrary observables [6–8]. These results first emerged in general Noether form and were then re-formulated in a fully gaugeinvariant way, with the shifting symmetry’s generators closing a noncommutative Lie algebra. Notably, even finite yet non-infinitesimal shifts in phase space have been shown to leave all equilibrium observables invariant, establishing a genuine statistical-mechanical gauge symmetry that generates systematic sum rules and correlation constraints, including those that apply to multi-species mixtures [9]. This form of thermal gauge invariance provides a new route for uncovering conservation-law analogues in equilibrium many-body systems. It yields not only the standard force-balance relations but also hyperforce sum rules that apply to arbitrary phase-space observables, in a manner reminiscent of Hirschfelder’s classical hypervirial theorem [7, 10, 11].

Classical and quantum many-body systems host a rich spectrum of symmetries, including spatial translations, rotations, and scale transformations, as well as time translations, Galilean boosts, particle-permutation or speciesexchange operations, and symmetries associated with order parameters. Each of these symmetries, according to Noether’s principle, gives rise to a corresponding conservation law or exact identity within the system [1]. Combining distinct symmetries has historically yielded many central results in physics. For example, in equilibrium liquids the Euclidean spatial invariances (translations+rotations, and uniform scaling) link structural correlations to thermodynamic responses via the compressibility and virial routes [12–14]. The translational symmetry enforces the compressibility sum rule relating the static structure factor to the isothermal compressibility, while rotational and Galilean invariance underlie the Irving–Kirkwood virial theorem connecting microscopic forces to pressure [13–15].

Space–time translations combined with Galilean boosts enforce Ward or Slavnov–Taylor identities for hydrodynamics and Kubo relations for viscosity and conductivity [12, 16]. Adding scale to translations and rotations in conformal or Schrödinger settings fixes correlation functions and dispersion relations with little freedom [17]. Rigid motions together with additivity support morphometric thermodynamics, which reduces interfacial free energies to four geometric invariants. Translational invariance with charge neutrality gives the Stillinger–Lovett moment conditions for electrolytes [18]. Time translation with microreversibility yields fluctuation theorems such as Jarzynski and Crooks [19, 20]. What remains absent is a single unified symmetry group that simultaneously incorporates all conventional symmetries together with the additional phase-space gauge shifting. Earlier formulations addressed each symmetry, or at most pairs of 2 symmetries, in isolation, and no general framework has existed that integrates geometric, dynamical, internal, and gauge-like invariances within a coherent theoretical structure.

In this work, we embed spatial symmetries– translational, rotational, and scaling–together with dynamical, compositional, and phase-space gauge-shifting symmetries into a single Lie-group structure. By analyzing the noncommuting elements of its Lie algebra, we obtain a series of cross-symmetry identities that extend existing relations between different symmetry sectors. The resulting equivariant framework suggests several connections that may link distinct aspects of equilibrium physics. Among these, a Wigner–Eckart–Ward reduction appears to yield selection rules that simplify tensorial and higherrank hyperforce correlation functions into scalar forms under combined symmetry operations, in close analogy with the reduction of tensor operators in the Wigner–Eckart theorem [21, 22]. The same framework also provides exact scattering-to-stress relations, leading to sum-rule identities that connect static two-point structure factors with stress correlations and viscoelastic responses. These relations may allow estimation of elastic and viscous properties without relying on direct stress autocorrelation measurements.

In addition, we propose an equivariant gaugeconstrained-DFT variational principle, conceived as an extension of classical density functional theory, whose Euler–Lagrange equations are constructed to satisfy the derived Ward identities and cross-Ward constraints. Within this formulation, equilibrium densities and orderparameter profiles are expected to respect the full set of symmetry-induced sum rules, promoting internal consistency between structure and thermodynamics. Established results such as the compressibility relation, virial pressure expression, and known hyperforce balance laws are recovered as limiting cases of the general framework.

II. THEORETICAL CALCULATION

settings, the Ward–Takahashi identity expresses the invariance of the generating functional under an infinitesimal symmetry and thus relates correlation functions of currents and fields. In the present equilibrium context, the same logic yields exact relations among ensemble averages and correlators that follow solely from symmetry and measure invariance. For further reading see Refs. [23–27].

Consider a classical many-body system in n spatial dimensions. For a fixed particle number N , the phase space is M N = (R n × R n ) N , with microstate (phasespace point) Γ = (r N ; p N ) ∈ M N and Liouville mea-

p N ≡ (p 1 , . . . , p N ) collect the particle positions r i ∈ R n and canonical momenta p i ∈ R n .

The N -particle Hamiltonian is

If the particle number is fixed (canonical ensemble), the same formulas apply with the sum over N removed. The normalization is the canonical partition function

and thermal averages are

The guiding idea of what follows is simple: once we know how a symmetry acts on microscopic variables, Liouville-volume preservation together with integration by parts under the Gibbs weight produces identities that relate the symmetry variation of observables to that of the energy.

3

A. Unified symmetry group

We introduce a unified symmetry group G that collects geometric–time kinematics, internal symmetries, and discrete involutions:

(7) Here Gal(n) is the Galilei group generated by spatial translations, rotations, Galilean boosts, and time translations; the dilation factor R D implements uniform spatial scaling and acts by automorphisms on Gal(n). The internal sector combines the phase–space gauge–shift group G sh , the species symmetry G sp (permutations of components or continuous mixing, depending on the model), and the order–parameter symmetry G ord (for example O(n) for vector order, U (1) for phase fields, or a tensor P representation). The final factor Z T

2 and Z 2 collects the discrete symmetries of time reversal T and parity P; they act by conjugation on all other factors and constrain correlation functions even though they are not connected Lie symmetries. The symbol × denotes a direct product, while the symbol ⋉ denotes a semidirect product.

The infinitesimal counterpart of the group G is its Lie algebra g ≡ Lie(G) = T e G, i.e. the tangent space at the identity element e of the connected component. Each generator X ∈ g defines a one-parameter subgroup g s = exp(sX) ∈ G. We work primarily at this Lie-algebra level because the Stein identities below follow from integration by parts applied to the corresponding infinitesimal P phase-space flows. The discrete factors Z T

2 and Z 2 are not connected to e and therefore do not contribute generators to g, although they constrain correlators by conjugation. The semidirect products in Eq. (7) represent how different symmetry sectors act on each other; at the Lie-algebra level this yields nonvanishing commutators and is the origin of the cross-Ward identities derived below. The gauge-shift factor G sh may be regarded as a suitable (in general infinite-dimensional) subgroup of the diffeomorphism group of space; its Lie algebra consists of smooth vector fields ϵ(r) and forms a subalgebra g sh ⊂ g. For background on Lie groups, Lie derivatives, and related geometric notions, see e.g. Refs. [23–25].

We use a (left) representation of g by phase-space vector fields on M: to each generator X ∈ g we associate a vector field v X on M, whose flow is the infinitesimal transformation of the microscopic variables generated by

X. The induced action on an observable A(Γ) is the Lie

derivative

and the representation map is a Lie-algebra homomor-

Here ∇ Γ = (∇ r N , ∇ p N ) is the full phase-space gradient with respect to Γ = (r N ; p N ), with ∇ r N = (∇ r 1 , . . . , ∇ r N ) and ∇ p N = (∇ p 1 , . . . , ∇ p N ). We reserve ∇ for the spatial gradient with respect to its displayed argument, and write ∇ r i and ∇ p i for derivatives with respect to particle coordinates and momenta.

A central ingredient is the gauge-shift generator G[ϵ] ∈ g, parametrized by a smooth spatial vector field ϵ(r). Under the representation map X 7→ v X it is represented on phase space by the cotangent-lift vector field

where (∇ϵ) ij = ∂ r j ϵ i is the Jacobian of ϵ (so

Equation (9) is the standard cotangent-lift canonical form of the local phase-space shifting symmetry; it is equivalent to the global shifting generator in Ref. [8] [Eq. (7) therein], rewritten here as a phase-space vector field. Thus v G is simply the phase-space representative of the Lie-algebra element G[ϵ]. When no confusion arises, we identify a generator with its representative and write

A straightforward calculation shows that v G preserves the Liouville phase-space volume, i.e. ∇ Γ · v G = 0. However, when weighted by the Gibbs density P, one finds a nonzero weighted divergence

with H ≡ H N (Γ) (or H − µN in the grand ensemble) understood as the energy entering the Gibbs weight.

Assuming sufficient regularity of A and vanishing boundary terms, integration by parts under the Gibbs measure implies the finite-temperature Stein identity [28] for any generator X ∈ g whose representative flow is

Here ⟨·⟩ denotes the thermal ensemble average (canoni-

Equation (11) follows from integration by parts under the Gibbs weight and therefore assumes that the phasespace domain has no boundary (as for periodic boundary conditions) or that the associated surface term vanishes. More generally, let the accessible configuration space for N particles be a region Ω N ⊂ (R n ) N with boundary ∂Ω N (e.g. Ω N = ({z > 0}) N for a planar hard wall). For a volume-preserving generator, the same calculation yields the boundary Stein identity

where v X r denotes the configuration-space part of v X and n is the outward unit normal on ∂Ω N . Equation (11)

4 is recovered whenever this surface contribution vanishes, which holds for periodic boundaries, for sufficiently confining soft potentials, or for generators whose flows are tangent to the boundary (v X r ·n = 0 on ∂Ω N ). When a boundary breaks a symmetry, keeping the surface term in Eq. (12) produces boundary Ward identities; Eq. (42) is an example for a planar wall.

B. Lie algebra: nonzero commutators

We denote T (a) as the generator of spatial translations by the constant vector a; R(Θ) is the rotation generator with antisymmetric matrix Θ acting on R n ; D is the (spatial) dilation generator; B(u) is the Galilean boost generator with constant velocity u; t(τ ) is the time-translation generator by the amount τ ; and Π[χ] is the momentum–fiber translation generated by the smooth field χ(r). The symbols ϵ 1,2 denote two independent gauge fields. The composite field χ ϵ,u that appears in the boost–gauge bracket is defined particlewise, for a particle of mass m i , by

The explicit forms of the operators are given by

tors are (see App. A for a detailed derivation):

In what follows we use the Stein identity (11) to derive the constraints collected in Sec. II C. For any generators X, Y ∈ g, the phase-space representation X 7→ v X is a Lie-algebra homomorphism, so the induced Lie derivatives satisfy

with [X, Y ] ∈ g the Lie bracket. Consequently, each nonvanishing commutator in Eqs. (22)–(28) yields a crosscorrelation constraint by applying Eq. (11) to the bracket generator [X, Y ]. For generators parametrized by arbitrary test fields, this Ward identity localizes and produces the correlator relations given below, in direct analogy with Ward–Takahashi/Schwinger–Dyson mechanisms in field theory [4, 27].

C. Gauge master identity and cross-Ward

constraints

1. Local gauge shift × global U (1)

The Lie bracket of phase-space vector fields U, V representing the generators is

while for spatial vector fields ϵ 1 , ϵ 2 we use

The nonzero commutators involving the gauge generaP Introduce the microscopic density ρ̂(r) =

i δ(r − r i )

δ(r − r j )] with pair force f ij , where δ(·) is the Dirac delta function. For any configuration-only observable B(r N ) define the cross-correlators

In a translationally invariant bulk fluid (V ext = 0), the one-point averages ⟨ρ̂(r)⟩ = ρ and ⟨ F̂ int (r)⟩ = 0 are constants. The mixed correlators C Bρ (r) and C BF int (r) are nevertheless nontrivial for local or translation-covariant insertions B (e.g. B = δ ρ̂(0) or B = δ ρ̂ −k ): translational invariance then implies dependence only on the separation from the insertion point, i.e. on r (or k). Only for translationally invariant B (in particular B = 1) do these correlators collapse to one-point averages and the bulk identity reduces to 0 = 0.

P The local spatial gauge-shift symmetry G[ϵ] =

i ϵ(r i )· ∇ r i , together with global U(1) particle-number invariance, 5 implies an exact equilibrium identity. In bulk (V ext = 0) and for configuration observables B(r N ), and upon the standard separation of contact terms at the insertion where ρ ≡ ⟨ρ̂(r)⟩ is the bulk number density, the same gauge–U(1) identity yields the divergence form

6

with n ∥ the outward normal on the lateral surface ∂A × (0, L). For periodic boundary conditions parallel to the wall (or, equivalently, for large A with short-ranged correlations), the lateral surface integrals vanish. Taking L → ∞ removes the z = L terms because C δρF int (r) → 0 and ∇h(r) → 0 far from the wall. The only remaining contribution comes from the wall plane at z = 0, whose outward normal is −ẑ (hence the appearance of the difference between z = L and z = 0 + above). Rearranging then yields the wall contact relation [Π[χ], G[ϵ]] in Eq. (27) encodes the covariant way this statement is preserved under local gauge shifts.

dynamic identity Define the time–domain cross correlators

and their Fourier transforms

7 Combining Eq. (48) with the gauge master identity,

gauge–boost (or gauge–time) cross constraint

(49)

For single-mass systems, Galilean invariance identifies momentum and current densities, so Eq. (49) equivalently links C BF int and C Bj . In the equal-time limit, both sides vanish consistently with Eq. (45). The appearance of Π in Eq. (25) is the geometric reason why boosts couple the gauge sector to current correlations.

9. Gauge × gauge: integrability and Lie–compatibility

The gauge sector closes under commutation, [G[ϵ 1 ], G[ϵ 2 ]] = G[[ϵ 1 , ϵ 2 ]], consistent with the scalar–potential structure of the master identity. It is convenient to write the vector correlator as a differential one–form,

where d denotes the exterior differential. Thus ϖ B is the covariant representation of the same correlator, ϖ B,i = C BF int ,i dr i , which exposes its gradient origin and implies the integrability condition.

∇ × C BF int (r) = 0, k × C BF int (k) = 0, so Lie–dragging the one–form equals taking a gradient of the Lie–dragged scalar. The commutator of two such draggings satisfies

and by Eq. (54) this reduces to

The same reasoning applied to ϖ B = k B T dC Bρ shows that the gradient representation is preserved under any finite composition of local gauge flows. In particular, the transverse projector annihilates C BF int in every frame fixed by a gauge diffeomorphism, hence longitudinality

The integrability also constrains topology. On periodic domains the potential C Bρ is periodic, hence Eq. (52) holds and no harmonic one–form component

Rcan enter ϖ B . At boundaries Stokes’ theorem converts Σ dϖ B = 0 into a statement about the line integral of C BF int along ∂Σ, which is the geometric origin of the wall contact relations once the gauge master identity is invoked in the bulk. Altogether, closure of the gauge algebra and exactness of ϖ B ensure that the scalar potential C Bρ is the unique generator of the mixed force correlator, that the full vector is fixed by its longitudinal projection, and that nested local gauge transformations act compatibly on both scalar and vector sectors.

10. Dynamic cross identities and the stress channel

(51)

independent of isotropy or central forces. In simply connected domains every closed one–form is exact, hence there exists a scalar potential Ψ B with C BF int = ∇Ψ B . Equation (50) identifies Ψ B = k B T C Bρ , so the longitudinal representation holds globally and circulation vanishes along any closed loop,

The dynamic stress relation follows from two geometric inputs. The Irving–Kirkwood identity gives the divergence of the microscopic stress in terms of the internal force density, ∂ j σ̂ ij = F̂ int,i , where σ̂ ij is the microscopic Irving– Kirkwood stress tensor and F̂ int,i is the ith component of the microscopic internal–force density F̂ int . Hence in Fourier space i k j σ̂ ij (k, t) = F̂ int,i (k, t). Taking the mixed correlator with a configuration observable B and contracting with k i yields

In Fourier space a curl–free vector is necessarily longitudinal. Writing C BF int (k) = −i k ϕ B (k) and using the divergence identity i k·C BF int (k) = −k B T k 2 C Bρ (k) gives ϕ B (k) = −k B T C Bρ (k) and the full vector relation

The gauge master identity in the time–frequency domain,

8

dynamic density–density correlator, so that

The right-hand side is the dynamic density spectrum, proportional to the measured dynamic structure factor under the usual conventions. Nonnegativity of C δρδρ (k, ω) implies M L (k, ω) ≥ 0 at each frequency, and for any nonnegative weight w(ω) one obtains the sharp integral bound

(61)

Inverse transforming Eq. (59) gives the equal–time rela-

which at small k connects the static longitudinal stress to the compressibility via S(0) = ρ k B T κ T when S includes the self contribution.

The continuity identity in Eq. (48) combines with Eq. (59) into

(62)

which links stress to the longitudinal current channel. In generalized hydrodynamics this identity pins the part of the longitudinal memory function that enters sound and Brillouin peaks. In the hydrodynamic regime one has C δρδρ (k, ω) controlled by the isothermal sound speed and the longitudinal viscosity; Equation (60) therefore fixes the corresponding stress spectrum without direct stress measurements and provides a route from dynamic scattering to the combination ζ + 4η/3 at low k and small ω, where ζ is the bulk viscosity and η is the shear viscosity.

Causality of the retarded longitudinal stress susceptibility entails analyticity in the upper half-plane and Kramers– Kronig dispersion relations [30]. Through the fluctuation–dissipation theorem [31] these analyticity properties are inherited by the correlator representation. Equation (60) then implies that the dispersion constraints and sum rules for M L are

the same as those of C δρδρ , including

and the small–k compressibility limit.

D. Wigner–Eckart–Ward reduction

Rotational invariance and the gauge master identity combine to fix the tensorial structure of mixed correlations between an irreducible tensor T (ℓ) and the gauge–shift one–body field J G . Here we take a concrete microscopic realization. The rank–ℓ irreducible tensor is built from the local density via spherical harmonics,

(63)

with r̂ ′ = r ′ /r ′ and w ℓ a short–ranged radial weight and Y ℓm the usual harmonics (bond–orientational choices recover Q ℓm [32, 33]). We take J G as the gauge–shift one– body field [7, 34]

(64)

configuration–only B. We note that J G is different from the observable-dependent hyperforce density of Schmidt et al. [9], Ŝ A (r) = σ(r) Â (with σ the local shifting operator, Ref. [9], Eq. (20)). Here we define J G as the universal one-body gauge (Noether) field entering the configurationsector gauge master identity, ⟨B F̂ int (r)⟩ = −⟨B J G (r)⟩ (up to the standard contact term at the insertion point). With Eq. (64), J G (r) = −k B T ∇ρ̂(r), this object is simply the thermal/entropic (ideal) force density and may be written as the divergence of the ideal stress, J G = ∇· σ̂ id

The symbol ⊗ denotes the tensor product of irreducible tensors, so

and all directional content is carried by the angular dependence of the scalar correlator with δρ, while the vector character is supplied by k̂. See App. B for details.

Decompose the scalar spectrum using irreducible harmonics,

where A ℓm (k) are the scalar radial amplitudes obtained by angular projection,

which depend only on the wavenumber magnitude k = |k| by isotropy; re-expand the product k̂ Y ℓm into irreducible vector or tensor harmonics. The vector k̂ transforms as a spherical tensor of rank one, hence the tensor product 1 ⊗ ℓ reduces into total angular momentum L = ℓ − 1, ℓ, ℓ + 1. Longitudinality excludes the purely toroidal channel associated with L = ℓ, since the gradient of a scalar generates only radial and poloidal components. The mixed correlator therefore contains exactly two channels,

with direction carried by the irreducible vector or ten(L) sor harmonics Y

m ( k̂) and physics captured by the two 9 scalar radial spectra R ℓ±1 (k). The reduced coefficients C ℓ,1→L depend only on angular momentum algebra and are independent of k̂ and k.

components and C 2m,i (k) ≡ ⟨Q 2m (0) J G,i (k)⟩. Here i = −1, 0, 1 denotes a spherical vector component (conversion to Cartesian components is a fixed unitary map and does not affect the radial spectra). The selection rule gives a decomposition into L = 1 and L = 3 only,

The reduction has a direct operational meaning. Angle– resolved or polarized scattering that couples to Q 2m isolates the L = 1 and L = 3 channels through the kernels in Eqs. (69)–(70), which allows one to extract α 1 (k) and α 3 (k) from measured cross responses without reference to stress correlators. In the time domain the same selection rule applies to longitudinal spectra, once the continuity identity has been folded into the gauge relation, so that dynamic analogues of Eq. (68) follow with ω dependence carried entirely by the two radial functions. The Wigner– Eckart–Ward reduction thus provides a compact map from symmetry to measurable amplitudes, explains the vanishing of transverse components, and collapses tensor– hyperforce correlations to two scalar spectra in isotropic central–force fluids.

with α 1 and α 3 the two radial spectra and ⟨2 m; 1 i|L M ⟩ Clebsch–Gordan coefficients [22]. The absence of the L = 2 sector follows from the longitudinal/gradient structure in Eq. (65): the product k̂ Y 2m ( k̂) contains only L = 1 and L = 3 (equivalently, ⟨2 0; 1 0|2 0⟩ = 0), so the toroidal L = 2 vector–harmonic channel is forbidden. This expresses the purely longitudinal nature of J G , consistent with P T (k) C δρF int (k) = 0.

Orthogonality of the harmonics produces explicit projectors for the radial spectra. Introduce the angular kernels 1

E. Equivariant gauge-constrained DFT

The Ward and cross–Ward identities link structure and mechanics at equilibrium. To build the gauge constraint already at the level of stationary points, we extend the usual DFT field content to the multiplet Φ = (ρ, J G ), where ρ(r) is the one–body density and J G (r) is the gauge–shift one–body field. The role of the auxiliary onebody field J G is to represent the universal ideal/entropic force density that appears on the right-hand side of the configuration-sector gauge master relation. Using the microscopic definition already introduced in Eq. (64), one has for any configuration observable B(r N )

and let N L denote the corresponding normalization constants fixed by

read

Small wave number scaling follows directly from Eq. (68) (or Eq. (65)), which contains an explicit prefactor k: hence C 2m,i (k) = O(k) as k → 0 (at fixed reduced amplitudes), and α 1 (k) and α 3 (k) remain finite at long wavelength. The two spectra are therefore genuine radial response functions in the sense of Wigner–Eckart, with all angular dependence stripped away.

For the translation-covariant insertions used throughout the bulk discussion (e.g. B = δ ρ̂(0)), the last term in Eq. (75) is a local contact contribution at the insertion point and vanishes for r ̸ = 0. Thus, away from the contact D term one

Erecovers the compact bulk master relation B F̂ int (r) = − ⟨B J G (r)⟩, which is equivalent to the gauge–U(1) master identity Eqs. (31)–(32).

The hyperdensity functional theory of Ref. [34] constructs density functionals for general observables A[ρ] via the hyperdirect correlation functional c A (r, [ρ]) = δA[ρ]/δρ(r) together with a hyper–Ornstein–Zernike relation for the hyperfluctuation profile χ A . Our construction here has a different purpose: we keep the physical Hamiltonian fixed and introduce J G as an auxiliary field so that the stationary Euler–Lagrange equations enforce the gauge Ward identity and its correlated consequences.

In standard classical DFT the intrinsic free energy is decomposed as F [ρ] = F id [ρ] + F ex [ρ]. The ideal-gas 10 contribution is given explicitly by

where Λ is the constant thermal de Broglie wavelength; the additive constant ln Λ n can be absorbed into a shift of µ.

where F ex [ρ] is the interaction part of the intrinsic Helmholtz free-energy functional. We approximate it by a hard-sphere reference treated in fundamental measure

11

× 10 2

C F int r

Magnitude[

/ 8 ] 4 k

FIG. 1. Radial validation of the gauge–U(1) identity ∇ · C δρF int = k B T ρ 2 ∇ 2 h in a bulk Lennard–Jones fluid (T ∗ = 1.5, ρ ∗ = 0.85). Solid line: shell–averaged divergence ⟨∇· C δρF int ⟩ r computed from the density–internal–force cross–correlator measured in MD. Dashed line: k B T ρ 2 ∇ 2 h(r) obtained from the total correlation h(r) = g(r) − 1 of the same run.

relations in Sec. II C 10. The parameter λ G controls the cost of fluctuations transverse to the constraint (and can improve numerical conditioning), while the Ward relation for mixed ρ–J G correlators remains fixed by the covariant

Finally, for observables B(r N , p N ) with explicit momentum dependence, the microscopic gauge Ward identity contains additional kinematic source terms generated by the momentum part of the gauge generator. A field-level implementation of those identities would require extending the multiplet Φ beyond (ρ, J G ), for example by adding the current or kinetic-stress sector. The present minimal functional is therefore tailored to the configuration sector used throughout this work.

III. SIMULATION

We simulate a one–component Lennard–Jones (LJ) fluid in a cubic periodic box containing N = 6912 particles. The pair potential is

where ε is the characteristic energy scale corresponding to the depth of the potential well, while σ specifies the characteristic length scale that defines the effective particle diameter. The potential is truncated at r c = 2.5σ without applying long–range tail corrections.

Molecular dynamics trajectories are generated with HOOMD-blue [35, 36] using a velocity–Verlet integrator. Temperature is controlled with a Bussi stochastic velocity–rescaling thermostat (canonical sampling), and the simulations are performed at reduced temperature

FIG. 2. Shell-averaged magnitudes of the longitudinal (circles) and transverse (squares) projections of the density–internal-force cross correlator C δρF int (k) = V −1 δ ρ̂ −k F̂ int (k) for a bulk Lennard–Jones fluid at T ∗ = 1.5 and ρ ∗ = 0.85. The transverse signal is strongly suppressed relative to the longitudinal one, consistent with the gauge×rotation cross identity P T (k) C δρF int (k) = 0. [Eq. (38)], while the small-k growth of the longitudinal branch follows the gauge–U(1) prediction |P L C δρF int (k)| ∝ k |S(k)−1| from Eq. (35).

T ∗ = k B T /ε and reduced density ρ ∗ = ρσ 3 . Each system is equilibrated for 10 7 integration steps followed by a production run of 10 6 steps.

Figure 1 shows stringent real-space validation of the gauge–U(1) divergence identity in Eq. (35) for the bulk Lennard–Jones state point (T ∗ = 1.5, ρ ∗ = 0.85). The two sides of Eq. (35)–the shell-averaged divergence ⟨∇· C δρF int ⟩ r measured directly from the trajectory and the structural expression k B T ρ 2 ∇ 2 h(r) obtained independently from h(r) = g(r) − 1–are essentially indistinguishable on the scale of the plot over the full r range shown. In particular, the strong near-contact peak, the subsequent negative lobe at the first minimum, and the damped oscillations across successive coordination shells coincide in both phase and amplitude, demonstrating that the mixed density–internal-force response is completely fixed by the equilibrium pair structure as predicted by the Ward identity. At larger separations both signals rapidly decay toward zero and remain consistent with zero within the remaining statistical noise, as expected for a homogeneous bulk fluid with short-ranged interactions.

In Fig. 2 we resolve the Fourier–space structure of the density–internal–force cross correlator, C δρF int (k) = V −1 ⟨δ ρ̂ −k F int (k)⟩, by projecting onto the longitudinal and transverse sectors with P L (k) = k̂ k̂ ⊤ and P T (k) = 1 − P L (k), followed by thin |k|-shell averaging of the magnitudes. Across all accessible wave numbers the signal is overwhelmingly longitudinal: ⟨|P T C δρF int |⟩ shell ≪ ⟨|P L C δρF int |⟩ shell . This pure longitudinality is the expected consequence of rotational invariance combined with the gauge identity, Eq. (38), which forbids a trans-

IV. CONCLUSION

FIG. 3. Test of the momentum–fiber/gauge identity, Eq.

(45), using the

signed longitudinal current cross correlator

The exact prediction is zero (horizontal line); error bars are block-averaged standard errors. Units are Lennard–Jones re-

verse component in an isotropic, central-force fluid. The small but finite transverse baseline is attributable to finite statistics, discrete anisotropy of the reciprocal grid, and the use of shell averages of absolute values.

The longitudinal amplitude grows approximately linearly with k at small wave number, in agreement with the gauge–U(1) master identity, Eq. (35), which gives

Fourier-space evidence for the gauge×rotation cross-Ward constraint and illustrates the reduction of the vector hyperforce response to a single radial spectrum in isotropic liquids.

The isothermal compressibility κ T was extracted by two independent routes. Using the compressibility relation

contribution, we obtain the value κ T = 3.612 × 10 −2 . An independent estimate follows from the gauge–U(1) master identity [Eq. (35)] gives κ T = 3.595(70) × 10 −2 . The calculated isothermal compressibility values from both methods are in good agreement.

The two estimators provide a parameter-free check of the first-moment (gauge×dilation) identity, Eq. (39); the resulting values 3.664 × 10 −2 and 3.648 × 10 −2 are in close agreement, confirming that the same κ T governs both the density–density and density–hyperforce sectors.

We test the momentum–fiber/gauge cross identity, Eq. (45), which predicts vanishing equal-time current cross correlators with any configuration observable B. We evaluate instead

the signed longitudinal shell aver-

B = ρ̂ −k and B = δU . As shown in Fig. 3, both channels fluctuate about zero with no systematic drift over the accessible k range, and remain consistent with the exact We have formulated a single symmetry framework for equilibrium many-body systems by embedding standard geometric and dynamical invariances together with the phase-space gauge-shifting symmetry into the unified group G. Via Noether’s calculus for Liouville-preserving flows we derived a hierarchy of Ward identities and, crucially, new cross-Ward constraints that originate from the noncommutation of generators. The gauge–U(1) master identity fixes force–density cross correlators [Eqs. (31)–(35)]; in concert with rotations it enforces strict longitudinality and curl-freeness [Eqs. (38), (51)]; with dilations it yields a first-moment/compressibility sum rule [Eqs. (39), (40)]; with translations it implies a k-derivative identity [Eq. (44)]; with boundaries it produces a wall-contact relation [Eq. (42)]; and with boosts/time translations it links force–, current–, and stress–response channels, fixing the longitudinal stress from density dynamics and providing bounds [Eqs. (49)–(61)]. The Wigner–Eckart–Ward reduction further collapses tensor–hyperforce correlators to a small set of scalar radial spectra in isotropic fluids [Eqs. (67), (68)]. We also outlined an equivariant gauge-constrained DFT whose Euler–Lagrange equations automatically satisfy the full set of identities [Eqs. (77)–(81)].

The unified gauge–geometry perspective offers a systematic organizer for structure–mechanics relations in liquids, mixtures, and interfaces, and suggests practical routes to infer elastic/viscoelastic kernels from static or angle-resolved scattering. Extensions to multicomponent electrolytes [37], confined [38, 39] or anisotropic media [40] and quantum statistics [31, 41] are natural next steps.

Appendix A: Derivation of the nonzero commutators This Appendix derives Eqs. (22)–(28) by direct evaluation of Lie brackets of the phase-space vector fields that represent the generators. Throughout we work at fixed particle number N and write the phase-space point Γ = (r N ; p N ). All generators below are first-order differential operators acting on smooth phase-space functions F (Γ); for such operators the commutator [U, V ] ≡ U V − V U coincides with the Lie bracket of the associated vector fields.

13 1. 3. Conventions and one-particle reduction For each particle i and Cartesian component α =

Translation–gauge commutator Let T i (a) = a α ∂ iα and G i [ϵ] as in (A8). Using (A1) and noting that T i has constant coefficients, the r-component

equivalently [U, V ]F = U (V F ) − V (U F ) for any test function F (Γ). Here ∇ Γ denotes the full gradient with respect to all (r i , p i ).

All generators

Pused in Sec. II B are sums of one-particle pieces, U = i U i , where U i acts only on (r i , p i ). Hence [U i , V j ] = 0 for i ̸ = j, and therefore

Generators in components We use the phase-space representatives (Sec. II B):

14

5. Dilation–gauge commutator

is

The p-component follows similarly and collapses to the cotangent-lift form

Hence

(A15)

Boost–gauge commutator and the origin of the Π

term

which is the spatial Lie bracket (A9). A straightforward (but slightly longer) component calculation for the pcomponent yields

so the commutator retains the cotangent-lift structure. Summing over i gives

8. Momentum-fiber translation–gauge commutator

[ Π i [χ], G i [ϵ] ] vanishes because Π i differentiates only in