\section{O($3$) order parameter}
\label{sec:oop}

The vison condensation in the quantum loop model on the triangular lattice can be analyzed by mapping it to the ferromagnetic transverse-field Ising model on the dual honeycomb lattice~\cite{Krishanu2015,Wang2017,Yan2021,Huh2011}, where
\begin{equation}
H^{}_{dual}=J^{}_1\sum_{\langle I J \rangle}v_{I}^{z}v_{J}^{z}+J^{}_2\sum_{\langle\langle I J \rangle\rangle}v_{I}^{z}v_{J}^{z}+J^{}_3\sum_{\langle\langle\langle I J \rangle\rangle\rangle}v_{I}^{z}v_{J}^{z}-\Gamma\sum_{I}v_{I}^{x}.
\label{eq:h}
\end{equation}
Here, $I$ and $J$ are the sites on dual lattice while $J_{1,2, 3}$ is the first-, second-, and third-neighbor Ising coupling. In the large-$\Gamma$ limit, all spins point in the $x$ direction. A flipped spin generates a vison excitation, thus $v_{I}^{z}=+1(-1)$ is the vison creation (annihilation) operator. As we increase $J$, the visons acquire a dispersion and eventually condense at its minima. The Fourier transform of the Hamiltonian in Eq.~\eqref{eq:h} can be written as
\begin{equation}
H^{}_{dual}=-\sum_{\mathbf{q}}\psi^{\phantom{\dagger}}_{\mathbf{q}}\,J(\mathbf{q})\,\psi_{\mathbf{q}}^{\dagger},
\label{eq:h2}
\end{equation}
where $\psi^{\phantom{\dagger}}_\mathbf{q}=(v_{1q}^{z},\ v_{2q}^{z})$, and $J(\mathbf{q})$ is the interaction matrix, which takes the form of
\begin{equation}      
	J(\mathbf{q})=\left(                
	\begin{array}{cc}  
		a & b  \\  
		b^* & a \\  
	\end{array}
	\right)             
	\label{jq}
\end{equation}
with $J_1=J_2=J_3=1$. $a$ and $b$ in the above matrix are defined as
\begin{align}
a &= e^{i\mathbf{q}\cdot \mathbf{s}_1}+e^{-i\mathbf{q}\cdot \mathbf{s}_1}+e^{i\mathbf{q}\cdot \mathbf{s}_2}+e^{-i\mathbf{q}\cdot \mathbf{s}_2}+e^{i\mathbf{q}\cdot (\mathbf{s}_1+\mathbf{s}_2)}+e^{i\mathbf{q}\cdot (\mathbf{s}_1+\mathbf{s}_2)}, \\
b &= 1+e^{i\mathbf{q}\cdot \mathbf{s}_1}+e^{i\mathbf{q}\cdot \mathbf{s}_2}+e^{i\mathbf{q}\cdot (\mathbf{s}_1+\mathbf{s}_2)}+e^{-i\mathbf{q}\cdot (\mathbf{s}_1+\mathbf{s}_2)}+e^{i\mathbf{q}\cdot (\mathbf{s}_1-\mathbf{s}_2)},
\end{align}
where $\mathbf{s}_1=(\sqrt{3}/2,-3/2)$ and $\mathbf{s}_2=(\sqrt{3}/2,3/2)$ are the primitive vectors as shown in Fig.~\ref{fig:fig1}(a). The minima of the energy dispersion $J(\mathbf{q})$ occur at three momenta $\mathbf{M}_{1}=(\pi/\sqrt{3},\pi/3)$, $\mathbf{M}_{2}=(\pi/\sqrt{3},-\pi/3)$, and $\mathbf{M}_{3}=(0,2\pi/3)$ in the hexagonal Brilliouin zone, as demonstrated in Fig~\ref{fig:sfig1}. The eigenvectors at these momenta are
\begin{equation}
	\label{eq:u}
	\begin{split}
		u_1=u_2=(1,\ 1)^T,\quad
		u_3=(1,\ -1)^T.
	\end{split}
