\section{QMC methods}
\subsection{Finite temperature auxiliary field QMC method}
As mentioned in the main text, we investigate an half-filling $\pi$-flux extended Hubbard model on square lattice, its Hamiltonian is $H=H_0+H_U$ with non-interacting part $H_0=-\sum_{\langle i j\rangle, \alpha} t_{i j}\left(c_{i \alpha}^{\dagger} c_{j \alpha}+\text { H.c. }\right)$ and interacting part $H_U=U \sum_{\square}\left(n_{\square}-2\right)^{2}$. 
Since $H_0$ and $H_U$ do not commute, we use Trotter decomposition to separate $H_0$ and $H_U$ in the imaginary time propagation
\begin{equation}
	Z=\operatorname{Tr}\left[\mathrm{e}^{-\beta(H-\mu N)}\right]=\operatorname{Tr}\left[\left(\mathrm{e}^{-\Delta_{\tau} H_{U}} \mathrm{e}^{-\Delta_{\tau} H_{t}}\right)^{M}\right]+\mathcal{O}\left(\Delta_{\tau}^{2}\right),
\end{equation}
where $Z$ is the partition function, $\beta=M\Delta\tau$ is the inverse temperature, $\mu$ is chemical potential and equals $0$ at half-filling.
In QMC method, we can only deal with quadratic fermionic operator, while $H_U$ contains the quartic term, thus we should employ a Hubbard-Stratonovich (HS) decomposition as
\begin{equation}
e^{-\Delta\tau U(n_{\square}-2)^{2}}=\frac{1}{4}\sum_{\{s_{\square,\tau}\}}\gamma(s_{\square,\tau})e^{-2\alpha\eta(s_{\square,\tau})}e^{\alpha\eta(s_{\square,\tau})n_{\square}}
\label{eq:decompo}
\end{equation}
with $\alpha=\sqrt{-\Delta\tau U}$, $\gamma(\pm1)=1+\sqrt{6}/3$,
$\gamma(\pm2)=1-\sqrt{6}/3$, $\eta(\pm1)=\pm\sqrt{2(3-\sqrt{6})}$,
$\eta(\pm2)=\pm\sqrt{2(3+\sqrt{6})}$ and the sum symbol is taken over the auxiliary fields $s_{\square,\tau}$ on each $\tau$-th time slice square plaquette. Now, the interacting part is transformed into quadratic term but coupled with an auxiliary field.
Following simulations are based on the single-particle basis $\boldsymbol{c} = \{c_1, c_2 \cdots c_{N} \}$, $N=L\times L$ is the system size, so we can use the matrix notation $K$ and $V$ to represent $H_0$ and $H_U$ operators. 
We define the imaginary time propagators 
\begin{equation}
\begin{aligned}
&U_{s_{\square}}\left(\tau_{2}, \tau_{1}\right)=\prod_{m=m_{1}+1}^{m_{2}} e^{\boldsymbol{c}^{\dagger} V\left(s_{\square,m\Delta\tau}\right) \boldsymbol{c}} e^{-\Delta_{\tau} \boldsymbol{c}^{\dagger} K \boldsymbol{c}}, \\
&B_{s_{\square}}\left(\tau_{2}, \tau_{1}\right)=\prod_{m=m_{1}+1}^{m_{2}} e^{V\left(s_{\square,m\Delta\tau}\right)} e^{-\Delta_{\tau} K},
\end{aligned}
\end{equation}
where $m_1 \Delta\tau=\tau_1$ and $m_2 \Delta\tau=\tau_2$. 
Then partition funciton $Z$ can be rewritten as 
\begin{equation}
Z=\sum_{\{s_{\square,\tau} \}} \operatorname{Tr}\left[U_{s_{\square}}(\beta, 0)\right] \prod_{m=1}^{M} \gamma(s_{\square,m\Delta\tau})e^{-2\alpha\eta(s_{\square,m\Delta\tau})} = \sum_{\{s_{\square,\tau} \}} \operatorname{det}\left[1+B_{s_{\square}}(\beta, 0)\right] \prod_{m=1}^{M} \gamma(s_{\square,m\Delta\tau})e^{-2\alpha\eta(s_{\square,m\Delta\tau})}.
\end{equation}
Physical observables are measured according to 
$\langle O \rangle = \frac{\operatorname{Tr}\left[\mathrm{e}^{-\beta H} O\right]}{\operatorname{Tr}\left[\mathrm{e}^{-\beta H}\right]}$.
The equal-time single-particle Green function is given by
$
G_{i,j}(\tau,\tau)=\left\langle c_{i,\tau} c_{j,\tau}^{\dagger}\right\rangle=\left(1+B_{s_{\square}}(\tau, 0) B_{s_{\square}}(\beta, \tau)\right)_{i, j}^{-1}
$, and the dynamical single-particle Green function is given by
$G_{i,j}(\tau_1,\tau_2)=\left\langle c_{i,\tau_1} c_{j,\tau_2}^{\dagger}\right\rangle=-\left[\left(\mathbf{1}-\mathbf{G}\left(\tau_{1}, \tau_{1}\right)\right) B_{s_{\square}}^{-1}\left(\tau_{2}, \tau_{1}\right)\right]_{i, j}$.
Other physical observables can be calculated from single-particle Green function through Wick's theorem. More technical details of the finite-temperature QMC algorithms can be found in the reference book~\cite{assaadWorld-line2008}.

