# 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
## [Chart/Diagram Type]: Comparative Analysis of $V_{\phi^4}$ and $V_{XY}$ Models
### Overview
The image presents a comparative analysis of two physical models, $V_{\phi^4}$ (left column) and $V_{XY}$ (right column), across three different metrics (rows). Each column displays three vertically stacked plots showing the behavior of specific functions ($Y_g^{(2)}$, $\partial_\epsilon S$, and $\partial_\epsilon^2 S$) relative to a parameter $\epsilon$. The plots illustrate critical phenomena, likely phase transitions, with specific critical values ($\epsilon_c$) highlighted by vertical lines.
### Components/Axes
**Top Formulas:**
* **Left:** $V_{\phi^4} := \sum_i [\frac{\mu^2}{2}\phi_i^2 + \frac{\lambda}{4!}\phi_i^4] + \frac{J}{2}\sum_{\mu=1}^n (\phi_i - \phi_{i+e_\mu})^2$
* **Right:** $V_{XY} = \frac{J}{2N} \sum_{i,j} [1 - \cos(\theta_i - \theta_j)]$
**Plot Structure:**
* **Left Column (a.1, a.2, a.3):** X-axis represents $\epsilon$ ranging from 9 to 14.
* **Right Column (b.1, b.2, b.3):** X-axis represents $\epsilon$ ranging from 0.55 to 1.0.
* **Legends:**
* **a.1/b.1:** "x micro average" (represented by black 'x' markers).
* **a.2/a.3/b.2/b.3:** "PHT" (black filled circle) and "EFE" (red open circle).
* **Annotations:**
* **a.2:** $\epsilon_c^{MIPA} \approx 11.1$ (indicated by a red dashed vertical line).
* **b.2:** $\epsilon_c^{MIPA} \approx 0.74$ (indicated by a black dashed vertical line) and $\epsilon_c^\infty = \frac{3J}{4}$ (indicated by a green dashed vertical line).
### Detailed Analysis
#### Left Column: $V_{\phi^4}$
* **a.1 ($Y_g^{(2)}(\epsilon)$):**
* **Trend:** The data points (black 'x') form a sharp peak. The values increase from $\approx 1.4 \times 10^{-3}$ at $\epsilon=9$ to a maximum of $\approx 2.0 \times 10^{-3}$ at $\epsilon \approx 11.1$, then decrease sharply to $\approx 0.2 \times 10^{-3}$ at $\epsilon=14.5$.
* **a.2 ($\partial_\epsilon S(\epsilon)$):**
* **Trend:** The data (PHT/EFE) shows a smooth, monotonically decreasing curve.
* **Values:** Starts at $\approx 6.05 \times 10^{-2}$ at $\epsilon=9$ and decreases to $\approx 5.0 \times 10^{-2}$ at $\epsilon=14.5$.
* **Note:** The PHT and EFE data points are perfectly superimposed.
* **a.3 ($\partial_\epsilon^2 S(\epsilon)$):**
* **Trend:** The data (PHT/EFE) shows a sharp cusp/peak at $\epsilon \approx 11.1$.
* **Values:** Increases from $\approx -2.2 \times 10^{-3}$ at $\epsilon=9$ to a peak of $\approx -1.2 \times 10^{-3}$ at $\epsilon \approx 11.1$, then decreases back to $\approx -2.2 \times 10^{-3}$ at $\epsilon=14.5$.
#### Right Column: $V_{XY}$
* **b.1 ($Y_g^{(2)}(\epsilon)$):**
* **Trend:** The data (black 'x') shows a discontinuity at $\epsilon \approx 0.74$.
* **Values:** Decreases from $\approx 3.0$ at $\epsilon=0.55$ to $\approx 2.5$ at $\epsilon \approx 0.74$. Immediately after the discontinuity, the value drops to $\approx -3.5$ and then increases towards $\approx -1.0$ at $\epsilon=1.0$.
* **b.2 ($\partial_\epsilon S(\epsilon)$):**
* **Trend:** The data (PHT/EFE) shows a discontinuity at $\epsilon \approx 0.74$.
* **Values:** Decreases from $\approx 2.5$ at $\epsilon=0.55$ to $\approx 2.0$ at $\epsilon \approx 0.74$. After the discontinuity, it drops to $\approx 1.5$ and continues to decrease to $\approx 1.0$ at $\epsilon=1.0$.
* **b.3 ($\partial_\epsilon^2 S(\epsilon)$):**
* **Trend:** The data (PHT/EFE) shows a sharp discontinuity at $\epsilon \approx 0.74$.
* **Values:** Increases from $\approx -3.2$ at $\epsilon=0.55$ to $\approx -1.5$ at $\epsilon \approx 0.74$. After the discontinuity, it drops sharply to $\approx -6.5$ and then increases towards $\approx -2.0$ at $\epsilon=1.0$.
### Key Observations
* **Data Superposition:** In all plots containing both PHT and EFE markers (a.2, a.3, b.2, b.3), the red open circles (EFE) and black filled circles (PHT) are perfectly aligned, indicating identical results from both methods.
* **Critical Behavior:** The vertical lines denote critical values ($\epsilon_c$).
* The left column ($V_{\phi^4}$) exhibits continuous behavior with a peak/cusp at $\epsilon_c$.
* The right column ($V_{XY}$) exhibits a clear discontinuity (jump) at $\epsilon_c$, characteristic of a first-order phase transition.
* **Consistency:** The "micro average" data points (black 'x') are consistent with the trends observed in the PHT/EFE data series.
### Interpretation
The data demonstrates the behavior of two distinct physical systems near critical points.
* **$V_{\phi^4}$ (Left):** The continuous nature of the derivatives suggests a second-order (continuous) phase transition. The peak in $Y_g^{(2)}$ and the cusp in $\partial_\epsilon^2 S$ at $\epsilon \approx 11.1$ are hallmarks of critical fluctuations.
* **$V_{XY}$ (Right):** The sharp discontinuities in all three metrics at $\epsilon \approx 0.74$ strongly suggest a first-order (discontinuous) phase transition. The system undergoes a sudden change in state at this critical value.
* **Methodology:** The perfect overlap between PHT (likely Perturbative Hartree) and EFE (likely Effective Field Theory/Equation) suggests that both theoretical frameworks are in excellent agreement for these specific models. The inclusion of $\epsilon_c^\infty = \frac{3J}{4}$ in plot b.2 provides an analytical reference point for the observed numerical discontinuity.
</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/Diagram: 2D $\phi^4$ nearest neighbours
### Overview
This image presents a 3x3 grid of plots analyzing the behavior of a 2D $\phi^4$ model across three different system sizes ($N = 64^2$, $N = 128^2$, and $N = 256^2$). The columns represent the system size, while the rows represent three different physical quantities: $\Upsilon_g^{(2)}(\epsilon)$, $\partial_\epsilon S(\epsilon)$, and $\partial_\epsilon S^2(\epsilon)$. The plots investigate the behavior of these quantities as a function of $\epsilon$ (ranging from 9 to 15), highlighting a critical value $\epsilon_c$.
### Components/Axes
* **Legend (Top Center):**
* **RICCATI:** Represented by black filled circles ($\bullet$).
* **PEARSON:** Represented by red open circles ($\circ$).
* **MICRO AVERAGE:** Represented by black crosses ($\times$).
* **Columns (System Size $N$):**
* **Left Column:** $N = 64^2$ (labeled a.1, a.2, a.3).
* **Middle Column:** $N = 128^2$ (labeled b.1, b.2, b.3).
* **Right Column:** $N = 256^2$ (labeled c.1, c.2, c.3).
* **Rows (Quantities):**
* **Top Row:** Y-axis labeled $\Upsilon_g^{(2)}(\epsilon)$. Scale: 0.0005 to 0.0020.
* **Middle Row:** Y-axis labeled $\partial_\epsilon S(\epsilon)$. Scale: 0.050 to 0.060.
* **Bottom Row:** Y-axis labeled $\partial_\epsilon S^2(\epsilon)$. Scale: -0.0022 to -0.0012.
* **X-Axis:** Labeled $\epsilon$ at the bottom of the grid. Range: 9 to 15.
* **Vertical Markers:** A red dashed vertical line is present in all plots, indicating the critical value $\epsilon_c$.
* Left column: $\epsilon_c \simeq 11.2$.
* Middle and Right columns: $\epsilon_c \simeq 11.1$.
### Detailed Analysis
#### Row 1: $\Upsilon_g^{(2)}(\epsilon)$ (Plots a.1, b.1, c.1)
* **Trend:** All plots show a sharp peak at $\epsilon_c$. The curve rises from $\epsilon=9$, peaks at $\epsilon_c$, and decays as $\epsilon$ increases toward 15.
* **Data Series:** Only "MICRO AVERAGE" (black crosses) is plotted.
* **Observation:** As the system size $N$ increases from $64^2$ to $256^2$, the peak becomes progressively sharper and higher, characteristic of finite-size scaling near a phase transition.
#### Row 2: $\partial_\epsilon S(\epsilon)$ (Plots a.2, b.2, c.2)
* **Trend:** All plots show a monotonic downward slope as $\epsilon$ increases.
* **Data Series:** Both "RICCATI" (black filled circles) and "PEARSON" (red open circles) are plotted.
* **Observation:** The two data series overlap almost perfectly across all three system sizes, indicating high consistency between the Riccati and Pearson methods. The red dashed line at $\epsilon_c$ intersects the curve at a point that appears to be an inflection point or a specific reference value.
#### Row 3: $\partial_\epsilon S^2(\epsilon)$ (Plots a.3, b.3, c.3)
* **Trend:** All plots show a peak (maximum value, closest to 0) at $\epsilon_c$. The curve rises from $\epsilon=9$ to the peak and then decreases (becomes more negative) as $\epsilon$ increases.
* **Data Series:** Both "RICCATI" (black filled circles) and "PEARSON" (red open circles) are plotted.
* **Observation:** Similar to Row 2, the two methods show excellent agreement. As $N$ increases, the peak at $\epsilon_c$ becomes sharper, consistent with the behavior observed in Row 1.
### Key Observations
* **Finite Size Scaling:** The sharpening of peaks in rows 1 and 3 as $N$ increases from $64^2$ to $256^2$ is a classic signature of a phase transition in a lattice model.
* **Methodological Agreement:** The "RICCATI" and "PEARSON" methods produce nearly identical results, suggesting that both are robust estimators for the quantities being measured.
* **Critical Point Shift:** The critical value $\epsilon_c$ is estimated at 11.2 for the smallest system size ($64^2$) and shifts to 11.1 for the larger system sizes ($128^2$ and $256^2$), suggesting the critical point is being refined as the system approaches the thermodynamic limit.
### Interpretation
This figure demonstrates the numerical analysis of a phase transition in a 2D $\phi^4$ model. The data suggests that the system undergoes a critical change at $\epsilon_c \approx 11.1$. The quantities $\Upsilon_g^{(2)}(\epsilon)$ and $\partial_\epsilon S^2(\epsilon)$ act as indicators of this transition, exhibiting singular behavior (peaks) at the critical point. The fact that these peaks sharpen with increasing system size ($N$) is strong evidence of a second-order phase transition. The high degree of overlap between the Riccati and Pearson data points validates the computational approach used to derive these thermodynamic quantities.
</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
## Multi-Panel Scatter Plot: 1D XY mean field
### Overview
The image presents a 3x3 grid of scatter plots analyzing the behavior of a "1D XY mean field" system. The grid is organized by system size ($N$) in columns and by different mathematical functions in rows. The plots illustrate a phase transition occurring at a critical value of the parameter $\epsilon$ (epsilon), consistently located between 0.7 and 0.8.
### Components/Axes
**Legend (Top Center):**
* **RICCATI**: Represented by black filled circles.
* **PEARSON**: Represented by red open circles.
* **MICRO AVERAGE**: Represented by black 'x' markers.
**Grid Structure:**
* **Columns (System Size $N$):**
* Left Column: $N = 8100$
* Center Column: $N = 14400$
* Right Column: $N = 40000$
* **Rows (Functions):**
* Top Row: $\Upsilon_g^{(2)}(\epsilon)$
* Middle Row: $\partial_\epsilon S(\epsilon)$
* Bottom Row: $\partial_\epsilon S^2(\epsilon)$
**Axes:**
* **X-axis:** $\epsilon$ (epsilon), ranging from approximately 0.6 to 1.0. This axis is shared across all nine plots.
* **Y-axes:**
* Top Row: Range approx. -3 to +3.
* Middle Row: Range approx. 1.0 to 2.5.
* Bottom Row: Range approx. -7 to -2.
**Annotations:**
* Vertical dashed lines (Red and Green) appear in every plot at the critical point $\epsilon_c$.
* Middle Row, Left/Center: $\epsilon_c^{\text{MIPA}} \simeq 0.73$
* Middle Row, Right: $\epsilon_c^{\text{MIPA}} \simeq 0.74$
* Middle Row, Center: $\epsilon_c \simeq 0.75$
### Detailed Analysis
#### Top Row: $\Upsilon_g^{(2)}(\epsilon)$
* **Trend:** The data (using "MICRO AVERAGE" 'x' markers) shows a linear downward slope from $\epsilon \approx 0.6$ to $\approx 0.73$. At the critical point, there is a sharp, discontinuous drop to negative values (approx. -3). Following the discontinuity, the values increase gradually, curving upward toward 0 as $\epsilon$ approaches 1.0.
* **Consistency:** The trend is identical across all three system sizes ($N=8100, 14400, 40000$).
#### Middle Row: $\partial_\epsilon S(\epsilon)$
* **Trend:** The data (using overlapping "RICCATI" black dots and "PEARSON" red circles) shows a linear downward slope from $\epsilon \approx 0.6$ to $\approx 0.73$. At the critical point, there is a sharp discontinuity, followed by a steeper downward slope as $\epsilon$ increases toward 1.0.
* **Consistency:** The overlap between RICCATI and PEARSON data is near-perfect, indicating these two methods yield identical results for this metric.
#### Bottom Row: $\partial_\epsilon S^2(\epsilon)$
* **Trend:** The data (using overlapping "RICCATI" black dots and "PEARSON" red circles) shows a gradual upward slope from $\epsilon \approx 0.6$ to $\approx 0.73$. At the critical point, there is a sharp, discontinuous drop to a minimum value (approx. -7). Following this, the values increase rapidly.
* **Consistency:** Similar to the middle row, the RICCATI and PEARSON data points overlap perfectly.
### Key Observations
* **Method Equivalence:** In the middle and bottom rows, the RICCATI and PEARSON methods produce indistinguishable results, suggesting they are mathematically equivalent or converge to the same solution in this context.
* **Phase Transition:** All three functions exhibit a clear discontinuity at $\epsilon_c \approx 0.73 - 0.75$. This is characteristic of a first-order phase transition.
* **Finite-Size Effects:** While the qualitative behavior is consistent across all $N$, the annotated critical point $\epsilon_c^{\text{MIPA}}$ shifts slightly from 0.73 to 0.74 as the system size $N$ increases from 8100 to 40000.
* **Spatial Grounding:** The legend is positioned at the top center of the entire figure. The labels a.1, a.2, a.3 (left column), b.1, b.2, b.3 (center column), and c.1, c.2, c.3 (right column) serve as identifiers for the sub-plots.
### Interpretation
The data demonstrates a robust phase transition in a 1D XY mean field model. The functions $\Upsilon_g^{(2)}$, $\partial_\epsilon S$, and $\partial_\epsilon S^2$ appear to be thermodynamic derivatives (likely related to entropy $S$). The sharp discontinuities at $\epsilon_c$ indicate that the system undergoes a sudden change in state at this critical parameter value.
The "MIPA" label likely refers to a "Mean Field" approximation method. The fact that the critical point shifts slightly with increasing $N$ suggests that while the transition is present in the mean field limit, finite-size scaling effects are still observable. The perfect overlap between RICCATI and PEARSON data points suggests that for this specific physical model, both analytical approaches are equally valid and provide the same numerical output.
</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.