\end{equation}
Therefore, the vison O($3$) order parameter $\bm{\phi}=(\phi_1, \phi_2, \phi_3)$ can be expressed as
\begin{equation}
	\label{eq:phi}
	\begin{split}
		\phi_1=\sum_\mathbf{r}(v_{1,\mathbf{r}}+v_{2,\mathbf{r}})e^{i\pi x},\quad
		\phi_2=\sum_\mathbf{r}(v_{1,\mathbf{r}}+v_{2,\mathbf{r}})e^{i\pi y},\quad
		\phi_3=\sum_\mathbf{r}(v_{1,\mathbf{r}}-v_{2,\mathbf{r}})e^{i\pi (y-x)},
	\end{split}
\end{equation}
where $x$ and $y$ are the coordinates of all plaquettes of the triangular lattice.

\begin{figure}[htp]
	\centering
	\includegraphics[width=0.6\columnwidth]{sfig1.pdf}
	\caption{(a) The dispersion of $J(\mathbf{q})$  and (b) the Brilliouin zone (BZ) of the honeycomb lattice. The minima of the energy dispersion $J(\mathbf{q})$ occur at $\mathbf{M}_{1}=(\pi/\sqrt{3},\pi/3)$, $\mathbf{M}_{2}=(\pi/\sqrt{3},-\pi/3)$, $\mathbf{M}_{3}=(0,2\pi/3)$ and the three other symmetry-equivalent $\mathbf{M}$ points. }
	\label{fig:sfig1}
\end{figure}

The O($3$) order parameter $\bm{\phi}=(\phi_1,\ \phi_2,\ \phi_3)$ can transform under various symmetries. The transformation matrices for $\bm{\phi}$ with respect to translations ($T_x, T_y$), bond inversion ($I$), rotation ($R_6$) , global $Z_2$, and mirror symmetry $M$ along the $y$ axis, are given by
\begin{alignat}{5}
	T_x&=\left(
	\begin{array}{ccc}
		-1 & 0 & 0 \\
		0 & 1 &  0\\
		0 &  0 & -1\\
	\end{array}
	\right),\quad
	&&T_y&&=\left(
	\begin{array}{ccc}
		1 & 0 & 0 \\
		0 & -1 &  0\\
		0 &  0 & -1\\
	\end{array}
	\right),\quad
	&&I&&=\left(
	\begin{array}{ccc}
		1 & 0 & 0 \\
		0 & 1 &  0\\
		0 &  0 & -1\\
	\end{array}
	\right),
	\label{I}\\
	R_6&=\left(
	\begin{array}{ccc}
		0 & 0 & -1 \\
		1 & 0 &  0\\
		0 &  1 & 0\\
	\end{array}
	\right),
	\label{r6}
\quad
	&&Z_2&&=\left(
	\begin{array}{ccc}
		-1 & 0 & 0 \\
		0 &-1 &  0\\
		0 &  0 & -1\\
	\end{array}
	\right),
\quad
	&&M&&=\left(
	\begin{array}{ccc}
		0 & 0 & -1 \\
		0 &1 &  0\\
		-1 &  0 & 0\\
	\end{array}
	\right).
\end{alignat}
These matrices generate a finite subgroup of O($3$) which contain 48 elements~\cite{Krishanu2015}. This group is isomorphic to $Z_2\times S_4$, with $S_4$ being the permutation group of four elements.

\section{Order parameter histograms}
\label{sec:histo}

Besides the order parameter histograms presented in the main text, here, we add more detailed results of the evolution from the LN to the VP solid phase. Figure~\ref{fig:sfig2} compiles the histograms for a $L=\beta=16$ system, as $V$ is varied from $-1$ to $1$. It is clear that for cases with $V<0.1$, the histograms exhibit the face-cubic distribution, and for the cases with $V>0.1$, the histograms exhibit the corner-cubic distribution.

