\section{Introduction}\label{sec1}

The acceleration of the universe expansion is one of the most active discoveries of modern cosmology. Observational studies include Type Ia supernovae \cite{Riess,Perlmutter}, large scale structure \cite{T.Koivisto,S.F.}, baryon acoustic oscillations \cite{D.J.,W.J.} and cosmic microwave background radiation \cite{R.R.,Z.Y.}. The reason behind the late-time acceleration is a mystery. Several models hypothesize the existence of a component called {\it dark energy}, which makes up around 70\% of the entire universe and could possess the feature of speeding up the universe expansion.

The cosmological constant $\Lambda$ in General Relativity (GR) field equations plays the role of dark energy, i.e., a fluid with constant energy density and high negative pressure. There are some issues with the cosmological constant model, like the cosmic coincidence problem \cite{dalal/2001}, which is the fact that the density of non-relativistic matter and dark energy are the same order today. A more delicate issue surrounding the cosmological constant is the so-called {\it cosmological constant problem}, which is the high discrepancy between the astronomically observed value of $\Lambda$ \cite{Riess,Perlmutter} and the Particle Physics theoretically predicted value of the quantum vacuum energy \cite{weinberg/1989}.

With the main purpose of solving the above cosmological issues, dynamical (time-varying) dark energy models such as Chaplygin gas model \cite{M.C.,A.Y.}, k-essence \cite{T.Chiba,C.Arm.}, quintessence \cite{Carroll,Y.Fujii} and decaying vacuum models \cite{xu/2011,tong/2011,freese/1987,abdel-rahman/1992} have been proposed for some time in the literature.

Modified theories of gravity have also been intensively investigated to understand the origin of the cosmic acceleration as well as to address the cosmological constant model problems. It is possible to predict late-time cosmic acceleration by modifying GR action. Some possibilities can be seen within $f(R)$ \cite{appleby/2007,amendola/2007,saffari/2008}, $f(G)$ \cite{cognola/2006,li/2007}, $f(R,\mathcal{T})$ \cite{moraes/2016} and $f(T)$ \cite{ren/2021,nashed/2015,setare/2013} theories of gravitation, with $R$, $G$, $\mathcal{T}$ and $T$ being respectively the Ricci, Gauss-Bonnet, energy-momentum and torsion scalars.

In the present article we will work with the recently introduced $f(Q)$ theory of gravity \cite{J.B.}, for which $Q$ is the non-metricity scalar, to be presented below. The non-metricity formulation has been discussed earlier by Hehl and Ne'eman (see References \cite{H1,H2,H3,H4}). The symmetric teleparallel gravity is proven in the so-called coincidence gauge by imposing that the connection is symmetric \cite{nester/1999}. Weyl geometry is also observed to be a particular example of Weyl Cartan geometry in which torsion disappears. The non-metricity is interpreted as a massless spin 3-field in the case of symmetric connections \cite{N1,N2}. Also, it is noted in the literature that due to the appearance of non-metricity, the light cone structure is not preserved during parallel transport \cite{H5}. Further, fermions are an issue in TEGR because they couple to the axial contorsion of the Weitzenbock connection. This difficulty is eliminated in STEGR since Dirac fermions only couple to the completely antisymmetric component of the affine connection and is unaffected by any disformation piece.

Although recently proposed, the $f(Q)$ gravity theory already presents some interesting and valuable applications in the literature. The first cosmological solutions in $f(Q)$ gravity appear in References \cite{jimenez/2020,khyllep/2021}, while $f(Q)$ cosmography and energy conditions can respectively be seen in \cite{mandal/2020,mandal/2020b}.

Here we are going to consider $f(Q)$ cosmology in the presence of a viscous fluid. When a cosmic fluid expands too fast, the recovering of thermodynamic equilibrium generates an effective pressure. The high viscosity in a cosmic fluid is the manifestation of such an effective pressure \cite{J.R.,H.O.}.

Basically, there are two viscosity coefficients, namely shear viscosity and bulk viscosity. Shear viscosity is related to velocity gradients in the fluid, and by considering the universe as described by homogeneous and isotropic Friedmann-Lem\^aitre-Robertson-Walker (FLRW) background, it can be omitted. Anyhow, by dropping the FLRW background assumption, several cosmological models with shear viscosity fluid have been constructed, as one can check, for instance, \cite{bali/1988,bali/1987,deng/1991,huang/1990}. On the other hand, bulk viscosity, which we are going to consider here, introduces damping associated with volumetric straining. To get in touch with bulk viscous fluid cosmological models, one can check References \cite{samanta/2017,satish/2016,beesham/1993,colistete/2007}.

Researchers examine dark energy reconstruction with numerous observations as data increases. The majority of studies has been concentrated on observable evidences from Type Ia supernovae, cosmic microwave background and baryon acoustic oscillations (BAO), which are known to be helpful in constraining cosmological models. The Hubble parameter dataset shows the intricate structure of the expansion of the universe. The ages of the most massive and slowly evolving galaxies offer direct measurements of the Hubble parameter $H(z)$ at various redshifts $z$, resulting in the development of a new form of standard cosmological probe \cite{Jimenez/2002}.

