2512.02242

# Phase Transitions as Emergent Geometric Phenomena A Deterministic Entropy Evolution Law **Authors**: Loris Di Cairano > l.di.cairano.92@gmail.com, loris.dicairano@uni.luDepartment 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}\!·s dq^{N}\,dp_{1}\!·s 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(∇^{\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}∈\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∈ T_{\bm{x}}\Lambda:\eta(Y,∇^{\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 $∇_{\eta}H$ and generates the motion of points on an energy hypersurface to others. In general, $∇_{\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}∈\mathrm{span}\{∇^{\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}→\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}\!·s 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∈[\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 ## Chart/Diagram Type: Dual-Axis Graphs with Subplots ### Overview The image contains two side-by-side graphs (labeled **a** and **b**), each divided into three subplots (a.1, a.2, a.3 and b.1, b.2, b.3). These graphs depict mathematical functions and their derivatives as functions of a parameter **ε** (epsilon). Key elements include legends, vertical dashed lines, and annotations for critical points. --- ### Components/Axes #### Left Graph (a): - **a.1**: - **X-axis**: **ε** (epsilon), ranging from 9 to 14. - **Y-axis**: **V_φ4** (top plot), **Γ_g^(2)(ε)** (middle plot), **δ_εS(ε)** (bottom plot). - **Legend**: - **PHT**: Black dots (bottom-left corner). - **EFE**: Red circles (bottom-left corner). - **Annotations**: - Vertical red dashed line at **ε_c^MIPA ≈ 11.1**. - Micro average marked with a cross (x) at **ε ≈ 10.5**. - **a.2**: - **X-axis**: **ε** (epsilon), same range as a.1. - **Y-axis**: **δ_εS(ε)** (top plot), **δ²_εS(ε)** (bottom plot). - **a.3**: - **X-axis**: **ε** (epsilon), same range as a.1. - **Y-axis**: **δ²_εS(ε)** (bottom plot). #### Right Graph (b): - **b.1**: - **X-axis**: **ε** (epsilon), ranging from 0.6 to 1.0. - **Y-axis**: **V_XY** (top plot), **Γ_g^(2)(ε)** (middle plot). - **Legend**: - **PHT**: Black dots (bottom-left corner). - **EFE**: Red circles (bottom-left corner). - **Annotations**: - Vertical green dashed line at **ε_c^∞ = 3J/4**. - Micro average marked with a cross (x) at **ε ≈ 0.75**. - **b.2**: - **X-axis**: **ε** (epsilon), same range as b.1. - **Y-axis**: **δ_εS(ε)** (top plot), **δ²_εS(ε)** (bottom plot). - **b.3**: - **X-axis**: **ε** (epsilon), same range as b.1. - **Y-axis**: **δ²_εS(ε)** (bottom plot). --- ### Detailed Analysis #### Left Graph (a): 1. **a.1**: - **V_φ4**: Peaks near **ε ≈ 10.5**, then declines. - **Γ_g^(2)(ε)**: Dips sharply at **ε ≈ 10.5**, then rises. - **δ_εS(ε)**: Decreases monotonically from **ε = 9** to **ε = 14**. - **δ²_εS(ε)**: Shows a sharp minimum at **ε ≈ 10.5**, with values dropping to **-2.2**. 2. **a.2**: - **δ_εS(ε)**: Decreases from **6.0** at **ε = 9** to **5.0** at **ε = 14**. - **δ²_εS(ε)**: Peaks at **ε ≈ 10.5** (value ~-1.6), then declines. 3. **a.3**: - **δ²_εS(ε)**: Minimum at **ε ≈ 10.5** (value ~-2.2), with a secondary minimum at **ε ≈ 14**. #### Right Graph (b): 1. **b.1**: - **V_XY**: Peaks at **ε ≈ 0.7**, then declines. - **Γ_g^(2)(ε)**: Rises sharply at **ε ≈ 0.7**, then plateaus. - **δ_εS(ε)**: Decreases from **2.5** at **ε = 0.6** to **1.0** at **ε = 1.0**. - **δ²_εS(ε)**: Peaks at **ε ≈ 0.7** (value ~-5), then declines. 2. **b.2**: - **δ_εS(ε)**: Decreases from **2.5** at **ε = 0.6** to **1.0** at **ε = 1.0**. - **δ²_εS(ε)**: Peaks at **ε ≈ 0.7** (value ~-5), then declines. 3. **b.3**: - **δ²_εS(ε)**: Minimum at **ε ≈ 0.7** (value ~-5), with a secondary minimum at **ε ≈ 1.0**. --- ### Key Observations 1. **Critical Points**: - Left graph: **ε_c^MIPA ≈ 11.1** (red dashed line) marks a phase transition or critical threshold. - Right graph: **ε_c^∞ = 3J/4** (green dashed line) indicates a theoretical limit. 2. **Micro Averages**: - Crosses (x) in both graphs highlight averaged values at **ε ≈ 10.5** (left) and **ε ≈ 0.75** (right). 3. **Function Behavior**: - **V_φ4** and **V_XY** exhibit peaks near their respective critical points. - **δ²_εS(ε)** in both graphs shows minima at critical points, suggesting stability or optimal conditions. 4. **Legend Consistency**: - **PHT** (black dots) and **EFE** (red circles) are consistently placed in the bottom-left corner of all subplots. --- ### Interpretation The graphs likely represent a physical or mathematical model where **ε** governs system behavior. The critical points (**ε_c^MIPA** and **ε_c^∞**) demarcate distinct regimes: - **Left Graph (a)**: Dominated by **V_φ4** and **Γ_g^(2)(ε)**, with **δ²_εS(ε)** minima indicating stability at **ε ≈ 10.5**. - **Right Graph (b)**: Governed by **V_XY** and **Γ_g^(2)(ε)**, with **δ²_εS(ε)** minima at **ε ≈ 0.7**, suggesting a phase transition or optimal parameter range. The equations embedded in the plots (e.g., **V_φ4 = Σ[μ²φ_i²/2 + λ/4!φ_i⁴] + J/2 Σ(φ_i - φ_{i+ε})²**) imply interactions between variables, possibly in a lattice or network model. The **PHT** and **EFE** markers may represent different experimental or theoretical datasets, with **EFE** showing sharper transitions near critical points. The **δ²_εS(ε)** minima in both graphs highlight regions of minimal variance, critical for understanding system stability or phase behavior. The **micro averages** (crosses) suggest averaged trends over a range of **ε**, contrasting with the detailed data points (dots and circles). This analysis underscores the interplay between **ε**, critical thresholds, and system stability, with potential applications in statistical mechanics, material science, or complex systems theory. </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^{∞}_{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}=∇^{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)=∂_{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∂_{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}[∇_{\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 $∂_{E}S_{g}$ and $∂_{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, $∂_{E}S_{g}$ displays an inflection point (panel a.3 and $∂_{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}^{∞}\!=\!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 $∂_{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→∞$ , 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}}∈\Sigma_{E}$ and consider the restriction $\omega|_{T_{\bm{x}}\Sigma_{E}}$ , which we denote by $\omega_{E}$ . For any vector $v∈ 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∈ 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}∈ T_{\bm{x}}\Sigma_{E}$ , we have $\bm{X}_{H}∈\ker(\omega_{E})$ . The kernel is non-trivial; hence, $\omega_{E}$ is degenerate. Moreover, $X_{H}(x)∈ 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}}⊂eq 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}∈(T_{\bm{x}}\Sigma_{E})^{\perp_{\omega}}$ . For any $\bm{v}∈ 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}∈(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≠ 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}}⊂eq T_{x}\Sigma_{E}$ . We have already shown that $X_{H}∈ 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}∈\Lambda$ , we seek vectors $\bm{\zeta}∈ 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}=∂_{p_{1}}/(∂ H/∂ p_{1})$ if $∂ H/∂ p_{1}≠ 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 $∇ 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 $∇ 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 $\|∇ H\|$ varies across $\Sigma_{E}$ , for two distinct points $\bm{x}≠\bm{y}$ on the same $\Sigma_{E}$ , we have $dE_{\bm{x}}≠ dE_{\bm{y}}$ , mapping them to different target hypersurfaces. To resolve this, it is sufficient to determine the vector field $\bm{\zeta}∈\text{span}\{∇ 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,...,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(∂_{\alpha}x)=0$ and $\bm{\zeta}=∇_{\eta}H/\|∇_{\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}:=∂_{\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(∂_{E})=1$ , $dE(∂_{\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⊂ M$ , the area form is determined by the unit normal vector field $\bm{n}$ : ( $\|n\|_{g}=1$ , $g(n,T)=0$ for all $T∈ 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},...,u^{m-1},u^{m})$ with $\Sigma=\{u^{m}=0\}$ it reduces to the usual expression $d\sigma^{g}=\sqrt{\det g_{→p}}\,du^{1}\wedge·s\wedge du^{m-1}$ , where $g_{→p}$ is the metric induced on $T\Sigma$ . Then, the area of a Borel set $A⊂\Sigma$ is $\mathrm{area}_{g}(A)=∈t_{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∈[1,2n-1]$ . In these new coordinates $(x^{0},x^{1},...,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∂_{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}=∇_{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 $∇_{g}H$ on $T_{\bm{x}}\Sigma_{E}$ must vanish: $∂ H/∂ 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 $∂_{E}=\mathscr{L}_{\bm{\xi}}=(\bm{\xi}\!·\!∇)=(u/G)\!·\!∇$ , we can compute $∂_{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,...,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→∞$ ), 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\%-60\%$ . 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}}∉[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∈[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}≤ E_{1}<...<E_{M}≤ 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 $→$ 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}∈[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 ## Chart/Diagram Type: 2D φ⁴ Nearest Neighbours Analysis ### Overview The image presents three panels (a, b, c) comparing three mathematical functions (Riccati, Pearson, Micro Average) across varying system sizes (N = 64², 128², 256²). Each panel contains three subplots (1, 2, 3) showing distinct metrics: 1. **Top subplot**: Y_g^(2)(ε) (second-order correlation function) 2. **Middle subplot**: δε(ε) (energy density derivative) 3. **Bottom subplot**: ∂_ε S²(ε) (entropy derivative) A vertical red dashed line marks the critical value ε_c ≈ 11.1 across all panels. ### Components/Axes - **X-axis**: ε (parameter, range: 9–14) - **Y-axes**: - Top: Y_g^(2)(ε) (0–0.002) - Middle: δε(ε) (0.05–0.06) - Bottom: ∂_ε S²(ε) (-0.0022–0.0005) - **Legends**: - **Riccati**: Black dots (dotted line) - **Pearson**: Red circles (dashed line) - **Micro Average**: Crosses (not explicitly plotted in visible data) ### Detailed Analysis #### Panel a (N = 64²) - **a.1 (Y_g^(2)(ε))**: - Riccati (dotted): Peaks at ε ≈ 11.1, then decays. - Pearson (dashed): Monotonically decreasing from ε = 9 to 14. - Critical line: ε_c ≈ 11.1 (red dashed). - **a.2 (δε(ε))**: - Pearson (dashed): Linear decline from ~0.06 (ε=9) to ~0.05 (ε=14). - Riccati (dotted): Not visible (likely negligible). - **a.3 (∂_ε S²(ε))**: - Pearson (dashed): Sharp peak at ε ≈ 11.1, then decays. #### Panel b (N = 128²) - **b.1 (Y_g^(2)(ε))**: - Riccati (dotted): Broader peak centered at ε ≈ 11.1. - Pearson (dashed): Slightly smoother decline than panel a. - **b.2 (δε(ε))**: - Pearson (dashed): Steeper slope than panel a, dropping from ~0.06 to ~0.052. - **b.3 (∂_ε S²(ε))**: - Pearson (dashed): Narrower peak at ε ≈ 11.1 compared to panel a. #### Panel c (N = 256²) - **c.1 (Y_g^(2)(ε))**: - Riccati (dotted): Sharper peak at ε ≈ 11.1, narrower than panel b. - Pearson (dashed): Minimal deviation from panel b. - **c.2 (δε(ε))**: - Pearson (dashed): Slightly reduced slope, ending at ~0.051. - **c.3 (∂_ε S²(ε))**: - Pearson (dashed): Peak height reduced compared to panel b. ### Key Observations 1. **Critical Value Consistency**: ε_c ≈ 11.1 is marked identically across all panels, suggesting it is a universal parameter. 2. **Riccati Behavior**: Peaks align with ε_c, indicating a phase transition or critical point. 3. **Pearson Trends**: Monotonic decline in all subplots, with magnitude decreasing as N increases. 4. **Micro Average**: Not explicitly plotted; inferred to lie between Riccati and Pearson curves. 5. **System Size Effects**: Larger N (256²) sharpens peaks in Riccati and reduces Pearson magnitudes. ### Interpretation - **Critical Phenomena**: The alignment of Riccati peaks with ε_c ≈ 11.1 suggests this parameter marks a critical threshold (e.g., phase transition) in the 2D φ⁴ model. - **Pearson Decay**: The linear decline of δε(ε) and ∂_ε S²(ε) implies diminishing energy density and entropy sensitivity as ε increases, possibly reflecting system stabilization. - **Micro Average Ambiguity**: The absence of explicit Micro Average data (crosses) leaves its role unclear, but its placement between Riccati and Pearson suggests it may represent an intermediate regime. - **Scaling Effects**: Larger system sizes (N = 256²) enhance precision in locating ε_c and refine the behavior of Riccati/Pearson functions, indicating improved convergence with system size. This analysis highlights the interplay between mathematical functionals and critical parameters in statistical physics models, with implications for understanding phase transitions and criticality. </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 $∂_{\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 $∂_{\epsilon}^{2}S$ together with an inflection point in $∂_{\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, $∂_{E}S$ manifests a stable inflection point (it does not move for increasing size) while the peak in $∂_{E}^{2}S$ does. This originates the transition. <details> <summary>x3.png Details</summary>  ### Visual Description ## 1D XY Mean Field Graphs: Critical Behavior Analysis ### Overview The image presents three panels (a, b, c) corresponding to system sizes *N* = 8100, 14400, and 40000, respectively. Each panel contains three sub-panels (1, 2, 3) plotting distinct thermodynamic quantities against the coupling parameter ε. The graphs compare three theoretical models (RICCATI, PEARSON, MICRO AVERAGE) across critical thresholds marked by vertical lines at ε = 0.75 (green) and ε = 0.8 (red). Key critical exponents and MIPA values are annotated. --- ### Components/Axes - **X-axis**: Coupling parameter ε (ranging from 0.6 to 1.0 across panels). - **Y-axes**: - **Panel a**: - Subpanel 1: Γ_g^(2)(ε) (normalized pair correlation function). - Subpanel 2: ε_c S_c(ε) (scaled susceptibility). - Subpanel 3: ∂_ε S_c^2(ε) (derivative of squared susceptibility). - **Panels b and c**: Identical y-axis labels as panel a. - **Legends**: - **RICCATI**: Black dots (•). - **PEARSON**: Red circles (○). - **MICRO AVERAGE**: Black crosses (×). - **Vertical Lines**: - Green dashed line at ε = 0.75 (critical threshold). - Red dashed line at ε = 0.8 (secondary threshold). --- ### Detailed Analysis #### Panel a (N = 8100) - **a.1 (Γ_g^(2)(ε))**: - RICCATI (black dots) shows a flat trend near ε = 0.8. - PEARSON (red circles) exhibits a sharp drop at ε ≈ 0.75, then rises. - MICRO AVERAGE (black crosses) remains constant. - **a.2 (ε_c S_c(ε))**: - All models converge to ε_c^MIPA ≈ 0.73 (annotated). - PEARSON deviates slightly below ε = 0.75. - **a.3 (∂_ε S_c^2(ε))**: - RICCATI (black dots) shows a sharp peak at ε ≈ 0.75. - PEARSON (red circles) has a broad peak spanning ε = 0.7–0.8. - MICRO AVERAGE (black crosses) remains flat. #### Panel b (N = 14400) - **b.1 (Γ_g^(2)(ε))**: - RICCATI (black dots) flattens earlier (ε ≈ 0.78). - PEARSON (red circles) peaks at ε ≈ 0.75, then declines. - MICRO AVERAGE (black crosses) remains stable. - **b.2 (ε_c S_c(ε))**: - ε_c^MIPA ≈ 0.73 (same as panel a). - PEARSON aligns closely with critical threshold. - **b.3 (∂_ε S_c^2(ε))**: - RICCATI (black dots) peaks sharply at ε ≈ 0.75. - PEARSON (red circles) shows a broader peak (ε = 0.7–0.8). - MICRO AVERAGE (black crosses) remains flat. #### Panel c (N = 40000) - **c.1 (Γ_g^(2)(ε))**: - RICCATI (black dots) flattens at ε ≈ 0.76. - PEARSON (red circles) peaks at ε ≈ 0.75, then declines. - MICRO AVERAGE (black crosses) remains constant. - **c.2 (ε_c S_c(ε))**: - ε_c^MIPA ≈ 0.74 (slightly higher than smaller N). - PEARSON aligns with critical threshold. - **c.3 (∂_ε S_c^2(ε))**: - RICCATI (black dots) peaks sharply at ε ≈ 0.75. - PEARSON (red circles) shows a broad peak (ε = 0.7–0.8). - MICRO AVERAGE (black crosses) remains flat. --- ### Key Observations 1. **Critical Thresholds**: - Green line (ε = 0.75) marks the primary phase transition. - Red line (ε = 0.8) indicates a secondary threshold where RICCATI models flatten. 2. **Model Behavior**: - **RICCATI**: Peaks sharply at ε ≈ 0.75 (a.3, b.3, c.3) and flattens near ε = 0.8. - **PEARSON**: Exhibits broad peaks (ε = 0.7–0.8) in susceptibility derivatives (a.3, b.3, c.3), suggesting slower convergence. - **MICRO AVERAGE**: Remains flat across all panels, indicating mean-field stability. 3. **MIPA Values**: - ε_c^MIPA ≈ 0.73 (panels a, b) and 0.74 (panel c), showing slight N-dependence. 4. **Anomalies**: - PEARSON data points deviate from critical thresholds in a.1 and b.1, possibly due to model limitations or noise. --- ### Interpretation - **Phase Transition Dynamics**: The convergence of ε_c^MIPA values (0.73–0.74) across N suggests a universal critical point, with larger N refining the estimate. - **Model Discrepancies**: PEARSON’s broad peaks in susceptibility derivatives (a.3, b.3, c.3) contrast with RICCATI’s sharp peaks, highlighting differences in handling long-range correlations. - **N-Scaling Effects**: Larger N (panel c) shifts ε_c^MIPA slightly upward, indicating improved accuracy in critical exponent estimation. - **Thermodynamic Stability**: MICRO AVERAGE’s flat trends suggest mean-field theory remains valid across ε ranges, while RICCATI and PEARSON capture finer critical behavior. This analysis underscores the importance of model selection in capturing critical phenomena, with RICCATI and PEARSON offering complementary insights into phase transitions. </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 $∂_{\epsilon}^{2}S$ show a jump around $\epsilon_{c}≈ 0.74$ closed to the infinite-size energy value $\epsilon_{c}^{∞}=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 $∂_{\epsilon}^{2}S$ admits a negative-valued maximum in $∂_{\epsilon}^{2}S$ but a visible jump appears around the infinity-size critical energy $\epsilon^{∞}_{c}=3J/4=0.75$ . The first-derivative $∂_{\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^{∞}_{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.