\section{Computation of the geometric observables}

We consider a generic autonomous Hamiltonian system, described by the Hamiltonian function $H:\Lambda\subset\mathbb{R}^{N}\times\mathbb{R}^{N}\rightarrow\mathbb{R}$ where $\Lambda$ is the phase space. It is well-known that, picking a fixed energy value $H(\bm{x})=E$, the dynamics of the representative point of the system, $\bm{x}=\{\bm{p},\bm{q}\}\in\Lambda$, lies on the energy level set \cite{pettini2007geometry} ($n=2N$):
\begin{equation}\label{def:energy_level_sets}
    \Sigma^{H}_{E}:=\{\bm{x}\in \Lambda\,|\, H(\bm{x})=E\}\subset\mathbb{R}^n.
\end{equation} 
By allowing the energy to vary over a set of accessible energies denoted by $\mathcal{E}:=[E_0,E_1]$, we obtain a collection $\{\Sigma_{E}^{H}\}_{E\in\mathcal{E}}$ of energy level sets which defines a foliation of the phase space, i.e., $\Lambda=\bigcup_{E\in\mathcal{E}}\;\Sigma_{E}^{H}$.
We thus need to introduce an appropriate metric tensor on the phase space which, in turn, induces a metric tensor, $h_{E}$, on each hypersurface. Now, from a purely mathematical viewpoint, we have no restriction about the specific structure of the metric tensor and this led to introduce in Ref. \cite{di2021topology} the metric tensor
\begin{equation}\label{def:metric_tensor_gromov}
    g = dE\otimes dE+h_E,
\end{equation}
where $E$ is interpreted as an energy parameter and $h_E$ is the Euclidean metric tensor induced on the hypersurface $\Sigma_{E}^{H}$.\\
By adopting the metric tensor \eqref{def:metric_tensor_gromov}, we can study the extrinsic geometry of each energy level set through the introduction of the Weingarten operator
\begin{equation}\label{def:weingarten_g}
    \mathcal{W}^{g}_{\bm{\nu}}(\bm{X}):=\nabla_{\bm{X}}\bm{\nu}\equiv\nabla\left(\frac{\nabla H}{\|\nabla H\|}\right)\cdot\bm{X},
\end{equation}
where $\bm{X}$ is any vector field lying on the tangent space to the energy level set $\Sigma_{E}^{H}$, $\nabla$ and $\|\cdot\|$ are, respectively, the gradient operator and the norm defined on the phase space.
The phase space can be equipped with a more appropriate metric tensor \cite{gori_thesis,gori2018topological} defined by
\begin{equation}\label{def:conformal_metric}
    \widetilde{g}=dE\otimes dE+\chi^{\frac{2}{n-1}}(h_E)_{ij}~dx^{i}\otimes dx^{j},
\end{equation}
where $\chi:=1/\|\nabla H\|$. The collection of energy level sets can be generated evolving the energy level set $\Sigma_{E_0}^{H}$ through the vector field $\partial_{E}$ such that $dE(\partial_{E})=1$ \cite{hirsch2012differential}. In particular, it is easy to show that \cite{di2021topology}:
\begin{equation}\label{def:zeta_vector}
    \bm{\xi}:=\partial_{E}=\chi\bm{\nu}.
\end{equation}
We note that the vector field $\bm{\xi}$ has the following properties: \textit{(i)} it generates the flow so that $dE(\partial_{E})=1$, \textit{(ii)} it is normalized to the unity in the metric tensor \eqref{def:conformal_metric}. This allows us to introduce a new Weingarten operator defined by
\begin{equation}\label{def:weingarten_g_tilde}
    \mathcal{W}^{\widetilde{g}}_{\bm{\xi}}(\bm{X}):=\nabla^{\widetilde{g}}_{\bm{X}}\bm{\xi}\cdot \bm{X}\equiv\nabla\left(\frac{\nabla H}{\|\nabla H\|^{2}}\right)\cdot \bm{X}