In our present work we include 31 measurements of Hubble expansion spanned using differential age method \cite{Sharov/2017} and BAO data consisting of six points \cite{Blake/2011}. Scolnic et al. recently published a large Type Ia supernovae sample named Pantheon, with 1048 points and covering the redshift range $0.01< z< 2.3$ \cite{Scolnic/2018}. Our analysis uses the $H(z)$, BAO and Pantheon samples to constrain the cosmological model.

This work aims to describe the recently observed late-time acceleration with the help of bulk viscosity of cosmic fluid (without including any dark energy component) in the framework of  $f(Q)$ theory of gravity. The manuscript is organized as follows: in Sec. \ref{sec2} we discuss the $f(Q)$ gravity formalism. In Sec. \ref{sec3} we describe the FLRW universe dominated by bulk viscous non-relativistic matter and derive the expression for the Hubble parameter and the deceleration parameter. Further, in Sec. \ref{sec4}, we analyze the observational data to find the best-fit ranges for the parameters using the Hubble data set containing 31 points, BAO sample and the Pantheon data set of 1048 samples. Moreover, we analyze the behavior of different cosmological parameters such as density, effective pressure, deceleration parameter and effective equation of state (EoS) parameter. In Sec. \ref{sec5}, we analyze the behavior of statefinder parameters on the values constrained by the observational data to differentiate between dark energy models. Lastly, we discuss our results in Sec. \ref{sec6}.

\section{Fundamental formulations in $f(Q)$ gravity}\label{sec2}

In $f(Q)$ gravity theory, the spacetime is established with the help of non-metricity and symmetric teleparallelism condition, i.e. $\nabla_\alpha g_{\mu\nu} \neq 0$ and $R^\rho_{\sigma\mu\nu}=0$. The associated affine connection is given by

\begin{equation}\label{fq1}
	\Upsilon^\alpha_{\ \mu\nu}=\Gamma^\alpha_{\ \mu\nu}+L^\alpha_{\ \mu\nu},
\end{equation}
with

\begin{eqnarray}
	\Gamma^\alpha_{\ \mu\nu}\equiv\frac{1}{2}g^{\alpha\lambda}(g_{\mu\lambda,\nu}+g_{\lambda\nu,\mu}-g_{\mu\nu,\lambda}),\label{fq2}\\
	L^\alpha_{\ \mu\nu}\equiv\frac{1}{2}(Q^{\alpha}_{\ \mu\nu}-Q_{\mu \ \nu}^{\ \alpha}-Q_{\nu \ \mu}^{\ \alpha}),\label{fq3}
\end{eqnarray}
being the Christoffel symbols and the distortion tensor, respectively, where

\begin{eqnarray}
	Q_{\alpha\mu\nu}\equiv\nabla_\alpha g_{\mu\nu},\label{fq4}
\end{eqnarray}
is the non-metricity tensors.

The non-metricity tensor given by Eq.\eqref{fq4} has following two traces,

\begin{equation}\label{2a}
Q_\alpha = Q_\alpha\:^\mu\:_\mu,
\end{equation}

\begin{equation}\label{2b}
	\tilde{Q}_\alpha = Q^\mu\:_{\alpha\mu}.
\end{equation}

In addition, the superpotential tensor or non-metricity conjugate is given by

\begin{eqnarray}\label{2c}
4P^\lambda\:_{\mu\nu} &=& -Q^\lambda_{\ \mu\nu}+Q_{\mu\ \nu}^{\ \lambda}+Q_{\nu\ \mu}^{\ \lambda}+(Q^\lambda-\tilde{Q}^\lambda)g_{\mu\nu}\nonumber\\
&-&\frac{1}{2}(\delta^\lambda_{\ \mu}Q_\nu+\delta^\lambda_{\ \nu}Q_\mu).
\end{eqnarray}

Then the trace of the non-metricity tensor can be acquired as \cite{Ruth}

\begin{equation}\label{2d}
Q = -Q_{\lambda\mu\nu}P^{\lambda\mu\nu}.
\end{equation}

The definition of the energy-momentum tensor for matter is
\begin{equation}\label{2e}
\mathcal{T}_{\mu\nu} = \frac{-2}{\sqrt{-g}} \frac{\delta(\sqrt{-g}L_m)}{\delta g^{\mu\nu}}.
\end{equation}

Furthermore, one can obtain the following relation between the curvature tensor $R^\rho_{\sigma\mu\nu}$ and $\mathring{R}^\rho_{\sigma\mu\nu}$ corresponding to the connection $\Upsilon$  and $\Gamma$ as

\begin{equation}\label{2f}
R^\rho_{\sigma\mu\nu} = \mathring{R}^\rho_{\sigma\mu\nu} + \mathring{\nabla}_\mu L^\rho_{\nu\sigma} - \mathring{\nabla}_\nu L^\rho_{\mu\sigma} + L^\rho_{\mu\lambda} L^\lambda_{\nu\sigma} - L^\rho_{\nu\lambda} L^\lambda_{\mu\sigma}
\end{equation}

and so

\begin{equation}\label{2g}
R_{\sigma\nu} = \mathring{R}_{\sigma\nu} + \frac{1}{2} \mathring{\nabla}_\nu Q_\sigma + \mathring{\nabla}_\rho L^\rho_{\nu\sigma} - \frac{1}{2} Q_\lambda L^\lambda_{\nu\sigma} - L^\rho_{\sigma\lambda} L^\lambda_{\rho\sigma}
\end{equation}