\subsection{Projection QMC method}
Since we also want to investigate the ground state properties, the projection quantum Monte Carlo (PQMC) method is a good choice. 
We can calculate the ground state wave function $\vert \Psi_0 \rangle$ through the projection of a trial wave function $\vert \Psi_T \rangle$ as $\vert \Psi_0 \rangle = \lim\limits_{\Theta \to \infty} e^{-\frac{\Theta}{2} \mathbf{H}} \vert \Psi_T \rangle$, where $\Theta$ is the projection length. 
Actually, in PQMC method, we use $\langle\Psi_{T}|e^{-\Theta H}|\Psi_{T}\rangle$ to replace the role of partition function in finite-temperature version.
Physical observables $O$ can be measured according to
\begin{equation}
\label{eq:observablepqmc}
\langle O \rangle = \frac{\langle \Psi_0 \vert O \vert \Psi_0 \rangle}{\langle \Psi_0 \vert \Psi_0 \rangle} 
						= \lim\limits_{\Theta \to \infty} \frac{\langle \Psi_T \vert  e^{-\frac{\Theta}{2} \mathbf{H}} O  e^{-\frac{\Theta}{2} \mathbf{H}} \vert \Psi_T \rangle}{\langle \Psi_T \vert  e^{-\Theta \mathbf{H}} \vert \Psi_T \rangle} .
\end{equation}

After the Trotter and HS decomposition, we have
\begin{equation}
\langle\Psi_{T}|e^{-\Theta H}|\Psi_{T}\rangle=\sum_{\{s_{\square,\tau}\}}\left[\left(\prod_{\tau}\prod_{\square}\gamma(s_{\square,\tau})e^{-2\alpha\eta(s_{\square,\tau})}\right)\det\left[P^{\dagger}B(\Theta,0)P\right]\right]
\label{eq:mcweight}
\end{equation}
here $P$ is the coefficient matrix of trial wave function $|\Psi_T\rangle$; $B$ matrix is defined as
$B(\tau+1,\tau)=e^{V_\tau}e^{-\Delta_\tau K}$
and has a property $B(\tau_3,\tau_1)=B(\tau_3,\tau_2)B(\tau_2,\tau_1)$.
In the practice,  we choose the ground state wavefunction of the half-filled non-interacting parts of Hamiltonian as the trial wave function. 
The Monte Carlo sampling of auxiliary fields are further performed based on the weight defined in the sum of  Eq.~\eqref{eq:mcweight}. Single particle observables are measured by Green's function directly and many body correlation functions are measured from the products of single-particle Green's function based on their corresponding form after Wick-decomposition. 
The equal time Green's function are calculated as
\begin{equation}
G(\tau,\tau)=1-R(\tau)\left(L(\tau)R(\tau)\right)^{-1}L(\tau)
\end{equation}
with $R(\tau)=B(\tau,0)P$, $L(\tau)=P^{\dagger}B(\Theta,\tau)$. 
The imaginary-time displaced Green's function  
$G(\tau,0) \equiv \left\langle \mathbf{c} \left(\frac{\tau}{2}\right) \mathbf{c}^{\dagger}\left(-\frac{\tau}{2}\right)\right\rangle$
are calculated as
\begin{equation}
\left\langle \mathbf{c} \left(\frac{\tau}{2}\right) \boldsymbol{c}^{\dagger}\left(-\frac{\tau}{2}\right)\right\rangle=\frac{\left\langle\Psi_{T}\left|e^{-\left(\frac{\Theta}{2}+\frac{\tau}{2}\right) H} \boldsymbol{c} e^{-\tau H} \boldsymbol{c}^{\dagger} e^{-\left(\frac{\Theta}{2}-\frac{\tau}{2}\right) H}\right| \Psi_{T}\right\rangle}{\left\langle\Psi_{T}\left|e^{-\Theta H}\right| \Psi_{T}\right\rangle}
\end{equation}
More technique details of PQMC method, please also refer to Refs~\cite{assaadWorld-line2008}. 
\subsection{Infinite $U$ QMC algorithm}
For infinite $U$ calculation, we use following formula to perform the $SU(N)$ infinite $U$ projection~\cite{ouyangProjection2021}
\begin{equation}
 \left.e^{-\Delta_\tau U (n_{\square}- 2)^2}\right|_{U\rightarrow+\infty} 