\begin{figure}[b]
	\centering
	\includegraphics[width=1.0\columnwidth]{sfig2.pdf}
	\caption{Order parameter histograms for a $L=16$ system. (a) Projection of the O($3$) order parameter on the ($\phi_1,\phi_2$) plane and ($\phi_2,\phi_3$) plane. The histograms in the three phases (LN with $V=-0.9$, VP with $V=0.3$, and QSL with $V=0.9$) are shown. (b) Evolution of the histogram by tuning $V$, from LN, through VP, to QSL.}
	\label{fig:sfig2}
\end{figure}

Close to the LN-to-VP transition, i.e., the panel with $V=0.1$ in Fig.~\ref{fig:sfig2}(b), the histogram develops asymptotic O($3$) symmetry with very weak corner-cubic anisotropy (this is because $V=0.1$ lies on the VP side of the transition). Inside the VP phase, and as $V$ increases towards the VP--QSL continuous transition point $V_c = 0.6$, the histograms [in panels with $V=0.3, 0.5$ in Fig.~\ref{fig:sfig2}(b)] begin to shrink in diameters, suggesting the reduction of the vison order parameter. Eventually, inside the QSL phase [panels with $V=0.7,0.9$ in Fig.~\ref{fig:sfig2}(b)], the histogram becomes a point at $(0,0,0)$, indicating the complete vanishing of the vison order parameter in the topologically ordered QSL phase.


\section{Real-space vison correlations}
\label{sec:vison}

To investigate the symmetry-breaking pattern of the VP phase, the correlations of the triangle-$T$ operators are calculated. As described in the main text, this operator is defined on a triangular plaquette from three related kinetic ($t$) terms, $T_i=t_{1,i}+t_{2,i}+t_{3,i}$, [see inset of Fig.~\ref{fig:sfig3}(a)]. Here, $i$ is the triangle label, and $1,2,3$ connote the three kinds of $t$-operators (rhombi) in the Hamiltonian acting on this triangle. The correlation of the operators between the triangular plaquettes $i$ and $j$ is $\langle T_iT_j\rangle$.

We calculate the real-space vison correlations with selecting the leftmost and lowest red triangle below as a reference (its correlation is identically 1)
as illustrated in Fig.~\ref{fig:sfig3}(b).
Two visons living on the centers of related triangular plaquettes $i$ and $j$ can be connected with a open string; their correlation is $\langle v_i v_j\rangle=\langle(-1)^{N_{P_{ij}}} \rangle$ and ${N_{P_{ij}}}$ is the number of dimers cut along the path $P$ between plaquettes $i$ and $j$. In the two-dimer-per-site case, $\langle v_i v_j \rangle$ is independent on the chosen path, i.e., the visons here are gauge invariant and measurable.

We found that the vison correlations clearly exhibit the translational symmetry breaking of the VP solid, but preserve $C_3$ rotational symmetry. The triangles with positive (negative) values of vison correlations gather into big triangles. In addition, the central triangle of the thus-formed big triangle always has a larger absolute value, which is about three times that of the other triangles in the same big triangle. The translation periods along the $\mathbf{s}_{1}$ and $\mathbf{s}_{2}$ directions are both $2$, which means that the unit cell is a $2\times 2$ rhomboid.

\begin{figure}[htp]
	\centering
	\includegraphics[width=1\columnwidth]{sfig3.pdf}
	\caption{(a) Real-space $T$-operator correlations of the $L=8$ system at $V=0.3$. The inset shows the summation of the three $t$ terms in the Hamiltonian.  (b) Real-space vison correlations of the $L=8$ system at $V=0.3$; $\mathbf{s}_{1}$ and $\mathbf{s}_{2}$ are the primitive vectors. For both figures, we set the  leftmost and lowest triangle as the reference. The red (grey) color describes the positive (negative) value of the correlation.}
	\label{fig:sfig3}
\end{figure}