\begin{equation}\label{2h}
R = \mathring{R} + \mathring{\nabla}_\lambda Q^\lambda - \mathring{\nabla}_\lambda \tilde{Q}^\lambda - \frac{1}{4} Q_\lambda Q^\lambda + \frac{1}{2} Q_\lambda \tilde{Q}^\lambda - L_{\rho\nu\lambda} L^{\lambda\rho\nu}
\end{equation}

Now, the connection \eqref{fq1} can be parameterized as \cite{J.B.}

\begin{equation}
\Upsilon^\alpha \,_{\mu \beta} = \frac{\partial x^\alpha}{\partial \xi^\rho} \partial_ \mu \partial_\beta \xi^\rho.
\end{equation}
Here, $\xi^\alpha = \xi^\alpha (x^\mu)$ is an invertible relation. Hence, it is always possible to find a coordinate system so that the connection $ \Upsilon^\alpha_{\ \mu\nu} $ vanishes. This situation is called coincident gauge and the covariant derivative $\nabla_{\alpha}$ reduces to the partial one $\partial_{\alpha} $. But in any other coordinate system in which this affine connection does not vanish, the metric evolution will be affected and result in a completely different theory \cite{AP,Avv}. Thus in the coincident gauge coordinate , we have
\begin{equation}\label{2i}
Q_{\alpha\mu\nu} = \partial_\alpha g_{\mu\nu}
\end{equation}
while in an arbitrary coordinate system,
\begin{equation}\label{2j}
Q_{\alpha\mu\nu}= \partial_\alpha g_{\mu\nu} - 2 \Upsilon^\lambda_{\alpha (\mu} g_{\nu)\lambda}.
\end{equation}

Locally, GR does not distinguish between gravitational and inertial effects, however, by invoking frame fields, it is possible to covariantly define gravitational energy in the teleparallel approach \cite{maluf/2013}. The canonical frame is identified in the absence of both curvature and torsion and the canonical coordinates in the absence of inertial effects. New physics can emerge from such a formalism.

The $f(Q)$ gravity is described by the action \cite{J.B.}:
\begin{equation}\label{2k}
S= \int{\frac{1}{2}f(Q)\sqrt{-g}d^4x} + \int{L_m\sqrt{-g}d^4x},
\end{equation}
in which $f(Q)$ is an arbitrary function of the non-metricity, $g=\text{det}\  g_{\mu\nu}$ and $L_m$ is the lagrangian density of matter. One can check in Reference \cite{J.B.} that the  functional form $f(Q)=-Q $ corresponds to the STEGR (symmetric teleparallel equivalent to General Relativity) limit.

The gravitational field equations obtained by varying action \eqref{2k} with respect to the metric are
\begin{widetext}
\begin{equation}\label{2l}
\frac{2}{\sqrt{-g}}\nabla_\lambda (\sqrt{-g}f_Q P^\lambda\:_{\mu\nu}) + \frac{1}{2}g_{\mu\nu}f+f_Q(P_{\mu\lambda\beta}Q_\nu\:^{\lambda\beta} - 2Q_{\lambda\beta\mu}P^{\lambda\beta}\:_\nu) = -\mathcal{T}_{\mu\nu},
\end{equation}
\end{widetext}
in which we defined $ f_Q =df/dQ $.

Moreover, by varying the action with respect to the connection, one obtains

\begin{equation}\label{2m}
\nabla_\mu \nabla_\nu (\sqrt{-g}f_Q P^{\mu\nu}\:_\lambda) =  0.
\end{equation}

\section{The cosmological model}\label{sec3}

We consider the FLRW metric for our analysis \cite{ryden/2003}, such that:
\begin{equation}\label{3a}
ds^2= -dt^2 + a^2(t)(dx^2+dy^2+dz^2).
\end{equation}
In Equation (\ref{3a}) above, $ a(t) $ is the cosmic scale factor, as usually.

From now, we will fix the coincident gauge so that connection becomes trivial and the metric is only a fundamental variable.

We are going to assume a bulk viscous fluid and below we present some considerations favoring such an assumption.

Firstly, to consider bulk viscosity in a fluid can be seen as an attempt to refine its  description, minimizing its ideal properties. This can be checked, for instance, in the stellar astrophysics realistic models in References  \cite{gusakov/2007,gusakov/2008,haensel/2002}.

Under conditions of spatial homogeneity and isotropy (which refers to the cosmological principle, as one can check, for instance, Reference \cite{ryden/2003}), the bulk viscous pressure is the unique admissible dissipative phenomenon. In a gas dynamical model, the existence of an effective bulk pressure can be traced back to a non-standard self-interacting force on the particles of the gas \cite{colistete/2007}. The bulk viscosity contributes negatively to the total pressure, as one can chegk, for instance References \cite{odintsov/2020,fabris/2006,meng/2009}.

Due to spatial isotropy, the bulk viscous pressure is the same in all spatial directions and hence proportional to the volume expansion $\theta=3H$, with $H=\dot{a}/a$ being the Hubble parameter and a dot represents time derivative.