=  \frac{1}{M}\sum_{s_{\square,\tau}=1}^M e^{\frac{\text{i} 8\pi s_{\square,\tau}}{M} (n_{\square}- \frac{N}{2})},
\end{equation}
where we can set $M=2N+1$.
As the infinite $U$ term has now been replaced by fermion bilinears coupled to auxiliary fields, we can further use above finite temperature and projection QMC scheme to perform the calculation.

\begin{figure}[htb!]
\includegraphics[width=1.0\columnwidth]{figS1}
\caption{(a) $1/L$ extrapolation of single-particle gap. (b) $1/L$ extrapolation of VBS structure factor. The dash lines are obtained by cubic polynomial fitting with least squares method. (c) The AF structure factors as function of interaction strength. Error bars are smaller than the symbols.}
\label{fig:AF-order}
\end{figure}

\section{Other QMC data}
\subsection{The single-particle gap and VBS order parameter near QCP of SM-VBS }
The Dirac SM to VBS phase transition is a transition from massless Dirac fermion to insulator, thus single-particle gap $\Delta_{\text{sp}}$ will open at Dirac cones, i.e. at the momentum point $\mathbf{K}_0$.
We can extract $\Delta_{\mathrm{sp}}(\mathbf{K}_0,L)$ from a fit to the asymptotic long imaginary time behavior of the single-particle Green's function
$G(\mathbf{k}, \tau, L)=\left(1/L^{4}\right) \sum_{i, j, \sigma} e^{i \mathbf{k} \cdot\left(\mathbf{r}_{i}-\mathbf{r}_{j}\right)}\left\langle c_{i, \sigma}^{\dagger}(\tau) c_{j, \sigma}(0)\right\rangle \propto e^{-\Delta_{\mathrm{sp}}(\mathbf{k},L) \tau}$, then extrapolate it to the thermodynamic limit (TDL).
As shown in Fig.~\ref{fig:AF-order}~(a), it is clear that $\Delta_{\mathrm{sp}}(\mathbf{K}_0,L\rightarrow\infty) \rightarrow 0$ at $U/t=32$ in the SM phase, $\Delta_{\mathrm{sp}}(\mathbf{K}_0,L\rightarrow\infty) > 0$ at $U/t=38$ in the VBS phase, and  $\Delta_{\mathrm{sp}}(\mathbf{K}_0,L\rightarrow\infty) \sim 0$ at $U_c/t=36$.
The $1/L$ extrapolation of $C_{\text{VBS}}$ has the similar behavior as $\Delta_{\mathrm{sp}}(\mathbf{K}_0,L)$ near $U_c$, as shown in Fig.~\ref{fig:AF-order}~(b).

\subsection{The absence of AF order}
In Hubbard-like model, the AF order usually dominate when VBS order become weak in strong coupling region~\cite{xuKekule2018,liaoGross2022}. 
Here, we demonstrate that the AF order is absent in our model. 
We define the structure factor of AF order as 
$
C_\text{AF} (\mathbf{k}, L) \equiv \frac{1}{L^4}\sum_{\mathbf{r}} e^{-i\mathbf{q}\cdot\mathbf{r}}\sum_{\alpha,\beta}\langle S^{\alpha}_{\beta}(i) S^{\beta}_{\alpha}(j)\rangle
$, $C_\text{AF}$ will peak at momentum point $\mathbf{M}$ for AF order.
We can perform $1/L$ extrapolation of $C_\text{AF}(\mathbf{M},L)$ to get the AF structure factor at TDL, markded as $C_\text{AF}(\mathbf{M},\infty)$.
As shown in Fig.~\ref{fig:AF-order} (c), $C_\text{AF}(\mathbf{M},\infty)$ vanish at whole interaction strength range $U/t \in [0,160]$, which mean the AF order is absent in our model.


