\section{Adiabatic evolution of eigenstates in an ESQPT}
In Fig. \ref{supp-1} we show the adiabatic evolution of an eigenstate  
when it goes through an ESQPT. In the top panel of \ref{supp-1}, we plotted a region of the correlation diagram of Fig.~1 (of the main text) with eigenvalues $\phi_i$ that have overlap square $|c_i|^2>10^{-6}$. We have highlighted with red dots the adiabatic evolution of one of the eigenstates. In the bottom panels of Fig. \ref{supp-1}, we show the Husimi function of the eigenstates marked in the top panel with the corresponding label. It can be clearly seen that before the non-isolated avoided crossing, the eigenstates are localized in tori of rotational motion (panels 1, 2 and 3).
As we get closer to the middle of the avoided crossing, the state is located on the periodic orbit $z_0$ (panels 4 and 5) and then the eigenstates are located in tori  of libration motion (panels 6, 7 and 8).

\begin{figure}[h] 
    \centering
    \includegraphics[width=0.65\linewidth]{fig1_supp}
    \caption{Top panel: A region of the correlation diagram of Fig.~1(main text) with overlap square $|c_i|^2>10^{-6}$ with the resonance of the 
    PO $(q,p)=(0,0)$.  
    Bottom panel: Eigenstates adibatically crossing an non-isolated avoided crossing marked in the top panel. 
    \label{supp-1}}
\end{figure}

\section{Resonance $\ket{z_0}$ of the fixed point $z_0$}
A resonance of an unstable PO is a Gaussian beam constructed
in the neighborhood of the PO; see equation (124) in \cite{Vergini_2020}. Here, such construction reduces
to a wave packet centered at the fixed point $z_0=(0,0)$. We first
write the wave packet in the plane $q-p$
\begin{equation}
    W(q)=\left( \frac{\sinh \lambda}{\pi \hbar}\right)^{1/4}
    \exp \left[-\frac{q^2}{2 \hbar} \left(\sinh \lambda +
    {\rm i} \frac{k}{2} \right) \right].
\end{equation}
Later, one projects $W(q)$ on the unit torus
\begin{equation}
 \langle j | z_0 \rangle=W(j/N)+W(j/N-1) 
 \quad {\rm for} \quad j=0,\dots,N-1.
\end{equation}
\par
The mean value of the evolution operator in the state $\ket{z_0}$ 
can be estimated 
by the simple semiclassical expression
\begin{equation}
    \langle z_0 | \hat{U} | z_0 \rangle \simeq 
    \frac{{\rm e}^{{\rm i} \phi_{BS}}}{\sqrt{\cosh \lambda}},
\end{equation}
with the Bohr-Sommerfeld phase $\phi_{BS}$ being the action of the 
map at the fixed point $z_0$, divided by $\hbar$.

\section{Semiclassical smoothing of the spectral decomposition of $\ket{z_0}$}
We rewrite equation (5) of reference \cite{Vergini_2017} for the 
case of an are-preserving map in the unit torus (we use the 
notation introduced in the letter)
\begin{equation}
    \sum |c_i|^2 g\left(\frac{2 (\phi-\tilde{\phi}_i)}{\lambda \beta}\right)\simeq K
    \tilde{F}(x) \left( 1 + \frac{1}{\lambda \sqrt{N}} 
  \sum \frac{\cos(\psi_j)}{\sqrt{A_j |L_j|}} \right),
\end{equation}
where the sum runs over the set of HOs with relevance $A_j< A_{max}$,
with $\beta=2/\ln(2 \pi N A_{max} )$ and $L_j$ the Lazutkin canonical
invariant. Moreover, $K=\sqrt{\pi/8} \beta/(1+\beta)$ and  
\begin{equation}
 g(y)=\frac{\sin(y)/y+\beta \cos(y)}{(1+\beta^2 y^2) (1+\beta)}
\quad {\rm for}\quad y\ne 0\quad {\rm and} \quad g(0)=1.
\end{equation}
\par
As $A_{max}$ increases, one has to include more HOs into the sum.
Then the fuction $g(2 (\phi-\tilde{\phi}_i)/\lambda \beta)$ is
closer to a delta function
and the eigenphases are better defined. 
\par
\begin{figure}[h] 
    \centering
    \includegraphics[width=0.65\linewidth]{fig2_supp}
    \caption{The points are the calculated values of $\Delta S $ and the solid line corresponds to the estimated function of Eq.~(6) in the main text.
    \label{figDeltaS}}
\end{figure}
\section{Approximation of $\Delta S$ }
In Fig.~\ref{figDeltaS} we show the estimation Eq.~(6) (main text) of the canonical invariant $\Delta S=S_2-S_1$ (black shaded area in  Fig.~1(f)) as a function of $k$  
\section{The functions $\eta (x)$ and $\tilde{F}(x)$  }
\begin{figure}[h] 
    \centering
    \includegraphics[width=0.45\linewidth]{fig3_supp_a}
    \includegraphics[width=0.45\linewidth]{fig3_supp_b}
    \caption{Error in the approximation of $\varphi(x)$ and 
    $\tilde{F}(x)$.
   
    \label{varphiFourier}}