The effective pressure of the cosmic fluid becomes \cite{brevik/2005,gron/1990,C.E./1940}
\begin{equation}
\bar{p}= p-\zeta\theta=p-3\zeta H,
\end{equation}
in which $p$ is the usual pressure and $ \zeta >0 $ is the bulk viscosity coefficient, which we will assume as a free parameter of the model.

The corresponding energy-momentum tensor is given by
\begin{equation}\label{3c}
 \mathcal{T}_{\mu\nu} = (\rho+\bar{p})u_\mu u_\nu+ \bar{p}g_{\mu\nu},
\end{equation}
in which $\rho$ is the matter-energy density and the four-velocity $ u^\mu $ is such that its components are $ u^\mu = (1,0,0,0)$.

The relation between normal pressure and matter-energy density follows \cite{J/2006} $p=(\gamma-1)\rho$, with $\gamma$ being a constant lying in the range $0\leq \gamma \leq 2$. Then the effective equation of state for the bulk viscous fluid is given by the following:
\begin{equation}\label{3i}
\bar{p}= (\gamma -1)\rho - 3\zeta H.
\end{equation}

The Friedmann-like equations for our $f(Q)$ gravitational model are obtained from the substitution of Equations (\ref{3a})-(\ref{3i}) into Equation (\ref{2g}) and read  as follows (check, for instance, References \cite{Ruth,T/2018})
\begin{equation}\label{3d}
3H^2= \frac{1}{2f_Q}\left(-\rho + \frac{f}{2}\right),
\end{equation}
\begin{equation}\label{3e}
\dot{H}+\left(3H+\frac{\dot{f_Q}}{f_Q}\right)H = \frac{1}{2f_Q} \left(\bar{p}+\frac{f}{2}\right).
\end{equation}
In particular, for $f(Q)=-Q$ we retrieve the usual Friedmann equations \cite{T/2018}, as expected, since as we have mentioned above, this particular choice for the functional form of the function $f(Q)$ is the STEGR limit of the theory.

For our investigation of bulk viscosity fluid cosmological model, we consider the following $f(Q)$ functional form:
\begin{equation}\label{3f}
f(Q)= \alpha Q^n,
\end{equation}
with $\alpha\neq0 $ and constant $n$. This particular functional form for $f(Q)$ was motivated by a polinomial form applied, for instance, in Reference \cite{mandal/2020b}.

For the above choice of the $f(Q)$ function (Equation (\ref{3f})), we rewrite Equations (\ref{3d})-(\ref{3e}) as follows

\begin{equation}\label{3g}
\rho = \alpha6^n \left(\frac{1}{2}-n \right) H^{2n},
\end{equation}
\begin{equation}\label{3h}
\dot{H}+\frac{3}{2n} H^2= \frac{6^{1-n}\bar{p}}{2\alpha n(2n-1)} H^{2(1-n)},
\end{equation}
in which for the former we have isolated $\rho$.

From Equation \eqref{3h} and Equation \eqref{3i}, we have the following:

\begin{equation}\label{3j}
\dot{H}+ \frac{3\gamma}{2n} H^2 = - \frac{6^{2-n}\zeta}{4\alpha n(2n-1)} H^{3-2n}.
\end{equation}

Now, we replace the term $d/dt$ by $d/dln a$ via the expression $d/dt=Hd/dln a$, such that Equation \eqref{3j} becomes
\begin{equation}\label{3k}
\frac{dH}{dln a} + \frac{3\gamma}{2n}H = -\frac{6^{2-n}\zeta}{4\alpha n(2n-1) } H^{2(1-n)}.
\end{equation}

The integration of Equation \eqref{3k} yields the following solution

\begin{equation}\label{3l}
H(a) = \left\{( H_0 a^{-\frac{3\gamma}{2n}} )^{2n-1} + \frac{6^{1-n}\zeta}{\gamma \alpha (2n-1)^2}[ a^{-\frac{3\gamma(2n-1)}{2n}}-1 ] \right\}^{\frac{1}{2n-1}},
\end{equation}
 with $H_0$ being a constant of integration to be found below.

We obtain the Hubble parameter in terms of redshift by using relation \cite{ryden/2003} $ a(t)= 1/(1+z)$ as
\begin{widetext}
\begin{equation}\label{3m}
H(z) = \left\{[ H_0 (1+z)^{\frac{3\gamma}{2n}} ]^{2n-1} + \frac{6^{1-n}\zeta}{\gamma \alpha (2n-1)^2} [ (1+z)^{\frac{3\gamma(2n-1)}{2n}}-1 ] \right\}^{\frac{1}{2n-1}}.
\end{equation}
\end{widetext}
 By making $z=0$ in (\ref{3m}) we find that $H(0)=H_0$ so that $H_0$ is the Hubble parameter present value.

The deceleration parameter is defined as $q= -\ddot{a}a/\dot{a}^2 = -\ddot{a}/(H^2a) $. From Equation \eqref{3j} we then have

\begin{widetext}
\begin{equation}\label{3l}
q(z)=\frac{3}{2n}\left\{\frac{\zeta}{\alpha 6^{n-1}(2n-1) \left\{ [ H_0(1+z)^{\frac{3\gamma}{2n}} ]^{2n-1} + \frac{6^{1-n}\zeta}{\gamma\alpha (2n-1)^2} [ (1+z)^{\frac{3\gamma(2n-1)}{2n}} -1 ] \right\}} + \gamma\right\} -1.
\end{equation}
\end{widetext}

