2512.02242

# Phase Transitions as Emergent Geometric Phenomena A Deterministic Entropy Evolution Law **Authors**: Loris Di Cairano > Department of Physics and Materials Science, University of Luxembourg, L-1511 Luxembourg City, Luxembourg ## Abstract We show that thermodynamics can be formulated naturally from the intrinsic geometry of phase space alone—without postulating an ensemble, which instead emerges from the geometric structure itself. Within this formulation, phase transitions are encoded in the geometry of constant-energy manifold: entropy and its derivatives follow from a deterministic equation whose source is built from curvature invariants. As energy increases, geometric transformations in energy-manifold structure drive thermodynamic responses and characterize criticality. We validate this framework through explicit analysis of paradigmatic systems—the 1D XY mean-field model and 2D $\phi^{4}$ theory—showing that geometric transformations in energy-manifold structure characterize criticality quantitatively. The framework applies universally to long-range interacting systems and in ensemble-inequivalence regimes. In physics, predictable theories are usually defined by a coherent triad: state variables, equations of motion, and intrinsic sources. State variables describe the system’s configuration; equations of motion govern their evolution. In classical and quantum mechanics, this structure is self-contained: forces in classical mechanics and Hamiltonian operators in quantum mechanics provide intrinsic, deterministic dynamics. Statistical mechanics occupies a singular position. Its scope is to bridge the deterministic microscopic world with macroscopic thermodynamics. Fundamental state variables—entropy, for instance—do not obey intrinsic evolution laws but rather condense the global properties of microscopic configurations. This is not a defect but the essential feature of thermodynamics. It stems naturally from the symplectic structure $(\Lambda,H,\omega)$ of phase space, which furnishes Hamilton’s equations and the Liouville measure. However, this structure alone does not determine which microstates are accessible at a given macroscopic configuration. Bridging this gap has required postulating statistical ensembles $\rho_{H}$ , making the theory predictive through integration rather than by solving the equations of motion. Although this asymmetry between statistical mechanics and other theories does not constitute a scientific crisis, it motivates a foundational question: Is the ensemble-based structure of statistical mechanics the only one that predicts experimentally verifiable thermodynamics? If not, can statistical mechanics admit a deterministic evolution law for entropy intrinsic to the theory itself? In this Letter, we show that thermodynamics emerges naturally from the intrinsic geometry of phase space. We formulate the theory axiomatically from $(\Lambda,H,\omega,\eta)$ , where the metric tensor is introduced alongside the Hamiltonian and the symplectic form. Crucially, the ensemble is not postulated but emerges necessarily from geometry alone. Uncovering a sort of thermodynamic covariance within this formalism, namely an equivalence class of metrics all inducing the same thermodynamic description, we show that the microcanonical measure coincides exactly with the area of constant-energy hypersurfaces. Then, we derive an exact deterministic evolution law for entropy—the Entropy Flow Equation (EFE)—whose source is built entirely from the curvature invariants of energy hypersurfaces. This geometric identification reveals that phase transitions (PTs) correspond to geometric changes in the energy hypersurfaces. The notion that geometry encodes thermodynamic information has a long history, explored through both geometric and topological perspectives. Rugh showed that inverse temperature is connected to Gauss curvature [1, 2, 3]. Pettini, Franzosi, and collaborators developed topological approaches linking topological changes of the energy-manifold to PTs [4, 5, 6, 7], stimulating an debate on the relation between topology and thermodynamic [8, 9, 10, 11, 12, 13, 14, 15, 16, 17] Later work established direct identifications between geometric properties and microcanonical observables [18, 19, 20, 21]. However, these contributions treat geometry as an additional structure imposed on a prescribed ensemble. Our work reveals that thermodynamics arises from phase-space geometry itself, without ensemble assumptions. The framework is intrinsically microcanonical, applying naturally to finite systems, long-range interactions, and ensemble-inequivalent regimes where conventional methods—Lee–Yang theory, renormalization group—may be ill-defined [22, 23, 24]. We solve the EFE exactly and validate it against independent thermodynamic methods on paradigmatic systems—the long-range 1D XY mean-field model and the 2D $\phi^{4}$ model—demonstrating that geometric transformations in energy-manifold curvature quantitatively characterize PTs. This geometric foundation unifies the treatment of criticality across disparate physical regimes. Classical statistical mechanics is grounded in the Hamiltonian quadruple $(\Lambda,H,\omega,\rho(H))$ . Here $\Lambda$ is the phase space with coordinates $\bm{x}=(\bm{p},\bm{q})$ and dimension $2n$ , $H$ is the Hamiltonian, $\omega=dp_{i}\wedge dq^{i}$ is the canonical symplectic form, and $\rho(H)$ is the statistical weight. We will focus on the microcanonical statistical ensemble. The dynamics is generated by the Hamiltonian vector field $\bm{X}_{H}$ via [25] $$ \iota_{X_{H}}\omega=dH, $$ so that trajectories lie on energy-manifold $\Sigma_{E}:=H^{-1}(E)$ and energy is conserved, $$ \dot{H}=dH(\bm{X}_{H})=\omega(\bm{X}_{H},\bm{X}_{H})=0. $$ While $\omega$ furnishes the Liouville measure $d\mu_{\Lambda}=\omega^{\wedge n}/n!=dq^{1}\!\cdots dq^{N}\,dp_{1}\!\cdots dp_{N}$ (for $\dim\Lambda=2n$ ), it does not provide a canonical measure on $\Sigma_{E}$ because of three structural obstructions; see Sec. A in Supplemental Materials (SM): (i) dimensional mismatch (energy shells have odd dimension $2n\!-\!1$ ); (ii) non-uniqueness of a transverse projection; (iii) coisotropy/degeneracy of $\omega|_{\Sigma_{E}}$ . In the standard framework, this gap is bridged by postulating an ensemble; in the microcanonical case, one assumes $$ d\mu_{\rho}=\rho(H)\,d\mu_{\Lambda},\qquad\rho(H)\propto\delta(H-E), $$ which makes the theory predictive by connecting the microscopic Hamiltonian flow to macroscopic thermodynamics. In this work, we take a step back and construct the bridge geometrically. Retaining $(\Lambda,H,\omega)$ , we introduce a natural metric tensor, the Euclidean one: $\eta=\delta^{ij}dp_{i}\!\otimes dp_{j}+\delta_{ij}dq^{i}\!\otimes dq^{j}$ on $\Lambda$ . The metric tensor allows the connection $\eta(\nabla^{\eta}H,Y)=dH(Y)$ , due to the induced isomorphism between covectors and vectors. As a consequence, the energy conservation law is converted into the orthogonality relation $$ \eta(\nabla^{\eta}H,\bm{X}_{H})=dH(\bm{X}_{H})=\omega(\bm{X}_{H},\bm{X}_{H})=0, $$ which identifies, at each $\bm{x}\in\Lambda$ , $$ T_{\bm{x}}\Lambda=\mathrm{span}\{\nabla_{\eta}H(\bm{x})\}\oplus T_{\bm{x}}\Sigma_{E}, $$ with $T_{\bm{x}}\Sigma_{E}=\{Y\in T_{\bm{x}}\Lambda:\eta(Y,\nabla^{\eta}H)=0\}$ . In other words, $\Lambda$ decomposes into two orthogonal directions: one is identified by $\bm{X}_{H}$ , which always lies on $\Sigma_{E}$ and generates the Hamiltonian flow; the second is identified by $\nabla_{\eta}H$ and generates the motion of points on an energy hypersurface to others. In general, $\nabla_{\eta}H$ moves points on the same $\Sigma_{E}$ to different shells (see Sec. B in SM). The only way to generate an energy motion that advances all points of $\Sigma_{E}$ to the same nearby shell $\Sigma_{E+\epsilon}$ consists of introducing an energy clock, i.e., a $\bm{\zeta}\in\mathrm{span}\{\nabla^{\eta}H\}$ such that [26] (see Sec. B in SM) $$ dH(\bm{\zeta})=1\quad\Longrightarrow\quad\bm{\zeta}=\frac{\nabla^{\eta}H}{\|\nabla^{\eta}H\|^{2}}. $$ This structure is not imposed by geometry but is clearly dictated by physics; that is, the presence of a Hamiltonian function which selects the right way to move from a value of energy to another. This vector field generates the thermodynamic dynamics: the flow $\Phi^{\mathrm{diff}}_{\epsilon}:\Sigma_{E}\to\Sigma_{E+\epsilon}$ defined by $$ \frac{d}{d\epsilon}\,\Phi^{\mathrm{diff}}_{\epsilon}(\bm{x})=\bm{\zeta}(\bm{x}(\epsilon)),\qquad H(\bm{x}(\epsilon))=E+\epsilon, $$ which yields diffeomorphisms between neighboring energy hypersurfaces. This motion naturally defines the adapted coordinates $(E,y^{\alpha})$ with $E=H(\bm{x})$ and $y^{\alpha}$ coordinates on $\Lambda$ . In these coordinates, the metric is block diagonal (see Sec. C in the SM for more details): $$ \eta=\frac{dE\otimes dE}{\|\nabla^{\eta}H\|_{\eta}^{2}}+(\sigma^{\eta}_{E})_{\alpha\beta}\,dy^{\alpha}\!\otimes dy^{\beta}, $$ where $\sigma^{\eta}_{E}$ is the metric induced on $\Sigma_{E}$ . The natural measure induced on $\Sigma_{E}$ is $$ \begin{split}\begin{split}d\mu^{\eta}_{E}&=d\mu_{\Lambda}\Big|_{\Sigma_{E}}:=\frac{d\sigma^{\eta}_{E}}{\|\nabla^{\eta}H\|_{\eta}},\end{split}\end{split} $$ with $d\sigma^{\eta}_{E}=\sqrt{\det\sigma^{\eta}_{E}}\,dy^{1}\!\cdots dy^{2n-1}$ and this yields the (geometric) density of states $$ \begin{split}\Omega_{\eta}(E):=\int_{\Sigma_{E}}d\mu_{E}^{\eta}=\int_{\Sigma_{E}}\frac{d\sigma^{\eta}_{E}}{\|\nabla^{\eta}H\|_{\eta}}=\int_{\Lambda}\delta\!\big(H(\bm{x})-E\big)\,d\mu_{\Lambda}.\end{split} $$ where the coarea formula has been used [27]. The last integral on the right coincides with the Boltzmann density of state $\Omega_{B}(E)$ . This result embodies a conceptual inversion of the standard statistical mechanics framework. Ordinarily, one postulates an ensemble $\rho_{H}$ first, and then—if desired—adds geometric structure and passes from right to left. Here, we reverse the logical order: the metric tensor is introduced axiomatically alongside the fundamental Hamiltonian and symplectic structure, and the microcanonical ensemble emerges necessarily as a purely geometric consequence. The metric supplies the missing transverse geometry, internalizing the ensemble not as an external postulate but as an intrinsic property of phase space itself. The only physical constraint imposed by the symplectic structure on the metric tensor is the conservation of Liouville measure, $\mathscr{L}_{\bm{X}_{H}}(d\mu_{\Lambda})=0$ , and then $\mathscr{L}_{\bm{X}_{H}}d\mu^{\eta}_{E}=0$ . This constraint fixes the form of the measure but not the metric tensor itself. Hence, there exists an equivalence class $[\eta]$ of metrics such that if $g$ is such that $$ d\mu^{g}_{E}=\frac{d\sigma^{g}_{E}}{\|\nabla^{g}H\|_{g}}\overset{\text{area}}{\simeq}\frac{d\sigma^{\eta}_{E}}{\|\nabla^{\eta}H\|_{\eta}}=d\mu^{\eta}_{E}\,, $$ then $g\in[\eta]$ and it is thermodynamically equivalent to $\eta$ . This is a geometric gauge symmetry of thermodynamics that we formulate as follows: Thermodynamic covariance. The thermodynamic content of a Hamiltonian system is independent of the geometric representation of phase space, provided that the microcanonical weight is preserved. The equivalence class $[\eta]=\{g:g\sim\eta\}$ consists of all metrics that induce the same microcanonical measure $\Omega_{\eta}$ . Therefore, all microcanonical observables, namely, $S(E)$ and its energy derivatives, depend only on the equivalence class $[\eta]$ . <details> <summary>x1.png Details</summary>  ### Visual Description ## Multi-Panel Scientific Figure: Critical Behavior Analysis ### Overview The image is a composite scientific figure containing two mathematical equations at the top and six plots arranged in a 3x2 grid below. The left column (labeled a.1, a.2, a.3) and the right column (labeled b.1, b.2, b.3) present related analyses for two different systems or parameter regimes. The plots show the behavior of various functions with respect to a parameter ε (epsilon). Key features include vertical dashed lines indicating critical points and annotations for specific values. ### Components/Axes **Equations (Top):** 1. **Left Equation:** `V_φψ := Σ_i [ (μ²/2)φ_i² + (λ/4!)φ_i⁴ ] + (J/2) Σ_{μ=1}^n (φ_i - φ_{i+e_μ})²` 2. **Right Equation:** `V_XY = (J/2N) Σ_{i,j} [1 - cos(θ_i - θ_j)]` **Plot Grid Structure:** * **X-axis (All plots):** Labeled `ε` (epsilon). * **Column a (Left):** Range approximately 9 to 14.5. * **Column b (Right):** Range approximately 0.55 to 1.0. * **Y-axes (By Row):** * **Row 1 (a.1, b.1):** Labeled `γ_g^(2)(ε)`. Scale is linear. * **Row 2 (a.2, b.2):** Labeled `∂_ε S(ε)`. Scale is linear. * **Row 3 (a.3, b.3):** Labeled `∂_ε^2 S(ε)`. Scale is linear. * **Legends & Annotations:** * **"micro average"**: Marked with a black 'x' symbol in plots a.1 and b.1. * **"PHT"**: Represented by a solid black circle (●). Appears in the legend of plots a.2, a.3, b.2, b.3. * **"EFE"**: Represented by a red open circle (○). Appears in the legend of plots a.2, a.3, b.2, b.3. * **Vertical Dashed Lines:** * **Column a:** A red vertical dashed line at `ε ≈ 11.1`. Annotated in plot a.2 as `ε_c^MIPA ≈ 11.1`. * **Column b:** A green vertical dashed line at `ε ≈ 0.74`. Annotated in plot b.2 as `ε_c^MIPA ≈ 0.74`. An additional green text annotation in plot b.2 states `ε_c^∞ = 3J/4`. * **Data Series:** Plots in rows 2 and 3 contain two data series: black dots (presumably for "PHT") and red circles (presumably for "EFE"). ### Detailed Analysis **Column a (Left Side):** * **a.1 (`γ_g^(2)(ε)`):** The curve (gray 'x' markers) shows a sharp peak at the critical point `ε ≈ 11.1`. The value rises from ~1.4x10⁻³ at ε=9 to a peak of ~2.0x10⁻³, then decays. * **a.2 (`∂_ε S(ε)`):** Both data series (black dots/PHT and red circles/EFE) show a nearly linear, decreasing trend. The values drop from ~6.2x10⁻² at ε=9 to ~5.0x10⁻² at ε=14.5. The two series are nearly indistinguishable. The critical point `ε_c^MIPA ≈ 11.1` is marked. * **a.3 (`∂_ε^2 S(ε)`):** The data shows a pronounced, sharp peak (a cusp) at the critical point `ε ≈ 11.1`. The value rises from ~-2.2x10⁻³ at ε=9 to a peak of ~-1.2x10⁻³, then decays back down. The black (PHT) and red (EFE) series overlap closely. **Column b (Right Side):** * **b.1 (`γ_g^(2)(ε)`):** The curve (gray 'x' markers) shows a discontinuity or very sharp change at the critical point `ε ≈ 0.74`. The value is high (~2.8) for ε < 0.74, drops abruptly, and then increases from ~-3.5 to ~-1.0 for ε > 0.74. * **b.2 (`∂_ε S(ε)`):** Both data series show a decreasing trend. The values drop from ~2.5 at ε=0.55 to ~1.0 at ε=1.0. The slope appears to change near the critical point `ε_c^MIPA ≈ 0.74`. The annotation `ε_c^∞ = 3J/4` is present. * **b.3 (`∂_ε^2 S(ε)`):** The data shows a sharp, downward cusp (a minimum) at the critical point `ε ≈ 0.74`. The value falls from ~-3.0 at ε=0.55 to a minimum of ~-7.0, then rises back to ~-2.0 at ε=1.0. The black (PHT) and red (EFE) series overlap closely. ### Key Observations 1. **Critical Points:** Both systems exhibit critical behavior at specific ε values (`~11.1` for column a, `~0.74` for column b), marked by vertical lines and singularities in the plotted functions. 2. **Nature of Singularity:** The singularity manifests differently. In column a, `γ_g^(2)` shows a peak and `∂_ε^2 S` shows a cusp. In column b, `γ_g^(2)` shows a discontinuity and `∂_ε^2 S` shows a sharp minimum. 3. **Data Series Agreement:** The "PHT" (black dots) and "EFE" (red circles) data series are in very close agreement across all plots where they appear (rows 2 and 3), suggesting two different methods or approximations yield nearly identical results for these quantities. 4. **Equation Context:** The equations at the top define potentials (`V_φψ` and `V_XY`) for what appear to be field theories (possibly φ⁴ theory and an XY model), providing the theoretical context for the analyzed data. ### Interpretation This figure presents a comparative analysis of critical phenomena in two different statistical mechanical or field-theoretic models. The left column (a) likely corresponds to a system with a continuous (second-order) phase transition, evidenced by the continuous but non-differentiable cusp in the second derivative of entropy (`∂_ε^2 S`). The right column (b) suggests a system with a more abrupt transition, possibly first-order or of a different universality class, indicated by the apparent discontinuity in `γ_g^(2)`. The parameter ε is the control parameter driving the system through its critical point (`ε_c`). The functions plotted—`γ_g^(2)` (likely a correlation function or susceptibility), the first derivative of entropy `∂_ε S` (related to specific heat or a response function), and its second derivative `∂_ε^2 S`—are standard quantities for characterizing phase transitions. The sharp features at `ε_c` are hallmarks of critical behavior, where correlation lengths diverge and the system becomes scale-invariant. The close match between "PHT" and "EFE" data suggests the figure is validating a theoretical or computational method (perhaps "EFE" stands for an Effective Field Theory approach) against another standard ("PHT"). The annotation `ε_c^∞ = 3J/4` in plot b.2 provides an analytical prediction for the critical point in a certain limit (∞), which the numerical/data-driven `ε_c^MIPA ≈ 0.74` approximates. The overall purpose is to demonstrate and quantify the critical behavior of these models, locate their critical points, and compare different calculational approaches. </details> Figure 1: Comparison of entropy derivatives obtained as the solution of the entropy flow equation and through thermodynamic methods. Empty red circles are associated with the solutions (derivatives of entropy) of the Entropy Flow Equation (EFE). Black disks are the numerical estimations of entropy derivatives obtained from microcanonical simulations with Pearson-Halicioglu-Tiller method [28]. The geometric function $\Upsilon^{(2)}_{g}$ represented by black crosses is computed through numerical simulation according with Eq. (S-27) in SM. In panels a.1-a.3, we report results for the 2D $\phi^{4}$ -model with nearest interactions and coupling parameters $\lambda=3/5$ , $\mu=\sqrt{2}$ and $J=1$ and size $N=256^{2}$ . The vertical red dashed line indicates the critical energy detected by MIPA and in agreement with the literature. Panels b.1-b.3 are associated with the comparison of entropy derivatives in the 1D XY mean-field model with coupling parameter $J=1$ and size $N=40000$ . Notice that the infinite-size critical energy $\epsilon^{\infty}_{c}=3J/4$ is closed to the finite-size critical energy $\epsilon_{c}^{\text{MIPA}}=0.74$ estimated by MIPA. See SM for more details. We are then allowed to choose the most convenient element $g$ within the class $[\eta]$ . Such a selection, corresponding to fix the gauge, follows by the request that the generator $\bm{\xi}$ of the diffeomorphisms $\Phi^{\rm diff}_{\bm{\xi}}:\Sigma_{E}\mapsto\Sigma_{E+\delta E}$ has unit length, $\|\bm{\xi}\|_{g}=1$ . A straightforward calculation (see Sec. D in SM) shows that the unit-length implies $\bm{\xi}=\nabla^{g}H$ and $\text{area}_{g}=\Omega_{\eta}$ . We call this metric the unit-normal gauge. This leads to the following results: Geometric Microcanonical Equivalence. Within the thermodynamically equivalent metric class, it is possible to find a geometric representation $g$ such that counting accessible states is equivalent to measuring the occupied geometric area [29]. Entropy follows as an intrinsic geometric property: $$ S_{g}(E):=\ln\mathrm{area}^{g}(E)=\ln\Omega_{B}(E)=S_{B}(E)\,, \tag{1} $$ This result suggests that geometry itself dictates the full thermodynamic behavior. Intuitively, the thermodynamic properties of a physical system and the emergence of PTs are encoded within the geometry of the energy hypersurfaces and detected through the energy flow $\Phi^{\rm diff}_{\bm{\xi}}$ . We show how an evolution flow with an intrinsic geometric source naturally emerges. Differentiating Eq. (1) twice: $$ \displaystyle\partial_{E}S_{g}(E) \displaystyle=\frac{\partial_{E}\text{area}^{g}(E)}{\text{area}^{g}(E)}, \displaystyle\partial_{E}^{2}S_{g}(E) \displaystyle=\frac{\partial_{E}^{2}\text{area}^{g}(E)}{\text{area}^{g}(E)}-\left(\frac{\partial_{E}\text{area}^{g}(E)}{\text{area}^{g}(E)}\right)^{2}. \tag{2} $$ A recursive structure naturally appears, $\Upsilon^{(k)}(E)=\partial_{E}^{k}\text{area}^{g}(E)/\text{area}^{g}(E)$ , whose members can be expressed in terms of geometric invariants. The first variation of volume involves the trace of the Weingarten operator since $\mathscr{L}_{\bm{\xi}}d\mu_{g}=(\mathrm{div}_{g}\bm{\xi})\,d\mu_{g}$ with $$ \mathrm{div}_{g}\bm{\xi}=\mathrm{Tr}^{g}[W_{g}]=\Delta_{g}H\,. $$ Note that this can be written in terms of the Laplacian, gradient, and Hessian associated with the metric $\eta$ , and it reads $$ \begin{split}\text{Tr}^{g}[W_{g}]=\frac{\Delta_{\eta}H}{\|\nabla_{\eta}H\|^{2}_{\eta}}-2\frac{\langle\nabla_{\eta}H,{\rm Hess}_{\eta}\,H\,\nabla_{\eta}H\rangle_{\eta}}{\|\nabla_{\eta}H\|^{4}_{\eta}}\;.\end{split} \tag{4} $$ Therefore: $$ \displaystyle\partial_{E}S_{g}(E)\equiv\Upsilon^{(1)}(E)=\int_{\Sigma_{E}}\frac{\mathscr{L}_{\bm{\xi}}(d\mu_{g})}{\text{area}^{g}(E)}=\int_{\Sigma_{E}}\text{Tr}^{g}[W_{g}]~d\rho_{g}\,, \tag{1} $$ where $d\rho_{g}:=d\mu^{g}_{E}/\text{area}^{g}(E)$ . The microcanonical inverse temperature $\beta(E)\equiv\partial_{E}S_{g}(E)$ coincides with the geometric measure of the Weingarten operator, i.e., $$ \beta(E)=\langle\text{Tr}^{g}W_{g}\rangle_{\Sigma_{E}}\,. $$ This is the functional found by Rugh [1, 3], expressed here in a fully geometric way and already established in Ref. [19]. Taking the second variation, we have $\mathscr{L}_{\bm{\xi}}\!\left(\text{Tr}^{g}[W_{g}]\,d\mu_{g}\right)=\text{Tr}^{g}[\nabla_{\bm{\xi}}W_{g}]+\text{Tr}^{g}[W_{g}]^{2}d\mu_{g}$ . Then, using [30]: $$ \text{Tr}^{g}[\nabla_{\bm{\xi}}W_{g}]=-\text{Tr}^{g}[W_{g}^{2}]-Ric(\bm{\xi},\bm{\xi})\,, $$ where $Ric$ is the Ricci curvature along $\bm{\xi}$ and the Gauss-Codazzi equation $2\,Ric(\bm{\xi},\bm{\xi})=R^{g}_{\Lambda}-R^{g}_{\Sigma_{E}}-\text{Tr}^{g}[W_{g}^{2}]+\text{Tr}^{g}[W_{g}]^{2}$ , we obtain $$ \Upsilon^{(2)}_{g}(E)=\frac{1}{2}\!\int_{\Sigma_{E}}\!\!\Big\{\text{Tr}^{g}[W_{g}]^{2}-\text{Tr}^{g}[W_{g}^{2}]+R^{g}_{\Sigma_{E}}-R^{g}_{\Lambda}\Big\}d\rho_{g}, \tag{2} $$ where $R^{g}_{\Sigma_{E}}$ and $R^{g}_{\Lambda}$ are, respectively, the scalar curvatures of $\Sigma_{E}$ and of the phase space [31]. The explicit expression for $\Upsilon^{(2)}_{g}$ computed in simulations is given in Eq. (S-27) in the SM. Combining Eqs. (2)–(3), we obtain $$ \partial_{E}^{2}S_{g}(E)+[\partial_{E}S_{g}(E)]^{2}=\Upsilon^{(2)}_{g}(E), \tag{2} $$ a Riccati-type differential law for the entropy, hereafter called the Entropy Flow Equation (EFE). Equation (7) reveals that entropy is governed by a deterministic geometric flow: once the initial conditions are specified, its evolution with energy is uniquely determined by the curvature invariants of $\Sigma_{E}$ . This marks a change in perspective—entropy, traditionally introduced as a statistical construct, becomes a dynamical variable whose evolution is dictated by geometry itself. Variations of $\Upsilon^{(2)}_{g}$ , which encode changes in the curvature of $\Sigma_{E}$ , translate directly into thermodynamic responses. In this sense, geometry encodes thermodynamics. To validate the predictive power of the EFE, we applied it to two paradigmatic models: the 2D $\phi^{4}$ model with nearest-neighbor interactions and the long-range 1D XY mean-field model. All technical details about simulations are reported in Sec. H of the SM. In both cases, we evaluate the geometric function $\Upsilon^{(2)}(E)$ (Eq. (S-27) in the SM) through microcanonical simulations over a dense grid of energies, as reported in panels a.1 and b.1 of Fig. S2). The resulting $\Upsilon^{(2)}(E)$ was then inserted into Eq. (7), which we solved numerically using a high-order Runge–Kutta integrator, yielding the entropy derivatives $\partial_{E}S_{g}$ and $\partial_{E}^{2}S_{g}$ . These observables are compared with those obtained through the Pearson-Halicioglu-Tiller (PHT) method [28], which is a fully thermodynamic approach and geometry-independent. In Fig. S2, panels a.2-a.3 for the $\phi^{4}$ model and b.2-b.3 for the XY mean-field model, we report the comparison between geometric and thermodynamic entropy derivatives. The excellent quantitative agreement between them demonstrates that the EFE faithfully reproduces both the qualitative and quantitative energy dependence of entropy and its derivatives. Geometry contains information about thermodynamics and PTs. Both systems undergo a second-order PTs, as is well-known by the literature. Here, we show that the microcanonical inflection point analysis (MIPA) [32] is able to coherently detect the PT even at a finite-size. A deeper analysis with increasing system sizes, especially for the 1D XY mean-field model, is reported in Ref. [52]. For the $\phi^{4}$ model, $\partial_{E}S_{g}$ displays an inflection point (panel a.3 and $\partial_{E}^{2}S_{g}$ a negative-valued maximum at $\epsilon_{c}\simeq 11.1$ (panel a.2, in agreement with MIPA for a second-order PT. The critical point is also in agreement with the literature [11, 8, 34]. For the XY mean-field model, the same signatures appear at $\epsilon_{c}^{\mathrm{MIPA}}\!\simeq\!0.73$ ; see panels b.3 and b.2. This value corresponds to the finite-size critical energy, while in the thermodynamic limit, we have $\epsilon_{c}^{\infty}\!=\!3J/4$ [35, 36]. We will show below, through an analytical approach within the geometric framework, how to determine the infinite-size critical energy. Notably, the curvature of $\partial_{E}^{2}S_{g}$ changes more sharply in the XY case, reflecting the stronger collective character of its transition. Apart from the exact, that is, the non-perturbative solution provided by solving the EFE (7), it is crucial to determine the microscopic geometric mechanism that entails a PT. To this scope, we perform an explicit analytical expansion of $\mathrm{Tr}\,W_{g}$ in the paradigmatic 1D XY mean-field model, whose Hamiltonian is $$ H_{\text{XY}}(\theta,p)=\sum_{i=1}^{N}\frac{p_{i}^{2}}{2}+\frac{J}{2N}\sum_{i,j}\bigl[1-\cos(\theta_{i}-\theta_{j})\bigr]. \tag{8} $$ We remark that $\mathrm{Tr}\,W_{g}$ plays a central thermodynamic role, being related to the inverse temperature. The strategy consists of expanding $\mathrm{Tr}\,W_{g}$ in powers of $M$ , i.e., $$ \mathrm{Tr}[W_{g}]=C^{\infty}_{0}(\epsilon)+C^{\infty}_{2}(\epsilon)\,M^{2}+O(M^{4}), $$ where $M$ is the magnetization (order parameter), defined by $\bm{M}=N^{-1}\sum_{i}(\cos\theta_{i},\sin\theta_{i})=(M\cos\phi,M\sin\phi)$ , and where $C_{0}$ and $C_{2}$ are $\epsilon$ -dependent coefficients that survive in the thermodynamic limit. We will find that the critical point $\epsilon_{c}$ corresponds to the energy density at which the concavity of the energy manifolds changes sign. All details are reported in the SM. We begin by noting that $\mathrm{Tr}\,W_{g}$ in Eq. (4) can be exactly decomposed into three explicit blocks associated with momenta, $\kappa_{p}$ , diagonal angular curvature, $\kappa^{\rm d}$ , and angular interactions, $\kappa^{\rm int}$ : $$ \text{Tr}[W_{g}]=\sum_{i}\bigg(\kappa_{p_{i}}+\kappa^{\text{d}}_{\theta_{i}}+\sum_{j\neq i}\kappa^{\mathrm{int}}_{ij}\bigg) \tag{9} $$ and, setting $\Theta_{i}=\theta_{i}-\phi$ , they read $$ \displaystyle K:=\sum_{i=1}^{N}p_{i}^{2},\qquad\qquad G:=\|\nabla H\|^{2}=K+J^{2}M^{2}S_{2} \displaystyle\sum_{i}\kappa_{p_{i}}=\frac{N}{G}-\frac{2K}{G^{2}},\qquad\qquad\qquad S_{2}=\sum_{i}\sin^{2}\Theta_{i} \displaystyle\sum_{i}\kappa^{\text{d}}_{\theta_{i}}=\frac{JNM^{2}}{G}\!-\!\frac{2J^{3}M^{3}S_{21}}{G^{2}},\quad S_{21}=\sum_{i=1}^{N}\sin^{2}\Theta_{i}\,\cos\Theta_{i}, \displaystyle\sum_{i\neq j}\kappa^{\mathrm{int}}_{ij}=\frac{2J^{3}M^{2}}{N\,G^{2}}S_{\rm int},\;\;\;S_{\mathrm{int}}=\sum_{i\neq j}\sin\Theta_{i}\sin\Theta_{j}\cos(\theta_{i}-\theta_{j}). \tag{10} $$ The next step is to expand each building block above up to order $O(M^{2})$ and retain only terms of order $O(1)$ in $N$ , since they are the only ones that survive in the thermodynamic limit. The angular functions are written as follows (see Sec. F in SM): $$ \displaystyle S_{2} \displaystyle\simeq N\Big(\frac{1}{2}-\frac{1}{4}M^{2}\Big),\quad S_{21}\simeq\frac{MN}{4},\quad S_{\mathrm{int}}\simeq\frac{N^{2}}{4}. \tag{14} $$ Then, introducing the (twice) kinetic density $c:=K/N$ , the norm squared in Eq. (10) becomes $G\simeq N(c+J^{2}M^{2}/2)$ . This expansion is then plugged into the denominator of $\kappa_{p_{i}}$ , $\kappa^{\rm d}_{\theta_{i}}$ , and $\kappa^{\rm int}_{ij}$ , together with Eq. (14). Then, we use the binomial series $(1+x)^{\alpha}\simeq 1+\alpha x+O(x^{2})$ to expand $G$ in the denominators. In $\kappa_{p_{i}}$ , only the first term $N/G$ is leading in $M^{2}$ and of order $O(1)$ in $N$ ; it reads $$ \begin{split}\sum_{i}\kappa_{p_{i}}\simeq\frac{1}{c}-\frac{J^{2}}{2c^{2}}M^{2}+O(M^{4})\,,\end{split} $$ since $2K/G^{2}\simeq 2(1-2J^{2}M^{2}/c+O(M^{4}))/Nc$ is therefore negligible in the limit $N\to\infty$ , as it is of order $O\left(1/N\right)\left(1+O(M^{2})\right)$ . Similarly, in $\kappa^{\rm d}_{\theta_{i}}$ , only the first term, namely, $JNM^{2}/G$ , survives, and it reads $$ \sum_{i}\kappa^{\rm d}_{\theta_{i}}\simeq\frac{JM^{2}}{c}\;. $$ In fact, the second term gives $2J^{3}M^{3}S_{12}/G^{2}=J^{3}M^{4}/2c^{2}N+O(M^{6}/N)$ since $S_{21}=O(NM)$ and is thus negligible at order $M^{2}$ in the thermodynamic limit. Finally, the $\kappa^{\rm int}_{ij}$ -term is entirely sub-leading at order $M^{2}$ and vanishes in the large– $N$ limit; indeed $\sum_{i}\kappa^{\rm int}_{ij}\simeq J^{3}M^{2}/2Nc^{2}\simeq O(M^{2}/N)$ . Collecting all terms, the trace reads: $$ \mathrm{Tr}[W_{g}]^{\infty}\approx\frac{1}{c}+\Big(\frac{J}{c}-\frac{J^{2}}{2c^{2}}\Big)\,M^{2}, \tag{15} $$ and expressing it in terms of the energy density $\varepsilon$ yields $$ C^{\infty}_{2}(\varepsilon)=\frac{J}{2\varepsilon-J}-\frac{J^{2}}{2(2\varepsilon-J)^{2}}=\frac{2J}{(2\varepsilon-J)^{2}}\!\left(\varepsilon-\frac{3J}{4}\right)\!. \tag{16} $$ This analytical result makes the geometric nature of the transition explicit. The change of sign of $\mathrm{Tr}\,W_{g}$ at $\varepsilon_{c}=3J/4$ corresponds to a reversal of curvature in the energy manifolds—an instability of the $M=0$ branch where the geometry itself selects the ordered phase. The agreement with the known thermodynamic critical point confirms that the curvature of the energy manifolds encodes the full thermodynamic response. The 1D XY mean-field model, therefore, provides a transparent illustration of how a PT emerges as a purely geometric change. In summary, we have shown that thermodynamics can be formulated in a fully geometric fashion in terms of $(\Lambda,H,\omega,\eta)$ without postulating a priori the concept of (microcanonical) ensemble; rather, the latter arises naturally, and entropy becomes a geometric quantity whose evolution is governed by the deterministic entropy flow equation. This equation connects the derivatives of entropy to the curvature invariants of the energy manifold, thus providing a self-contained relation between dynamics and thermodynamics. The explicit analysis of the 1D XY mean-field model demonstrates that PTs correspond to geometric changes in the curvature of the energy manifolds—which reproduce known critical points. Combining the EFE with Bachmann’s microcanonical inflection-point analysis (MIPA) [32] we can establish a direct geometric foundation for PTs: each critical behavior classified by MIPA corresponds to a specific geometric profile of the source term $\Upsilon^{(2)}_{g}$ and recovers PT in the thermodynamic limit. The EFE thus provides a deterministic bridge between geometry and criticality, revealing the geometric mechanisms that generate PTs. Beyond its conceptual implications, the present framework has a universal scope and opens several directions for future investigation. It applies naturally to finite systems such as proteins [37, 38, 39, 40], where the microcanonical ensemble provides the most direct thermodynamic description [24]. It also offers tools and insights for addressing long-range interacting models [22, 24, 36, 41, 23, 42, 43], in which ensemble equivalence may break down and canonical methods or renormalization-group approaches may be limited. In such cases, the geometric formulation provides a natural language to investigate the origin and structure of PTs. When extended to include discrete configuration spaces, the framework can address spin systems where ensemble inequivalence also emerges [44, 45] and connects naturally with large-deviation formulations of non-equivalence [46, 47]. Acknowledgements. The calculations presented in this paper were carried out using the HPC facilities of the University of Luxembourg [48] (see hpc.uni.lu). ## References - [1] H. H. Rugh, “Dynamical approach to temperature,” Phys. Rev. Lett. 78, 772 (1997). - [2] H. H. Rugh, “A geometric, dynamical approach to thermodynamics,” J. Phys. A: Math. Gen. 31, 7761 (1998). - [3] H. H. Rugh, “Microthermodynamic formalism,” Phys. Rev. E 64, 055101 (2001). - [4] R. Franzosi, M. Pettini, and L. Spinelli, “Topology and phase transitions: Paradigmatic evidence,” Phys. Rev. Lett. 84, 2774 (2000). - [5] R. Franzosi and M. Pettini, “Topology and phase transitions II. Theorem on a necessary relation,” Nucl. Phys. B 782, 219 (2007). - [6] L. Casetti, M. Pettini, and E. Cohen, “Geometric approach to Hamiltonian dynamics and statistical mechanics,” Phys. Rep. 337, 237 (2000). - [7] R. Franzosi, L. Casetti, L. Spinelli, and M. Pettini, “Topological aspects of geometrical signatures of phase transitions,” Phys. Rev. E 60, R5009 (1999). - [8] M. Kastner and D. Mehta, “Phase transitions detached from stationary points of the energy landscape,” Phys. Rev. Lett. 107, 160602 (2011). - [9] L. Casetti, M. Kastner, and R. Nerattini, “Kinetic energy and microcanonical nonanalyticities in finite and infinite systems,” J. Stat. Mech. 2009, P07036 (2009). - [10] L. Casetti and M. Kastner, “Nonanalyticities of entropy functions of finite and infinite systems,” Phys. Rev. Lett. 97, 100602 (2006). - [11] D. Mehta, J. D. Hauenstein, and M. Kastner, “Energy-landscape analysis of the two-dimensional nearest-neighbor $\phi^{4}$ model,” Phys. Rev. E 85, 061103 (2012). - [12] F. Baroni, “The simplified energy landscape of the $\phi^{4}$ model and the phase transition,” J. Stat. Mech.: Theory Exp. 2024, 073201 (2024). - [13] L. Angelani, L. Casetti, M. Pettini, G. Ruocco, and F. Zamponi, “Topological signature of first-order phase transitions in a mean-field model,” Europhys. Lett. 62, 775 (2003). - [14] L. Angelani, L. Casetti, M. Pettini, G. Ruocco, and F. Zamponi, “Topology and phase transitions: From an exactly solvable model to a relation between topology and thermodynamics,” Phys. Rev. E 71, 036152 (2005). - [15] F. Baroni and L. Casetti, “Topological conditions for discrete symmetry breaking and phase transitions,” J. Phys. A: Math. Gen. 39, 529 (2005). - [16] M. Kastner, “Phase transitions and configuration space topology,” Rev. Mod. Phys. 80, 167 (2008). - [17] S. Risau-Gusman, A. C. Ribeiro-Teixeira, and D. A. Stariolo, “Topology, phase transitions, and the spherical model,” Phys. Rev. Lett. 95, 145702 (2005). - [18] M. Gori, R. Franzosi, and M. Pettini, “Topological origin of phase transitions in the absence of critical points of the energy landscape,” J. Stat. Mech. Theory Exp. 2018, 093204 (2018). - [19] M. Gori, “Configurational microcanonical statistical mechanics from Riemannian geometry of equipotential level sets,” arXiv:2205.14536 (2022). - [20] L. Di Cairano, “The geometric theory of phase transitions,” J. Phys. A: Math. Theor. 55, 27LT01 (2022). - [21] L. Di Cairano, M. Gori, and M. Pettini, “Topology and phase transitions: A first analytical step towards the definition of sufficient conditions,” Entropy 23, 1414 (2021). - [22] A. Campa, T. Dauxois, and S. Ruffo, “Statistical mechanics and dynamics of solvable models with long-range interactions,” Phys. Rep. 480, 57 (2009). - [23] F. Bouchet, S. Gupta, and D. Mukamel, “Thermodynamics and dynamics of systems with long-range interactions,” Physica A 389, 4389 (2010). - [24] J. Dunkel and S. Hilbert, “Phase transitions in small systems: Microcanonical vs. canonical ensembles,” Physica A 370, 390 (2006). - [25] V. I. Arnol’d, Mathematical Methods of Classical Mechanics, Vol. 60 (Springer Science & Business Media, 2013). - [26] M. W. Hirsch, Differential Topology, Vol. 33 (Springer Science & Business Media, 2012). - [27] H. Federer, Geometric Measure Theory (Springer, 2014). - [28] E. Pearson, T. Halicioglu, and W. Tiller, “Laplace-transform technique for deriving thermodynamic equations from the classical microcanonical ensemble,” Phys. Rev. A 32, 3030 (1985). - [29] M. Gori, R. Franzosi, G. Pettini, and M. Pettini, “Topological theory of phase transitions,” J. Phys. A: Math. 55, 375002 (2022). - [30] M. Gromov, “Four lectures on scalar curvature,” arXiv:1908.10612 (2019). - [31] P. Petersen, Riemannian Geometry, Vol. 171 (Springer, 2006). - [32] K. Qi and M. Bachmann, “Classification of phase transitions by microcanonical inflection-point analysis,” Phys. Rev. Lett. 120, 180601 (2018). - [33] L. Di Cairano, “A geometric foundation of thermodynamics: from statistical ensemble to entropy flow equation,” in preparation. - [34] G. Bel-Hadj-Aissa, M. Gori, V. Penna, G. Pettini, and R. Franzosi, “Geometrical aspects in the analysis of microcanonical phase transitions,” Entropy 22, 380 (2020). - [35] M. Antoni and S. Ruffo, “Clustering and relaxation in Hamiltonian long-range dynamics,” Phys. Rev. E 52, 2361 (1995). - [36] A. Campa, T. Dauxois, D. Fanelli, and S. Ruffo, Physics of Long-Range Interacting Systems (OUP Oxford, 2014). - [37] T. Koci and M. Bachmann, “Subphase transitions in first-order aggregation processes,” Phys. Rev. E 95, 032502 (2017). - [38] S. Schnabel, D. T. Seaton, D. P. Landau, and M. Bachmann, “Microcanonical entropy inflection points: Key to systematic understanding of transitions in finite systems,” Phys. Rev. E 84, 011127 (2011). - [39] M. Bachmann, Thermodynamics and Statistical Mechanics of Macromolecular Systems (Cambridge University Press, 2014). - [40] D. Aierken and M. Bachmann, “Comparison of conformational phase behavior for flexible and semiflexible polymers,” Polymers 12, 3013 (2020). - [41] T. Dauxois, S. Ruffo, E. Arimondo, and M. Wilkens, “Dynamics and thermodynamics of systems with long-range interactions: An introduction,” in Dynamics and Thermodynamics of Systems with Long-Range Interactions (Springer, 2002), pp. 1–19. - [42] N. Defenu, T. Donner, T. Macrì, G. Pagano, S. Ruffo, and A. Trombettoni, “Long-range interacting quantum systems,” Rev. Mod. Phys. 95, 035002 (2023). - [43] G. Giachetti, N. Defenu, S. Ruffo, and A. Trombettoni, “Villain model with long-range couplings,” J. Stat. Mech.: Theory Exp. 2023, 033201 (2023), see also arXiv:2209.11810. - [44] J. Barré, D. Mukamel, and S. Ruffo, “Inequivalence of ensembles in a system with long-range interactions,” Phys. Rev. Lett. 87, 030601 (2001). - [45] T. Dauxois, P. Holdsworth, and S. Ruffo, “Violation of ensemble equivalence in the antiferromagnetic mean-field XY model,” Eur. Phys. J. B 16, 659 (2000). - [46] R. S. Ellis, K. Haven, and B. Turkington, “Nonequivalent statistical equilibrium ensembles and refined stability theorems for most probable flows,” Nonlinearity 15, 239 (2002). - [47] R. S. Ellis, K. Haven, and B. Turkington, “Large deviation principles and complete equivalence and nonequivalence results for pure and mixed ensembles,” J. Stat. Phys. 101, 999 (2000). - [48] S. Varrette, P. Bouvry, H. Cartiaux, and F. Georgatos, “Management of an academic HPC cluster: The UL experience,” in Proc. of the 2014 Intl. Conf. on High Performance Computing & Simulation (HPCS 2014) (IEEE, Bologna, Italy, 2014), pp. 959–967. - [49] R. Abraham and J. E. Marsden, Foundations of Mechanics, 2nd ed. (Addison–Wesley, 1978). - [50] D. McDuff and D. Salamon, Introduction to Symplectic Topology, 3rd ed. (Oxford University Press, 2017). - [51] J. M. Lee, Introduction to Smooth Manifolds, 2nd ed. (Springer, 2013). - [52] L. Di Cairano, A Geometric Foundation of Thermodynamics: from statistical ensemble to Entropy Flow Equation, in preparation. - [53] K. V. Mardia and P. E. Jupp, Directional Statistics (Wiley, 2000). - [54] N. I. Fisher, Statistical Analysis of Circular Data (Cambridge University Press, 1993). - [55] S. R. Jammalamadaka and A. SenGupta, Topics in Circular Statistics (World Scientific, 2001). - [56] M. Abramowitz and I. A. Stegun (eds.), Handbook of Mathematical Functions (Dover, 1965). - [57] J. R. Ray, “Microcanonical ensemble Monte Carlo method,” Phys. Rev. A 44, 4061 (1991). - [58] J. R. Ray and C. Freléchoz, “Microcanonical ensemble Monte Carlo method for discrete systems,” Phys. Rev. E 53, 3402 (1996). Supplemental Material Phase Transitions as Emergent Geometric Phenomena A Deterministic Entropy Evolution Law ## A Obstructions In this section, we show and discuss the obstruction of the symplectic form in naturally inducing a measure on the energy hypersurfaces. For further details about proofs and tools used in this section, we refer to Refs. [25, 49, 51, 50]. ### A.1 Degeneracy of the symplectic form on the energy hypersurfaces We show that the restriction of the symplectic form $\omega$ to the tangent bundle of $\Sigma_{E}$ is degenerate, with the Hamiltonian vector field $\bm{X}_{H}$ in its kernel. Let ${\bm{x}}\in\Sigma_{E}$ and consider the restriction $\omega|_{T_{\bm{x}}\Sigma_{E}}$ , which we denote by $\omega_{E}$ . For any vector $v\in T_{x}\Sigma_{E}$ (that is, $\bm{v}$ tangent to the level set $H=E$ ), we have, by definition: $$ dH(\bm{v})=0. $$ By the defining equation of the Hamiltonian vector field, $\iota_{\bm{X}_{H}}\omega=dH$ , we have for any vector $v$ : $$ \omega(\bm{X}_{H},\bm{v})=dH(\bm{v}). $$ In particular, for $v\in T_{x}\Sigma_{E}$ : $$ \omega_{E}(\bm{X}_{H},\bm{v})=\omega(\bm{X}_{H},\bm{v})=dH(\bm{v})=0. $$ Since this holds for all $\bm{v}\in T_{\bm{x}}\Sigma_{E}$ , we have $\bm{X}_{H}\in\ker(\omega_{E})$ . The kernel is non-trivial; hence, $\omega_{E}$ is degenerate. Moreover, $X_{H}(x)\in T_{x}\Sigma_{E}$ because $$ dH(\bm{X}_{H})=\omega(\bm{X}_{H},\bm{X}_{H})=0 $$ by the antisymmetry of $\omega$ . Thus $X_{H}$ is both tangent to $\Sigma_{E}$ and in the kernel of $\omega$ restricted to $\Sigma_{E}$ . ### A.2 The Symplectic Complement and Coisotropic Submanifolds Here, we show that the symplectic complement of $T_{\bm{x}}\Sigma_{E}$ is one-dimensional and is generated by $\bm{X}_{H}({\bm{x}})$ . Furthermore, $(T_{\bm{x}}\Sigma_{E})^{\perp_{\omega}}\subseteq T_{\bm{x}}\Sigma_{E}$ , making $\Sigma_{E}$ a coisotropic submanifold. Consider the symplectic complement, which is defined as $$ (T_{\bm{x}}\Sigma_{E})^{\perp_{\omega}}=\{\bm{w}\in T_{\bm{x}}\Lambda:\omega(\bm{w},\bm{v})=0\;\forall\bm{v}\in T_{\bm{x}}\Sigma_{E}\}. $$ We first show that $\bm{X}_{H}\in(T_{\bm{x}}\Sigma_{E})^{\perp_{\omega}}$ . For any $\bm{v}\in T_{\bm{x}}\Sigma_{E}$ : $$ \omega(\bm{X}_{H},\bm{v})=dH(\bm{v})=0, $$ as shown in Sec. A.1. Thus $\bm{X}_{H}\in(T_{\bm{x}}\Sigma_{E})^{\perp_{\omega}}$ . To show that $\bm{X}_{H}$ generates the entire complement, we use a dimensional argument. For a subspace $W$ of a symplectic vector space $(V,\omega)$ of dimension $2n$ , the symplectic complement satisfies $$ \dim(W)+\dim(W^{\perp_{\omega}})=2n. $$ In our case, $\dim(T_{\bm{x}}\Sigma_{E})=2n-1$ (since $\Sigma_{E}$ is a codimension-1 hypersurface), therefore: $$ \dim((T_{\bm{x}}\Sigma_{E})^{\perp_{\omega}})=2n-(2n-1)=1. $$ Since $\bm{X}_{H}$ is non-zero (assuming $H$ is a regular function, i.e., $dH\neq 0$ on $\Sigma_{E}$ ) and belongs to $(T_{\bm{x}}\Sigma_{E})^{\perp_{\omega}}$ , and this complement has dimension $1$ , we conclude: $$ (T_{\bm{x}}\Sigma_{E})^{\perp_{\omega}}=\mathrm{span}\{X_{H}\}. $$ Finally, we verify that $(T_{x}\Sigma_{E})^{\perp_{\omega}}\subseteq T_{x}\Sigma_{E}$ . We have already shown that $X_{H}\in T_{x}\Sigma_{E}$ (since $dH(X_{H})=0$ ). Thus, the symplectic complement is contained in the tangent space itself. This is precisely the defining property of a coisotropic submanifold [50]. ### A.3 Non-uniqueness of Transverse Vectors Finally, we show that the condition $dH(\bm{\zeta})=1$ does not uniquely determine a vector field $\bm{\zeta}$ transverse to $\Sigma_{E}$ . The solution space has dimension $2n-1$ at each point. At each point $\bm{x}\in\Lambda$ , we seek vectors $\bm{\zeta}\in T_{\bm{x}}\Lambda$ satisfying the linear equation $$ dH(\bm{\zeta})=1. $$ This is a single linear equation in a $2n$ -dimensional vector space. The solution set is an affine hyperplane: $$ S=\{\bm{\zeta}\in T_{\bm{x}}\Lambda:dH(\bm{\zeta})=1\}. $$ To find the dimension, we first note that the solutions to the homogeneous equation $dH(\bm{\zeta})=0$ form the tangent space $T_{\bm{x}}\Sigma_{E}$ , which has dimension $2n-1$ [51]. The solution space $S$ to the inhomogeneous equation $dH(\bm{\zeta})=1$ is a translation of this vector space (an affine subspace); hence, it also has dimension $2n-1$ . Explicitly, if $\bm{\zeta}_{0}$ is any particular solution (for example, $\bm{\zeta}_{0}=\partial_{p_{1}}/(\partial H/\partial p_{1})$ if $\partial H/\partial p_{1}\neq 0$ ), then the general solution is: $$ \bm{\zeta}=\bm{\zeta}_{0}+\bm{v},\quad\text{where }\bm{v}\in T_{\bm{x}}\Sigma_{E}. $$ Since $T_{\bm{x}}\Sigma_{E}$ has dimension $2n-1$ , we have a $(2n-1)$ -parameter family of transverse vectors. Different choices of $\bm{\zeta}$ yield different $(2n-1)$ -forms $\iota_{\bm{\zeta}}d\mu_{\Lambda}$ when restricted to $\Sigma_{E}$ . Without additional structure (such as a metric to define a canonical normal direction), there is no geometric principle within the symplectic framework to select one element of $S$ over another. We remark that these three obstructions establish that symplectic geometry, while sufficient to identify the energy shell $\Sigma_{E}$ and govern the dynamics upon it, does not provide the geometric data necessary to measure the volume of $\Sigma_{E}$ . The introduction of a metric tensor resolves all three obstructions simultaneously: it eliminates the dimensional mismatch by providing the gradient $\nabla H$ as a canonical transverse direction; it removes the non-uniqueness by defining orthogonality; and it repairs the degeneracy by inducing a non-degenerate metric on $\Sigma_{E}$ itself. ## B Necessity of the energy clock The gradient of $H$ , while providing the transverse direction to $\Sigma_{E}$ , does not provide a uniform parametrization of the energy increments across $\Sigma_{E}$ . In order to see that consider a small displacement along $\nabla H$ under parametrization $\epsilon$ : $$ \bm{x}(\epsilon)=\bm{x}_{0}+\nabla H(\bm{x}_{0})\,d\epsilon\,. $$ The target energy is $$ E(\bm{x}(\epsilon))=H(\bm{x}_{0})+\nabla H(\bm{x}_{0})\cdot\nabla H(\bm{x}_{0})\,d\epsilon. $$ Therefore, the energy step is $$ dE_{\bm{x}}=E(\bm{x}(\epsilon))-E_{0}=\|\nabla H(\bm{x}_{0})\|^{2}\;d\epsilon $$ Since $\|\nabla H\|$ varies across $\Sigma_{E}$ , for two distinct points $\bm{x}\neq\bm{y}$ on the same $\Sigma_{E}$ , we have $dE_{\bm{x}}\neq dE_{\bm{y}}$ , mapping them to different target hypersurfaces. To resolve this, it is sufficient to determine the vector field $\bm{\zeta}\in\text{span}\{\nabla H\}$ such that $$ dH(\bm{\zeta})=1\implies\bm{\zeta}=\frac{\nabla_{\eta}H}{\|\nabla_{\eta}H\|_{\eta}^{2}}\,. $$ This vector field generates the motion in energy that we call thermodynamic dynamics, namely, the diffeomorphisms between energy hypersurfaces, $$ \Phi^{\rm diff}_{\bm{\xi}}:\Sigma_{E}\to\Sigma_{E^{\prime}},\qquad\frac{d\Phi^{\rm diff}}{d\epsilon}=\bm{\zeta}(\bm{x}(\epsilon)) $$ It thus provides the transverse direction that allows the decomposing $$ T_{\bm{x}}\Lambda=\text{span}\{\bm{\zeta}\}\oplus T_{\bm{x}}\Sigma_{E}\,, $$ and ensures that the energy motion maps all the points on $\Sigma_{E}$ into points of $\Sigma_{E^{\prime}}$ simultaneously. ## C Euclidean metric tensor in adapted frame The $\bm{\zeta}$ -decomposition of $T_{\bm{x}}\Sigma_{E}$ naturally defines adapted coordinates $(E,y^{\alpha})$ where $E=H(\bm{x})$ parametrizes the energy and $\{y^{\alpha}\}_{\alpha=1}^{2n-1}$ on $\Sigma_{E}$ . Let $E:=H(\bm{x})$ and let $y^{\alpha}$ ( $\alpha=1,\dots,2n-1$ ) be local coordinates on the level set $\Sigma_{E}=\{H=E\}$ . Choose the embedding $x=x(E,y)$ so that its coordinate vectors are $$ \partial_{E}x=\bm{\zeta},\qquad\partial_{\alpha}x\in T\Sigma_{E},\ \ \ dH(\partial_{\alpha}x)=0. $$ In canonical coordinates $(q^{i},p_{i})$ , the phase–space metric is $$ \eta=\delta_{ij}\,dq^{i}\!\otimes dq^{j}+\delta^{ij}\,dp_{i}\!\otimes dp_{j}. $$ and in adapted coordinates $(E,y)$ , its components become $$ \eta_{AB}(E,y)\;=\;\eta\big(\partial_{A}x,\partial_{B}x\big),\qquad A,B\in\{E,\alpha\}. $$ Then, using $dH(\partial_{\alpha}x)=0$ and $\bm{\zeta}=\nabla_{\eta}H/\|\nabla_{\eta}H\|_{\eta}^{2}$ , $$ \eta_{E\alpha}=\eta(\partial_{E}x,\partial_{\alpha}x)=\eta(\bm{\zeta},\partial_{\alpha}x)=\frac{dH(\partial_{\alpha}x)}{\|\nabla_{\eta}H\|_{\eta}^{2}}=0,\qquad\eta_{EE}=\eta(\bm{\zeta},\bm{\zeta})=\frac{1}{\|\nabla_{\eta}H\|_{\eta}^{2}}. $$ For the tangential components, we obtain the induced metric $$ \sigma^{\eta}_{\alpha\beta}(E,y):=\eta(\partial_{\alpha}x,\partial_{\beta}x), $$ where $\{\bm{e}_{\alpha}\}_{\alpha=1}^{2n-1}$ with $\bm{e}_{\alpha}:=\partial_{\alpha}x$ is the basis of $T_{\bm{x}}\Sigma_{E}$ . Therefore, the Euclidean metric in adapted coordinates takes the block–diagonal form $$ \eta\;=\;\frac{dE\otimes dE}{\|\nabla_{\eta}H\|_{\eta}^{2}}\;+\;\sigma^{\eta}_{\alpha\beta}(E,y)\,dy^{\alpha}\!\otimes dy^{\beta}\,, $$ with $dE(\partial_{E})=1$ , $dE(\partial_{\alpha})=0$ , and $\sigma^{\eta}_{E}$ the metric induced on $\Sigma_{E}$ . Note that the microcanonical density of states naturally emerges as the induced measure on $\Sigma_{E}$ . Indeed, the Liouville measure in adapted coordinates is expressed as $$ d\mu^{\eta}_{\Lambda}=\frac{dEd\sigma_{E}}{\|\nabla_{\eta}H\|_{\eta}}=dp_{1}\ldots dp_{n}dq^{1}\ldots dq^{n},\quad\implies d\mu^{\eta}_{E}=d\mu^{\eta}_{\Lambda}\Big|_{\Sigma_{E}}=\frac{d\sigma^{\eta}_{E}}{\|\nabla_{\eta}H\|_{\eta}} $$ which means $$ \int_{\Sigma_{E}}d\mu_{E}=\int_{\Sigma_{E}}\frac{d\sigma_{E}}{\|\nabla_{\eta}H\|_{\eta}}=\int_{\Lambda}\delta(H-E)~dp_{1}\ldots dp_{n}dq^{1}\ldots dq^{n}=\Omega_{B}(E). $$ where the coarea formula is used [27]. We should remark on a subtle consequence due to the physical request behind the energy clock discussed in Sec. B. From a pure mathematical viewpoint, given a generic oriented Riemannian manifold $(M,g)$ with volume form $d\mu_{X}^{g}$ , and smooth hypersurface $\Sigma\subset M$ , the area form is determined by the unit normal vector field $\bm{n}$ : ( $\|n\|_{g}=1$ , $g(n,T)=0$ for all $T\in T\Sigma$ ). This reads: $$ d\sigma^{g}\;:=\;\iota_{n}\,d\mu_{M}^{g}\Big|_{\Sigma}. $$ This is a coordinate-free definition (reversing $n$ flips the orientation and the sign of the form), and in adapted coordinates $(u^{1},\dots,u^{m-1},u^{m})$ with $\Sigma=\{u^{m}=0\}$ it reduces to the usual expression $d\sigma^{g}=\sqrt{\det g_{\top}}\,du^{1}\wedge\cdots\wedge du^{m-1}$ , where $g_{\top}$ is the metric induced on $T\Sigma$ . Then, the area of a Borel set $A\subset\Sigma$ is $\mathrm{area}_{g}(A)=\int_{A}d\sigma^{g}$ . In our setting $(\Lambda,H,\omega,\eta)$ , the energy shells are $\Sigma_{E}=\{H=E\}$ , and the geometric area is $$ \mathrm{area}_{\eta}(\Sigma_{E})=\int_{\Sigma_{E}}d\sigma_{E}^{\eta},\qquad d\sigma_{E}^{\eta}:=\sqrt{\det\sigma^{\eta}_{E}}\,dx^{1}\cdots dx^{2n-1}. $$ However, thermodynamic quantities are described as function of energy, and therefore defined on $\Sigma_{E}$ . Energy variation of a thermodynamic quantities must then be expressed in units of energy. Therefore, moving from one energy value to another of a quantity $\delta E$ implies that we are following the energy motion identified by $\bm{\zeta}$ . This vector indeed “moves” every point of $\Sigma_{E}$ to $\Sigma_{E+\delta E}$ as guaranteed by the condition $dH(\bm{\zeta})=1$ . Therefore, decomposing Liouville’s volume along this clock yields the canonical microcanonical surface measure on $\Sigma_{E}$ which is given by $$ d\mu_{E}^{\eta}:=\iota_{\zeta}\,d\mu_{\Lambda}\Big|_{\Sigma_{E}}=\frac{\iota_{n}d\mu_{\Lambda}}{\|\nabla_{\eta}H\|_{\eta}}=\frac{d\sigma_{E}^{\eta}}{\|\nabla_{\eta}H\|_{\eta}}. $$ ## D Mapping from Euclidean metric tensor to unit-norm gauge tensor Let us consider the Euclidean metric tensor in adapted coordinates in Eq. (S14). We define the transformation [19]: $$ dx^{0}=\chi\,dE,\quad dx^{i}=\chi^{-\frac{1}{n-1}}\,dy^{i},\qquad\chi:=\frac{1}{\|\nabla H\|_{\eta}} $$ where $i\in[1,2n-1]$ . In these new coordinates $(x^{0},x^{1},\ldots,x^{2n-1})$ , the metric tensor is (bi-)conformally related to $\eta$ : $$ g_{00}=\chi^{-2}\,\eta_{EE},\qquad(\sigma^{g}_{E})_{ij}=\chi^{\frac{2}{n-1}}\,(\sigma^{\eta}_{E})_{ij}, $$ yielding the equivalent form $$ g=dx^{0}\otimes dx^{0}+(\sigma^{g}_{E})_{ij}\,dx^{i}\otimes dx^{j}. $$ This defines a new metric tensor $g$ that we call the unit-norm gauge. This metric tensor has peculiar properties. ### D.1 The vector field generating the energy flow in the unit-norm gauge First, in this metric, the vector field generating the flow of the diffeomorphism $\bm{\xi}$ must satisfy $dx^{0}(\bm{\xi})=1$ , and therefore $\bm{\xi}\equiv\partial_{x^{0}}$ . From Eq. (S16), we have $$ \begin{split}dx^{0}(\bm{\xi})=\chi dE(\bm{\xi})=1,\end{split} $$ on the other side, $dE(\bm{\zeta})=1$ therefore $$ dE(\bm{\zeta})=dx^{0}(\bm{\xi})=1\quad\implies\bm{\xi}=\chi^{-1}\bm{\zeta}. $$ Finally, $\bm{\xi}$ is such that $\|\bm{\xi}\|_{g}=1$ . This is easy to check $$ \begin{split}g(\bm{\xi},\bm{\xi})=g(\partial_{x^{0}},\partial_{x^{0}})=g_{00}=\chi^{-2}\eta_{EE}=1,\end{split} $$ as follows from Eq. (S13). Finally, we observe that $\bm{\xi}=\nabla_{g}H$ . Indeed: $$ \nabla_{g}H=g^{00}\frac{\partial H}{\partial x^{0}}\partial_{x^{0}}+g^{ij}\frac{\partial H}{\partial x^{i}}\partial_{x^{j}} $$ However, independently by the metric tensor ( $g$ remains compatible with $\omega$ ): $$ g(\nabla_{g}H,\bm{X}_{H})=dH(\bm{X}_{H})=0 $$ therefore, any components of $\nabla_{g}H$ on $T_{\bm{x}}\Sigma_{E}$ must vanish: $\partial H/\partial x^{i}=0$ which leads to $$ \nabla_{g}H=g^{00}\frac{\partial H}{\partial x^{0}}\partial_{x^{0}} $$ Finally, since $g^{00}=g_{00}=1$ , we notice that $$ \frac{\partial H}{\partial x^{0}}=dH(\partial_{x^{0}})=1\;\;\implies\nabla_{g}H\equiv\bm{\xi}\equiv\partial_{x^{0}}. $$ Finally, we notice that in the unit-norm gauge the induced measure on $\Sigma_{E}$ coincides with the area form $d\mu_{E}^{g}=d\sigma_{E}^{g}$ , being $\|\bm{\xi}\|_{g}=1$ . Then, the Boltzmann density of states coincides with the geometric area: $$ \mathrm{area}_{g}(\Sigma_{E})=\int_{\Sigma_{E}}d\sigma_{E}^{g}=\int_{\Sigma_{E}}\frac{d\sigma^{\eta}_{E}}{\|\nabla_{\eta}H\|_{\eta}^{2}}=\int_{\Lambda}\delta(H-E)\,d\mu_{\Lambda}=\Omega_{B}(E). $$ ## E Calculation of second-order GCF In this section, we compute the explicit form of $\Upsilon^{(2)}_{g}$ in terms of gradient, Hessian and Laplacian associated to the metric tensor $\eta$ . The expression that we will obtain is the one used in numerical simulations, namely, evaluated during the dynamics. We consider $$ \mathrm{Tr}[W_{\bm{\xi}}]=\frac{\Delta H}{\|\nabla H\|^{2}}-2\,\frac{\langle\nabla H,(\text{Hess}\,H)\,\nabla H\rangle}{\|\nabla H\|^{4}},\qquad\bm{\xi}=\frac{\nabla H}{\|\nabla H\|^{2}}, $$ we introduce the shorthand notation: $$ u:=\nabla H,\qquad G:=\|\nabla H\|^{2},\qquad A:=\Delta H,\qquad B:=\langle u,(\text{Hess}\,H)u\rangle. $$ Then, $$ \mathrm{Tr}[W_{\bm{\xi}}]=\frac{A}{G}-\frac{2B}{G^{2}}. $$ Since $\partial_{E}=\mathscr{L}_{\bm{\xi}}=(\bm{\xi}\!\cdot\!\nabla)=(u/G)\!\cdot\!\nabla$ , we can compute $\partial_{E}\mathrm{Tr}[W_{\bm{\xi}}]$ , i.e.: $$ \displaystyle\partial_{E}\mathrm{Tr}[W_{\bm{\xi}}] \displaystyle=\frac{\langle u,\nabla A\rangle}{G^{2}}-\frac{2A\langle u,\nabla G\rangle}{G^{3}}-\frac{2\langle u,\nabla B\rangle}{G^{3}}+\frac{4B\langle u,\nabla G\rangle}{G^{4}}. $$ Using the auxiliary relations $$ \langle u,\nabla G\rangle=2B,\qquad\langle u,\nabla B\rangle=2\,\langle u,(\text{Hess}\,H)^{2}u\rangle+\nabla^{3}H(u,u,u), $$ we obtain $$ \partial_{E}\mathrm{Tr}[W_{\bm{\xi}}]=\frac{\langle u,\nabla A\rangle}{G^{2}}-\frac{2AB}{G^{3}}-\frac{4\,\langle u,(\text{Hess}\,H)^{2}u\rangle+2\,\nabla^{3}H(u,u,u)}{G^{3}}+\frac{8B^{2}}{G^{4}}. $$ Furthermore, $$ \big(\mathrm{Tr}[W_{\bm{\xi}}]\big)^{2}=\frac{A^{2}}{G^{2}}-\frac{4AB}{G^{3}}+\frac{4B^{2}}{G^{4}}\,. $$ Therefore: $$ \partial_{E}\mathrm{Tr}[W_{\bm{\xi}}]+\big(\mathrm{Tr}[W_{\bm{\xi}}]\big)^{2}=\frac{A^{2}+\langle u,\nabla A\rangle}{G^{2}}-\frac{6AB+4\,\langle u,(\text{Hess}\,H)^{2}u\rangle+2\,\nabla^{3}H(u,u,u)}{G^{3}}+\frac{12B^{2}}{G^{4}} $$ Collecting terms gives the final result: $$ \displaystyle\partial_{E}^{2}{\rm vol}^{g}(E) \displaystyle=\int_{\Sigma_{E}}\!\!\Bigg[\frac{A^{2}+\langle u,\nabla A\rangle}{G^{2}}-\frac{6AB+4\,\langle u,(\text{Hess}\,H)^{2}u\rangle+2\,\nabla^{3}H(u,u,u)}{G^{3}}+\frac{12B^{2}}{G^{4}}\Bigg]d\mu_{g} \displaystyle=\int_{\Sigma_{E}}\!\!\Bigg[\frac{(\Delta H)^{2}+\langle\nabla H,\nabla(\Delta H)\rangle}{\|\nabla H\|^{4}}-\frac{6\,(\Delta H)\langle\nabla H,(\text{Hess}\,H)\nabla H\rangle}{\|\nabla H\|^{6}}-\frac{4\,\langle\nabla H,(\text{Hess}\,H)^{2}\nabla H\rangle}{\|\nabla H\|^{6}} \displaystyle\hskip 142.26378pt-\frac{2\,\nabla^{3}H(\nabla H,\nabla H,\nabla H)}{\|\nabla H\|^{6}}+\frac{12\,\langle\nabla H,(\text{Hess}\,H)\nabla H\rangle^{2}}{\|\nabla H\|^{8}}\Bigg]d\mu_{g}. $$ Finally, the second-order GCF reads $$ \begin{split}\Upsilon^{(2)}_{g}(E)=\frac{\partial_{E}^{2}{\rm vol}^{g}(E)}{\text{area}^{g}(E)}&=\int_{\Sigma_{E}}\!\!\Bigg[\frac{(\Delta H)^{2}+\langle\nabla H,\nabla(\Delta H)\rangle}{\|\nabla H\|^{4}}-\frac{6\,(\Delta H)\langle\nabla H,(\text{Hess}\,H)\nabla H\rangle}{\|\nabla H\|^{6}}-\frac{4\,\langle\nabla H,(\text{Hess}\,H)^{2}\nabla H\rangle}{\|\nabla H\|^{6}}\\ &\hskip 142.26378pt-\frac{2\,\nabla^{3}H(\nabla H,\nabla H,\nabla H)}{\|\nabla H\|^{6}}+\frac{12\,\langle\nabla H,(\text{Hess}\,H)\nabla H\rangle^{2}}{\|\nabla H\|^{8}}\Bigg]d\rho_{g}.\end{split} \tag{2} $$ with $d\rho_{g}=d\mu^{g}_{E}/\text{area}^{g}(E)$ . Considering the microcanonical average $$ \langle O\rangle_{E}=\int O(\pi,\phi)\delta(H(\pi,\phi)-E)~D\pi\,D\phi $$ then this object is estimated by computing $$ \begin{split}\Upsilon^{(2)}_{g}(E)=\left\langle\frac{(\Delta H)^{2}+\langle\nabla H,\nabla(\Delta H)\rangle}{\|\nabla H\|^{4}}\right\rangle_{E}&-\left\langle\frac{6\,(\Delta H)\langle\nabla H,(\text{Hess}\,H)\nabla H\rangle}{\|\nabla H\|^{6}}\right\rangle_{E}-\left\langle\frac{4\,\langle\nabla H,(\text{Hess}\,H)^{2}\nabla H\rangle}{\|\nabla H\|^{6}}\right\rangle_{E}\\ &-\left\langle\frac{2\,\nabla^{3}H(\nabla H,\nabla H,\nabla H)}{\|\nabla H\|^{6}}\right\rangle_{E}+\left\langle\frac{12\,\langle\nabla H,(\text{Hess}\,H)\nabla H\rangle^{2}}{\|\nabla H\|^{8}}\right\rangle_{E}\end{split} \tag{2} $$ ## F Expansion of the Weingarten operator ### F.1 Angular statistics and small- $M$ expansions Let $\Theta_{i}=\theta_{i}-\phi$ and $M:=\big\|\frac{1}{N}\sum_{i}(\cos\theta_{i},\sin\theta_{i})\big\|\ll 1$ . On the disordered branch we model the angular fluctuations by the von Mises density $$ f(\Delta)=\frac{e^{h\cos\Delta}}{2\pi I_{0}(h)},\qquad h\ll 1, $$ whose moments satisfy (see standard Refs. [53, 54, 55]) $$ \langle\cos(k\Delta)\rangle=I_{k}(h)/I_{0}(h) $$ Another standard identity is $\langle\sin(k\Delta)\rangle=0$ by symmetry, whence $$ \langle\cos\Delta\rangle=\frac{I_{1}}{I_{0}}=\frac{h}{2}-\frac{h^{3}}{16}+O(h^{5})=:M\quad\Rightarrow\quad h=2M+O(M^{3}). $$ Small– $h$ expansions for the modified Bessel functions $I_{k}$ follow from classical handbooks [56]. Using $\cos^{2}\Delta=\tfrac{1}{2}(1+\cos 2\Delta)$ and the small- $h$ expansions of $I_{k}$ , $$ \langle\cos^{2}\Delta\rangle=\frac{1}{2}+\frac{1}{2}\frac{I_{2}}{I_{0}}=\frac{1}{2}+\frac{h^{2}}{16}+O(h^{4})=\frac{1}{2}+\frac{M^{2}}{4}+O(M^{4}), $$ hence $$ \langle\sin^{2}\Delta\rangle=\frac{1}{2}-\frac{M^{2}}{4}+O(M^{4}). $$ Moreover, using $\cos^{3}\Delta=\tfrac{1}{4}(3\cos\Delta+\cos 3\Delta)$ , $$ \langle\cos^{3}\Delta\rangle=\frac{3}{4}M+\frac{1}{4}\frac{I_{3}}{I_{0}}=\frac{3}{4}M+\frac{h^{3}}{192}+O(h^{5})=\frac{3}{4}M+\frac{M^{3}}{24}+O(M^{5}), $$ so that $$ \langle\sin^{2}\Delta\,\cos\Delta\rangle=\langle\cos\Delta-\cos^{3}\Delta\rangle=\frac{1}{4}M-\frac{1}{24}M^{3}+O(M^{5}). $$ By the law of large numbers, sums over $i=1,\dots,N$ are replaced by $N$ times the expectations up to $O(\sqrt{N})$ fluctuations. Therefore, $$ S_{2}:=\sum_{i=1}^{N}\sin^{2}\Theta_{i}=N\,\langle\sin^{2}\Delta\rangle=N\Big(\frac{1}{2}-\frac{1}{4}M^{2}\Big)+O(NM^{4})+O(\sqrt{N}), $$ $$ S_{21}:=\sum_{i=1}^{N}\sin^{2}\Theta_{i}\,\cos\Theta_{i}=N\,\langle\sin^{2}\Delta\,\cos\Delta\rangle=\frac{MN}{4}+O(NM^{3}). $$ For the scaling of $S_{\mathrm{int}}$ we use $\cos(\Theta_{i}-\Theta_{j})=\cos\Theta_{i}\cos\Theta_{j}+\sin\Theta_{i}\sin\Theta_{j}$ , and we get $$ \displaystyle S_{\mathrm{int}} \displaystyle=\sum_{i\neq j}\sin\Theta_{i}\sin\Theta_{j}\cos(\Theta_{i}-\Theta_{j}) \displaystyle=\sum_{i\neq j}\big(\sin\Theta_{i}\cos\Theta_{i}\big)\big(\sin\Theta_{j}\cos\Theta_{j}\big)+\sum_{i\neq j}\sin^{2}\Theta_{i}\,\sin^{2}\Theta_{j} \displaystyle=\Big(\sum_{i}\sin\Theta_{i}\cos\Theta_{i}\Big)^{2}-\sum_{i}\sin^{2}\Theta_{i}\cos^{2}\Theta_{i}+\Big(\sum_{i}\sin^{2}\Theta_{i}\Big)^{2}-\sum_{i}\sin^{4}\Theta_{i}. $$ By symmetry $\sum_{i}\sin\Theta_{i}\cos\Theta_{i}=O(\sqrt{N})$ , hence its square is $O(N)$ and subleading versus $N^{2}$ terms. Set $A:=\sum_{i}\sin^{2}\Theta_{i}$ ; from (S30), $$ A=\frac{N}{2}-\frac{N}{4}M^{2}+O(NM^{4}) $$ Using $$ \sin^{4}\Delta=\frac{3-4\cos 2\Delta+\cos 4\Delta}{8} $$ and $I_{4}/I_{0}=O(h^{4})$ , we get $$ \langle\sin^{4}\Delta\rangle=\frac{3}{8}-\frac{1}{2}\frac{I_{2}}{I_{0}}+\frac{1}{8}\frac{I_{4}}{I_{0}}=\frac{3}{8}-\frac{h^{2}}{16}+O(h^{4})=\frac{3}{8}-\frac{M^{2}}{4}+O(M^{4}), $$ so that $\sum_{i}\sin^{4}\Theta_{i}=N\big(\tfrac{3}{8}-\tfrac{M^{2}}{4}\big)+O(NM^{4})+O(\sqrt{N})$ . Plugging into (S32) and retaining $O(N^{2})$ terms, $$ S_{\mathrm{int}}=\Big(\tfrac{N}{2}-\tfrac{N}{4}M^{2}\Big)^{2}-\sum_{i}\sin^{4}\Theta_{i}+O(N)=\frac{N^{2}}{4}+O(N^{2}M^{2})-\frac{3N}{8}+O(N), $$ which yields, to leading order in $N$ , $$ S_{\mathrm{int}}\simeq\frac{N^{2}}{4}\qquad(N\to\infty,\ M\to 0). $$ ### F.2 Denominator expansions and block sums Recall $K:=\sum_{i}p_{i}^{2}$ and $$ G:=\|\nabla H\|^{2}=K+J^{2}M^{2}S_{2}. $$ Write $c:=K/N$ (finite in the thermodynamic limit). Using (S30), $$ G=N\Big(c+\frac{J^{2}}{2}M^{2}\Big)+O(NM^{4}). $$ Therefore, $$ \displaystyle\frac{1}{G} \displaystyle=\frac{1}{Nc}\Big(1+\frac{J^{2}}{2c}M^{2}\Big)^{-1}=\frac{1}{Nc}\Big(1-\frac{J^{2}}{2c}M^{2}+O(M^{4})\Big), \displaystyle\frac{1}{G^{2}} \displaystyle=\frac{1}{N^{2}c^{2}}\Big(1+\frac{J^{2}}{2c}M^{2}\Big)^{-2}=\frac{1}{N^{2}c^{2}}\Big(1-\frac{J^{2}}{c}M^{2}+O(M^{4})\Big). $$ Inserting (S30), (S31), (S33), (S35)–(S36) into the blocks $$ \sum_{i}\kappa_{p_{i}}=\frac{N}{G}-\frac{2K}{G^{2}},\qquad\sum_{i}\kappa^{\mathrm{d}}_{\theta_{i}}=\frac{JNM^{2}}{G}-\frac{2J^{3}M^{3}S_{21}}{G^{2}},\qquad\sum_{i\neq j}\kappa^{\mathrm{int}}_{ij}=\frac{2J^{3}M^{2}}{N\,G^{2}}S_{\mathrm{int}}, $$ and keeping only $O(N^{0})$ terms (the ones that survive as $N\to\infty$ ), reproduces the expressions for $\bm{P}$ and $\bm{I}$ of the main text. Summary (small- $M$ , thermodynamic limit). $$ S_{2}\simeq N\Big(\frac{1}{2}-\frac{1}{4}M^{2}\Big),\qquad S_{21}\simeq\frac{MN}{4},\qquad S_{\mathrm{int}}\simeq\frac{N^{2}}{4}, $$ $$ \frac{1}{G}\simeq\frac{1}{Nc}\Big(1-\frac{J^{2}}{2c}M^{2}\Big),\qquad\frac{1}{G^{2}}\simeq\frac{1}{N^{2}c^{2}}\Big(1-\frac{J^{2}}{c}M^{2}\Big), $$ ## G Numerical methods and integration schemes Consider $\phi$ as a generalized coordinates: this can be the real field configuration in the 2D $\phi^{4}$ model or the angle variable $\theta$ in the XY model. ### G.1 Configurational microcanonical Monte Carlo algorithm The microcanonical ensemble has been reproduced by adopting the method proposed in Refs. [58, 57] and consisting of a microcanonical Monte Carlo (MICROMC) algorithm in the configuration space. We have chosen a method based on a random proposal to enforce ergodicity. Note that the numerical procedure that we adopted and reported below is system-independent. ### G.2 Microcanonical distribution function The microcanonical algorithm is based on the identification of a suitable density probability function that is used to produce a reliable sampling. To find such a probability function, we start with the microcanonical partition function: $$ \Omega(E)=\int\delta(H[\pi,\phi]-E)D\pi\,D\phi\,, $$ where $$ D\pi=\prod_{\bm{n}\in\mathbb{L}}d\pi_{\bm{n}},\qquad D\phi=\prod_{\bm{n}\in\mathbb{L}}d\phi_{\bm{n}}\,, $$ and where the Hamiltonian is $$ H(\pi,\phi)=\sum_{\bm{n}\in\mathbb{L}}\frac{\pi_{\bm{n}}^{2}}{2}+V(\phi)\,, $$ for any potential function. Finally, $\delta$ represents the Dirac delta function. From now on, we introduce $L=N^{2}$ to denote the total number of degrees of freedom. We focus on the momentum integral that separates as $$ I_{p}=\int\delta\left(\sum_{\bm{n}\in\mathbb{L}}\frac{\pi_{\bm{n}}^{2}}{2}+V(\phi)-E\right)D\pi\,. $$ Thus, we notice that the kinetic energy variable, $$ K=\sum_{\bm{n}\in\mathbb{L}}\frac{p_{\bm{n}}^{2}}{2}\;, $$ can be interpreted as the equation for a $L$ -dimensional sphere of radius $r^{2}=2K$ . Then, exploiting the spherical coordinates in $L$ dimensions, we have $$ \prod_{\bm{n}\in\mathbb{L}}d\pi_{\bm{n}}=\rho^{L-1}\,d\rho\,d\Omega_{L-1}\,, $$ where $d\Omega_{L-1}$ is the differential solid angle element on the unit sphere $S^{L-1}$ , expressed as: $$ d\Omega_{L-1}=\prod_{i=1}^{L-1}d\theta_{i}\sin^{L-1-i}\theta_{i},\qquad\Omega_{L}=\int d\Omega_{L-1}\;. $$ Then, we define $$ \rho=\sqrt{2K},\quad d\rho=\frac{dK}{\sqrt{2K}}\implies\rho^{L-1}d\rho=(2K)^{L/2-1}\,dK\,. $$ Finally, integral (S37) rewrites $$ \begin{split}I_{p}&=\Omega_{L}\int\delta\left[K-(E-V(\phi))\right](2K)^{L/2-1}\;dK\\ &=2^{L-1}\Omega_{L}(E-V(\phi))^{L/2-1}\Theta(E-V(\phi))\,.\end{split} $$ Reintroducing above the integral over field configurations, we recover the configurational microcanonical partition function that reads $$ \Omega(E)=\int(E-V(\phi))^{L/2-1}\Theta(E-V(\phi))\;D\phi. $$ Notice that we have dropped irrelevant constants that play no role in the calculation of expectation values. The numerical algorithm then exploits the probability distribution function arising from Eq. (S39) $$ W_{E}[\phi]=\left(E-V(\phi)\right)^{L/2-1}\Theta[E-V(\phi)]\,. $$ The sampling is carried out by randomly picking a lattice site, $\bm{n}$ , and proposing a new random configuration $\phi^{old}_{\bm{n}}\mapsto\phi^{new}_{\bm{n}}=\phi^{old}_{\bm{n}}+\eta\Delta\phi$ where $\eta$ is a random number uniformly sampled from the interval $[-1,1]$ and $\Delta\phi$ is adjusted to achieve a target acceptance rate of $50\$ . The new configuration is accepted according to the Monte Carlo acceptance probability $$ W(\phi^{old}\to\phi^{new})=\min\left(1,\frac{W_{E}[\phi^{new}]}{W_{E}[\phi^{old}]}\right)\,. $$ Notice that an efficient way to evaluate this acceptance probability is to rewrite the acceptance ratio as follows: $$ \frac{W_{E}[\phi^{new}]}{W_{E}[\phi^{old}]}=\exp\left[\bigg(\frac{L}{2}-1\bigg)\log\left(\frac{E-V(\phi^{new})}{E-V(\phi^{old})}\right)\right]\,. $$ once we have checked that $E-V(\phi^{new})>0$ . ### G.3 Initial configuration Initial conditions are randomly proposed in order to start with a high-energy configuration, $E_{rand}$ , larger than the desired input energy $E_{inp}$ and such that $E_{rand}-V(\phi_{ini})>0$ . Then, we perform $10^{4}$ steps of equilibration with the MICROMC algorithm. At this stage, we search for the correct value for $\Delta\phi$ . To do that, we start with a small value for $\Delta\phi$ , say $0.001$ and run a few MICROMC steps using the desired input energy $E_{inp}$ and computing the acceptance rate, $N_{\text{acc}}$ . If $N_{\text{acc}}\notin[0.5,0.6]$ , then we repeat the procedure by replacing $\Delta\phi$ with $\Delta\phi+0.001$ . It should be stressed that other strategies, such as the molecular dynamics algorithm, have also been used to select the initial conditions. Comparison of the thermodynamic observables obtained with both methods did not yield appreciable differences. Finally, once the suitable tune parameter’s value has been obtained, the configuration is equilibrated for $10^{4}$ steps, and the trajectory is evolved with the MICROMC method and used for computing averages. For each energy value, we have produced $N_{trj}=32$ realizations for each system size. Each thermodynamic observable has been evaluated through $N_{\text{avg}}=10^{6}$ measurements, performed in every $N_{\text{step}}=100$ MICROMC step. The microcanonical average of a given observable, $f$ , is computed by $$ \langle f\rangle_{\varepsilon}=\frac{1}{N_{\text{trj}}\cdot N_{\text{avg}}}\sum_{i=1}^{N_{trj}}\sum_{\alpha=1}^{N_{\text{avg}}}f(\phi^{(i)}_{\alpha})\,, $$ where $f$ is the observable evaluated on the configuration $\phi^{(i)}_{\alpha}$ at the $\alpha$ -th MC step for the $i$ -th realization. ### G.4 Numerical calculation of microcanonical entropy’s derivatives The MIPA method requires the estimation of higher-order derivatives of the microcanonical entropy from a microcanonical sampling of the phase space. To this purpose, we adopt the Pearson-Halicioglu-Tiller (PHT) method [28] which allows a direct calculation of the first- and second-order derivatives of the microcanonical entropy as follows. Any energy derivative of the entropy function can be rewritten in terms of averages of power of the kinetic energy $k=(E-V(\phi)/L$ . In this framework, the first- and second-order derivatives of the microcanonical entropy can be written as $$ \displaystyle\partial_{\varepsilon}S(\varepsilon) \displaystyle=\left(\frac{1}{2}-\frac{1}{L}\right)\langle k^{-1}\rangle_{\varepsilon}\,, \displaystyle\partial^{2}_{\varepsilon}S(\varepsilon) \displaystyle=L\bigg[\bigg(\dfrac{1}{2}-\dfrac{1}{L}\bigg)\bigg(\dfrac{1}{2}-\dfrac{2}{L}\bigg)\langle k^{-2}\rangle_{\varepsilon}-\bigg(\dfrac{1}{2}-\dfrac{1}{L}\bigg)^{2}\langle k^{-1}\rangle^{2}_{\varepsilon}\bigg]\,, $$ where the microcanonical averages have been estimated by Eq. (S40) over the realizations. ## H Numerical solution of the Entropy Flow Equation (EFE) #### Equation and unknowns. We solve the EFE $$ S_{g}^{\prime\prime}(E)+\big(S_{g}^{\prime}(E)\big)^{2}\;=\;\Upsilon^{(2)}_{g}(E), \tag{2} $$ on an interval $E\in[E_{\min},E_{\max}]$ . Set $$ g(E):=S_{g}^{\prime}(E),\qquad S_{g}^{\prime}(E)=g(E),\qquad S_{g}^{\prime\prime}(E)=g^{\prime}(E), $$ so that (S43) reduces to the Riccati ODE $$ g^{\prime}(E)+g(E)^{2}\;=\;\Upsilon^{(2)}_{g}(E). \tag{2} $$ Once $g$ is known, the entropy follows by a single quadrature, $$ S_{g}(E)=S_{g}(E_{0})+\int_{E_{0}}^{E}g(\tilde{E})\,d\tilde{E}. $$ #### Input data: geometric source. The source $\Upsilon^{(2)}_{g}$ is computed numerically at a discrete set of energies $\{(E_{i},\Upsilon_{i})\}_{i=1}^{M}$ , with $E_{\min}\leq E_{1}<\dots<E_{M}\leq E_{\max}$ . To obtain a smooth right–hand side in (S44), we construct a high–order interpolant $$ \widehat{\Upsilon}(E)\;:=\;\mathcal{I}_{\!p}\big[\{(E_{i},\Upsilon_{i})\}\big](E), $$ using a polynomial spline of order $p=9$ (Mathematica’s Interpolation with InterpolationOrder $\to$ 9). This yields a $C^{p-1}$ function $\widehat{\Upsilon}$ on $[E_{\min},E_{\max}]$ that we use in place of $\Upsilon^{(2)}_{g}$ in (S44). #### Initial data. We pick a left endpoint $E_{0}\in[E_{\min},E_{\max})$ and prescribe: $$ g(E_{0})=\beta_{0},\qquad S_{g}(E_{0})=0, $$ where $\beta_{0}$ is the microcanonical inverse temperature at $E_{0}$ obtained from numerical simulation through Eq. (S41) and $S_{g}(E_{0})$ sets the additive entropy, since it cancels in derivatives, it can be fixed to zero. #### Time–integration (energy stepping). We integrate the first–order system $$ \begin{cases}S^{\prime}(E)=g(E),\\[2.0pt] g^{\prime}(E)=\widehat{\Upsilon}(E)-g(E)^{2},\end{cases}\qquad E\in[E_{0},E_{\max}], $$ with an explicit high–order Runge–Kutta method. In practice we used Mathematica’s NDSolve with $$ \texttt{Method}\;=\;\{\texttt{"TimeIntegration"}\to\{\texttt{"ExplicitRungeKutta"},\ \texttt{"DifferenceOrder"}\to 9\}\}, $$ and uniform step $\Delta E$ on $[E_{0},E_{\max}]$ . The numerical solution returns $S(E)$ , $g(E)$ , and, where needed, $g^{\prime}(E)$ . #### Discretization and sampling. For visualization and post–processing we evaluate $(g,g^{\prime})$ on the uniform grid $$ E_{j}\;=\;E_{0}+j\,\Delta E,\qquad j=0,1,\dots,J,\quad E_{J}\leq E_{\max}. $$ This produces the arrays $\{(E_{j},g(E_{j}))\}$ and $\{(E_{j},g^{\prime}(E_{j}))\}$ used in the figures. #### Stability and cross–checks. The Riccati form (S44) can be validated against the equivalent second–order linear problem via the substitution $g=\psi^{\prime}/\psi$ , namely $$ \psi^{\prime\prime}(E)\;=\;\widehat{\Upsilon}(E)\,\psi(E),\qquad g(E)=\frac{\psi^{\prime}(E)}{\psi(E)},\qquad S_{g}(E)=S_{0}+\ln\frac{\psi(E)}{\psi(E_{0})}. $$ We solved (S49) with the same RK scheme and verified that the reconstructed $g$ matches the direct Riccati integration within the stated tolerances. Convergence was monitored by halving $\Delta E$ and by lowering the interpolant order $p$ (from 9 to 7) to ensure insensitivity to the smoothing. #### Reproducibility notes. (i) $\widehat{\Upsilon}(E)$ is evaluated only within the range $\{E_{i}\}$ used in numerical simulations (no extrapolation). (ii) All ODE solves use a fixed step and order (RK9); decreasing $\Delta E$ by a factor of 2 changes $g$ and $g^{\prime}$ by less than the line width. (iii) The additive constant $S_{0}$ is fixed by matching $S_{g}$ to the reference at $E_{0}$ (or by setting $S_{0}=0$ ); derivatives are independent of $S_{0}$ . #### Output. The procedure returns smooth functions $g(E)=S_{g}^{\prime}(E)$ and $g^{\prime}(E)=S_{g}^{\prime\prime}(E)$ on $[E_{0},E_{\max}]$ , together with the entropy $S_{g}(E)$ from (S45). These are the curves reported in the figures and used to compare with microcanonical observables. ## I Investigation of phase transitions The Hamiltonian for the 2D $\phi^{4}$ model with nearest neighbor interactions is $$ \begin{split}H:=\sum_{\bm{i}}\Bigg[\frac{\pi^{2}_{\bm{i}}}{2}+\frac{\lambda}{4!}\phi^{4}_{\bm{i}}&-\frac{\mu^{2}}{2}\phi^{2}_{\bm{i}}+\frac{J}{4}\sum_{\bm{k}\in N(\bm{i})}(\phi_{\bm{i}}-\phi_{\bm{i}})^{2}\Bigg]\end{split} $$ and we choose values $$ \lambda=3/5,\qquad\mu=2,\qquad J=1. $$ The label $\bm{i}:=(i_{1},i_{2})$ is the two-dimensional array of integer numbers used for labeling the sites; finally, we denote with $N(\bm{i})$ the set of the nearest neighbors of the $\bm{i}$ -th site. The Hamiltonian of the 1D XY mean-field is given by $$ H(\theta,p)=\sum_{i=1}^{N}\frac{p_{i}^{2}}{2}+\frac{J}{2N}\sum_{i,j}\bigl[1-\cos(\theta_{i}-\theta_{j})\bigr], $$ where $J=1$ . For the $\phi^{4}$ -model we chose the the energy range $[9,16]$ with energy resolution $\Delta E=0.1$ and system sizes $N=64^{4},\,128^{2},\,256^{2}$ . For XY model, the energy range is $[0.5,1]$ with energy resolution $\Delta E=0.01$ and studied three different system sizes $N=8100,\,14400,40000$ . In both cases, we compute the energy derivatives of entropy using the PHT method explained in Sec. G.4. Notice that the PHT is our reference for the microcanonical observables. This estimation is actually non-perturbative and exact (within statistical precision). In order words, the PHT method expresses thermodynamic derivatives as exact averages of microscopic quantities without approximations, model assumptions, or fitting parameters. Then, we also solved the EFE equation following the procedure of Sec. H. Then, for each system size, we computed $\beta^{0}_{N}(E_{\text{min}})$ and use it as an initial condition for the differential equation. This analysis is reported, respectively, in Fig. S1 and S2. <details> <summary>x2.png Details</summary>  ### Visual Description ## Multi-Panel Scientific Chart: 2D φ⁴ Nearest Neighbours Analysis ### Overview The image is a multi-panel scientific figure containing nine subplots arranged in a 3x3 grid. It presents a comparative analysis of three statistical quantities (γ_g^(2)(ε), ∂_ε S(ε), and ∂_ε S²(ε)) as a function of a parameter ε, for three different system sizes (N = 64², 128², 256²). The analysis compares results from two methods, "RICCATI" and "PEARSON," alongside a "MICRO AVERAGE." ### Components/Axes * **Main Title:** "2D φ⁴ nearest neighbours" (located at the top center). * **Legend:** Positioned below the main title. * **RICCATI:** Represented by a black filled circle (●). * **PEARSON:** Represented by a red open circle (○). * **MICRO AVERAGE:** Represented by a black cross (×). * **Grid Structure:** * **Columns (System Size):** * **Left Column (a):** N = 64² * **Middle Column (b):** N = 128² * **Right Column (c):** N = 256² * **Rows (Measured Quantity):** * **Top Row (1):** Y-axis label: γ_g^(2)(ε) * **Middle Row (2):** Y-axis label: ∂_ε S(ε) * **Bottom Row (3):** Y-axis label: ∂_ε S²(ε) * **X-Axis:** Common to all subplots. Label: ε (epsilon). Scale: Linear, ranging from 9 to 14, with major ticks at 9, 10, 11, 12, 13, 14. * **Subplot Labels:** Each subplot is labeled in its top-right corner (e.g., a.1, b.2, c.3). * **Critical Point Annotations:** Vertical, red, dashed lines in the middle and bottom rows mark a critical value ε_c. The annotated value is written vertically next to the line. * Subplot a.2: ε_c ≈ 11.2 * Subplot b.2: ε_c ≈ 11.1 * Subplot c.2: ε_c ≈ 11.1 * Subplots a.3, b.3, c.3: The line is present but not explicitly re-annotated; it aligns with the ε_c from the row above. ### Detailed Analysis **Row 1: γ_g^(2)(ε) vs. ε (Top Row - Subplots a.1, b.1, c.1)** * **Trend:** All three plots show a similar shape: a peak followed by a decay. The curve rises from ε=9, reaches a maximum near ε=11, then decreases as ε increases to 14. * **Data Series:** The plotted points are black crosses (×), corresponding to the "MICRO AVERAGE" from the legend. No distinct RICCATI (●) or PEARSON (○) series are visible in this row. * **Y-Axis Scale:** Ranges from 0.0000 to 0.0020. The peak value is approximately 0.0017 for N=64², increasing slightly to ~0.0018 for N=128² and N=256². * **Spatial Grounding:** The peak occurs just to the left of the vertical grid line at ε=11. **Row 2: ∂_ε S(ε) vs. ε (Middle Row - Subplots a.2, b.2, c.2)** * **Trend:** All three plots show a nearly linear, decreasing trend. The value starts high at ε=9 and decreases steadily as ε increases to 14. * **Data Series:** The plotted points are red open circles (○), corresponding to "PEARSON." A red dashed line connects these points. No RICCATI or MICRO AVERAGE series are visible. * **Y-Axis Scale:** Ranges from 0.050 to 0.060. The value at ε=9 is ~0.060, decreasing to ~0.050 at ε=14. * **Critical Point:** A vertical red dashed line is present at ε_c. The annotation reads "ε_c ≈ 11.2" for N=64² and "ε_c ≈ 11.1" for N=128² and N=256². **Row 3: ∂_ε S²(ε) vs. ε (Bottom Row - Subplots a.3, b.3, c.3)** * **Trend:** All three plots show a pronounced peak. The curve rises sharply from ε=9, peaks near ε=11, then falls off, approaching zero or a small negative value by ε=14. * **Data Series:** The plotted points are red open circles (○), corresponding to "PEARSON." A red dashed line connects these points. No RICCATI or MICRO AVERAGE series are visible. * **Y-Axis Scale:** Ranges from -0.0022 to -0.0012. The peak value is approximately -0.0013 for all system sizes. * **Spatial Grounding:** The peak aligns precisely with the vertical red dashed line marking ε_c from the row above. ### Key Observations 1. **Consistency Across Scales:** The qualitative behavior of all three measured quantities (γ_g^(2), ∂_ε S, ∂_ε S²) is remarkably consistent across the three system sizes (N=64², 128², 256²). This suggests the observed phenomena are robust and not finite-size artifacts. 2. **Methodological Focus:** The "PEARSON" method (red circles) is used for the derivative quantities (∂_ε S and ∂_ε S²), while the "MICRO AVERAGE" (black crosses) is used for γ_g^(2). The "RICCATI" method (black circles) is defined in the legend but does not appear to be plotted in any of these nine subplots. 3. **Critical Point Signature:** The vertical dashed line at ε_c ≈ 11.1-11.2 marks a significant point. It coincides with the peak in ∂_ε S²(ε) and lies within the decaying region of γ_g^(2)(ε) and on the linear slope of ∂_ε S(ε). 4. **Quantitative Shift:** There is a subtle but clear shift in the annotated critical value ε_c from ~11.2 for the smallest system (N=64²) to ~11.1 for the larger systems (N=128², 256²), indicating a possible finite-size effect on the critical point location. ### Interpretation This figure presents results from a numerical study of a 2D φ⁴ model, likely investigating a phase transition or critical phenomenon. The parameter ε is probably a control parameter like temperature or coupling strength. * **γ_g^(2)(ε)** likely represents a two-point correlation function or a related susceptibility. Its peak near ε=11 suggests enhanced fluctuations, a hallmark of approaching a critical point. * **∂_ε S(ε)** is the first derivative of an entropy-like quantity S with respect to ε. Its linear decrease indicates a smooth, continuous change in the system's disorder or information content across the parameter range. * **∂_ε S²(ε)**, the derivative of the squared entropy, shows a sharp peak at ε_c. This is a strong indicator of a critical point, as such derivatives often diverge or show non-analytic behavior at phase transitions. The peak's alignment with the annotated ε_c confirms this interpretation. * The consistency across system sizes (N) implies the study is examining scaling behavior. The slight shift in ε_c with increasing N is typical for numerical analyses of critical phenomena, where the true critical point is approached in the thermodynamic limit (N → ∞). **In summary, the data strongly suggests the system undergoes a continuous phase transition at a critical value ε_c ≈ 11.1.** The different plotted quantities provide complementary views of this transition: γ_g^(2) shows critical fluctuations, ∂_ε S shows the smooth background change, and ∂_ε S² provides a clear, sharp signature for pinpointing the critical point. The absence of the RICCATI data in these plots may mean it was used for a different analysis or is presented in another figure. </details> Figure S1: Comparison of entropy derivatives at different system sizes obtained as the solution of the entropy flow equation and through thermodynamic methods. The quantitative and qualitative agreement is excellent. As the system size increases, the peak in $\Upsilon^{(2)}_{g}$ and $\partial_{\epsilon}^{2}S$ becomes sharper and sharper indicating that in the thermodynamic limit the peak converges to zero thus giving rise to the divergence of the specific heat. The solutions of the EFE (black disks) are then plotted in comparison with the PHT-estimation (red circles) for different systems sizes. For the $\phi^{4}$ -model, we see a peak (negative-valued maximum) in $\partial_{\epsilon}^{2}S$ together with an inflection point in $\partial_{\epsilon}S$ at $\epsilon_{c}=11.2$ for $N=64^{2}$ and $\epsilon_{c}=11.1$ for $N=128^{2}$ and $N=256^{2}$ . Moreover, this peak increases with increasing system sizes and approaches zero. According to MIPA [32], this behavior is in agreement with the second-order PT expected in the thermodynamic limit for this system. Moreover, the specific heat $$ C_{v}(E)=\frac{(\partial_{E}S)^{2}}{\partial_{E}^{2}S}\, $$ must diverge in the thermodynamic limit (consistently with the definition of second-order transition). This divergence originates from this mechanism. Indeed, $\partial_{E}S$ manifests a stable inflection point (it does not move for increasing size) while the peak in $\partial_{E}^{2}S$ does. This originates the transition. <details> <summary>x3.png Details</summary>  ### Visual Description \n ## Multi-Panel Chart: 1D XY Mean Field Analysis ### Overview The image is a 3x3 grid of scientific plots analyzing the 1D XY mean field model across three different system sizes (N). The charts display the behavior of three different physical quantities as a function of the parameter ε (epsilon), likely representing an inverse temperature or coupling strength. The analysis compares results from three different methods: RICCATI, PEARSON, and MICRO AVERAGE. ### Components/Axes * **Title:** "1D XY mean field" (centered at top). * **Legend:** Positioned at the top center, below the title. * **RICCATI:** Represented by a black filled circle (●). * **PEARSON:** Represented by a red open circle (○). * **MICRO AVERAGE:** Represented by a black cross (×). * **Grid Structure:** * **Columns:** Correspond to different system sizes (N). From left to right: N = 8100, N = 14400, N = 40000. * **Rows:** Correspond to different measured quantities. From top to bottom: 1. **Top Row (a.1, b.1, c.1):** Y-axis label: `Y_g^(2)(ε)`. Y-axis scale: approximately -3 to 3. 2. **Middle Row (a.2, b.2, c.2):** Y-axis label: `∂_ε S(ε)`. Y-axis scale: approximately 1.0 to 2.5. 3. **Bottom Row (a.3, b.3, c.3):** Y-axis label: `∂_ε S^2(ε)`. Y-axis scale: approximately -7 to -2. * **Common X-Axis:** All nine plots share the same x-axis label `ε` (epsilon). The scale runs from approximately 0.6 to 1.0. * **Vertical Reference Lines:** Each column contains vertical dashed lines marking critical values of ε: * **Green dashed line:** Labeled `ε_c^MIPA ≈ 0.73` in column 1 (a.2), `ε_c^MIPA ≈ 0.73` in column 2 (b.2), and `ε_c^MIPA ≈ 0.74` in column 3 (c.2). * **Red dashed line:** Present in all columns. In column 2 (b.2), it is explicitly labeled `ε_c ≈ 0.75`. * **Black dashed line:** Present in all columns, typically to the right of the red line. ### Detailed Analysis **Row 1: Y_g^(2)(ε) vs. ε** * **Trend:** For all system sizes (N), the quantity `Y_g^(2)(ε)` shows a generally decreasing trend as ε increases from 0.6 to 1.0. The data points (primarily black crosses for MICRO AVERAGE) form a curve that slopes downward. * **Data Points:** The curve starts near a value of 3 at ε=0.6 and decreases to approximately -1 at ε=1.0. There is a noticeable discontinuity or sharp change in slope near the vertical reference lines (ε ≈ 0.73-0.75). **Row 2: ∂_ε S(ε) vs. ε** * **Trend:** This derivative shows a clear, monotonic decreasing trend across all panels. The slope is negative and becomes steeper as ε increases. * **Data Points:** The series (red open circles for PEARSON) begins at a value of approximately 2.5 at ε=0.6 and falls to approximately 1.0 at ε=1.0. The decline is smooth but exhibits a change in curvature near the critical ε region. **Row 3: ∂_ε S^2(ε) vs. ε** * **Trend:** This second derivative shows a clear, monotonic increasing trend (becoming less negative) across all panels. * **Data Points:** The series (red open circles for PEARSON) starts at a highly negative value of approximately -7 at ε=0.6 and rises to approximately -2 at ε=1.0. The increase is steep initially and then gradually flattens. **Critical Values (ε_c):** * The green-dashed `ε_c^MIPA` is consistently marked around 0.73-0.74. * The red-dashed `ε_c` is marked at approximately 0.75 in the N=14400 column. * The black-dashed line appears at a slightly higher ε value than the red line in each column. ### Key Observations 1. **Consistency Across Scales:** The qualitative trends for all three quantities (`Y_g^(2)`, `∂_ε S`, `∂_ε S^2`) are remarkably consistent across the three different system sizes (N=8100, 14400, 40000), suggesting robust physical behavior. 2. **Method Agreement:** The data points from the PEARSON (red circles) and MICRO AVERAGE (black crosses) methods appear to lie on the same continuous curves in each plot, indicating good agreement between these computational approaches for the measured quantities. 3. **Critical Region Signatures:** All three plotted quantities show distinct changes in behavior (slope, curvature) in the vicinity of ε ≈ 0.73-0.75, which is highlighted by the vertical reference lines. This region is identified as containing critical points (`ε_c`). 4. **Finite-Size Effect:** The exact position of the critical point `ε_c^MIPA` shows a very slight shift to higher ε values as the system size N increases (from ~0.73 to ~0.74), which is a common signature of finite-size effects in statistical physics. ### Interpretation This chart presents a finite-size scaling analysis of a phase transition in the 1D XY mean field model. The parameter ε likely controls the system's proximity to a critical point. * **What the data demonstrates:** The quantities plotted are likely related to the free energy or its derivatives. `∂_ε S(ε)` and `∂_ε S^2(ε)` resemble the first and second derivatives of an entropy (S) or free energy with respect to the control parameter ε. The smooth but rapidly changing behavior near ε_c is characteristic of a continuous (second-order) phase transition. The quantity `Y_g^(2)(ε)` may be a Binder cumulant or a similar dimensionless ratio used to locate critical points. * **Relationship between elements:** The vertical lines mark candidate critical points determined by different methods or criteria (MIPA, another denoted by `ε_c`). The close proximity of these lines (green, red, black) suggests they are different estimates for the same underlying phase transition. The consistency of trends across system sizes reinforces the physical reality of the transition. * **Notable patterns/anomalies:** The primary pattern is the clear, systematic evolution of all three thermodynamic quantities as a function of ε, with a focal point of change in the 0.73-0.75 range. There are no obvious data anomalies; the curves are smooth and well-defined. The slight rightward shift of `ε_c^MIPA` with increasing N is an expected finite-size effect, where the apparent critical point in a finite system converges to the true thermodynamic limit value as N → ∞. </details> Figure S2: Comparison of entropy derivatives at different system sizes obtained as the solution of the entropy flow equation and through thermodynamic methods. The quantitative and qualitative agreement is excellent also in the mean-field case. Here, the transition is sharper and the discontinuity is visible already at finite size. The geometric function $\Upsilon^{(2)}_{g}$ and the second derivative $\partial_{\epsilon}^{2}S$ show a jump around $\epsilon_{c}\approx 0.74$ closed to the infinite-size energy value $\epsilon_{c}^{\infty}=3J/4=0.75$ for our choice $J=1$ . Notice that the larger the size, the closer the critical point detected by MIPA. For the XY mean-field model, MIPA analysis works as well but the mechanism is different. This transition is related to a symmetry breaking associated to the magnetization order parameter and it manifests quite strongly already at finite size. The function $\partial_{\epsilon}^{2}S$ admits a negative-valued maximum in $\partial_{\epsilon}^{2}S$ but a visible jump appears around the infinity-size critical energy $\epsilon^{\infty}_{c}=3J/4=0.75$ . The first-derivative $\partial_{\epsilon}S$ admits a cuspid-like behavior (a drastic inflection point) at $\epsilon_{c}=0.73$ for $N=8100$ and $\epsilon_{c}=0.74$ for $N=14400$ and $N=40000$ . Interestingly, this maximum converges to the vertical green dashed line corresponding to $\epsilon^{\infty}_{c}$ for increasing size and the jump becomes a stronger and stronger discontinuity. This analysis will be reported more deeply in Ref. [52]. The finite-size MIPA analysis is in agreement with the infinite-size PT which is thus classified as a second-order PT.