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 ## Charts: Correlation Functions and Susceptibility ### Overview The image presents six sub-charts (a.1 through b.3) displaying correlation functions and susceptibility values as a function of a parameter denoted by 'ε' (epsilon). The charts appear to relate to a physical system, potentially magnetic or related to fluctuations, based on the labels. The charts are arranged in a 2x3 grid. ### Components/Axes * **X-axis (all charts):** ε (epsilon). The scales vary between charts. * **Y-axis:** Varies per chart, representing different correlation functions or susceptibility measures. * **Legends:** Each chart has a legend indicating the data series represented by different marker styles (dots, crosses). * **Labels:** * a.1: Y_g⁽²⁾(ε) * a.2: ∂_εS(ε) * a.3: ∂²_εS(ε) * b.1: Y_g⁽²⁾(ε) * b.2: ∂_εS(ε) * b.3: ∂²_εS(ε) * **Markers:** * 'x' - micro average * '•' - PHT * 'o' - EFE * **Vertical Lines:** Red dashed vertical lines are present in a.1 and b.1, indicating specific values of ε. * **Text Annotations:** * a.2: ε^MIPA ≈ 11.1 * b.2: ε^MIPA ≈ 0.74 * b.3: ε_c^∞ = 3J/4 ### Detailed Analysis or Content Details **a.1: Y_g⁽²⁾(ε)** * X-axis range: 9 to 14 * Y-axis range: approximately -0.002 to 0.002 * Data Series 1 (crosses): Starts at approximately (9, 0.0017), decreases to approximately (11, 0.0002), then increases to approximately (14, 0.0015). The line exhibits a minimum around ε = 11. * Vertical dashed line at approximately ε = 11.1. **a.2: ∂_εS(ε)** * X-axis range: 9 to 14 * Y-axis range: approximately 5.0 to 6.2 * Data Series 1 (dots): Starts at approximately (9, 5.8), decreases to approximately (11.1, 5.2), then increases to approximately (14, 5.4). * Markers: PHT at approximately (9, 5.7), EFE at approximately (9, 5.1). **a.3: ∂²_εS(ε)** * X-axis range: 9 to 14 * Y-axis range: approximately -2.2 to 0.001 * Data Series 1 (dots): Starts at approximately (9, -2.1), decreases to approximately (11.1, -2.2), then increases to approximately (14, -1.2). * Markers: PHT at approximately (9, -1.9), EFE at approximately (9, -2.0). **b.1: Y_g⁽²⁾(ε)** * X-axis range: 0.6 to 1.0 * Y-axis range: approximately -3 to 3 * Data Series 1 (crosses): Starts at approximately (0.6, 2.5), decreases to approximately (0.74, -0.5), then increases to approximately (1.0, 2.5). The line exhibits a minimum around ε = 0.74. * Vertical dashed line at approximately ε = 0.74. **b.2: ∂_εS(ε)** * X-axis range: 0.6 to 1.0 * Y-axis range: approximately 1.0 to 2.5 * Data Series 1 (dots): Starts at approximately (0.6, 2.2), decreases to approximately (0.74, 1.5), then increases to approximately (1.0, 2.0). * Markers: PHT at approximately (0.6, 2.1), EFE at approximately (0.6, 1.7). **b.3: ∂²_εS(ε)** * X-axis range: 0.6 to 1.0 * Y-axis range: approximately -7 to 1.0 * Data Series 1 (dots): Starts at approximately (0.6, -5.5), decreases to approximately (0.74, -6.5), then increases to approximately (1.0, -2.0). * Markers: PHT at approximately (0.6, -5.0), EFE at approximately (0.6, -6.0). * Text Annotation: ε_c^∞ = 3J/4 ### Key Observations * All charts show a similar trend: a peak or minimum around a specific value of ε, followed by a change in direction. * The vertical dashed lines in a.1 and b.1 highlight critical values of ε, potentially representing phase transitions or critical points. * The markers PHT and EFE appear to represent specific conditions or points of interest within the system. * The values of ε^MIPA are indicated in a.2 and b.2. * The value of ε_c^∞ is indicated in b.3. ### Interpretation The charts likely represent the behavior of a system's correlation functions and susceptibility as a function of a control parameter ε. The correlation functions (Y_g⁽²⁾(ε)) indicate how correlated fluctuations are at different scales. The susceptibility (∂_εS(ε) and ∂²_εS(ε)) measures the system's response to changes in ε. The peaks and minima in the charts suggest the presence of critical points or phase transitions. The vertical lines at ε ≈ 11.1 and ε ≈ 0.74 likely mark these critical points. The different data series (PHT and EFE) may represent different conditions or configurations of the system, and their positions relative to each other provide insights into how these conditions affect the system's behavior. The annotation ε_c^∞ = 3J/4 in b.3 suggests a connection to an infinite-range model, where J represents an interaction strength. The MIPA values (11.1 and 0.74) likely represent critical exponents or other parameters related to the system's critical behavior. The overall pattern suggests a system exhibiting a transition from one state to another as ε is varied, with the specific details depending on the system's parameters and conditions. </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 ## Chart: 2D φ⁴ nearest neighbours ### Overview The image presents a series of six subplots (a.1-a.3, b.1-b.3, c.1-c.3) displaying plots of different quantities related to a 2D φ⁴ nearest neighbours model. Each subplot corresponds to a different value of N (64², 128², 256²). The plots show relationships between variables ε (on the x-axis) and various functions of ε (on the y-axis). ### Components/Axes * **Title:** "2D φ⁴ nearest neighbours" (top-center) * **X-axis Label:** "ε" (appears on all subplots) * **Y-axis Labels:** * a.1, b.1, c.1: "Y(2ε)" * a.2, b.2, c.2: "∂₂S(ε)" * a.3, b.3, c.3: "∂₂S(ε)" * **N Values:** N = 64², N = 128², N = 256² (displayed above each set of three subplots) * **Legend:** * Black circles: "RICCATI" * Red circles: "PEARSON" * Black crosses: "MICRO AVERAGE" * **Vertical Lines:** Lines labeled "ε ≈ 11.2" (in a.2) and "ε ≈ 11.1" (in b.2 and c.2) are present in the subplots. * **X-axis Scale:** The x-axis (ε) ranges approximately from 9 to 15 in all subplots. * **Y-axis Scale:** The y-axis scales vary for each subplot, ranging from approximately -0.022 to 0.020 for a.3, -0.0022 to 0.022 for b.3, and -0.0018 to 0.002 for c.3. ### Detailed Analysis Let's analyze each subplot and data series: **a.1 (N = 64²):** * **RICCATI (Black circles):** The curve starts at approximately Y(2ε) = 0.002 at ε ≈ 9, reaches a maximum of approximately 0.014 at ε ≈ 11.5, and then decreases to approximately 0.001 at ε ≈ 15. * **PEARSON (Red circles):** The curve starts at approximately Y(2ε) = 0.018 at ε ≈ 9, decreases monotonically to approximately 0.002 at ε ≈ 15. * **MICRO AVERAGE (Black crosses):** The curve starts at approximately Y(2ε) = 0.001 at ε ≈ 9, reaches a maximum of approximately 0.008 at ε ≈ 11.5, and then decreases to approximately 0.0005 at ε ≈ 15. **a.2 (N = 64²):** * **RICCATI (Black circles):** The curve starts at approximately ∂₂S(ε) = -0.055 at ε ≈ 9, increases to a maximum of approximately 0.060 at ε ≈ 11.2, and then decreases to approximately -0.050 at ε ≈ 15. * **PEARSON (Red circles):** The curve starts at approximately ∂₂S(ε) = -0.050 at ε ≈ 9, increases to a maximum of approximately 0.055 at ε ≈ 11.2, and then decreases to approximately -0.045 at ε ≈ 15. * **MICRO AVERAGE (Black crosses):** The curve starts at approximately ∂₂S(ε) = -0.052 at ε ≈ 9, increases to a maximum of approximately 0.058 at ε ≈ 11.2, and then decreases to approximately -0.048 at ε ≈ 15. **a.3 (N = 64²):** * **RICCATI (Black circles):** The curve starts at approximately ∂₂S(ε) = -0.022 at ε ≈ 9, increases to a maximum of approximately 0.002 at ε ≈ 11.2, and then decreases to approximately -0.020 at ε ≈ 15. * **PEARSON (Red circles):** The curve starts at approximately ∂₂S(ε) = -0.020 at ε ≈ 9, increases to a maximum of approximately 0.001 at ε ≈ 11.2, and then decreases to approximately -0.018 at ε ≈ 15. * **MICRO AVERAGE (Black crosses):** The curve starts at approximately ∂₂S(ε) = -0.021 at ε ≈ 9, increases to a maximum of approximately 0.0015 at ε ≈ 11.2, and then decreases to approximately -0.019 at ε ≈ 15. **b.1, b.2, b.3 (N = 128²):** The trends are similar to the N = 64² case, but the curves appear to be shifted slightly and have different amplitudes. **c.1, c.2, c.3 (N = 256²):** The trends are again similar, with further shifts and amplitude changes. The y-axis scales are smaller, indicating smaller values for the plotted functions. ### Key Observations * All three methods (RICCATI, PEARSON, MICRO AVERAGE) show similar trends across all values of N. * As N increases (from 64² to 256²), the magnitude of the plotted functions generally decreases. * The vertical lines at ε ≈ 11.1 and ε ≈ 11.2 seem to mark points of interest in the curves, potentially indicating critical points or singularities. * The PEARSON method consistently shows a monotonic decrease in Y(2ε) in the a.1, b.1, and c.1 subplots. ### Interpretation The plots likely represent the behavior of different approximations or solutions to the 2D φ⁴ model as a function of a parameter ε. The different methods (RICCATI, PEARSON, MICRO AVERAGE) represent different ways of solving or approximating the model. The parameter N likely relates to the number of samples or grid points used in the approximation. The decreasing magnitude of the functions as N increases suggests that the approximations become more accurate or converge as the number of samples increases. The vertical lines at ε ≈ 11.1 and ε ≈ 11.2 may indicate a phase transition or a critical point in the system, where the behavior of the model changes significantly. The consistent monotonic decrease of the PEARSON method in the Y(2ε) plots suggests that it may be a more stable or conservative approximation. The data suggests that the 2D φ⁴ nearest neighbours model exhibits complex behavior as a function of ε, and that different approximation methods can provide insights into this behavior. The choice of approximation method and the number of samples used can significantly impact the results. </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 ## Chart: ID XY mean field ### Overview The image presents a comparative analysis of three different methods – Riccati, Pearson, and Micro Average – within an "ID XY mean field" context. The analysis is visualized through three sets of subplots (a, b, c), each containing three smaller plots (numbered .1, .2, .3) displaying relationships between variables. The plots appear to explore the behavior of these methods across a range of epsilon (ε) values. ### Components/Axes * **Title:** "ID XY mean field" (centered at the top) * **Legend:** Located at the top-right corner, identifying the data series: * Riccati (represented by black circles) * Pearson (represented by orange circles) * Micro Average (represented by black crosses) * **X-axis:** ε (epsilon), ranging from approximately 0.6 to 1.0, with tick marks at 0.6, 0.7, 0.8, 0.9, and 1.0. * **Y-axis (subplots .1):** γ₂(ε) (gamma 2 of epsilon), ranging from approximately -3 to 3, with tick marks at -3, -2, -1, 0, 1, 2, and 3. * **Y-axis (subplots .2):** δ₂S(ε) (delta 2 S of epsilon), ranging from approximately -7 to 2.5, with tick marks at -7, -6, -5, -4, -3, -2, -1, 0, 1, 2, and 2.5. * **Y-axis (subplots .3):** δS²(ε) (delta S squared of epsilon), ranging from approximately -7 to 3, with tick marks at -7, -6, -5, -4, -3, -2, -1, 0, 1, 2, and 3. * **Vertical Dashed Lines:** Three vertical dashed lines are present, marking ε values of approximately 0.75 (in subplot b), 0.73 (in subplot a), and 0.74 (in subplot c). * **Text Annotations:** * "N = 8100" (top-left, above subplot a) * "N = 14400" (center-top, above subplot b) * "N = 40000" (top-right, above subplot c) * "ε_C ≈ 0.75" (within subplot b) * "ε_MIPA ≈ 0.73" (within subplot a) * "ε_MIPA ≈ 0.74" (within subplot c) * "a.1", "a.2", "a.3", "b.1", "b.2", "b.3", "c.1", "c.2", "c.3" (labels for each subplot) ### Detailed Analysis or Content Details **Subplot a (N = 8100):** * **Riccati (black circles):** The line in a.1 slopes downward, starting at approximately γ₂(ε) = 2.7 at ε = 0.6 and decreasing to approximately γ₂(ε) = -2.5 at ε = 1.0. In a.2, the line starts at approximately δ₂S(ε) = 1.2 at ε = 0.6 and decreases to approximately δ₂S(ε) = -1.5 at ε = 1.0. In a.3, the line starts at approximately δS²(ε) = -1.5 at ε = 0.6 and decreases to approximately δS²(ε) = -6.5 at ε = 1.0. * **Pearson (orange circles):** The line in a.1 is relatively flat, fluctuating around γ₂(ε) = 0. In a.2, the line is also relatively flat, fluctuating around δ₂S(ε) = 0. In a.3, the line is relatively flat, fluctuating around δS²(ε) = 0. * **Micro Average (black crosses):** The line in a.1 slopes downward, starting at approximately γ₂(ε) = 2.5 at ε = 0.6 and decreasing to approximately γ₂(ε) = -2.5 at ε = 1.0. In a.2, the line starts at approximately δ₂S(ε) = 1.0 at ε = 0.6 and decreases to approximately δ₂S(ε) = -2.0 at ε = 1.0. In a.3, the line starts at approximately δS²(ε) = -2.0 at ε = 0.6 and decreases to approximately δS²(ε) = -6.0 at ε = 1.0. **Subplot b (N = 14400):** * **Riccati (black circles):** The line in b.1 slopes downward, starting at approximately γ₂(ε) = 2.5 at ε = 0.6 and decreasing to approximately γ₂(ε) = -2.5 at ε = 1.0. In b.2, the line starts at approximately δ₂S(ε) = 1.2 at ε = 0.6 and decreases to approximately δ₂S(ε) = -1.5 at ε = 1.0. In b.3, the line starts at approximately δS²(ε) = -1.5 at ε = 0.6 and decreases to approximately δS²(ε) = -6.5 at ε = 1.0. * **Pearson (orange circles):** The line in b.1 is relatively flat, fluctuating around γ₂(ε) = 0. In b.2, the line is also relatively flat, fluctuating around δ₂S(ε) = 0. In b.3, the line is relatively flat, fluctuating around δS²(ε) = 0. * **Micro Average (black crosses):** The line in b.1 slopes downward, starting at approximately γ₂(ε) = 2.5 at ε = 0.6 and decreasing to approximately γ₂(ε) = -2.5 at ε = 1.0. In b.2, the line starts at approximately δ₂S(ε) = 1.0 at ε = 0.6 and decreases to approximately δ₂S(ε) = -2.0 at ε = 1.0. In b.3, the line starts at approximately δS²(ε) = -2.0 at ε = 0.6 and decreases to approximately δS²(ε) = -6.0 at ε = 1.0. **Subplot c (N = 40000):** * **Riccati (black circles):** The line in c.1 slopes downward, starting at approximately γ₂(ε) = 2.7 at ε = 0.6 and decreasing to approximately γ₂(ε) = -2.5 at ε = 1.0. In c.2, the line starts at approximately δ₂S(ε) = 1.2 at ε = 0.6 and decreases to approximately δ₂S(ε) = -1.5 at ε = 1.0. In c.3, the line starts at approximately δS²(ε) = -1.5 at ε = 0.6 and decreases to approximately δS²(ε) = -6.5 at ε = 1.0. * **Pearson (orange circles):** The line in c.1 is relatively flat, fluctuating around γ₂(ε) = 0. In c.2, the line is also relatively flat, fluctuating around δ₂S(ε) = 0. In c.3, the line is relatively flat, fluctuating around δS²(ε) = 0. * **Micro Average (black crosses):** The line in c.1 slopes downward, starting at approximately γ₂(ε) = 2.5 at ε = 0.6 and decreasing to approximately γ₂(ε) = -2.5 at ε = 1.0. In c.2, the line starts at approximately δ₂S(ε) = 1.0 at ε = 0.6 and decreases to approximately δ₂S(ε) = -2.0 at ε = 1.0. In c.3, the line starts at approximately δS²(ε) = -2.0 at ε = 0.6 and decreases to approximately δS²(ε) = -6.0 at ε = 1.0. ### Key Observations * The Riccati and Micro Average methods exhibit similar behavior across all subplots, with downward-sloping lines in each plot. * The Pearson method consistently shows minimal variation, with lines fluctuating around zero in all plots. * The vertical dashed lines indicate critical epsilon values for each method, with slight variations between them. * The value of N (sample size) increases from subplot a to c, potentially influencing the observed trends. ### Interpretation The chart suggests that the Riccati and Micro Average methods are sensitive to changes in epsilon (ε), exhibiting a clear relationship between ε and the plotted variables (γ₂, δ₂S, δS²). The Pearson method, however, appears to be largely unaffected by variations in ε. The increasing value of N across the subplots might indicate a convergence of the Riccati and Micro Average methods as the sample size grows. The vertical lines highlight specific epsilon values where the methods may exhibit unique behavior or transitions. The consistent downward trend in the Riccati and Micro Average methods suggests a potential inverse relationship between epsilon and the plotted variables. The differences in the epsilon values marked by the vertical lines (0.73, 0.74, 0.75) could represent critical points or phase transitions in the system being modeled. The data suggests that the Pearson method is a more stable estimator, while the Riccati and Micro Average methods are more sensitive to the underlying data and parameter variations. </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.