\section{Observational constraints}\label{sec4}

To examine the observational features of our cosmological model, we use the most recent cosmic Hubble and Supernovae observations. We use 31 points of the Hubble data sets, 6 points of the BAO data sets and 1048 points from the Pantheon supernovae samples. We apply the Bayesian analysis and likelihood function along with the Markov Chain Monte Carlo (MCMC) method in \texttt{emcee} python library \cite{Mackey/2013}.

\subsection{Hubble datasets}

The Hubble parameter can be expressed as $H(z)=-dz/[dt(1+z)]$. As $dz$ is derived from a spectroscopic survey, the model-independent value of the Hubble parameter may be calculated by measuring the quantity $dt$.

{We incorporate the set of 31 data points that are measured from the differential age approach \cite{Sharov/2018} to avoid extra correlation with BAO data. The mean values of the model parameters $\zeta$, $\alpha$, $\gamma$ and $n$ are calculated using the chi-square function  as follows:

\begin{equation}
\chi _{H}^{2}(\zeta,\alpha,\gamma,n)=\sum\limits_{k=1}^{31}
\frac{[H_{th}(\zeta,\alpha,\gamma,n, z_{k})-H_{obs}(z_{k})]^{2}}{
\sigma _{H(z_{k})}^{2}}.  \label{4a}
\end{equation}

Here, $H_{th}$ is the Hubble parameter value predicted by the model, $H_{obs}$ represents its observed value and the standard error in the observed value of $H$ is $\sigma _{H(z_{k})}$.

From the Hubble dataset, we obtain the best fit values for $\zeta$, $\alpha$, $\gamma$, $n$ as the $1-\sigma$ and $2-\sigma$ contour plots in Fig.\ref{ContourHub}.

\begin{widetext}

\begin{figure}[H]
\centering
\includegraphics[scale=0.67]{HUBC}
\caption{The $1-\sigma$ and $2-\sigma$ likelihood contours for the model parameters using the Hubble datasets.}
\label{ContourHub}
\end{figure}

\begin{figure}[H]
\includegraphics[scale=0.52]{ErrorHubble}
\caption{The error bar plot of $H$ versus $z$ for the considered $f(Q)$ model. The solid red line is the curve for $f(Q)$ model whereas the black dotted line represents the $\Lambda$CDM model. The blue dots depict the 31 points of the Hubble data.} \label{ErrorHubble}
\end{figure}

\end{widetext}

The best fit values of the model parameters are $\zeta= 0.65^{+0.10}_{-0.10} $, $\alpha=-0.0166^{+0.0079}_{-0.0053}$, $\gamma=1.36^{+0.14}_{-0.12}$ and $n=0.974^{+0.019}_{-0.027}$.

Fig.\ref{ErrorHubble} shows the error bar plot of the considered model and $\Lambda$CDM or  standard cosmological model, with  cosmological constant density parameter $\Omega_{\Lambda_{0}}=0.7$, matter density parameter $\Omega_{m_{0}}=0.3$ and $H_{0}= 69$ km/s/Mpc.

\subsection{BAO datasets}
The BAO distance dataset, which includes the 6dFGS, SDSS and WiggleZ surveys, comprise BAO measurements at six different redshifts in Table 1. The characteristic scale of BAO is ruled by the sound horizon $r_{s}$ at the epoch of photon decoupling $z_{\ast }$ that is given by the following relation:

\begin{equation}\label{4b}
r_{s}(z_{\ast })=\frac{c}{\sqrt{3}}\int_{0}^{\frac{1}{1+z_{\ast }}}\frac{da}{
a^{2}H(a)\sqrt{1+(3\Omega _{b0}/4\Omega _{\gamma 0})a}}.
\end{equation}
Here, $\Omega _{b0}$ and $\Omega _{\gamma 0}$ correspond to the present densities of baryons and photons respectively.

The following relations are used in BAO measurements

\begin{equation}\label{4c}
 \triangle \theta =\frac{r_{s}}{d_{A}(z)},
\end{equation}

\begin{eqnarray}\label{4d}
d_{A}(z)=\int_{0}^{z}\frac{dz^{\prime }}{H(z^{\prime })},
\end{eqnarray}

\begin{equation}
	\triangle z=H(z)r_{s},
\end{equation}
where $\triangle \theta $ represents the measured angular separation, $d_{A}$ is the angular diameter distance and $\triangle z$ represents the measured redshift separation of the BAO feature in the 2 point correlation function of the galaxy distribution on the sky along the line of sight.

In this work, BAO datasets of six points for $d_{A}(z_{\ast })/D_{V}(z_{BAO})$ is taken from the References \cite{BAO1, BAO2, BAO3, BAO4, BAO5, BAO6}, where the redshift at the epoch of photon decoupling is taken as $z_{\ast }\approx 1091$ and  $d_{A}(z)$ is the co-moving angular diameter
distance together with the dilation scale $D_{V}(z)=\left[
d_{A}(z)^{2}z/H(z)\right] ^{1/3}$. The chi-square function for the BAO datasets is taken to be \cite{BAO6}

\begin{equation}\label{4e}
\chi _{BAO}^{2}=X^{T}C^{-1}X\,,
\end{equation}

\begin{widetext}

\begin{table}[H]

\begin{center}

\begin{tabular}{|c|c|c|c|c|c|c|}
\hline
$z_{BAO}$ & $0.106$ & $0.2$ & $0.35$ & $0.44$ & $0.6$ & $0.73$ \\ \hline
$\frac{d_{A}(z_{\ast })}{D_{V}(z_{BAO})}$ & $30.95\pm 1.46$ & $17.55\pm 0.60$
& $10.11\pm 0.37$ & $8.44\pm 0.67$ & $6.69\pm 0.33$ & $5.45\pm 0.31$ \\
\hline
\end{tabular}
\caption{Values of $d_{A}(z_{\ast })/D_{V}(z_{BAO})$ for distinct values of $z_{BAO}$.}
\end{center}
\label{Table-1}
\end{table}
where

\begin{equation*}
X=\left(
\begin{array}{c}
\frac{d_{A}(z_{\star })}{D_{V}(0.106)}-30.95 \\
\frac{d_{A}(z_{\star })}{D_{V}(0.2)}-17.55 \\
\frac{d_{A}(z_{\star })}{D_{V}(0.35)}-10.11 \\
\frac{d_{A}(z_{\star })}{D_{V}(0.44)}-8.44 \\
\frac{d_{A}(z_{\star })}{D_{V}(0.6)}-6.69 \\
\frac{d_{A}(z_{\star })}{D_{V}(0.73)}-5.45%
\end{array}%
\right) \,,
\end{equation*}

and the inverse covariance matrix $C^{-1}$ is defined in
\cite{BAO6} as

\begin{equation*}
C^{-1}=\left(
\begin{array}{cccccc}
0.48435 & -0.101383 & -0.164945 & -0.0305703 & -0.097874 & -0.106738 \\
-0.101383 & 3.2882 & -2.45497 & -0.0787898 & -0.252254 & -0.2751 \\
-0.164945 & -2.454987 & 9.55916 & -0.128187 & -0.410404 & -0.447574 \\
-0.0305703 & -0.0787898 & -0.128187 & 2.78728 & -2.75632 & 1.16437 \\
-0.097874 & -0.252254 & -0.410404 & -2.75632 & 14.9245 & -7.32441 \\
-0.106738 & -0.2751 & -0.447574 & 1.16437 & -7.32441 & 14.5022%
\end{array}%
\right) \,.
\end{equation*}

\begin{figure}[H]
\centering
\includegraphics[scale=0.75]{BAOC}
\caption{The $1-\sigma$ and $2-\sigma$ likelihood contours for the model parameters using the BAO datasets.}
\label{ContourBAO}
\end{figure}

\end{widetext}

The values that fit observations are } $\zeta=0.66^{+0.12}_{-0.12}$, $\alpha=-0.0125^{+0.0047}_{-0.0023}$, $\gamma=1.191^{+0.070}_{-0.070}$ and $n=0.996^{+0.010}_{-0.012}$ in Fig.\ref{ContourBAO}.