\section{Sweeping cluster quantum Monte Carlo with two dimers per site}
\label{sec:qmc}


The sweeping cluster QMC approach employed in this work is a new method developed by us, which works well in constrained quantum lattice models~\cite{Yan2019,ZY2020improved,Yan2021,ZYan2022}. 
Prior to sweeping cluster QMC, to solve the QDM or QLM types of constrained models, one had to rely on either exact diagonalization of small systems, or variational approaches such as DMRG that suffer from finite-size effects on the cylindrical geometry~\cite{Krishanu2015}, or the projector Monte Carlo approaches, which include the Green's function~\cite{Ivanov2004,Ralko2005,Ralko2006,Vernay2006,Ralko2007,Plat2015} and diffusion Monte Carlo schemes~\cite{OFS2005,OFS2006}. These projector Monte Carlo methods obey the geometric constraints but are not efficient away from the RK point~\cite{OFS2005walk} and only work at $T=0$. Furthermore, there does not exist any cluster update for the projector methods. On the contrary, the sweeping cluster algorithm is based on the world-line Monte Carlo scheme~\cite{OFS2002,Alet2005a,Alet2005b} to sweep and update layer by layer along the imaginary time direction so that the local constraints (gauge field) are recorded by update lines. In this way, all the samplings are done in the restricted Hilbert space and it contains the cluster update scheme for constrained systems~\cite{ZY2020improved} and works at all temperatures. Proper finite size scaling analysis can then be carried out to explore phase transitions and critical phenomena. 


To implement the constraint of having two dimer per site in our QMC simulation, we set the initial state as one of the three LN patters with the same probability~\cite{Plat2015}. In the region of the LN phase, the equilibrium process will not change the orientation of the initial state; therefore, the randomly selected initial state ensures each pattern of the LN phase has an equal probability of $\frac{1}{3}$  to emerge. In the VP and QSL regions, one can also use such initialization schemes, and we have tested that in these regions too, they produce the same QMC results as a completely random initialization. We simulate $L\times L$ triangular lattices with system sizes $L=8,12,16$ with the inverse temperature set to $\beta=L$ and $10^4$ Monte Carlo samplings used to obtain average values of the observables in all calculations.

In this article, we can measure the information about single visons because in a strictly constrained space, the energy gap of other quasiparticles such as spinons, becomes infinitely large and thus these quasiparticles do not exist in the restricted Hilbert space. Therefore, the vison is well-defined here.