\end{equation}
where $\nabla^{\widetilde{g}}$ is defined with respect to the metric tensor $\widetilde{g}$ whereas $\nabla$ is the gradient operator on the phase space as defined before.\\
The Weingarten operator \eqref{def:weingarten_g_tilde} is related to the previous one (see Eq. \eqref{def:weingarten_g}) through the relation
\begin{equation}\label{def:transf_wing}
    \mathcal{W}^{g}_{\bm{\nu}}\mapsto   \mathcal{W}^{\widetilde{g}}_{\bm{\xi}}=\chi~\mathcal{W}^{g}_{\bm{\nu}}+\chi^{-1}\partial_{E}\chi\frac{1\!\!\!1_{\Sigma^{H}_{E}}}{n-1}, 
\end{equation}
where 
\begin{equation}
    \partial_{E}\chi=\langle \bm{\xi},\nabla\chi\rangle_{\widetilde{g}}=\chi\langle \bm{\nu},\nabla\chi\rangle_{\widetilde{g}}.
\end{equation}
and $1\!\!\!1_{\Sigma_{E}^{H}}$ is the identity operator defined on the tangent space to $\Sigma_{E}^{H}$.


\subsection{Expression of the geometric observables}





We now derive the explicit expression for the geometric curvature functions (GCFs) adopting the metric tensor $\widetilde{g}$ together with the expression \eqref{def:transf_wing} for the Weingarten.\\
We recall below the expression for the first and second order GCFs as reported in the main text:  
\begin{equation}
\begin{split}\label{def:Omega_1_2}
    \Omega_{\widetilde{g}}^{(1)}(E)&:=\int_{\Sigma_{E}^{H}}Tr[\mathcal{W}^{\widetilde{g}}_{\bm{\xi}}]~d\rho\\
    \Omega^{(2)}_{\widetilde{g}}(E)&:=\frac{1}{2}\int_{\Sigma_{E}^{H}}\bigg\{Tr[\mathcal{W}^{\widetilde{g}}_{\bm{\xi}}]^{2}-Tr[(\mathcal{W}^{\widetilde{g}}_{\bm{\xi}})^{2}]
    +R^{\widetilde{g}}(\Sigma_{E}^{H})-R^{\widetilde{g}}(\Lambda)\bigg\}\;d\rho,
\end{split}
\end{equation}
Thus, the Weingarten operator is
\begin{equation}
    \mathcal{W}^{\widetilde{g}}_{\bm{\xi}}=\frac{Hess\, H}{\|\nabla H\|^{2}}-2\frac{\nabla H\otimes Hess\, H\cdot\nabla H}{\|\nabla H\|^{4}},
\end{equation}
so that the trace reads:
\begin{equation}
    Tr[\mathcal{W}^{\widetilde{g}}_{\bm{\xi}}]=\frac{\Delta H}{\|\nabla H\|^{2}}-2\frac{\langle\nabla H, Hess\, H\cdot\nabla H\rangle}{\|\nabla H\|^{4}}.
\end{equation}
For what concerns the second order GCF, this can be explicitly computed from its definition:
\begin{equation}
    \Omega_{\widetilde{g}}^{(2)}(E):=\frac{\partial^2_{E}vol^{\widetilde{g}}(E)}{vol^{\widetilde{g}}(E)}.
\end{equation}
In other words, we compute
\begin{equation}
    \partial^2_{E}vol^{\widetilde{g}}(E)=\int_{\Sigma_{E}^{H}}\partial_{E}\left(Tr[\mathcal{W}^{\widetilde{g}}_{\bm{\xi}}]d\mu^{\widetilde{g}}\right),
\end{equation}
which yields
\begin{equation}
\begin{split}\label{def:omega_function}
    \Omega^{(2)}_{\widetilde{g}}(E):=\int_{\Sigma_{E}}\bigg\{\frac{Tr[\mathcal{W}^{g}_{\bm{\nu}}]^2-Tr[(\mathcal{W}^{g}_{\bm{\nu}})^2]}{\|\nabla H\|^{2}}&-3\frac{Tr[\mathcal{W}^{g}_{\bm{\nu}}]}{\|\nabla H\|^{5}}\langle\nabla H,HessH\cdot\nabla H\rangle_{g^{\mathbb{E}}}-\frac{\nabla\nabla\nabla H(\nabla H,\nabla H,\nabla H)}{\|\nabla H\|^{6}}\\
    &+\frac{5}{\|\nabla H\|^{8}}\langle\nabla H,HessH\cdot\nabla H\rangle^{2}_{g^{\mathbb{E}}}-2\frac{\|HessH\cdot\nabla H\|^{2}}{\|\nabla H\|^{6}} \bigg\}d\rho^{\widetilde{g}},