\subsection{Pantheon datasets}

Scolnic et al. \cite{Scolnic/2018} put together the Pantheon samples consisting of 1048 type Ia supernovae in the redshift range $0.01 < z < 2.3$. The PanSTARSS1 Medium Deep Survey, SDSS, SNLS and numerous low-z and HST samples contribute to it. The empirical relation used to calculate the distance modulus of SNeIa from the observation of light curves is given by $\mu= m_{B}^{*}+\alpha X_{1}-\beta C-M_{B} + \Delta_{M}+ \Delta_{B}$, where $X_{1}$ and $C$ denote the stretch and color correction parameters, respectively \cite{Scolnic/2018,Mukherjee/2021}, $m_{B}^*$ represents the observed apparent magnitude and $M_{B}$ is the absolute magnitude in the B-band for SNeIa. The parameters $\alpha$ and $\beta$ are the two nuisance parameters describing the luminosity stretch and luminosity color relations, respectively. Further, the distance correction factor is $\Delta_{M}$ and $\Delta_{B}$ is a distance correction based on predicted biases from simulations.

The nuisance parameters in the Tripp formula \cite{Tripp/1998} were reconstructed using a novel technique called BEAMS with Bias Corrections \cite{Kessler/2017, Fotios/2021} and the observed distance modulus was reduced to the difference between the corrected apparent magnitude $m_{B}$ and the absolute magnitude $M_{B}$, which is $\mu= m_{B}-M_{B}$. We shall avoid marginalizing the over nuisance parameters $\alpha$ and $\beta$ but marginalize over the Pantheon data for $M_{B}$. Hence, we ignore the values of $\alpha$ and $\beta$ for the present investigation of the model.

The luminosity distance read as

\begin{equation}\label{4f}
D_{L}(z)= (1+z) \int_{0}^{z} \frac{c dz'}{H(z')},
\end{equation}
 with $c$ being the speed of light.

The $\chi^{2}$ function for type Ia supernovae is obtained by correlating the theoretical distance modulus

\begin{equation}\label{4g}
\mu(z)= 5log_{10}D_{L}(z)+\mu_{0},
\end{equation}
with
\begin{equation}\label{4h}
\mu_{0} =  5log(1/H_{0}Mpc) + 25,
\end{equation}
such that

