\subsection*{Emergent Rokhsar-Kivelson point in realistic quantum Ising models}

\begin{center}
    \noindent\large{\textbf{Supplemental Materials}}
\end{center}

\vspace{1mm}

\begin{enumerate}[label=\Roman*.]
    \item Height field representation and the topological flux
    \item Exact correspondence between TFIM and the extended QDM
    \item Finite size scaling of the multicritical point
    \item Spin-spin correlation function
    \item Emergent $U(1)$ symmetry
    \item Histogram of the clock order parameter across RK-QCP
    \item Linear-quadratic crossover
    \item Curvatures of quadratic dispersion
    \item The intermediate regime
    \item Stochastic series expansion (SSE)
    \item Stochastic analytical continuation (SAC)
\end{enumerate}

\section{Height field representation and the topological flux}

In this section, we give an introduction to the height field representation and the topological flux of the quantum dimer model. 

In the height representation, each hexagonal plaquette is assigned an integer number $h(\mathbf r)$. Turning clockwise around a site of the even, the height $h(\mathbf r)$ changes by $+1$ when crossing an empty link, and by $-2$ when crossing an occupied link. The height field ressembles the vector potential $\mathbf A$ in the $U(1)$ gauge theory. Correspondingly, the magnetic field $\mathbf B=\nabla\times\mathbf A$ is perpendicular to the slope. 

When defined on a torus, the total change of the height field after going a period around the boundary determines the topological sector. Such sectors can be labelled by a pair of flux quantum numbers $(F_x,F_y)$. Going around a chosen loop~(Fig.~\ref{fig:s0}), we add $2$ when crossing an occupied dimer and $-1$ when crossing an empty link. The flux quantum numbers divided by the system size $L_x$, $L_y$ are defined as the magnetic flux. 
\begin{equation}
    f_x=F_x/L_x,\quad f_y=F_y/L_y