\end{split}
\end{equation}
where $\mathcal{W}_{\bm{\nu}}^{g}$ is the Weingarten operator defined in Eq. \eqref{def:weingarten_g} whose traces are
\begin{equation}
\begin{split}
    Tr[\mathcal{W}^{g}_{\bm{\nu}}]&=\frac{\Delta H}{\|\nabla H\|}-\frac{\langle\nabla H, Hess\, H\cdot\nabla H\rangle}{\|\nabla H\|^{3}},\\
    Tr[(\mathcal{W}^{g}_{\bm{\nu}})^2]&=\frac{Tr[Hess\,H^2]}{\|\nabla H\|^2}+\frac{\langle \nabla H,Hess\,H\cdot\nabla H\rangle^2}{\|\nabla H\|^6}-2\frac{\|Hess\,H\cdot\nabla H\|^2}{\|\nabla H\|^4}.
\end{split}
\end{equation}
Note that the trace of the square of the Weingarten operator $\mathcal{W}^{\widetilde{g}}_{\bm{\xi}}$ is
\begin{equation}
    Tr[(\mathcal{W}^{\widetilde{g}}_{\bm{\xi}})^2]=\frac{Tr[Hess\;H^2]}{\|\nabla H\|^4}
    + 4 \frac{\langle \nabla H,Hess\;H\nabla H\rangle^2}{\|\nabla H\|^8} - 4\frac{\|Hess\;H\nabla H\|^2}{\|\nabla H\|^6}.
\end{equation}






\section{Numerical details for the $\phi^{4}$-model}






\subsection{Integration of the Hamilton's equation and computation of the thermodynamic observables}
\label{app:thermodynamics}

We performed molecular dynamics simulations in the microcanonical ensemble for the 2-dimensional $\phi^{4}$ Hamiltonian system with $N=30^{2}=900$ particles. The Hamiltonian is
\begin{equation}
\begin{split}\label{def:hamiltonian_phi4}
    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\!\!N(\bm{i})}(\phi_{\bm{i}}-\phi_{\bm{i}})^{2}\Bigg]
\end{split}
\end{equation}
with $\lambda=3/5$, $\mu=2$ and $J=1$, then, $\bm{i}:=(i_1,i_2)$ is the two-dimensional index for labeling the sites, i.e., $1\leq i_1\leq n=30$ and $1\leq i_2\leq n=30$; finally $N\!\!N(\bm{i})$ is the set of the nearest neighbour lattice sites associated to $\bm{i}$.\\
The numerical integration of the Hamilton's equations of motion derived from the Hamiltonian function \eqref{def:hamiltonian_phi4} has been performed, with periodic boundary conditions, choosing random initial conditions and using a bilateral symplectic integration scheme \cite{casetti1995efficient,gori2018topological} with a time step such that the energy is conserved within the relative precision of $\Delta E/E\approx 10^{-6}$. The $2D$-$\phi^{4}$ model undergoes a phase transition at the critical energy density value $\epsilon_{t}=E_{t}/N\approx 11.1$ as already observed in Refs. \cite{bel2020geometrical,gori2018topological}. Note that by working with the potential function, instead of the Hamiltonian one, the critical potential density value is $v_{t}\approx 2.2$ as found in Refs. \cite{mehta2012energy,kastner2011phase,caiani1998hamiltonian}. Finally, since the study of the thermodynamic observables' behaviors as function of the number of the degrees of freedom, $N$, has been already done in Ref. \cite{bel2020geometrical}, we restricted our analysis to a system with $900$ particles. We evaluated the geometric observables in Eq. \eqref{def:Omega_1_2} exploiting the ergodic hypothesis, that is, we converted the microcanonical averages into time averages \cite{bel2021geometrical,gori2018topological}. More precisely, given a generic phase-space-valued function, $f:\Lambda\rightarrow\mathbb{R}$, we have
\begin{equation}
\begin{split}\label{def:time_average}
    \langle f(\Sigma_{E}^{H})\rangle_{M}:=\int_{\Sigma_{E}^{H}}f\;d\rho\equiv \lim_{T\to\infty}\frac{1}{T}\int_{0}^{T} f(\bm{X}(\tau))\;d\tau=:\overline{f(\Sigma_{E}^{H})},