\section{Dynamical spin spectrum of free $\pi$-flux model}
We have a free $SU(4)$ $\pi$-flux model on square lattice with Hamiltonian
\begin{equation}
	H=- \sum_{\langle i j\rangle, \alpha} t_{i j}\left(c_{i \alpha}^{\dagger} c_{j \alpha}+\text { H.c. }\right).
\end{equation}
, and consider two-site unit cell, with inner cell coordinates $\mathbf{v}_{1}=(0,0)$,
$\mathbf{v}_{2}=(1,0)$.
We transform the Hamiltonian to the momentum space, 
\begin{equation}
H=-t\sum_{\mathbf{k}}c_{\mathbf{k}}^{\dagger}\left(\begin{array}{cc}
2\cos k_y & 1+e^{-2ik_{x}}\\
1+e^{2ik_{x}} & -2\cos k_y
\end{array}\right)c_{\mathbf{k}}
\end{equation}
where $c_{\mathbf{k}}=(c_{1,\mathbf{k}},c_{2,\mathbf{k}})^{T}$. We further write
$H$ in the diagonal basis
\begin{align}
H & =\sum_{\mathbf{k},a,b,m,n}c_{a,\mathbf{k}}^{\dagger}(U_{\mathbf{k}})_{am}(D_{\mathbf{k}})_{mn}(U_{\mathbf{k}}^{-1})_{nb}c_{b,\mathbf{k}}\nonumber \\
 & \equiv\sum_{\mathbf{k},a}\epsilon_{a,\mathbf{k}}f_{a,\mathbf{k}}^{\dagger}f_{a,\mathbf{k}}
\end{align}
where we use $a,b,m,n\in [1,2] $ to denote sublattice index, $(D_{\mathbf{k}})_{mn}=\delta_{mn}\epsilon_{m,\mathbf{k}}$, $\epsilon_{m,\mathbf{k}}=\pm 2 \sqrt{\cos^2 k_x + \cos^2 k_y}$ is the dispersion,
and we define $f_{a,\mathbf{k}}\equiv U_{\mathbf{k}}^{-1}c_{a,\mathbf{k}}$. Here the dispersion
is flavor degenerate, and we have omitted the flavor index, . The sum
is over $\mathbf{k}$ in the small (rectangle) BZ, i.e. $k_{x}\in(-\frac{\pi}{2},\frac{\pi}{2})$ and $k_{y}\in(-\pi,\pi$).

To calculate spin spectrum, we can use fluctuation-dissipation theorem,
which relates the spin spectrum to the spin susceptibility. The spin
susceptibility has following form