\end{equation}
Such flux also corresponds to the tilt of the height field $h$. For the ground states, the `tilt' is locked along one of the three dimer directions $\hat{\mathbf{e}}_i$~($i=1,2,3$). Therefore, without loss of generality, we can set $f_y=0$ only work with $f\equiv f_x$. The clock (plaquette) phase corresponds to zero-tilt $f=0$, and the stripe (staggered) phase corresponds to the maximum tilt $f=2$.

\begin{figure}[t!]
    \centering
    \includegraphics[width=0.3\linewidth]{fig_s0.pdf}
    \caption{An illustration of the height field. The red and blue points denote the even and odd lattice site on the dual lattice. The numbers in the center of hexagonal plaquettes illustrate the height field of the corresponding dimer configuration. The dark blue arrows denote the three unit vectors defined on the triangular lattice and its dual honeycomb lattice of dimers. The light orange and green paths denote the loop along which the winding numbers $F_x$ and $F_y$ are defined. }
    \label{fig:s0}
\end{figure}

\section{Exact correspondence between TFIM and the extended QDM}

For the $J_1$-$J_2$ TFIM, the transition from the clock to the stripe state is trivially first-order, suggesting $g_3>0$ at the transition. In the following we will establish the exact correspondence between our $J_1$-$J_2$-$J_3$ TFIM at $J_3=J_2/2$, and the extended QDM proposed in Ref. \cite{qdm_honey} at $v_0=v_3$. For the extended QDM, along the $v_0=v_3$ line the clock and stripe phases are sandwiched by a finite area of incommensurate states with intermediate flux. This indicates that  $g_3<0$ in the $J_3=J_2/2$ case. From the above we conclude that we can always find $J_3$ in the range $0<J_3<J_2/2$ where $g_3$ is fine-tuned to be zero, so that RK-QCP can be obtained at the clock-to-stripe transition.

The extended QDM proposed in Ref. \cite{qdm_honey} takes the form
\begin{equation}
H=v_0\hat{n}_0+v_3\hat{n}_3-t\sum_{\hexagon}(|\ha\rangle\langle\hb|+\mathrm{h.c.})
\end{equation}
where $\hat{n}_{\alpha}$ ($\alpha=0,1,2,3$) counts the total number of honeycomb plaquettes with $\alpha$ dimers, respectively. 

To to establish the mapping, let's consider a perturbation term above the nearest-neighbor TFIM on the triangular lattice
\begin{equation}
H'=J'\sum_i\left(\sum_{\langle ij \rangle _1}S_{j}^{z}\right)^{2}.
\end{equation}
On the one hand, this Hamiltonian can be decomposed into pairwise Ising interactions
\begin{equation}
    \begin{split}
        H' & =J'\sum_i\left(\sum_{\langle ij \rangle _1}(S_{i}^{z})^{2}+2\sum_{\langle ij \rangle _1 \langle ij' \rangle _1, j\neq j'}S_{i}^{z}S_{j}^{z}\right)\\
        & =2J'\sum_{\langle ij\rangle_{1}}S_{i}^{z}S_{j}^{z}+4J'\sum_{\langle ij\rangle_{2}}S_{i}^{z}S_{j}^{z}+2J'\sum_{\langle ij\rangle_{3}}S_{i}^{z}S_{j}^{z}+\mathrm{const.}
    \end{split}
    \label{eq:hpspin}
\end{equation}
The nearest neighbor terms is trivial that can be absorbed in $J_{1}$ and is ignored here. The main consequence of $H'$ is that it induces next-nearest- and third-nearest-neighbor interactions with $J_{2}=2J_{3}=4J'$. On the other hand, the term $\left(\sum_{\langle ij \rangle _1}S_{j}^{z}\right)^{2}=\left(\sum_{\langle ij \rangle _1}S_{i}^{z}S_{j}^{z}\right)^{2}$ is related to density of dimers in the plaquette centered at $\mathbf{r}_i$: for given $i$, $|\sum_{\langle ij \rangle _1}S_{j}^{z}|=3,2,1$ and $0$ corresponds to 0,1,2,3-dimer in that plaquette respectively. We have
\begin{equation}
    \begin{split}
        H'&=J'(0^2\hat{n}_3+1^2\hat{n}_2+2^2\hat{n}_1+3^2\hat{n}_0)\\
        &=2J'(\hat{n}_0+\hat{n}_3)+\mathrm{const.}
    \end{split}
    \label{eq:hpqdm}
\end{equation}
In the simplification we have used the dimer sum rule identity $\hat{n}_{0}+\hat{n}_{1}+\hat{n}_{2}+\hat{n}_{3}=N$ and $2\hat{n}_{0}+\hat{n}_{1}-\hat{n}_{3}=0$. Comparing \eqref{eq:hpspin} and \eqref{eq:hpqdm} we conclude that the our Ising spin model with $J_{2}=2J_{3}$ correspond to the extended QDM with $v_{0}=v_{3}=J_{2}/2$. 

A quantum Monte Carlo calculation verifies this correspondence. On the $J_2=2J_3$ line, an incommensurate phase region with the flux $0<f<2$ emerges at roughly $0.14<J_2<0.17$, or $0.69<J_2/\Delta<0.85$. On the other hand, the devil's stailcase on $V_0=V_3$ line locates at roughly $0.7<V_3/t<1$. The lower boundaries coincides very well, whereas the deviation of the upper boundary results from finite size effect. 

\begin{figure}[ht]
    \centering
    \includegraphics[height=0.3\linewidth]{fig_s1a.pdf}
    \includegraphics[height=0.3\linewidth]{fig_s1b.pdf}
    \caption{(a) The flux $f$ as a function of $J_2$ on the $J_2=2J_3$ line. $L=24$, $\beta=L$, $\Delta=0.2$ are taken. (b) The phase diagram of extendend dimer model taken from Ref. \cite{qdm_honey}. The $v_0=v_3$ line is highlighted. }
    \label{fig:s1}
\end{figure}

\section{Finite size scaling of the multicritical point}

The phase diagram in the main text is calculated at finite size $L=24$. To calculate the movement of the RK point when $L\rightarrow\infty$, we carried out a finite size scaling. The RK point is determined by determining the intersection of clock-stripe transition line and the stripe-incommensurate critical line, \textit{i.e.} $E(f=0)=E(f=2)=E(f=2-3/L_x)$. Through an extrapolation, we determine the position of the RK point in the thermaldynamic limit, which is $J_2=0.1093(11)$, $J_3=0.0429(7)$ (Fig. \ref{fig:s5}). 

\begin{figure}[ht]
    \centering
    \includegraphics[width=0.4\linewidth]{fig_s5a.pdf}
    \includegraphics[width=0.4\linewidth]{fig_s5b.pdf}
    \caption{The finite size scaling which gives the position of the RK point in the thermodynamic limit}
    \label{fig:s5}
\end{figure}

\section{Spin-spin correlation function}

We calculate the $S^z$-$S^z$ correlation at the multicritical point (Fig. \ref{fig:s2}). The correlations are measured along $x$-direction. Due to the fact that the RK wavefunction is the superposition of all the possible classical configurations, we can simultate the result of RK wavefunction by carrying out a classical Monte Carlo simulation of honeycomb-lattice dimer model at $\beta=0$. As cross-sector update is inaccesible in TFIM QMC, non-contractable loop updates are also banned in classical Monte Carlo simulation. These two correlation functions agrees considerably well, which confirms the RK property of the Ising tricritical point. 

However, as the correlation function oscillates as well as decays with the distance $r$, it's hard to determine the decaying behavior. So we need to turn to a correlation which decays monotonously, namely, the correlation of $U(1)$ spins. The $U(1)$ spin is defined on each triangle as a superposition of the $\mathbb{Z}_2$ spins on three sublattices:
\begin{equation}
    \psi=S_A^z+S_B^z\mathrm e^{\mathrm i2\pi/3}+S_C^z\mathrm e^{\mathrm i4\pi/3}
    \label{eq:clock}
\end{equation}
Such spin remains the same in single clock domain. Its correlation shows a algebraic decay behavior, which further confirms the RK nature of the Ising tricritical point. 

\begin{figure}[ht]
    \includegraphics[width=.6\linewidth]{fig_s2.pdf}
    \caption{The spin-spin correlation function along $x$-axis, calculated for both $J$-$J_2$-$J_3$ model at the RK point and for tricritical wavefunction. Inset: the $U(1)$ spin correlation along $x$-axis in $J$-$J_2$-$J_3$ model at the RK point plotted in double-logarithm scale. Parameters $L=48$, $\beta=L^2$, $\Delta=0.2$ are taken for $J$-$J_2$-$J_3$ model and $L=48$ is taken for RK wavefunction}
    \label{fig:s2}
\end{figure}

\section{Assymptotic behaviour of the correlators of RK wavefunction}

Here we evaluate the behaviour of the correlation functions $\langle n^x(0)n^x(\mathbf{r})\rangle$ and $\langle S^z(0)S^z(\mathbf{r})\rangle$ of RKQCP at large distance $|r|\rightarrow\infty$. The RK wave function is the equal weight superposition of classical dimer coverings $|\psi_\mathrm{RK}\rangle\sim\sum_{c}|c\rangle$. The observable with respect to this wavefunction $\langle\psi_\mathrm{RK}|\hat{\mathscr{O}}|\psi_\mathrm{RK}\rangle$ is equivalent to the statistical problem in classical dimer model at infinite temperature $T=\infty$, where the statistical weight of all dimer coverings are identical. 

For the dimer-dimer correlation function $G_d(\mathbf{R})=\langle n^x(0)n^x(\mathbf{R})\rangle$, such statistical problem has already been solved in Ref.~\cite{fisher1963statistical}, and the assymptotic behaviour is predicted to be $G_d(\mathbf{R})\sim|\mathbf{R}|^{-2}$. To measure the spin-spin correlation function $\langle S_i^zS_j^z\rangle$, we can construct a non-intersecting string $\mathscr{C}_{ij}=\{i,i_2,...,i_{N-1},j\}$ along nearest-neighbor triangular lattice bonds that terminates at $\mathbf{r}_i$ and $\mathbf{r}_j$. The spin-spin correlator $G_s(\mathbf{r}_i-\mathbf{r}_j)=\langle{S}_i^z{S}_j^z\rangle=2^{2-N}\langle({S}_i^z{S}_{i_2}^z)({S}_{i_2}^z{S}_{i_3}^z)\dots({S}_{i_{N-1}}^z{S}_j^z)\rangle$ can then be expressed in terms of density of dimers across the string, i.e., $\langle{S}_i^z{S}_j^z\rangle\sim\langle n^x_{ii_2}n^x_{i_2i_3}\dots n^x_{i_{N-1}j}\rangle$, where $n^x_{kl}$ denotes the dimer density on the link crossing the bond $\langle kl\rangle$. The assymptotic behaviour of such correlator is then predicted to be $G_s(\mathbf{r}_i-\mathbf{r}_j)\sim|\mathbf{r}_i-\mathbf{r}_j|^{-0.5}$.

\section{Emergent $U(1)$ symmetry}

To further confirm the emergent $U(1)$ symmetry, we measure the angular distribution of the VBS order parameter at the RK point and in the clock phase. We find that whereas in the clock phase, such distribution has three peaks with a $2\pi/3$ angle, at RK point, such distribution has no observable angular dependence, thus confirming the emergent $U(1)$ symmetry~(Fig.~\ref{fig:s8}). 

\begin{figure}
    \centering
    \includegraphics[width=.4\linewidth]{fig_s8.pdf}
    \caption{The angular probability density function (PDF) at point (i) in clock phase and (ii) at RK-DQCP. The system size of measurement is $L=24$}
    \label{fig:s8}
\end{figure}

\section{Histogram of the clock order parameter across RK-QCP}

\begin{figure}[ht]
    \centering
    \includegraphics[width=.8\linewidth]{fig_s3.jpg}
    \caption{(a) An illustration of the three lines on the phase diagram we scan. (b) The histogram with $U(1)$ symmetry at RK point. $L=24$, $\beta=L^2/2=288$, $\rho_s=2/3$ is taken throughout the calculation of histograms. (c) The histogram of $U(1)$ order parameter crossing the RK point, $J_2=0.10356$, $J_3=0.080,0.060,0.039(\mathrm{RK}),0.030,0.020$  respectively. (d) The histogram of $U(1)$ order parameter along the transition line, $J_2=0.08,0.09,0.10,0.10356(\mathrm{RK}),0.11$ respectively. (e) The histogram of $U(1)$ order parameter in the first order region, $J_2=0.07$, $J_3=0.014,0.008,0.006,0.005,0.004$ respectively.}
    \label{fig:s3}
\end{figure}

Another important quantity to evaluate the phase transition is the $U(1)$ order parameter defined in \ref{eq:clock}. Since QMC update cannot cross topological sectors, we restrict ourselves to the $f=0$ sector, in which stripe order is substituted by a quasi-long range order. 

The clock phase, due to its 6-fold degeneracy, is characterised by six distinct peaks connected by a $C_6$ rotational symmetry, whereas a centralized peak is seen in the disordered phase. When we go from the clock phase to the disordered phase, the six peaks shrinks into one centralized peak (Fig. \ref{fig:s3}, line 2). 

The distinction between first order transition and RK-type transition is that at the first order transition point, the $\mathbb Z_6$ angular dependence pertains, whereas at the RK point, the angular dependence of the $\mathbb Z_6$ peaks gradually fade away, and a $U(1)$ symmetry emerges, characterised by a ring in the histogram (Fig. \ref{fig:s3}b) \cite{yan_sweeping}. This can be seen clearly numerically that when we go along the transition line from the first order end to the RK point, the $\mathbb Z_6$ feature gradually disappears and the six peaks connect into one ring, and after crossing the RK point, the $U(1)$ symmetry is broken (Fig. \ref{fig:s3}, line 3). 

By contrast, in the first order region, when we enter the disordered region, central peak and $\mathbb Z_6$ peaks are detected simultaneously (Fig. \ref{fig:s3}, line 4), and the weight transfers from the $C_6$ peaks to the central peak. When we cross the transition line, central peak isn't detected immediately, as the clock phase is still a metastable state. Double peak only occurs when such metastable state collapses, in a region deep inside the disordered phase. 

\section{Linear-quadratic crossover}

When far away from the RK point, two linear modes are found at $\Gamma$ and $K$. When approaching the RK point, the nature of the effective gauge field changes from $(2+1)$d to $(2+2)$d. In the linear-quadratic crossover process (Fig. \ref{fig:s4}), the linear mode at $K$ softens into a quadratic mode, whereas the mode at $\Gamma$ gradually vanishes. 

\begin{figure}[ht]
    \centering
    \includegraphics[width=0.5\linewidth]{fig_s4.pdf}
    \caption{The linear-quadratic crossover in spin spectrum, measured at (a) $J_2=0.05$, (b) $J_2=0.08$, (c) $J_2=0.1035$ (RK point) and $J_3=0.0391$. $L=24$, $\beta=L^2/2=288$, $\rho_s=2/3$ is taken. }
    \label{fig:s4}
\end{figure}

\section{Curvatures of quadratic dispersions}

To explore the relation between the mode in spin correlation spectrum and the two modes in dimer correlation spectrum, we extracted the peak of intensity at different momentum and fit them to the form 
\begin{equation}
    \omega=\frac{1}{2}C_2(\mathbf q-\mathbf q_0)^2
    \label{eq:c2}
\end{equation}
in vicinity of the gapless point $\mathbf q_0$. The results are shown in Table~\ref{tbl:1}. We find the curvature of the pi0n* mode in spin spectrum roughly twice the curvature of pi0n in dimer spectrum.

\begin{table}[ht]
    \caption{The curvature defined in Eq. \ref{eq:c2} for different excitations measured in $J_1$-$J_2$-$J_3$ TFIM and RK-QDM. }
    \begin{tabular}{cc|cc}
        \hline\hline
        &&\multicolumn{2}{c}{Curvature $C_2$}\\
        \multicolumn{2}{c|}{Mode}&$J_1$-$J_2$-$J_3$ TFIM&RK-QDM\\
        \hline
        (\emph{Dimer spectrum})&Resonon&$0.080(3)$&$0.61(3)$\\
        &Pi0n&$0.057(4)$&$0.36(4)$\\
        (\emph{Spin spectrum})&Pi0n*&$0.095(2)$&$0.78(8)$\\
        \hline\hline
    \end{tabular}
    \label{tbl:1}
\end{table}

\section{The intermediate regime}

We measure the spin structure factor
\begin{equation}
    S(\mathbf{q})=\sum_{ij}\langle S_i^zS_j^z\rangle\mathrm{e}^{\mathrm{i}\mathbf{q}\cdot(\mathbf{r}_i-\mathbf{r}_j)}
\end{equation}
in the incommensurate phase. We fix $J_2=0.20$ and vary $J_3$ so that the system can be ordered at different momentum points. At sufficiently large $J_3$, the clock phase is stabilized. When gradually decrease $J_3$, we find the peak at $\mathrm{K}$ point splits into three peaks and gradually move towards $\mathrm{M}$ point~(Fig.~\ref{fig:s7}). The location of the peaks is connected with the flux of the corresponding ordered phase by 
\begin{equation}
    \mathbf{q}=\left(\pm\frac{2(1-2f)\pi}{3},\pm\frac{2\pi}{\sqrt{3}}\right),\left(\pm\frac{(4+f)\pi}{6},\pm\frac{(2-f)\pi}{\sqrt{3}}\right),\left(\pm\frac{(8-f)\pi}{6},\pm\frac{f\pi}{\sqrt{3}}\right).
\end{equation}
This is a clear hallmark of an intermediate ordered phase with intermediate tilt $0<f<2$. 

Such intermediate phase can be characterised by a $U(1)$ order parameter
\begin{equation}
    \psi_\mathbf{q}=\sum_{i}S_i\mathrm{e}^{\mathrm{i}\mathbf{q}\cdot\mathbf{r}_i},
\end{equation}
where $\mathbf{q}$ is chosen to be the momentum point at which the system is ordered. We find such order parameter distributes evenly on a ring~(Fig.~\ref{fig:s7}b). This is an evidence for the incommensurate nature of the order. The $U(1)$ symmetry emerges because when one translate the system by one unit along the $\mathbf{x}$ direction, the order parameter gains an extra incommensurate phase factor,
\begin{equation}
    \mathbf{r}\to\mathbf{r}+\hat{\mathbf{x}},\quad \psi\to\psi\mathrm{e}^{\mathrm{i}q_x}
\end{equation}

\begin{figure}
    \centering
    \includegraphics[width=0.6\linewidth]{fig_s7.pdf}
    \caption{(a) The spin structure factor and (b) the complex order parameter distribution in the incommensurate phase, measured at $J_2=0.20$ and (i) $J_3=0.113$ where $f=0.25$, (ii) $J_3=0.107$ where $f=0.75$ and (iii) $J_3=0.102$ where $f=1.25$. (c) An illustration of the structure factor peaks in the Brillouin zone.}
    \label{fig:s7}
\end{figure}

\section{Stochastic series expansion (SSE)}

For the numerical works in this paper, we use a quantum Monte Carlo (QMC) method with stochastic series expansion (SSE) algorithm \cite{qmc_1,qmc_2,qmc_3} to calculate the ground state properties and imaginary time Green function. This method will be briefly introduced below. 

In quantum statistics, the measurement of observables is closely related to the calculation of partition function $Z$
\begin{equation}
    \langle\mathscr{O}\rangle=\mathrm{tr}\:\left(\mathscr{O}\exp(-\beta H)\right)/Z,\quad Z=\mathrm{tr}\exp(-\beta H)
\end{equation}
where $\beta=1/T$ is the inverse temperature, $H$ is the Hamiltonian of the system and $\mathscr{O}$ is an arbitrary observable. Typically, in order to evaluate the ground state property, one takes a sufficiently large $\beta$ such that $\beta\sim L^z$, where $L$ is the system scale and $z$ is the dynamical exponent. In SSE, such evaluation of $Z$ is done by a Taylor expansion of the exponential and the trace is taken by summing over a complete set of suitably-chosen basis. 
\begin{equation}
    Z=\sum_\alpha\sum_{n=0}^\infty\frac{\beta^n}{n!}\langle\alpha|(-H)^n|\alpha\rangle
\end{equation}
We then write the Hamiltonian as the sum of a set of operators whose matrix elements are easy to calculate. 
\begin{equation}
    H=-\sum_iH_i
\end{equation}
In practice we truncate the Taylor expansion at a sufficiently large cutoff $M$ and it is convenient to fix the sequence length by introducing in identity operator $H_0=1$ to fill in all the empty positions despite it is not part of the Hamiltonian. 
\begin{equation}
    (-H)^n=\sum_{\{i_p\}}\prod_{p=1}^nH_{i_p}=\sum_{\{i_p\}}\frac{(M-n)!n!}{M!}\prod_{p=1}^nH_{i_p}
\end{equation}
and 
\begin{equation}
    Z=\sum_\alpha\sum_{\{i_p\}}\beta^n\frac{(M-n)!}{M!}\langle\alpha|\prod_{p=1}^nH_{i_p}|\alpha\rangle
\end{equation}

To the carry out the summation, a Monte Carlo procedure can be used to sample the operator sequence $\{i_p\}$ and the trial state $\alpha$ with according to their relative weight
\begin{equation}
    W(\alpha,\{i_p\})=\beta^n\frac{(M-n)!}{M!}\langle\alpha|\prod_{p=1}^nH_{i_p}|\alpha\rangle
\end{equation}
For sampling we adoopt a Metropolis algorithm where the configuration of one step is generated based on updating the configuration of the former step and the update is accepted at a probability
\begin{equation}
    P(\alpha,\{i_p\}\rightarrow\alpha',\{i'_p\})=\min\left(1,\frac{W(\alpha',\{i'_p\})}{W(\alpha,\{i_p\})}\right)
\end{equation}
Diagonal update, where diagonal operators are inserted into and removed from the operator sequence, and cluster update, where diagonal and off-diagonal operates convert into each other, are adopted in update strategy. 

In transverse field Ising model $H=J\sum_bS_{i_b}^zS_{j_b}^z-h\sum_i\sigma_i^x$, we write the Hamiltonian as the sum of following operators
\begin{equation}
    \begin{aligned}
        H_0&=1\\
        H_i&=h(S_i^++S_i^-)/2\\
        H_{i+n}&=h/2\\
        H_{b+2n}&=J(1/4-S_{i_b}^zS_{j_b}^z)
    \end{aligned}
\end{equation}
where a constant is added into the Hamiltonian for convenience. For the non-local update, a branching cluster update strategy is constructed \cite{qmc_2}, where a cluster is formed in $(D+1)$-dimensional by grouping spins and operators togather. Each cluster terminates on site operators and includes bond operators (Fig. \ref{fig:s6}a). All the spins in each cluster is flipped together at a probability $1/2$ agter all clusters are identified. 

\begin{figure}[b]
    \centering
    \includegraphics[width=0.6\linewidth]{fig_s6.pdf}    
    \caption{(a) An illustration of the SSE cluster update process where the cluster marked in red is identified and filpped as a whole. A vertical line shows a spin expanded. The solid and empty squares show the diagonal and off-diagonal site operators. The solid bars show the diagonal bond operators. (b) An illustration of the parametrization of a continuous spectral function into discrete $\delta$-functions. }
    \label{fig:s6}
\end{figure}

\section{Stochastic analytical continuation (SAC)}

For the spectra in this paper we adopted a stochastic analytical continuation (SAC) \cite{sac_1,sac_2,sac_3} method to obtain the spectral function $S(\omega)$ from the imaginary time correlation $G(\tau)$ measured from QMC, which is generally believed a numerically unstable problem. This method will be briefly introduced below. 

The spectral function $S(\omega)$ is connected to the imaginary time Green's function $G(\tau)$ through an integral equation 
\begin{equation}
    G(\tau)=\int_{-\infty}^{\infty}\mathrm{d}\omega S(\omega)K(\tau,\omega)
\end{equation}
where $K(\tau,\omega)$ is the kernal function depending on the temperature and the statistics of the particles. We restrict ourselves to the case of spin systems and with only positive frequencies in the spectral, where $K(\tau,\omega)=(e^{-\tau\omega}+e^{-(\beta-\tau)\omega})/\pi$. To ensure the normalization of spectral function, we further modity the transformation and come to the following equation :
\begin{equation}
    G(\tau)=\int_0^\infty\frac{\mathrm{d}\omega}{\pi}\frac{e^{-\tau\omega}+e^{-(\beta-\tau)\omega}}{1+e^{-\beta\omega}}B(\omega)
    \label{eq:kernal}
\end{equation}
where $B(\omega)=S(\omega)(1+e^{-\beta\omega})$ is the renormalized spectral function. 

In practice, $G(\tau)$ for a set of imaginary time $\tau_i(i=1,\cdots N_\tau)$ is measured in QMC simultation together with the statistical errors. The renormalized spectral function is parametrized into large number of equal-amplitude $\delta$-functions whose positions are sampled (Fig. \ref{fig:s6}b)
\begin{equation}
    B(\omega)=\sum_{i=0}^{N_\omega}a_i\delta(\omega-\omega_i)
\end{equation}
Then the fitted Green's functions $\tilde{G}_i$ from Eq. \ref{eq:kernal} and the measured Greens functions $\bar{G}_i$ are compared by the fitting goodness
\begin{equation}
    \chi^2=\sum_{i,j=1}^{N_\tau}(\tilde{G}_i-\bar{G}_i)(C^{-1})_{ij}(\tilde{G}_j-\bar{G}_j)
\end{equation}
where the covariance matrix is defined as 
\begin{equation}
    C_{ij}=\frac{1}{N_B(N_B-1)}\sum_{b=1}^{N_B}(G_i^b-\bar{G}_i)(G_j^b-\bar{G}_j),
\end{equation}
with $N_B$ the number of bins, the measured Green's functions of each $G_i^b$. 

A Metropolis process is utilized to update the series in sampling. The weight for a given spectrum is taken to follow a Boltzmann distribution
\begin{equation}
    W(\{a_i,\omega_i\})\sim\exp\left(-\frac{\chi^2}{2\Theta}\right)
\end{equation}
with $\Theta$ a virtue temperature to balance the goodness of fitting $\chi^2$ and the smoothness of the spectral function. All the spectral functions of sampled series $\{a_i,\omega_i\}$ is then averaged to obtain the spectrum as the final result. 

\end{document}