\end{split}
\end{equation}
where we recall that $\bm{x}=\{\bm{q}^{1},\ldots,\bm{q}^{N},\bm{p}_{1},\ldots,\bm{p}_{N}\}$ whereas $\bm{X}(\tau)$ is the phase space-trajectory---solution of the Hamilton's equations of motion---which has been computed numerically.\\
By means of Eq. \eqref{def:time_average}, we evaluated $\Omega_{\widetilde{g}}^{(1)}$ and $\Omega_{\widetilde{g}}^{(2)}$ along the Hamiltonian dynamics for each energy density value within the set $[0.5,40]$. The results are plotted in the main text. We observe a least sensitive inflection point in $\Omega_{\widetilde{g}}^{(1)}$ at the transition point $\epsilon_t=11.1$ and a negative-valued peak in $\partial_{E}^2 S_{\widetilde{g}}=\Omega_{\widetilde{g}}^{(2)}-(\Omega_{\widetilde{g}}^{(1)})^2$ still in $\epsilon_t=11.1$. Note that the energy-behaviors of both $\partial_{E} S_{\widetilde{g}}=\Omega_{\widetilde{g}}^{(1)}$ and $\partial_{E}^2 S_{\widetilde{g}}$ are qualitatively the same as those of $\beta(E)$ and $\gamma(E)$ in Fig. 4 (main text) predicted by Bachmann. Finally, we note that these results are in agreement with those obtained in Ref. \cite{bel2020geometrical}.

In order to get information about the phase transition and the critical energy point, we numerically computed time averages of the relevant thermodynamics observables such as the caloric curve, magnetization and specific heat. The magnetization is defined by \cite{pettini2007geometry}
\begin{equation}\label{def:magnetization}
    \mathcal{M}:=\sum_{\bm{i}}|\phi_{\bm{i}}|,
\end{equation}
whose average is computed in agreement with Eq. \eqref{def:time_average}. The caloric curve is give, instead, by ($k_B=1$) \cite{pettini2007geometry}
\begin{equation}
    T(E):=\frac{2}{N}\langle K\rangle
\end{equation}
where $\langle K\rangle$ is the time average of the kinetic energy defined through the momenta degrees of freedom $\pi_{i}$.\\
Finally, the specific heat can be computed through averaged of the kinetic energy as follows \cite{pettini2007geometry}
\begin{equation}\label{def:specific_heat}
    C_{v}(E):=\left(N-(N-2)\langle K\rangle\langle K^{-1}\rangle\right)^{-1},
\end{equation}
The numerical results are reported below. The magnetization as defined in Eq. \eqref{def:specific_heat} has been plotted in Fig. \ref{SFig:magnetization}, the caloric curve in Fig. \eqref{SFig:caloric_curve} and, finally, the specific heat defined by \eqref{def:specific_heat} is plotted in Fig. \ref{SFig:specific_heat}.
We remark that the results obtained in this work are in agreement with those published in Ref. \cite{bel2021geometrical}.



\begin{figure}[!htb]
    \begin{minipage}[t]{.3\textwidth}
        \centering
        \includegraphics[width=\textwidth]{plots/magnetization_30.pdf}
       \caption{Magnetization.}
        \label{SFig:magnetization}
    \end{minipage}
        \begin{minipage}[t]{.3\textwidth}
        \includegraphics[width=\textwidth]{plots/caloric_curve_30.pdf}
        \caption{Caloric curve. }
        \label{SFig:caloric_curve}
    \end{minipage}
    \begin{minipage}[t]{.3\textwidth}
        \includegraphics[width=\textwidth]{plots/specific_heat_30.pdf}
        \caption{Specific heat.}
        \label{SFig:specific_heat}
    \end{minipage}
    \caption*{Plots of the standard thermodynamics observables for the $\phi^{4}$-model. The vertical dashed lines mark the transition energy point, $\epsilon_t=E_t/N\approx 11.1$.}
    \label{}
\end{figure}




\subsection{Geometric observables in $\phi^4$-model}