\begin{equation}\label{4i}
\chi^{2}_{SN}(\zeta, \alpha, \gamma, n)= \sum _{k=1}^{1048} \dfrac{\left[ \mu_{obs}(z_{k})-\mu_{th}(\zeta, \alpha, \gamma, n, z_{k})\right] ^{2}}{\sigma^{2}(z_{k})},
\end{equation}
where $\mu_{th}$ is the theoretical value of distance modulus,  $\mu_{obs}$ is the observed value whereas $\sigma^{2}(z_{k})$ is the standard error in the observed value.

Using the Pantheon supernovae datasets, we obtain the best fit values for $\zeta$, $\alpha$, $\gamma$ and $n$ as the $1-\sigma$ and $2-\sigma$ contour plots in Fig.\ref{ContourPan}. The values that fit the model are $\zeta=0.67^{+0.12}_{-0.12} $, $\alpha=-0.00999^{+0.0047}_{-0.0024}$, $\gamma=1.34^{+0.15}_{-0.12}$ and $n=1.001^{+0.024}_{-0.024}$.

\begin{widetext}

\begin{figure}[H]
\centering
\includegraphics[scale=0.67]{PANC}
\caption{The $1-\sigma$ and $2-\sigma$ likelihood contours for the model parameters using the Pantheon datasets.}
\label{ContourPan}
\end{figure}

\end{widetext}

\subsection{Cosmological Parameters}

The evolution of deceleration parameter, energy density, pressure with bulk viscosity and the effective EoS parameter in redshift are presented below. In order to do so we use the set of values constrained by Hubble, BAO and the Pantheon data sets for the model parameters.

\begin{figure}[H]
{\includegraphics[scale=0.5]{q-z}}
\caption{Profile of the deceleration parameter for the given model corresponding to the values of the parameters constrained by the Hubble, BAO and the Pantheon data point sets.}
\label{q}
\end{figure}

\begin{figure}[H]
\centering
{\includegraphics[scale=0.5]{rho}}
\caption{Profile of the density parameter for the given model corresponding to the values of the parameters constrained by the Hubble, BAO and the Pantheon data point sets.}
\label{rho}
\end{figure}

\begin{figure}[H]
\centering
{\includegraphics[scale=0.5]{peff}}
\caption{Profile of the pressure for the given model corresponding to the values of the parameters constrained by the Hubble, BAO and the Pantheon data point sets.}
\label{p}
\end{figure}

\begin{figure}[H]
\centering
{\includegraphics[scale=0.5]{weff}}
\caption{Profile of the EoS parameter for the given model corresponding to the values of the parameters constrained by the Hubble, BAO and the Pantheon data point sets.}
\label{w}
\end{figure}

From Fig.\ref{q}, it is clear that the deceleration parameter shows the transition from a decelerated ($q>0$) to an accelerated ($q<0$) phase of the universe expansion for the constrained values of the model parameters. The transition redshift is $z_t \thickapprox 0.142$,   $ z_t \thickapprox 0.776 $ and $ z_t \thickapprox 0.622 $ corresponding to the Hubble, BAO and the Pantheon data sets, respectively. The present value of the deceleration parameter is, respectively,  $q_0 = -0.127$ , $q_0= -0.436$ and  $q_0= -0.454$.

From Fig.\ref{rho} it is clear that the density parameter shows the positive behavior for all the constrained values of the model parameters, which is expected.

Fig \ref{p} indicates that the bulk viscous cosmic fluid exhibits, for lower redshifts, the negative pressure that make bulk viscosity to be a viable candidate to drive the cosmic acceleration. Furthermore, the effective EoS parameter presented in Fig.\ref{w} indicates that the cosmic viscous fluid behaves like quintessence dark energy. The present values of EoS parameter corresponding to the Hubble, BAO and the Pantheon samples are  $w_0= -0.433$ , $w_0= -0.625$, and  $w_0= -0.635$.

\section{Statefinder analysis}\label{sec5}

The cosmological constant $\Lambda$ suffers from two major drawbacks, namely the aforementioned cosmological constant and cosmic coincidence problems. To surmount these problems, dynamic models of dark energy have been introduced in the literature as we have also mentioned in Introduction. To discriminate between these time-varying dark energy  models, an appropriate tool was required. In this direction, V. Sahni et al. introduced a new pair of geometrical parameters known as {\it statefinder parameters} $(r,s)$ \cite{VS}. The statefinder parameters are defined as
\begin{equation}
r=\frac{\dddot{a}}{aH^3},
\end{equation}
\begin{equation}
s=\frac{(r-1)}{3(q-\frac{1}{2})}
\end{equation}
The parameter $r$ can be rewriten as
\begin{equation}
r=2q^2+q-\frac{\dot{q}}{H}.
\end{equation}

For different values of the statefinder pair $(r,s)$, it represents the following dark energy models:
\\ \\
${r=1,s=0}$ represents $\Lambda$CDM model,\\
${r>1,s<0}$ represents Chaplygin gas model,\\
${r<1,s>0}$ represents quintessence model.\\

In Figures \ref{rs} and \ref{qr}, we plot the $s-r$ and $q-r$ diagrams for our cosmological model by taking the values of the parameters constrained by the Hubble, BAO and the Pantheon data sets.