\end{figure}

The functions $\eta(x)$ and $\tilde{F}(x)$ are defined
in terms of the Fourier transform of
other functions \cite{Vergini_2017}. 
For this reason, in order to provide simple expressions here we  present accurate interpolation formulas.
\par
The function $\varphi(x)= x \eta(x)$ is the phase of the complex function
\begin{equation}
 \tilde{f}(x)= \frac{1}{\sqrt{\pi}}\int_{-\infty}^{\infty} {\rm e}^{-y/2} K_0({\rm e}^{-y}) {\rm e}^{{\rm i} x y} {\rm d}y,
 \end{equation}
with $K_0$ the modified Bessel function of zeroth order.
That is $ \tilde{f}(x)=  |\tilde{f}(x)| {\rm e}^{{\rm i} \varphi(x)}$,
where we have selected the branch of the phase satisfying $\varphi(0)=0$.
The function $\eta(x)$ verifies the following asymptotic behaviours
for $x \to 0$ and $|x| \to \infty$
\begin{equation}
   \eta(x)=\eta(0) -a x^2 + b x^4 + \order{x^6} \quad {\rm and} \quad
   \eta(x)= - \ln |x|+1+\frac{\pi}{4 |x|}+\order{1/x^2},
\end{equation}
with $\eta(0)\approx 3.53430\;63528$, $a\approx 5.38865\;58307$
and $b\approx 12.80436\;182$.
Later, we propose an interpolation formula satisfying 
the previous behaviours
\begin{equation}
    \eta(x)\approx -\ln \sqrt{x^2+1/2a}+1+A 
    \left(1+B x^4\right)^{-1/4}+ C 
    \left(1+z_1 x^6 \right)^{z_2 \ln(z_3+x^2)},
\end{equation}
with $z_1=423$, $z_2=-0.4337$ and $z_3=1.78$ free parameters fitted
to minimize the error of the approximation for $\varphi(x)$ in 
the intermediate region. 
In contrast, the constants $A$, $B$ and $C$ are related 
to $\eta(0)$, $a$ and $b$ as follows
\begin{equation}
    B=\left[ \frac{16}{\pi} (a^2-b) \right] ^{4/5}\approx 34.190,
\end{equation}
$A=\pi B^{1/4} /4 \approx 1.8992$, and $C=\eta(0)-\ln\sqrt{2a}-A-1\approx-0.5536$.
\par
The function $\tilde{F}(x)$ is defined by the integral equation
\begin{equation}
   \tilde{F}(x)= \sqrt{\frac{2}{\pi}} \int_0^{\infty} \frac{\cos(x y)}{\sqrt{\cosh(y)}} {\rm d}y,
\end{equation}
with asymptotic behaviours for $x \to 0$ and $|x| \to \infty$
\begin{equation}
 \tilde{F}(x)= \tilde{F}(0)- c x^2+\order{x^4} \quad {\rm and}  
 \quad  \tilde{F}(x)= \sqrt{\frac{2}{|x|}} {\rm e}^{-\pi |x| /2} 
 \left(1+\order{1/x^2}\right),
\end{equation}
with $\tilde{F}(0)\approx 2.09209\;92401$ and 
$c\approx 8.99462\;98154$. Later, we propose an interpolation 
formula satisfying the previous behaviors
\begin{equation}
 \tilde{F}(x)\approx  \frac{(D^2+x^2)^{-1/4}}{\sqrt{\cosh(\pi x)}}
 +E (1+z_4 x^4)^{- \ln \sqrt{z_5+x^2}},
\end{equation}
with $z_4=103$ and $z_5=1.8$ free parameters fitted
to minimize the error of the approximation for $\tilde{F}(x)$ 
in the intermediate region. 
In contrast, the constants $D\approx 0.31247$ and 
$E\approx 0.30316$ were derived from the relationships
$\tilde{F}(0)=E+1/\sqrt{D}$ and $c=(1+D^2\pi^2)/(4 D^{5/2})$.
\par
Figure~\ref{varphiFourier} shows the error of the previous
approximations.
\par
\section{Detailed phase space portrait for $k=1.447$}
The infinite family of satellite POs of a given HO approximate
more and more the motion of the HO as the period of the PO goes
to infinity. This family consists of unstable POs with Bohr-Sommerfeld
phase well approximated by the condition $\psi_j(\phi)=2 \pi n$, with
$\psi_j(\phi)$ the pase of the given HO.
We have found that the shortest satellite PO of the first HO
generates a strong scar on the eigenstate number $3$ of 
Fig.~3 (of the main text), and the
shortest satellite PO of the second HO generates a strong scar
on the eigenstate number $4$. Fig.~\ref{funcsupmat} shows clearly 
these satellite POs, and their strong correlation with the 
Hussimis of states $3$ and $4$. 