In order to numerically evaluate $\Omega^{(1)}_{\widetilde{g}}$ and $\Omega^{(2)}_{\widetilde{g}}$, we have to provide an explicit mathematical expression for the Hessian and the gradient of the Hamiltonian function \eqref{def:hamiltonian_phi4}. Let us define $\bm{x}_{\bm{k}}:=(\pi_{\bm{k}},\phi_{\bm{k}})$, then, we have:
\begin{equation}
    \nabla_{\pi_{\bm{k}}}H=\pi_{\bm{k}},\qquad \nabla_{\phi_{\bm{k}}}H=\frac{\lambda}{3!}\phi^{3}_{\bm{k}}+(4J-\mu^2)\phi_{\bm{k}}-J\sum_{\bm{i}\in N\!\!N(\bm{k})}\phi_{\bm{i}},
\end{equation}
then, the components of the gradient of $H$ are given by
\begin{equation}
    \nabla_{\bm{x}_{\bm{k}}} H=(\nabla_{\pi_{\bm{k}}}H,\nabla_{\phi_{\bm{k}}}H).
\end{equation}
The Hessian of $H$ can be obtained by differentiating $\nabla_{\bm{k}}H$, we have
\begin{equation}
    Hess\,H=\left(\begin{matrix}
    1\!\!\!1_{\mathbb{R}^{N}} & \bm{0}_{\mathbb{R}^N}\\
    \bm{0}_{\mathbb{R}^N} & Hess\,V
    \end{matrix}\right),
\end{equation}
where
\begin{equation}
    Hess\,V_{\bm{kj}}=\nabla_{\phi_{\bm{k}}}\nabla_{\phi_{\bm{j}}}V=\left(\frac{\lambda}{2}\phi^{2}_{\bm{j}}+4J-\mu^2\right)\delta_{\bm{kj}}-J\delta_{\bm{j}, N\!\!N(\bm{k})}.
\end{equation}
Moreover, we have
\begin{equation}
\begin{split}
    \nabla\nabla\nabla\,V_{\bm{kjl}}&=\nabla_{\phi_{\bm{k}}}\nabla_{\phi_{\bm{j}}}\nabla_{\phi_{\bm{l}}}V=\lambda\phi_{\bm{j}}\delta_{\bm{kj}}\delta_{\bm{jl}},\\
    \nabla\nabla\nabla\nabla\,V_{\bm{kjli}}&=\nabla_{\phi_{\bm{k}}}\nabla_{\phi_{\bm{j}}}\nabla_{\phi_{\bm{l}}}\nabla_{\phi_{\bm{i}}}V=\lambda\delta_{\bm{ji}}\delta_{\bm{kj}}\delta_{\bm{jl}},\\
    \nabla\nabla\nabla\nabla\nabla\,V_{\bm{kjlim}}&=\nabla_{\phi_{\bm{k}}}\nabla_{\phi_{\bm{j}}}\nabla_{\phi_{\bm{l}}}\nabla_{\phi_{\bm{i}}}\nabla_{\phi_{\bm{m}}}V=0.
\end{split}
\end{equation}
The laplacian of $H$ is given by taking the trace of the Hessian:
\begin{equation}
    \Delta H=N+Tr[Hess\,V]=N+\frac{\lambda}{2}\sum_{\bm{i}}\phi^{2}_{\bm{i}}+N(4J-\mu^2).
\end{equation}



\section{Numerical details for the Ginzburg-Landau-like model}


This Hamiltonian system has been already studied in Ref.~\cite{di2021topology}, we report here just those numerical details which are the strictly necessary for the analysis of the geometric observables.\\

The Ginzburg-Landau-like model is defined by the Hamiltonian 
\begin{equation}\label{def:GL_hamiltonian}
    H_{GL}:=\sum_{i=1}^{n}\frac{p^{2}_{i}}{2}-\frac{\alpha}{2}\sum_{i=1}^{n}(q^{i})^{2}+\frac{\beta}{4}\left(\sum_{i=1}^{n}(q^{i})^{2}\right)^{2},
\end{equation}
where $\alpha,\;\beta\in\mathbb{R}^{+}$.
Depending on the value of $E$, the energy level sets show a change in their topology. In fact, by defining the \emph{order parameter}
\begin{equation}\label{def:order_parameter}
    r^{2}=\sum_{i=1}^{n}(q^{i})^{2},