\begin{figure}[H]
\includegraphics[scale=0.5]{r-s-fig}
\caption{Plot of the trajectories in the $r-s$ plane for the given cosmological model corresponding to the parameters values constrained by the Hubble, BAO and the Pantheon data sets.}
\label{rs}
\end{figure}

\begin{figure}[H]
\includegraphics[scale=0.5]{q-r-fig}
\caption{Plot of the trajectories in the $q-r$ plane for the given cosmological model corresponding to the parameters values constrained by the Hubble, BAO, and the Pantheon data sets.}
\label{qr}
\end{figure}

Figures \ref{rs} and \ref{qr} show that our bulk viscous model lies in the quintessence region. Also, the evolutionary trajectories of our model departure from the $\Lambda$CDM point. The present values of the statefinder parameters corresponding to the values of the model parameters constrained by the Hubble, BAO and the Pantheon samples are
$ r_0= 0.837$ and $s_0= 0.086$, $r_0= 0.62$ and $s_0= 0.135$,  $r_0= 0.762$ and $s_0=  0.083$ respectively.

\section{Discussions and conclusions}\label{sec6}

In the present section we will discuss the results obtained in Sections IV and V for the bulk viscous symmetric teleparallel cosmological model here developed and presented.

Cosmology has been on the agenda mainly for two reasons: dark energy and dark matter. While dark energy has been deeply discussed throughout the paper, dark matter is predicted within $\Lambda$CDM model as a sort of matter that does not interact electromagnetically, so that it cannot be seen, but its gravitational effects otherwise can well be detected. Still we have not yet detected or even associated dark matter to a particle of standard model or beyond \cite{cdmsii_collaboration/2010,akerib/2014,essig/2012}. Modified (or alternative) theories of gravity have also been used to describe dark matter effects \cite{bohmer/2008,mannheim/2012}. In these cases, dark matter is simply an effect of modification of gravity.

Returning to the dark energy question, it is highly counter-intuitive that the expansion of the universe is actually accelerating. Although the vacuum quantum energy can well explain this  dynamical effect via the cosmological constant in Einstein's field equations of GR, the aforementioned important and persistent problems related to $\Lambda$ supply the search for alternative explanations.

In the present article, as an attempt to describe dark energy we have assumed the symmetric teleparallel gravity as the underlying gravity theory.

The $f(Q)$ gravity was recently proposed by Jim\'enez et al. in \cite{J.B.} as a ramification of the geometric trinity, that says that space-time manifold can be described by curvature, torsion or non-metricity. Particularly, the symmetric teleparallel gravity describes gravitational interactions via the non-metricity scalar, with null curvature and torsion.

Our $f(Q)$ cosmological model was based on a spatially homogeneous and isotropic flat metric and an energy-momentum tensor describing a bulk viscous fluid. We let the $f(Q)$ function be a power of $n$ as $f(Q)\sim Q^n$, with $n$ a free parameter.

In Section IV we started testing our cosmological solutions. We started confronting the Hubble parameter \eqref{3m} with 31 data points that are measured from differential age approach. In Fig.\ref{ContourHub} we have obtained the best fit values for the model free parameters. We have then plotted Fig.\ref{ErrorHubble}, in which our $H(z)$ is confronted with cosmological data and compared with $\Lambda$CDM prediction. We can see the $f(Q)$ model describes observations with good agreement and specially for higher values of redshift it is clear that it provides a better fit when compared to $\Lambda$CDM model. Further, in Figs.\eqref{ContourBAO} and \eqref{ContourPan}, we have obtained the best fit values for the model free parameters.

Fig.\ref{p} indicates that the bulk viscous cosmic fluid exhibit the negative pressure that make bulk viscosity to be a viable candidate to drive the cosmic acceleration. This is also reflected in the deceleration parameter behavior in Fig.\ref{q}, which shows a transition from decelerated to accelerated phases of the universe expansion. Furthermore, the effective EoS parameter presented in Fig.\ref{w} indicates that the cosmic viscous fluid behaves like quintessence dark energy.

In Section V, Figs.\ref{rs} and \ref{qr} show that the evolutionary trajectories of our model is departed from $\Lambda$CDM fixed point {${r=1,s=0}$}. In the present epoch they lie in the quintessence region {${r<1,s>0}$}. The present model is, therefore, a good alternative to explain the universe dynamics, particularly with no necessity of invoking the cosmological constant.

 Finally, it is worth remarking that $f(Q)$ gravity still needs to be applied in several different regimes in order to be stablished as a viable gravitational formalism. In \cite{J.B.} it was shown that the $f(Q)$ gravity can be compatible with Solar System constraints. Recently, in \cite{frusciante/2021} it was shown that the $f(Q)$ gravity predicts a modified gravitational wave propagation and a different gravitational wave luminosity distance when compared to the standard electromagnetic one, although this investigation was done for a different functional form of the $f(Q)$ function.

{\bf ACKNOWLEDGMENTS}\\

RS acknowledges University Grants Commission (UGC), New Delhi, India for awarding Junior Research Fellowship (UGC-Ref. No.: 191620096030). SA  acknowledges  CSIR,  New  Delhi,  India  for  JRF. PKS   acknowledges   CSIR,   New   Delhi,   India   for   financial   support   to   carry   out   the   Research   project [No.03(1454)/19/EMR-II Dt.02/08/2019].


\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}