\section{Anisotropy parameters}
\label{sec:e}
The effective action of the O($3$) transition is given by,
\begin{align}
S=&\int dt dx^2 \sum_{i=1}^3(\partial^{}_{\mu} \phi^{}_i )^2+ r \sum_{i=1}^3 \phi_i^2 + \mu (\sum_{i=1}^3 \phi^{}_i\phi^{}_i)^2 +\nu^{}_4 \sum_{i=1}^3 (\phi^{}_i)^4+\mu^{}_6 (\sum_{i=1}^3 \phi_i^2)^3 + \nu^{}_6 (\phi^{}_1\phi^{}_2\phi^{}_3)^2+\nu_6' (\sum_{i=1}^3 \phi_i^2)(\sum_{i=1}^3 \phi_i^4).
\label{eq:seq8}
\end{align}
Assuming that the coupling constants $\nu_4$, $\nu_6$, $\nu_6' $ are all much smaller than one, the average of the spherical harmonics are
\begin{align}
\langle Y^{0}_{4}\rangle &= -\frac{1}{15\sqrt{\pi}}\nu^{}_4\langle\phi^4\rangle_0-\frac{1}{15\sqrt{\pi}}\nu_6'\langle\phi^6\rangle_0-\frac{1}{330\sqrt{\pi}}\nu^{}_6\langle\phi^6\rangle_0,\nonumber\\
\langle \phi^2Y^{0}_{4}\rangle &= -\frac{1}{15\sqrt{\pi}}\nu^{}_4\langle\phi^6\rangle_0-\frac{1}{15\sqrt{\pi}}\nu_6'\langle\phi^8\rangle_0-\frac{1}{330\sqrt{\pi}}\nu^{}_6\langle\phi^8\rangle_0,\nonumber\\
\langle Y^{0}_{6}\rangle &= -\frac{1}{231\sqrt{13\pi}}\nu^{}_6\langle\phi^6\rangle_0.
\label{eq:seq9}
\end{align}
Here, $\langle .\rangle_{0}$ denotes the averages under Boltzmann weights determined by $S_0$, where $S_0$ is the action~\eqref{eq:seq8} with the anisotropic couplings $\nu_4$, $\nu_6$ and $\nu_6'$ turned off. The symbol $\phi$ denote the length of $\phi_i$, i.e., $\phi=(\mu (\sum_{i=1}^3 \phi_i\phi_i)^2)^{1/2}$. The two spherical harmonics are
\begin{align}
Y^{0}_{4} &= \frac{3}{16\sqrt(\pi)}(3-30\cos^{2}{\theta}+35\cos^{4}{\theta}), \\
Y^{0}_{6} &= \frac{\sqrt{13}}{32\sqrt{\pi}}(-5+105\cos^{2}{\theta}-315\cos^{4}{\theta}+231\cos^{6}{\theta}).
\end{align}
As an example, we explain here how to derive the formula for $\langle Y^{0}_{4}\rangle$:
\begin{align}
\langle Y^{0}_{4}\rangle=&\nu^{}_4 \left\langle Y^{0}_{4} \sum_{i=1}^3 (\phi_i)^4\right\rangle_0+ \nu^{}_6 \left\langle Y^{0}_{4}  (\phi_1\phi_2\phi_3)^2\right\rangle_0+\nu_6' \left \langle Y^{0}_{4} (\sum_{i=1}^3 \phi_i^2)(\sum_{i=1}^3 \phi_i^4)\right\rangle_0\nonumber\\
=&\nu^{}_4 \left\langle \phi^4 \right\rangle_0 \left\{Y^{0}_{4} f(\theta,\psi)\right\}_{S^2}+\nu^{}_6 \left\langle \phi^6 \right\rangle_0 \left\{Y^{0}_{4}  g(\theta,\psi) \right\}_{S^2}
+\nu_6' \left\langle \phi^6 \right\rangle_0 \left\{Y^{0}_{4}  f(\theta,\psi)\right\}_{S^2}.\nonumber\\
=& -\frac{1}{15\sqrt{\pi}}\nu^{}_4\langle\phi^4\rangle_0-\frac{1}{15\sqrt{\pi}}\nu_6'\langle\phi^6\rangle_0-\frac{1}{330\sqrt{\pi}}\nu^{}_6\langle\phi^6\rangle_0,
\end{align}
with
\begin{align}
f(\theta,\psi)=\cos(\theta)^4+ \sin(\theta)^4 \cos (\psi)^4+  \sin(\theta)^4 \sin (\psi)^4, \quad g(\theta,\psi)=\cos(\theta)^2 \sin(\theta)^2 \cos (\psi)^2 \sin(\theta)^2 \sin (\psi)^2.
\end{align}
The bracket $\{ A \}_{S^2}=\frac{1}{4\pi}\int A \sin(\theta) d\theta d\psi $ is the average over the unit two-sphere. In the above derivation, we have used the parametrization $\phi_1$\,$=$\,$\phi \cos(\theta)$, $\phi_2$\,$=$\,$\phi \sin(\theta)\cos(\psi)$ and $\phi_3$\,$=$\,$\phi \sin(\theta)\sin(\psi)$.  We have also used the approximation $\langle\phi^n f(\theta,\psi)\rangle_0=\langle\phi^n\rangle_0 \{f(\theta,\psi)\}_{S_2}$. This crude approximation, strictly speaking, is valid in the semiclassical (approximate mean-field theory) region. The results in Fig.~\ref{fig:fig2} of the main text, however, show that around the first-order phase transition at $V/t=0.05(5)$, these formulae give reasonable qualitative results for the couplings $\nu_4,\nu_6$ and $\nu_6'$. In particular, $\nu_4$ changes sign at the phase transition as expected.