\end{equation}
together with the total momentum:
\begin{equation}\label{}
    P^{2}=\sum_{i=1}^{n}p_{i}^{2},
\end{equation}
the Hamiltonian function \eqref{def:GL_hamiltonian} rewrites
\begin{equation}\label{def:hamiltonian_function_GL}
    H_{GL}(P,r)=\frac{P^{2}}{2}-\frac{\alpha}{2}r^{2}+\frac{\beta}{4}r^{4}.
\end{equation}
We note that it admits three classes of stationary points given by the condition $\nabla H_{GL}=\bm{0}$, namely:
\begin{equation}
    \frac{\partial H_{GL}}{\partial P}=P=0,\qquad \frac{\partial H_{GL}}{\partial r}=(-\alpha +\beta r^{2})r=0.
\end{equation}
This implies that
\begin{equation}
    r^{\pm}_{m}=\pm\sqrt{\frac{\alpha }{\beta }} \quad {\rm and} \quad P_{m}=0,\qquad r_{M}=0\quad{\rm and}\quad P_{M}=0,
\end{equation}
where the subscript ${m}$ stands for minima and ${M}$ for maximum and, finally, in terms of particle coordinates, we have 
\begin{equation}
    \left(\sum_{i=1}^{n}(q_{\pm}^{i})^{2}\right)^{1/2}=\pm\sqrt{\frac{\alpha }{\beta }}\quad {\rm and}\quad  p_{j}=0,\qquad q^{i}=0\quad {\rm and}\quad  p_{j}=0,\qquad\qquad\forall \; j\in[1,n].
\end{equation}
Therefore, the energy level set are defined by 
\begin{equation}\label{def:energy_level_sets}
    H_{GL}(P,r)=\frac{P}{2}-\frac{\alpha}{2}r^{2}+\frac{\beta}{4}r^{4}=E_{GL}.
\end{equation}
Not all the energy values $E_{GL}$ are associated with an accessible energy level set, in fact:
\begin{equation} 
    \forall \; E_{GL}<H_{GL}(P_{m},r^{\pm}_{m})=-\frac{\alpha ^2}{4 \beta }\implies \Sigma_{E}^{H_{GL}}= \emptyset \ .
\end{equation}
Therefore, the lowest energy value is $E_{GL}^{m}=-\alpha^{2}/4 \beta$ and the accessible level sets are defined by the following range of energies 
\begin{equation}\label{def:energy_range_quartic}
    E\in[-\alpha^{2}/4 \beta,\infty).
\end{equation}
Thus, the energy level sets corresponding to the energy values $H(P_{m},r^{\pm}_{m})\leq E_{GL}< H_{GL}(P_{M},r_{M})$ are homeomorphic to two disjoint hyperspheres:
\begin{equation}
    \Sigma_{E}^{H_{GL}}\simeq \mathbb{S}^{n-1}\cup \mathbb{S}^{n-1},
\end{equation}
As $E_{GL}\to 0$, the energy level set $\Sigma_{0}^{H_{GL}}$ is homeomorphic to the one-point-union of the two previous hyperspheres and this can be written through the wedge sum:
\begin{equation}
    \Sigma_{0}^{H_{GL}}\simeq \mathbb{S}^{n-1}\wedge \mathbb{S}^{n-1}= (\mathbb{S}^{n-1}\cup \mathbb{S}^{n-1})/\sim,
\end{equation}
where $\sim$ is the equivalence relation which identifies a point $x_{1}$ on the first hypersphere with the point $x_{2}$ on the second hypersphere.\\
Finally, for $0\leq E_{GL}< \infty$, the energy level sets are homeomorphic to hyperspheres; hence, we have
\begin{equation}
    \Sigma_{E}^{H_{GL}}\simeq \mathbb{S}^{n-1}.