\begin{figure}[h] 
    \centering
    \includegraphics[width=0.65\linewidth]{fig4_supp_light} 
    \caption{Eigenfunctions No. 3 and 4. from Fig.~3 of the main text, with the corresponding satellite orbit in black symbols and the unstable (solid red) and stable (solid blue) manifolds. The yellow dots represent the classical phase space for $k=1.447$.
   
    \label{funcsupmat}}
\end{figure}

\section{Scaling of the participation ratio with the system size}
In Fig.~3 of the main text using a rescaling of the IPR $\xi$ we could show that a transition effectively takes place at an estimated value $k_{\text{break}}$. For clarity, and completeness in Fig.~\ref{PR} (left) we show a smoothing of the data obtained for $\xi$ (with no normalization) as a function of the (unrescaled) perturbation $k$, for several values of $N$. We observe that as $N$ grows, the step in the curves become steeper (and for smaller $k$). In fact, the -- correspondingly colored -- dots mark the $k_{\text {break}}$ value, where a clear displacement is observed, as it was shown in the inset of Fig.~3 [main text].

A complementary view of the IPR is provided by its inverse, i.e. the participation ratio or participation number (PR). It is usually associated with the area of Hilbert space occupied. If the system is fully chaotic then this saturates to a constant value proportional to $N$. 

In our case, the area filled by the chaotic layer is much smaller than unity (the area of the phase space) in 
the considered range of $k$. One should use an effective size $N_{\rm effective}$ instead of $N$, with
$N_{\rm effective}\approx N A_{\rm chaotic}$, and  where $A_{\rm chaotic}$ is the classical area of the chaotic layer which depends on $k$.
An elementary estimation, for small $k$, gives  
 $A_{\rm chaotic}\propto\, \Delta S  /\lambda$, where $\Delta S$ is the area of each
lobe (as represented in Fig.~1f), and $1/\lambda$ proportional to the number of
lobes along the broken separatrix. To clarify this point, in Fig.~\ref{PR} (right) we show
 $PR/N$ as a function of $k$ for a wide range of values  of $N$. It can be clearly observed that
for small $k$ values the scaling of $PR/N$ collapses onto the curve $\alpha \Delta S /\lambda$ (where $\alpha$ is an fitting prefactor). Furthermore, we observe that
the $PR$ of all the curves saturates at $N_{\rm effective}$ in accordance with our previous prediction.
On the other hand, in the small $k$ limit  $A_{\rm chaotic}< \xi_N$ (where $\xi_N$ is the normalization used in Fig.~3), and so  we have $PR/N\sim \xi_N/N$.
\begin{figure}[h] 
    \centering
    \includegraphics[width=0.85\linewidth]{fig5_supp} 
    \caption{(Left) Smoothed IPR $\xi$ as a function of $k$ for $N=200,400,1000,3000,5000, 10000$ (darkest shade is larger $N$). The dots correspond to the value of $k_{\rm break}$ obtained from $\Delta S/\hbar=3\pi/2$ [see main text].  (Right) $PR/N$ as a function of k for $N=200, 400,1000,3000,5000,10000$ (darker color means larger $N$). The red line represents $\alpha\Delta S/\lambda$ (we took $\alpha=14$). For all the curves a running average was performed to smooth out the fluctuations.
    \label{PR}}
\end{figure}
\section{Complementary data }
In this section we present data that supports the discussion presented in the main text. In Fig.~\ref{peine1026} we show the  logarithmic singularity in the level spacings  characteristic of the ESQPT, for a larger $N=1026$ value and the same perturbation $k=0.5$ as shown in Fig.~2 of the main text. The spectrum considered corresponds to the largest overlaps of the initial state $\ket{z_0}$ with the eigenstates of the map $U$. To compute the phase differences $\Delta\phi=\phi_{i+1}-\phi_i$, an unfolding procedure is necessary to account for contributions from different { Demkov structures  (see Fig.~2, main text)}. We note that  the semiclassical calculation  improves appreciably. This is expected because for $k<< k_{\text break}$ the error is of order $1/N$. For the perturbation value considered here we can estimate, by solving for $N=3/(4\Delta S)$ that the characteristic (singular) structure  will break at $N\approx 62900$.
\begin{figure}[h] 
    \centering
    \includegraphics[width=0.65\linewidth]{fig6_supp} 
    \caption{ Level spacing singularity characteristic of ESQPT, for the standard map  with $N=1026$, $k=0.5$. Circles and crosses indicate quantum and 
    semiclassical calculation, respectively. The solid line is a cubic spline interpolation of the semiclassical calculation.
    \label{peine1026}}
\end{figure}

\end{document}	