Through some simple algebra, the three anisotropy parameters can be expressed as
\begin{align}
\nu^{}_{4} &= -\frac{15\sqrt{\pi}}{\langle \phi^{4}\rangle}\left(\langle Y^{0}_{4}\rangle+\frac{\langle \phi^2 Y^{0}_{4}\rangle\langle \phi^4\rangle\langle \phi^6\rangle-\langle Y^{0}_{4}\rangle\langle\phi^6\rangle^2}{\langle \phi^6\rangle^2-\langle \phi^4\rangle\langle \phi^8\rangle}\right), \\
\nu^{}_{6} &= -\frac{231\sqrt{13\pi}}{\langle \phi^{6}\rangle}\langle Y^{0}_{6}\rangle,\\
\nu_6' &= 15\sqrt{\pi}\left(\frac{\langle \phi^2 Y^{0}_{4}\rangle\langle \phi^4\rangle-\langle Y^{0}_{4}\rangle\langle\phi^6\rangle}{\langle\phi^6\rangle^2-\langle \phi^4\rangle\langle \phi^8\rangle}+\frac{7\sqrt{13}}{10\langle\phi^6\rangle}\langle Y^{0}_{6}\rangle\right).
\end{align}
We have also used the approximation $\langle \phi^n \rangle_0=\langle \phi^n \rangle$, which is valid since their difference is of linear order in the perturbation parameters $\nu_4,\nu_6$, or $\nu_6'$. These differences only introduce higher-order corrections in Eq.~\eqref{eq:seq9}.

\section{Exact diagonalization results}
\label{sec:f}
In this section, we apply the Lanczos exact diagonalization (ED) method to a $4 \times 4$ ($N_{bond}=48$) lattice within the VP phase. The number of configurations of the constrained Hilbert space is $586,695$. From the ground state wavefunction, we analyze the wavefunction symmetry of 8 points of cubic order parameters without statistical error. As shown in Fig.~\ref{fig:sfig2}, the 8 points stay in 8 octants, so each configuration of the wavefunction can be classified into a certain octant, except the ones on the boundaries between different octants. Then---akin to what is done in QMC---we classify configurations according to related octants and average the ones in the same class to obtain the dimer density of a certain cubic order parameter point. It is worth noting that the two classes of opposite octants have opposite sign for vison configurations but their dimer configurations are the same. Thus, we have to average the dimer configurations of the two opposite octants in fact. The result [Fig.\ref{fig:sfig4}(a)] is similar to that obtained with QMC: there is still no any dimer order to within the numerical precision of ED. We found that for each link, the real-space dimer density is $1/3$.

Similarly, we average the ED vison configurations in each octant. The histogram of the order parameter and the real-space vison density obtained via ED are consistent with our QMC results. The real-space vison density is shown in Fig.~\ref{fig:sfig4}(b). All these results strongly support the QMC conclusions that there is indeed a hidden vison order phase between the spin liquid and nematic phases.

\begin{figure}[htp]
	\centering
	\includegraphics[width=0.8\columnwidth]{sfig4.pdf}
	\caption{(a) Exact-diagonalization result for the real-space dimer density for the $L=4$ system at $V=0.3$; as previously, $\mathbf{s}_{1}$ and $\mathbf{s}_{2}$ are the primitive vectors. All bonds hold $1/3$ dimer density. (b) The real-space vison density, where we have set $v_1(0,0)=1$. Similar to the QMC results, four small triangles with the same sign of the vison density form bigger triangles, and the small triangles at the centers of the big ones have a larger absolute value of the density.}
	\label{fig:sfig4}
\end{figure}

\end{widetext}

\end{document}