\end{equation}
We have chosen the following allowed energy subset $I_{GL}=[-\alpha^{2}/4\beta,1]$ and we have sampled energy values from these sets with an energy-step $\Delta E=10^{-4}$.
Then, we have numerically solved the Hamilton's equations with $n=150$ particles adopting a second order bilateral symplectic algorithm \cite{casetti1995efficient}. We set $\alpha=0.5$ and $\beta=0.7$ which implies that $\alpha^{2}/4\beta=0.089$ and an integration time step $\Delta t=10^{-4}$. We note that the closest the energy values to the lowest one the larger the scalar and mean curvatures. In fact, since the energy level sets are always homeomorphic to a sphere or two disjoint spheres, in the $E\to-\alpha^{2}/4\beta$ limit, the level sets reduce their volume till becoming two points. Therefore, in order to avoid floating point overflow, we chose as minimum energy value $E_{GL}^{m}=-0.08$ which is slightly larger than $-\alpha^{2}/4\beta$. For any energy value $E$, random initial conditions have been chosen. Then, we have evaluated $\Omega_{\widetilde{g}}^{(1)}(E)$ along the dynamics at each energy value within $I_{GL}$ through the relation \eqref{def:time_average} and we used it for integrating Eq. (33) numerically obtaining $S_{\widetilde{g}}(E)$. The results are plotted in the main text (see Fig. 1). A least-sensitive inflection point in the entropy function, $S_{\widetilde{g}}$, appears at the transition point $\epsilon_t=0$. In particular, the first-order geometric function $\Omega_{\widetilde{g}}^{(1)}(E)$ has a positive-valued minimum to the left of $\epsilon_t$. The apparent peak in $\Omega_{\widetilde{g}}^{(1)}(E)$ is due to the fact that the GL-model undergoes a PT in correspondence of a critical point. In fact, this is an \emph{extreme} behavior of the function $\beta(E)$ provided by Bachmann (see Fig. 4 in the main text) where the maximum and the minimum points of $\beta$ lie on the black dashed vertical line.

\subsection{Geometric observables in Ginzburg-Landau-like model}

In the Ginzburg-Landau-like model, we define $x_i:=(p_i,q_i)$ and the gradient of the Hamiltonian \eqref{def:GL_hamiltonian} reads:
\begin{equation}
    \nabla_{p_i}H_{GL}=p_i, \qquad \nabla_{q_i}H_{GL}=-\alpha q_i+\beta q_i\left(\sum_{i=1}^{n}q_i^2\right),
\end{equation}
so that 
\begin{equation}
    \nabla_{x_i}H_{GL}=(\nabla_{p_i}H_{GL},\nabla_{q_i}H_{GL}).
\end{equation}
The Hessian matrix is given by 
\begin{equation}
    Hess\,H_{GL}=\left(\begin{matrix}
    1\!\!\!1_{\mathbb{R}^{n}} & \bm{0}_{\mathbb{R}^n}\\
    \bm{0}_{\mathbb{R}^n} & Hess\,V_{GL}
    \end{matrix}\right),
\end{equation}
with
\begin{equation}
    [Hess\,V_{GL}]_{ij}=\left[-\alpha+\beta\left(\sum_{i=1}^{n}q_i^2\right)\right]\delta_{ij}+2\beta\, q_i~q_j
\end{equation}
Moreover, we have
\begin{equation}
\begin{split}
    [\nabla\nabla\nabla\,V_{GL}]_{ijk}&=2\beta\,q_k \delta_{ij}+2\beta (\delta_{ik}q_{j}+\delta_{jk}q_{i}),\\
    [\nabla\nabla\nabla\nabla\,V_{GL}]_{ijkl}&=2\beta\,\delta_{lk} \delta_{ij}+2\beta (\delta_{ik}\delta_{lj}+\delta_{jk}\delta_{il}),\\
    [\nabla\nabla\nabla\nabla\nabla\,V_{GL}]_{ijklm}&=0.
\end{split}
\end{equation}
The laplacian of $H$ is given by taking the trace of the Hessian:
\begin{equation}
    \Delta H_{GL}=n+Tr[Hess\,V_{GL}]=n-\alpha n+\beta\left(\sum_{i=1}^{n}q_i^2\right)(n+2).
\end{equation}