\begin{equation}
\chi_{ab}^{R}(\mathbf{r}-\mathbf{r}',t-t')=-i\theta(t-t')\sum_{\mu\nu}\langle[S_{\nu}^{\mu}(\mathbf{r},a,t),S_{\mu}^{\nu}(\mathbf{r}',b,t')]\rangle .
\end{equation}
Fourier transform
to momentum space, we get
\begin{equation}
\chi_{ab}^{R}(\mathbf{q},t-t')=-i\theta(t-t')\frac{1}{V}\sum_{\mu\nu}\langle[S_{\nu}^{\mu}(\mathbf{q},a,t),S_{\mu}^{\nu}(-\mathbf{q},b,t')]\rangle
\end{equation}
, where $S^{\alpha}_{\beta}(\mathbf{k},a,t)=c^{\dagger}_{\mathbf{k},\alpha}c_{\mathbf{k},\beta}-\frac{\delta_{\alpha,\beta}}{4}\sum_{\gamma}c^{\dagger}_{\mathbf{k},\gamma}c_{\mathbf{k},\gamma}$ are the $SU(4)$ spin full operators with $\alpha,\beta,\gamma \in [1,4]$.

For non-interacting $SU(4)$ case, the time dependence of spin operator is
given by

\begin{align}
S_{\nu}^{\mu}(\mathbf{q},a,t) & =\sum_{\mathbf{k}}e^{iHt}(c_{a,\mathbf{k},\mu}^{\dagger}c_{a,\mathbf{k}+\mathbf{q},\nu}-\frac{\delta_{\mu\nu}}{N}\sum_{\lambda}c_{a,\mathbf{k},\lambda}^{\dagger}c_{a,\mathbf{k}+\mathbf{q},\lambda})e^{-iHt}\nonumber \\
 & =\sum_{\mathbf{k}}\sum_{mn}(U_{\mathbf{k}}^{-1})_{ma}(U_{\mathbf{k}+\mathbf{q}})_{an}(f_{m,\mathbf{k},\mu}^{\dagger}f_{n,\mathbf{k}+\mathbf{q},\nu}-\frac{\delta_{\mu\nu}}{N}\sum_{\rho}f_{m,\mathbf{k},\rho}^{\dagger}f_{n,\mathbf{k}+\mathbf{q},\rho})e^{i(\epsilon_{m,\mathbf{k}}-\epsilon_{n,\mathbf{k}+\mathbf{q}})t}
\end{align}
Therefore
\begin{align}
\chi_{ab}^{R}(\mathbf{q},t-t')= & -i\theta(t-t')\frac{N^{2}-1}{V}\sum_{\mathbf{k}}\sum_{mn}(U_{\mathbf{k}}^{-1})_{ma}(U_{\mathbf{k}+\mathbf{q}})_{an}(U_{\mathbf{k}+\mathbf{q}}^{-1})_{nb}(U_{\mathbf{k}})_{bm} [n_{F}(\epsilon_{m,\mathbf{k}})-n_{F}(\epsilon_{n,\mathbf{k}+\mathbf{q}})]e^{i(\epsilon_{m,\mathbf{k}}-\epsilon_{n,\mathbf{k}+\mathbf{q}})(t-t')}
\end{align}
Perform Fourier transformation to frequency space, we get
\begin{align}
\chi_{ab}^{R}(\mathbf{q},\omega)
= \frac{N^2-1}{V}\sum_{\mathbf{k}}\sum_{mn}(U_{\mathbf{k}}^{-1})_{ma}(U_{\mathbf{k}+\mathbf{q}})_{an}(U_{\mathbf{k}+\mathbf{q}}^{-1})_{nb}(U_{\mathbf{k}})_{bm}\frac{n_{F}(\epsilon_{m,\mathbf{k}})-n_{F}(\epsilon_{n,\mathbf{k}+\mathbf{q}})}{\epsilon_{m,\mathbf{k}}-\epsilon_{n,\mathbf{k}+\mathbf{q}}+\omega+i\eta}
\end{align}

To get the full spin spectrum, we should keep in mind that above derivation
is in small BZ. Now we consider the spin correlation on a square lattice
(one site unit cell), and relate it to above formula.

\begin{equation}
\chi^{R}(\mathbf{r}_{i}-\mathbf{r}_{j},t-t')=-i\theta(t-t')\sum_{\mu\nu}\langle[S_{\nu}^{\mu}(\mathbf{r}_{i},t),S_{\mu}^{\nu}(\mathbf{r}_{j},t')]\rangle
\end{equation}
denote $\mathbf{r}_{i}=\mathbf{r}+\mathbf{v}_{a}$, $\mathbf{r}_{j}=\mathbf{r}'+\mathbf{v}_{b}$, we have

\begin{equation}
\chi^{R}(\mathbf{r}+\mathbf{v}_{a}-\mathbf{r}'-\mathbf{v}_{b},t-t')=-i\theta(t-t')\sum_{\mu\nu}\langle[S_{\nu}^{\mu}(\mathbf{r},a,t),S_{\mu}^{\nu}(\mathbf{r}',b,t')]\rangle\equiv\chi_{ab}^{R}(\mathbf{r}-\mathbf{r}',t-t')
\end{equation}
Perform fourier transformation, we get following relation
\begin{equation}
\chi^{R}(\mathbf{q},\omega)=\sum_{ab}\chi_{ab}^{R}(\mathbf{q},\omega)e^{-i\mathbf{q}(\mathbf{v}_{a}-\mathbf{v}_{b})},
\end{equation}
and the full spin spectrum is
\begin{equation}
S(\mathbf{q},\omega)=-2\text{Im}\chi^{R}(\mathbf{q},\omega)
\end{equation}


\end{document}
