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 ## Combined Chart: Energy Landscape Analysis ### Overview The image presents six plots arranged in two columns (left and right), each containing three subplots (top, middle, bottom). The plots depict the energy landscape analysis, showing the relationships between energy (epsilon) and various functions related to the system's state. The left column (a.1, a.2, a.3) focuses on one set of parameters, while the right column (b.1, b.2, b.3) focuses on another. Each subplot displays data points from two methods, PHT (solid black circles) and EFE (open red circles), along with a "micro average" (black x markers). Vertical lines indicate critical energy values. ### Components/Axes **General Layout:** * Two columns of plots, labeled 'a' (left) and 'b' (right). * Each column has three subplots, labeled '1', '2', and '3' from top to bottom. * Each subplot contains data from two methods: PHT (black filled circles) and EFE (red open circles). * Each subplot contains data from "micro average" (black x markers). * Legend: PHT (black filled circle), EFE (red open circle), micro average (black x marker) **Left Column (a):** * X-axis (all subplots): ε (epsilon), ranging from approximately 9 to 14. * a.1: Y-axis: γg(2)(ε), ranging from 0 to 2.0 x 10^-3. * a.2: Y-axis: ∂εS(ε), ranging from 5.0 x 10^-2 to 6.2 x 10^-2. * a.3: Y-axis: ∂ε2S(ε), ranging from -2.2 x 10^-3 to -1.2 x 10^-3. * Vertical dashed red line at εcMIPA ≈ 11.1. **Right Column (b):** * X-axis (all subplots): ε (epsilon), ranging from approximately 0.6 to 1.0. * b.1: Y-axis: γg(2)(ε), ranging from -3 to 3. * b.2: Y-axis: ∂εS(ε), ranging from 1.0 to 2.5. * b.3: Y-axis: ∂ε2S(ε), ranging from -7 to -2. * Vertical dashed green line at εc∞ = 3J/4, approximately at ε ≈ 0.74. * Horizontal line at y = 0 in b.1 **Equations (Top of Image):** * Left: Vφ4 := Σi [μ2/2 φi2 + λ/4! φi4] + J/2 Σμ=1n (φi - φi+eμ)2 * Right: VXY = J/2N Σij [1 - cos(θi - θj)] ### Detailed Analysis **Left Column (a):** * **a.1: γg(2)(ε):** The "micro average" data (black x markers) initially decreases slightly from ε = 9 to approximately ε = 10.5, then sharply increases to a peak near the vertical red line (ε ≈ 11.1), and then decreases steadily as ε increases to 14. * **a.2: ∂εS(ε):** Both PHT (black filled circles) and EFE (red open circles) data points form a single line. The line slopes downward linearly from approximately (9, 0.062) to (14, 0.050). The value of εcMIPA is indicated as approximately 11.1. * **a.3: ∂ε2S(ε):** Both PHT (black filled circles) and EFE (red open circles) data points form a single line. The line starts at approximately -2.2 x 10^-3 at ε = 9, rises sharply to a peak near the vertical red line (ε ≈ 11.1), and then decreases steadily as ε increases to 14, approaching -2.2 x 10^-3. **Right Column (b):** * **b.1: γg(2)(ε):** The "micro average" data (black x markers) starts at approximately 3 at ε = 0.6, decreases to a minimum near ε = 0.8, and then increases as ε approaches 1.0. * **b.2: ∂εS(ε):** Both PHT (black filled circles) and EFE (red open circles) data points form a single line. The line slopes downward linearly from approximately (0.6, 2.5) to (1.0, 1.1). The value of εcMIPA is indicated as approximately 0.74, and εc∞ = 3J/4 is indicated by a vertical green line. * **b.3: ∂ε2S(ε):** Both PHT (black filled circles) and EFE (red open circles) data points form a single line. The line starts at approximately -7 at ε = 0.6, rises sharply to a peak near the vertical green line (ε ≈ 0.74), and then decreases steadily as ε increases to 1.0, approaching -2. ### Key Observations * The plots show the behavior of different functions with respect to energy (epsilon). * The vertical lines in both columns indicate critical energy values (εcMIPA and εc∞). * The PHT and EFE methods produce nearly identical results for ∂εS(ε) and ∂ε2S(ε). * The "micro average" data (black x markers) in γg(2)(ε) shows distinct trends in both columns. * The functions ∂εS(ε) show a linear decrease with increasing epsilon in both columns. * The functions ∂ε2S(ε) show a peak near the critical energy values in both columns. ### Interpretation The data presented in these plots likely represents a theoretical analysis of a physical system, possibly related to phase transitions or critical phenomena. The functions γg(2)(ε), ∂εS(ε), and ∂ε2S(ε) likely represent different physical quantities or derivatives of a physical quantity with respect to energy. The critical energy values (εcMIPA and εc∞) indicate points where the system undergoes a significant change in behavior. The agreement between the PHT and EFE methods suggests that these methods are consistent in their predictions. The trends observed in the plots provide insights into the system's behavior near the critical points. The peaks in ∂ε2S(ε) near the critical energy values suggest that the system's sensitivity to energy changes is highest at these points. </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: 2D φ⁴ Nearest Neighbours ### Overview The image presents a series of line plots illustrating the behavior of different functions with respect to a variable epsilon (ε). The plots are arranged in a 3x3 grid, with each column representing a different value of N (64², 128², 256²). Each row represents a different function: Yg^(2)(ε), ∂εS(ε), and ∂εS²(ε). The plots compare results from RICCATI, PEARSON, and MICRO AVERAGE methods. ### Components/Axes * **Title:** 2D φ⁴ nearest neighbours * **Legend:** Located at the top of the image. * Black filled circle: RICCATI * Red open circle: PEARSON * Black cross: MICRO AVERAGE * **X-axis:** ε (epsilon), ranging from approximately 9 to 14 in all subplots. * **Y-axes:** * Top row (a.1, b.1, c.1): Yg^(2)(ε), ranging from approximately -0.001 to 0.002. * Middle row (a.2, b.2, c.2): ∂εS(ε), ranging from approximately 0.050 to 0.060. * Bottom row (a.3, b.3, c.3): ∂εS²(ε), ranging from approximately -0.0022 to -0.0012. * **Vertical dashed red lines:** Indicate critical values εc, approximately equal to 11.2 for N=64², and 11.1 for N=128² and N=256². * **Subplot labels:** a.1, a.2, a.3, b.1, b.2, b.3, c.1, c.2, c.3, located at the right of each subplot. * **N values:** N = 64², N = 128², N = 256², located at the top of each column. ### Detailed Analysis **Column 1: N = 64²** * **a.1: Yg^(2)(ε):** The MICRO AVERAGE data (black crosses) starts at approximately 0.0013 at ε=9, increases to a peak of approximately 0.0018 at ε=11.2, and then decreases to approximately 0.0003 at ε=14. * **a.2: ∂εS(ε):** The RICCATI (black filled circles) and PEARSON (red open circles) data are almost perfectly overlapping. The data decreases linearly from approximately 0.060 at ε=9 to approximately 0.050 at ε=14. * **a.3: ∂εS²(ε):** The RICCATI (black filled circles) and PEARSON (red open circles) data are almost perfectly overlapping. The data forms a peak centered around ε=11.2, rising from approximately -0.0022 at ε=9 to approximately -0.0014 at ε=11.2, and then decreasing back to approximately -0.0022 at ε=14. **Column 2: N = 128²** * **b.1: Yg^(2)(ε):** The MICRO AVERAGE data (black crosses) starts at approximately 0.0013 at ε=9, increases to a peak of approximately 0.0019 at ε=11.1, and then decreases to approximately 0.0002 at ε=14. * **b.2: ∂εS(ε):** The RICCATI (black filled circles) and PEARSON (red open circles) data are almost perfectly overlapping. The data decreases linearly from approximately 0.060 at ε=9 to approximately 0.050 at ε=14. * **b.3: ∂εS²(ε):** The RICCATI (black filled circles) and PEARSON (red open circles) data are almost perfectly overlapping. The data forms a peak centered around ε=11.1, rising from approximately -0.0022 at ε=9 to approximately -0.0013 at ε=11.1, and then decreasing back to approximately -0.0022 at ε=14. **Column 3: N = 256²** * **c.1: Yg^(2)(ε):** The MICRO AVERAGE data (black crosses) starts at approximately 0.0014 at ε=9, increases to a peak of approximately 0.0020 at ε=11.1, and then decreases to approximately -0.0004 at ε=14. * **c.2: ∂εS(ε):** The RICCATI (black filled circles) and PEARSON (red open circles) data are almost perfectly overlapping. The data decreases linearly from approximately 0.060 at ε=9 to approximately 0.050 at ε=14. * **c.3: ∂εS²(ε):** The RICCATI (black filled circles) and PEARSON (red open circles) data are almost perfectly overlapping. The data forms a peak centered around ε=11.1, rising from approximately -0.0022 at ε=9 to approximately -0.0012 at ε=11.1, and then decreasing back to approximately -0.0022 at ε=14. ### Key Observations * The RICCATI and PEARSON methods produce nearly identical results for ∂εS(ε) and ∂εS²(ε) across all values of N. * The peak in Yg^(2)(ε) and ∂εS²(ε) becomes more pronounced and shifts slightly to the left as N increases. * The linear decrease in ∂εS(ε) is consistent across all values of N. * The vertical dashed lines indicate a critical value εc, which appears to be related to the peak in Yg^(2)(ε) and ∂εS²(ε). ### Interpretation The plots compare the behavior of different functions related to a 2D φ⁴ model using different computational methods (RICCATI, PEARSON, and MICRO AVERAGE). The parameter N seems to represent a system size or resolution, and epsilon (ε) is a variable that influences the behavior of the functions. The near-identical results from the RICCATI and PEARSON methods suggest that they are converging to the same solution for the given functions. The MICRO AVERAGE method provides a different perspective, particularly in the Yg^(2)(ε) plots, where it shows a distinct peak. The shift and sharpening of the peaks in Yg^(2)(ε) and ∂εS²(ε) as N increases suggest that the system's behavior becomes more defined and potentially approaches a critical point or phase transition. The vertical lines at εc indicate a critical value of epsilon, which may be related to this transition. The linear decrease in ∂εS(ε) suggests a consistent relationship between epsilon and this function, regardless of the value of N or the computational method used. </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 ## Multiple Line Charts: 1D XY mean field ### Overview The image contains nine line charts arranged in a 3x3 grid. Each column represents a different value of 'N' (8100, 14400, and 40000), and each row represents a different function of epsilon (€): Y_g^(2)(€), ∂_€S(€), and ∂_€S^(2)(€). The x-axis for all charts is epsilon (€), ranging from approximately 0.6 to 1.0. The charts display the relationship between epsilon and the respective functions for different values of N. The legend at the top indicates three data series: RICCATI (black circles), PEARSON (red circles), and MICRO AVERAGE (black crosses). Vertical dashed lines indicate critical values of epsilon. ### Components/Axes * **Title:** 1D XY mean field * **Legend:** Located at the top of the image. * RICCATI (black circles) * PEARSON (red circles) * MICRO AVERAGE (black crosses) * **X-axis (all charts):** € (epsilon), ranging from 0.6 to 1.0. * **Y-axis (row 1):** Y_g^(2)(€), ranging from -3 to 3. * **Y-axis (row 2):** ∂_€S(€), ranging from 1.0 to 2.5. * **Y-axis (row 3):** ∂_€S^(2)(€), ranging from -7 to -2. * **Vertical Lines:** Each chart has a vertical dashed black line and a vertical dashed green line. The black line represents the critical value of epsilon (epsilon_c^MIPA), and the green line represents epsilon_c. * **Column Titles:** * Column 1: N = 8100 * Column 2: N = 14400 * Column 3: N = 40000 * **Row Titles:** * Row 1: Y_g^(2)(€) * Row 2: ∂_€S(€) * Row 3: ∂_€S^(2)(€) ### Detailed Analysis **Column 1: N = 8100** * **Chart a.1: Y_g^(2)(€)** * MICRO AVERAGE (black crosses): The line starts at approximately 2.5 for € = 0.6, remains relatively constant until approximately € = 0.73, then sharply decreases to approximately -3 at € = 0.78, and then gradually increases to approximately -2 at € = 1.0. * **Chart a.2: ∂_€S(€)** * RICCATI (black circles): The line starts at approximately 2.4 at € = 0.6 and decreases linearly to approximately 1.9 at € = 0.73. After € = 0.73, the line curves downwards, reaching approximately 1.0 at € = 1.0. * epsilon_c^MIPA ≈ 0.73 * **Chart a.3: ∂_€S^(2)(€)** * RICCATI (black circles): The line starts at approximately -2.5 at € = 0.6 and remains relatively constant until approximately € = 0.73. After € = 0.73, the line curves upwards, reaching approximately -3 at € = 1.0. **Column 2: N = 14400** * **Chart b.1: Y_g^(2)(€)** * MICRO AVERAGE (black crosses): The line starts at approximately 2.5 for € = 0.6, remains relatively constant until approximately € = 0.75, then sharply decreases to approximately -3 at € = 0.8, and then gradually increases to approximately -2 at € = 1.0. * **Chart b.2: ∂_€S(€)** * RICCATI (black circles): The line starts at approximately 2.4 at € = 0.6 and decreases linearly to approximately 1.9 at € = 0.73. After € = 0.73, the line curves downwards, reaching approximately 1.0 at € = 1.0. * epsilon_c^MIPA ≈ 0.73 * epsilon_c ≈ 0.75 * **Chart b.3: ∂_€S^(2)(€)** * RICCATI (black circles): The line starts at approximately -2.5 at € = 0.6 and remains relatively constant until approximately € = 0.73. After € = 0.73, the line curves upwards, reaching approximately -3 at € = 1.0. **Column 3: N = 40000** * **Chart c.1: Y_g^(2)(€)** * MICRO AVERAGE (black crosses): The line starts at approximately 2.5 for € = 0.6, remains relatively constant until approximately € = 0.77, then sharply decreases to approximately -3 at € = 0.82, and then gradually increases to approximately -2 at € = 1.0. * **Chart c.2: ∂_€S(€)** * RICCATI (black circles): The line starts at approximately 2.4 at € = 0.6 and decreases linearly to approximately 1.9 at € = 0.74. After € = 0.74, the line curves downwards, reaching approximately 1.0 at € = 1.0. * epsilon_c^MIPA ≈ 0.74 * **Chart c.3: ∂_€S^(2)(€)** * RICCATI (black circles): The line starts at approximately -2.5 at € = 0.6 and remains relatively constant until approximately € = 0.74. After € = 0.74, the line curves upwards, reaching approximately -3 at € = 1.0. ### Key Observations * The critical value of epsilon (epsilon_c^MIPA) increases slightly with increasing N (0.73, 0.73, 0.74). * The function Y_g^(2)(€) exhibits a sharp transition around the critical value of epsilon. * The functions ∂_€S(€) and ∂_€S^(2)(€) show a change in slope around the critical value of epsilon. * The PEARSON data series is not visible in any of the charts. ### Interpretation The charts illustrate the behavior of different functions related to a 1D XY mean field model as a function of epsilon (€) for different system sizes (N). The sharp transitions observed in Y_g^(2)(€) and the changes in slope in ∂_€S(€) and ∂_€S^(2)(€) around the critical value of epsilon (epsilon_c^MIPA) suggest a phase transition or a significant change in the system's behavior at this point. The slight increase in epsilon_c^MIPA with increasing N may indicate a size-dependent effect on the critical behavior. The absence of the PEARSON data series suggests that it may not be relevant or available for this particular analysis. </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.