\section{Proof of Universality of the first and second GCFs' energy behaviors}



In this section, we prove that the first and second order GCFs manifest an universal, that is, system-independent energy behavior. In other words, all the Hamiltonian systems which undergo a first (second) order PT à la Bachmann, admit the same qualitative energy behavior in the first (second) order GCF.\\
To reach such a purpose, it is actually sufficient to compute the GCF $\Omega_{\widetilde{g}}^{(1)}$ ($\Omega_{\widetilde{g}}^{(2)}$) for a specific system. Note that this can be easily done numerically. Then, exploiting the Riccati equation (38) and the Cauchy's theorem of existence and uniqueness of the solution of the ordinary differential equations, we can conclude that a different behavior of $\Omega_{\widetilde{g}}^{(1)}$ ($\Omega_{\widetilde{g}}^{(2)}$) cannot reproduce the behavior of $\partial_E^2 S$ predicted by Bachmann.\\
Then, let us start with the Riccati equation for the entropy that we recall below
\begin{equation}\label{eqn:entropy_flow_equation}
    \partial_E^2 S(E)+(\partial_E S(E))^2=\Omega_{\widetilde{g}}^{(2)}.
\end{equation}
then, let us plug the following transformation $S(E)=\log(f(E)/f_0)$ with $f(E_0)=f_0\in\mathbb{R}$ into the equation above. We get
\begin{equation}
\begin{split}\label{def:system_riccati}
    \partial_{E}f(E)= g(E),\qquad
    \partial_{E}g(E)= \Omega_{\widetilde{g}}^{(2)},
\end{split}
\end{equation}
We have now a set of linear first order differential equations which has to be solved for the unknown function $X=(f,g)$. Hence, Eq. \eqref{def:system_riccati} can be rewritten as
\begin{equation}\label{def:compact_riccati}
    \partial_E X(E)=A_{2}(E)\,X(E),
\end{equation}
where
\begin{equation}\label{def:A_matrix}
    A_{2}=\left(\begin{matrix}
    0 & 1\\
    \Omega_{\widetilde{g}}^{(2)} & 0
    \end{matrix}\right)
\end{equation}
By introducing the notation $F(E,X(E)):=A_{2}(E)\,X(E)$, a formal solution of Eq. \eqref{def:compact_riccati} is given by
\begin{equation}
\begin{split}\label{def:solution_R}
    X_{2}(E)=X_{2}(E_0)+\int_{E_0}^{E}F(s,X_{2}(s))~ds,
\end{split}
\end{equation}
Now, let us suppose by contradiction that there exists a $\widehat{\Omega}_{\widetilde{g}}^{(2)}$ with an energy behavior qualitatively different from that of $\Omega_{\widetilde{g}}^{(2)}$ but such that $\partial S^2_E$ admits a negative-valued peak as before. This means that, given the matrix
\begin{equation}\label{def:A_matrix_hat}
    A_{2}=\left(\begin{matrix}
    0 & 1\\
    \widehat{\Omega}_{\widetilde{g}}^{(2)} & 0
    \end{matrix}\right)
\end{equation}
together with the respective Riccati equation
\begin{equation}
    \partial_E \widehat{X}_{2}(E)=\widehat{A}_{2}(E)\widehat{X}_{2}(E).
\end{equation}
we have
\begin{equation}\label{eqn:same_solutions}
    X_{2}=\widehat{X}_{2}\equiv X.
\end{equation}
Hence, let us write the formal solutions
\begin{equation}
\begin{split}
    X_{2}(E)=X_{2}(E_0)+\int_{E_0}^{E}F(s,X_{2}(s))~ds,\qquad
    \widehat{X}_{2}(E)=\widehat{X}_{2}(E_0)+\int_{E_0}^{E}\widehat{F}(s,\widehat{X}_{2}(s))~ds,
\end{split}
\end{equation}
where $\widehat{F}(s,\widehat{X}_{2}(s)):=\widehat{A}_{2}(E)\widehat{X}_{2}(E)$ and let us subtract the first equation with the second one imposing Eq. \eqref{eqn:same_solutions}, we obtain:
\begin{equation}
    \int_{E_0}^{E}(F(s,X_{2}(s))-\widehat{F}(s,X(s)))~ds=0
\end{equation}
which implies $(A_{2}(E)-\widehat{A}_{2}(E))X(E)=0$ then, assuming that $X(E)$ is a non trivial solution, we get $A_{2}(E)=\widehat{A}_{2}(E)$ concluding that $\widehat{\Omega}_{\widetilde{g}}^{(2)}(E)=\Omega_{\widetilde{g}}^{(2)}(E)$, at least, in a neighborhood of the transition point. A similar result for the first-order GCF can easily deduced from the previous one.






\section*{References}
