# Manifold embeddings by heat kernels of connection Laplacian

Chen-Yun Lin  Chen-Yun Lin
Department of Mathematics
Lehman College, CUNY<EMAIL_ADDRESS>
###### Abstract.

We show that any closed $n$-dimensional manifold $(M,g)$ can be embedded by a
map constructed using the heat kernels of the connection Laplacian as well as
a maps constructed using truncated heat kernel at a certain time $t$ from a
$\delta$-net $\\{q_{i}\\}_{i=1}^{N_{0}}$ via a rescaling trick. Both the time
$t$ and $N_{0}$ are bounded in terms of the dimension, bounds on the Ricci
curvature and its derivative, the injectivity radius, and the volume.
Moreover, both maps can be made arbitrarily close to an isometry.

##  1\. Introduction

Data collected for the purpose of machine learning is often in a high-
dimensional space, but yet is believed to satisfy certain low dimensional
structure, that is, the collected dataset can be well approximated by a low
dimensional manifold sitting inside a high dimensional Euclidean space. See
[3, 10, 21, 15] for a far from complete list of available literature. How to
analyze a dataset under this assumption is generally called the manifold
learning problem. One particular goal is to recover the nonlinear low
dimensional structure of the manifold and to reduce the dimensionality of the
space where the dataset lies inside. Mathematically, this problem is
formulated as asking if it is possible to have an embedding to put the
manifold (hence the dataset) into a finite dimensional Euclidean space that is
bi-Lipschitz, even isometric. Although the embedding problem was first
positively answered by Whitney [36], and the isometrically embedding problem
was first solved by Nash [27], the approaches are not canonical and are not
essentially feasible for data analysis. In Berard-Besson-Gallot’s breakthrough
paper [5], the spectral embedding idea is explored. They show that a manifold
can be embedded into the sequence space $\ell^{2}$ with all eigenfunctions.

The spectral embedding idea is directly related to many manifold learning
algorithms, like eigenmaps [3], local linear embeddings [29] and diffusion
maps [10]. There is also a rich theoretical literature describing how the
graph Laplacian converges to the Laplace-Beltrami operator; for example, the
pointwise convergence wast discussed in [3, 10, 34] and the spectral
convergence was discussed in [4, 32, 14, 35, 6, 11], where the $L^{\infty}$
spectral convergence with rate was recently explored in [11]. Note that
numerically we are able to obtain only finite eigenfunctions and eigenvalues
in practice since only finite number of points are available. Thus, the next
natural question is whether we can embed the manifold by finite eigenfunctions
and eigenvalues. This question was positively answered separately by Bates [2]
and Portegies [28] for the Laplace-Beltrami operator; that is, one is able to
construct bi-Lipschitz embeddings of the manifold with finite eigenfunctions
and eigenvalues of the Laplace–Beltrami operator. In Portegies [28], it is
further shown that the embedding is almost isometric with a prescribed error
bound. Moreover, Jones et al. [23] observed that in terms of analysis it may
be easier to think of embedding with heat kernels.

The spectral embedding mentioned above depends on the Laplace-Beltrami
operator. The vector diffusion map (VDM) [31, 32], on the other hand, depends
on the connection Laplacian associated with a possibly nontrivial bundle
structure. The VDM is originally motivated by studying the cryo-electron
microscope problem, in particular the class averaging algorithm [19, 33, 37].
Its general goal is to integrate local information and the relationship
between these pieces of local information in order to obtain the global
information of the dataset; for example, the ptychographic imaging problem
[26], the synchronization problem [1], the vector nonlocal mean/median [24],
the orientability problem [30], etc. Numerically, the VDM depends on the
spectral study of the graph connection Laplacian [8, 7, 12, 13, 31, 32] ,
which is a direct generalization of the graph Laplacian discussed in the
spectral graph theory [9]. In [25], we showed that one is able to construct an
embedding of a manifold with finite eigenvector fields of the connection
Laplacian that is bi-Lipschitz.

The main contribution in this paper is to show that one can construct an
embedding of a manifold with heat kernels of the connection Laplacian that is
an almost isometry. Our approach is inspired by Portegies’ rescaling
technique. Note that due to the lack of canonical isometry from $L^{2}(TM)$ to
$\ell^{2}$, we do not have a similar almost isometry result for embeddings
with eigenvector fields. In Section 3, we introduce background on the heat
kernels and state the main results. In Section 4, we review the
$C^{k,\alpha}$-harmonic radius and introduce the rescaled heat kernel
$\bar{K}_{TM}$ defined in (4.5) that satisfies the rescaled parabolic system
(4.6). Following [17, Section 4 in Chapter 9], we construct a fundamental
solution $\Gamma$ for (4.6) and show that $\Gamma$ is close to the standard
Euclidean heat kernel on a local harmonic chart. We also show that any
solution with sufficient decay is close to this fundamental solution $\Gamma$.
Thus, the rescaled heat kernel $\bar{K}_{TM}$ is close to the stand Euclidean
heat kernel. In Section 5, we prove our main Theorems.

##  2\. Acknowledgments

I would like to thank Hau-Tieng Wu for the support as well as the helpful
discussions and comments. I would also like to thank Duke University for the
research environment. I am also supported by PSC-CUNY 63043-00 51.

##  3\. Background and Main Results

In this section, we review some background that is required for the main
results, as well as the Gaussian bound for the heat kernel that is used in the
proofs of the main results.

###  3.1. Heat kernel of connection Laplacian

Let $M$ be an $n$-diimensional smooth closed manifold. We denote by
$K_{TM}(p,t;q)$ the heat kernel of the connection Laplacian $\Delta_{TM}$ on
the manifold $M$. The heat kernel can be expressed as

$K_{TM}\left(p,t;q\right)=\sum_{i=1}^{\infty}e^{-\lambda_{i}t}X_{i}(p)\otimes
X_{i}(q),$

where $t>0$, $p,q\in M$, $\lambda_{i}$ are eigenvalues so that
$0\leq\lambda_{1}\leq\lambda_{2}\leq\cdots$ and $\\{X_{i}\\}_{i=1}^{\infty}$
is an $L^{2}(TM)$-orthonormal basis formed by eigenvector fields of
$\Delta_{TM}$. We define the $m$-th order truncated heat kernel $K_{TM}^{(m)}$
by

$K_{TM}^{(m)}(p,t;q):=\sum_{i=1}^{m}e^{-\lambda_{i}t}X_{i}(p)\otimes
X_{i}(q).$

Let $\|\cdot\|_{HS}$ denote the Hilbert-Schmidt norm. We have

$\displaystyle\left\|K_{TM}\left(p,t;q\right)\right\|_{HS}^{2}$
$\displaystyle=\mathrm{Tr}\left(K_{TM}(p,t;q)^{*}K_{TM}(p,t;q)\right)$
$\displaystyle=\sum_{i,j=1}^{\infty}e^{-\left(\lambda_{i}+\lambda_{j}\right)t}\left\langle
X_{i}(p),X_{j}(p)\right\rangle\left\langle X_{i}(q),X_{j}(q)\right\rangle;$

and

$\left\|K_{TM}^{(m)}\left(p,t;q\right)\right\|_{HS}^{2}=\sum_{i,j=1}^{m}e^{-\left(\lambda_{i}+\lambda_{j}\right)t}\left\langle
X_{i}(p),X_{j}(p)\right\rangle\left\langle X_{i}(q),X_{j}(q)\right\rangle.$

###  3.2. Dilatation

Given two metric spaces $(\mathfrak{M}_{1},d_{1})$ and
($\mathfrak{M}_{2},d_{2})$, the local dilatation of a map
$f:\mathfrak{M}_{1}\rightarrow\mathfrak{M}_{2}$ at a point
$p\in\mathfrak{M}_{1}$ is defined as

$\mathrm{dil}_{p}(f):=\lim_{r\rightarrow 0}\sup_{x,y\in
B_{r}(p)}\frac{d_{2}(f(x),f(y))}{d_{1}(x,y)}\,,$

where $B_{r}(p)$ is a ball in $\mathfrak{M}_{1}$ of radius $r$ and center $p$.
In the special case that $M$ is a smooth Riemannian manifold mapped into a
normed, finite-dimensional vector space via a smooth map $f$, the dilatation
is given by

$\mathrm{dil}_{p}(f)=|(df)_{p}|\,,$

where the norm of the right-hand side is interpreted as the operation norm of
the map from $T_{p}M$ to $T_{f(p)}f(M)$.

###  3.3. Main results

Consider a set of closed smooth manifolds of dimension $n$

$\mathcal{M}_{n,\kappa,i_{0},V}:=\\{(M^{n},g):|Rc|,|\nabla
Rc|\leq\kappa,\,\mathrm{inj(M)}\geq i_{0},\,\mathrm{Vol}(M)\leq V\\},$

where $\mathrm{inj}(M)$ denotes the injectivity radius of $M$ and $Rc$ and
$\nabla Rc$ denote the Ricci curvature and its covariant derivative,
respectively.

We prove the following theorems:

######  Theorem 1.

Given $\epsilon>0$, there exists $t_{0}=t_{0}(n,\kappa,i_{0},\epsilon)$ so
that for all $t<t_{0}$ there exists $\delta=\delta(n,\kappa,i_{0},\epsilon,t)$
such that for all $(M,g)\in\mathcal{M}_{n,\kappa,i_{0},V}$ and every
$\delta$-net $\\{q_{1},\cdots,q_{N_{0}}\\}$ on $M$, the map

$\displaystyle H:M^{n}$ $\displaystyle\rightarrow\mathbb{R}^{N_{0}}$ (3.1)
$\displaystyle p$
$\displaystyle\mapsto\frac{(2t)^{(3n+2)/4}}{V_{e}}\left(|A_{1}|^{1/2}\|K_{TM}(p,t;q_{1})\|^{2}_{HS},\cdots,|A_{N_{0}}|^{1/2}\|K_{TM}(p,t;q_{N_{0}})\|^{2}_{HS}\right)$

is an embedding and the dilatation is controlled by a given $\epsilon$,

(3.2) $1-\epsilon<|(dH)_{p}|^{2}<1+\epsilon\quad\mbox{for all $p\in M$}.$

Here $\\{A_{i}\\}_{i=1}^{N_{0}}$ is a partition of $M$ so that $A_{i}\subset
B_{\delta}(q_{i})$ for all $i$ and

(3.3)
$V_{e}:=\left(\int_{\mathbb{R}^{n}}(\partial_{x_{1}}\|\Gamma_{E}(0,\frac{1}{2};y)\|^{2}_{HS})^{2}dy\right)^{1/2}\,,$

where $\Gamma_{E}$ is the standard Euclidean heat kernel on
$\mathbb{R}^{n}\times\mathbb{R}$ whose definition is reviewed in (4.7).

Moreover, there exists an $N=N(n,\kappa,i_{0},V,\epsilon,t)$ so that for all
$m\geq N$, the same statements hold for the map $H^{(m)}$ defined by truncated
heat kernels

(3.4)
$H^{(m)}(p):=\frac{(2t)^{\frac{(3n+2)}{4}}}{V_{e}}\left(|A_{1}|^{\frac{1}{2}}\|K^{(m)}_{TM}(p,t;q_{1})\|^{2}_{HS},\cdots,|A_{N_{0}}|^{\frac{1}{2}}\|K^{(m)}_{TM}(p,t;q_{N_{0}})\|^{2}_{HS}\right)$

We comment that the $N_{0}$ is not necessary universal. It is for every choice
of covering sets that gives a $\delta$-net, where $\delta$ depends on the
geometry, $\epsilon$ and $t$.

######  Theorem 2.

Given $\epsilon>0$, there exists $t_{0}=t_{0}(n,\kappa,i_{0},\epsilon)$ so
that for all $t<t_{0}$ there exists $N_{0}=N_{0}(n,\kappa,i_{0},\epsilon,t,V)$
such that for all $(M,g)\in\mathcal{M}_{n,\kappa,i_{0},V}$, there exist points
$\\{q_{1},\cdots,q_{N_{0}}\\}$ on $M$ such that the map

$\displaystyle\mathcal{H}:M^{n}$ $\displaystyle\rightarrow\mathbb{R}^{N_{0}}$
(3.5) $\displaystyle p$
$\displaystyle\mapsto\frac{(2t)^{(3n+2)/4}}{V_{e}}A\left(\|K_{TM}(p,t;q_{1})\|^{2}_{HS},\cdots,\|K_{TM}(p,t;q_{N_{0}})\|^{2}_{HS}\right)$

where $A=A(n,\kappa,i_{0},V,\epsilon,t)>0$ is a small constant, is an
embedding and the local dilatation satisfies

(3.6) $1-\epsilon<|dH_{p}|<1+\epsilon\quad\mbox{for all $p\in M$.}$

Moreover, there exists an $N=N(n,\kappa,i_{0},V,\epsilon,t)$ so that for all
$m\geq N$, the same statements hold with every heat kernel $K_{TM}$ replaced
by the truncated heat kernel $K^{(m)}_{TM}$,

$\displaystyle\mathcal{H}^{(m)}:M^{n}$
$\displaystyle\rightarrow\mathbb{R}^{N_{0}}$ (3.7) $\displaystyle p$
$\displaystyle\mapsto\frac{(2t)^{(3n+2)/4}}{V_{e}}A\left(\|K_{TM}^{(m)}(p,t;q_{1})\|^{2}_{HS},\cdots,\|K_{TM}^{(m)}(p,t;q_{N_{0}})\|^{2}_{HS}\right)$

Note that this $N_{0}$ is universal as we choose a partition that depends on
the geometry, $t$ and $\epsilon$.

In order to prove the main results, we heavily rely on the Gaussian bound of
the heat kernel of the connection Laplacian, which is stated below.

###  3.4. Gaussian bound of the heat kernel

Let $K_{M}(p,t;q)$ denote the heat kernel of the Laplace-Beltrami operator
$\Delta_{M}$ on $(M,g)$. By the Faber-Krahn inequality [18, Section 14.2] and
the exponential decay [18, Theorem 15.14] of the heat kernel $K_{M}(p,t;q)$,
we have for $p,q$ in a ball of radius less than the $C^{2,\alpha}$-harmonic
radius $r_{h}$, whose definition is reviewed in Subsection 4.1,

(3.8)
$K_{M}(p,t;q)\leq\frac{C(n)(1+d(p,q)^{2}/t)^{n/2}}{(a(n)\min(t,r_{h}^{2}))^{n/2}}\exp\left(-\frac{d(p,q)^{2}}{4t}\right)$

where $d(p,q)$ denotes the geodesic distance between $p,q\in M$ for constants
$a(n)>0$ and $C(n)>0$ depending only on $n$.

In addition, the semigroup domination theorem by Hess-Schrader-Uhlenbrock [22,
Theorem 3.1] states the following in our setting

$\textrm{Tr}\left(e^{t\Delta_{TM}}\right)\leq
n\textrm{Tr}\left(e^{t\Delta_{M}}\right),$

which implies that

(3.9) $\displaystyle\|K_{TM}(p,t;q)\|_{HS}\leq
nK_{M}(p,t;q)\leq\frac{C(n)(1+d(p,q)^{2}/t)^{n/2}}{(a(n)\min(t,r_{h}^{2}))^{n/2}}\exp\left(-\frac{d(p,q)^{2}}{4t}\right).$

##  4\. Foundamental Solutions on Charts

The goal in this section is to show that the rescaled heat kernel,
$\bar{K}_{TM}$ defined in (4.5), is close to the standard heat Kernel
$\Gamma_{E}$. We start with reviewing the harmonic radius and there is a
universal lower bound for the harmonic radius that is independent of the point
on the manifold. This allows us to work on a harmonic chart of certain size.
In Lemma 5, we construct a fundamental solution $\Gamma$ of the rescaled
parabolic system (4.6) that is close to the Euclidean heat kernel using the
parametrix method as presented in Friedman’s book [17, Section 4 in Chapter
9]. In Lemma 6, we show that the rescaled heat kernel $\bar{K}_{TM}$ and the
fundamental solution $\Gamma$ are close by applying the Schauder estimates
[16] and the semigroup domination theory. Throughout the paper, the constant
$C$ varies line by line.

###  4.1. The $C^{k,\alpha}$-Harmonic Radius

Recall that (see [20]) given $Q>1$ and $\alpha\in(0,1)$, the
$C^{k,\alpha}$-harmonic radius at $p$ in $M$ is defined as the largest number
$r_{h}=r_{h}(Q,k,\alpha)(p)$ such that on the geodesic ball $B_{r_{h}}(p)$ of
center $p$ and radius $r_{h}$, there is a harmonic coordinate chart
$u:B_{r_{h}}(p)\subset M\rightarrow U\subset\mathbb{R}^{n}$ such that
$u(p)=0$, $g_{ab}(p)=\delta_{ab}$ and the metric tensor is $C^{k,\alpha}$
controlled in these coordinates. Namely, if $g_{ab}$, $a,b=1,\cdots,n$ are the
components of $g$ in these coordinates, then

(4.1) $\displaystyle Q^{-1}\delta_{ab}\leq g_{ab}\leq Q\delta_{ab}\mbox{ as
bilinear forms }$ (4.2) $\displaystyle\displaystyle\sum_{1\leq|\beta|\leq
k}r_{h}^{|\beta|}\sup_{p}|\partial^{\beta}g(p)|+\sum_{|\beta|=k}r_{h}^{2+\alpha}\sup_{p\neq
q}\frac{|\partial^{\beta}g_{ij}(p)-\partial^{\beta}g_{ij}(q)|}{d(p,q)^{\alpha}}\leq
Q-1\,,$

where $d(p,q)$ denotes the geodesic distance on $(M,g)$. The harmonic radius
of $(M,g)$, denoted as $r_{h}(Q,k,\alpha)(M)$, is defined by

(4.3) $\displaystyle r_{h}(Q,k,\alpha)(M):=\inf_{p\in M}r_{h}(Q,k,\alpha)(p).$

For any manifold $(M,g)\in\mathcal{M}_{n,\kappa,i_{0},V}$, it is proved by
Hebey-Herzlich [20, Corollary of Theorem 6] that there exist
$C^{2,\alpha}$-harmonic coordinate charts whose size is independent of the
choice of points. Here, we restate the Corollary for our case.

######  Corollary 3.

Let $\alpha\in(0,1)$ and $Q>1$. Let $(M,g)$ be a smooth closed $n$-dimensional
Riemannian manifold. Suppose for some $\kappa>0$, $i_{0}>0$

$|Rc|,|\nabla Rc|\leq\kappa\quad\mbox{and}\quad\mathrm{inj}(M)\geq i_{0}\,.$

Then, there exists a constant $C=C(n,Q,\alpha,\kappa,i_{0})$ such that the
$C^{2,\alpha}$-harmonic radius

$r_{h}=r_{h}(Q,2,\alpha)(M)\geq C.$

That is, for any $M\in\mathcal{M}_{n,\kappa,i_{0},V}$, there is a universal
$C^{2,\alpha}$-harmonic radius $r_{h}$ that is independent of the choice of
point. Let $p\in M$ and let $u:B_{r_{h}}(p)\rightarrow\mathbb{R}^{n}$ be a
harmonic coordinate chart with $u(p)=0$ satisfying (4.1) and (4.2). For a
vector field $X=X^{a}\partial_{a}$, by a direct calculation, we have

(4.4)
$\displaystyle\Delta_{TM}\left(X^{c}\partial_{c}\right)=-g^{ab}\left(\partial_{a}\partial_{b}X^{c}+2\Gamma_{bd}^{c}\partial_{a}X^{d}+\Gamma_{ae}^{c}\Gamma_{bd}^{e}X^{d}+\partial_{a}\Gamma_{bd}^{c}X^{d}\right)\partial_{c}.$

Note that the coefficients of $\Delta_{TM}X$ are controlled in
$C^{2,\alpha}$-harmonic coordinates by (4.2).

######  Remark 4.

In [25] where we have similar results but without discussing the Lipschitz
control, we write the term $\partial_{a}\Gamma_{bd}^{c}X^{d}$ in terms of the
Ricci curvature and hence having $|Rc|<\kappa$ is enough to have
$C^{1,\alpha}$-harmonic radius and bounded coeffients. However, in this paper
we will start with a rougher guess in the parametrix process to build a
solution that is close to the stand Euclidean heat kernel and need subtler
controls of the coefficients. Therefore, the class of manifold
$\mathcal{M}_{n,\kappa,i_{0},V}$ requires the extra condition $|\nabla
Rc|\leq\kappa$ to guarantee a universal $C^{2,\alpha}$-harmonic radius bound.

###  4.2.  Rescaled Heat Kernel

Let $u:B_{r_{h}}(p)\subset M\rightarrow U\subset\mathbb{R}^{n}$ be a harmonic
coordinate chart with $u(p)=0$. Define a rescaled heat kernel on
$\mathbb{R}^{n}\times\mathbb{R}_{+}$ associated with the heat kernel $K_{TM}$
of the connection Laplacian as:

(4.5) $\bar{K}_{TM}(x,s;y):=r^{n}K_{TM}(u^{-1}(xr),sr^{2};u^{-1}(yr)),$

where $s>0$, $r>0$ is a scale factor to be chosen later, and $x,y$ in the
proper domain, $U/r:=\\{x\in\mathbb{R}^{d}|\,rx\in U\\}\subset\mathbb{R}^{n}$.
Since $u$ is a differomorphism, denote $(u^{-1})^{*}g$ to be the induced
metric on $U$, and denote the rescaled metric as
$\bar{g}(x)=(u^{-1})^{*}g(rx)$. Denote the Christophe symbol of
$(u^{-1})^{*}g$ as $\Gamma_{ab}^{c}$, we have the Christophe symbol for the
rescaled metric satisfying $\bar{\Gamma}_{ab}^{c}(x)=\Gamma_{ab}^{c}(rx)$. On
$(U,(u^{-1})^{*}g)$, denote $L=\Delta_{(u^{-1})^{*}g}-\partial_{s}$ to be the
parabolic system, and $\bar{L}$ to be the rescaled version. By a direct
calculation, it satisfies the following expansion:

(4.6)
$\bar{L}X^{c}:=\bar{g}^{ab}\left(\partial_{a}\partial_{b}X^{c}+2r\bar{\Gamma}_{bd}^{c}\partial_{a}X^{d}+r^{2}\bar{\Gamma}_{ae}^{c}\bar{\Gamma}_{bd}^{e}X^{d}+r^{2}\partial_{a}\bar{\Gamma}_{bd}^{c}X^{d}\right)-\partial_{s}X^{c}=0,$

where $c=1\cdots,n$ and $X^{c}$ is a function defined on $U/r$.

We consider the standard Euclidean heat kernel on
$\mathbb{R}^{n}\times\mathbb{R}$ given by

(4.7) $\Gamma_{E}(x,t;y)=\frac{1}{(4\pi
t)^{n/2}}e^{-\frac{|x-y|^{2}}{4t}}I_{n},$

where $t>0$, $x,y\in\mathbb{R}^{n}$, and $I_{n}$ is the $n\times n$ identity
matrix. Note that $\Gamma_{E}$ is the heat kernel of the connection Laplacian
associated with the trivial tangent bundle of the canonical Euclidean space
$\mathbb{R}^{n}$. Below, we show that $\bar{K}_{TM}$ is close to the standard
Euclidean heat kernel $\Gamma_{E}$.

###  4.3. Foundamental solutions on charts

Denote $P_{R,T}(x)$ to be the parabolic cylinder on
$\mathbb{R}^{n}\times\mathbb{R}$; that is,

$P_{R,T}(x):=\left\\{(y,s)\in\mathbb{R}^{n}\times\mathbb{R}^{+}:0<s<T,|y-x|<R\right\\}\,,$

where $x\in\mathbb{R}^{n}$, $R,T>0$.

Define $Z(x,s;y)$ on $\mathbb{R}^{n}$ as

(4.8) $Z(x,s;y):=\frac{\sqrt{|\bar{g}(y)|}}{(4\pi
s)^{n/2}}\exp\left(-\frac{\bar{g}_{ab}(y)(x^{a}-y^{a})(x^{b}-y^{b})}{4s}\right)I_{n},$

where $|\bar{g}(y)|:=\det(\bar{g}_{ab}(y))$. Note that $Z(x,s;y)$ is a
solution of the principal part of the rescaled system
$\partial_{s}X^{c}=\bar{g}^{ab}(y)\partial_{a}\partial_{b}X^{c}$.

######  Lemma 5.

Let $\epsilon>0$, $0<\alpha<1$, $1<Q<\sqrt{2}$ and $R>0$ be given. Assume that
$r$ is small so that $rR<r_{h}/\sqrt{2}$, where $r_{h}$ is the
$C^{2,\alpha}$-harmonic radius. Consider the rescaled parabolic system (4.6)
on a domain $\Omega\times I\subset\mathbb{R}^{n}\times\mathbb{R}_{+}$ that
contains the rescaled domain $P_{R,4}(0)$. There exists a fundamental solution
$\Gamma$ for (4.6) on the domain $\Omega\times I$ such that

(4.9)
$|\Gamma(x,s;y)|\leq\frac{C(n)}{s^{n/2}}e^{-\frac{|x-y|^{2}}{8s}}\quad\mbox{
for $0<s<2$}.$

In particular, for all $y\in B_{R}(0)\subset\Omega$ and $(x,s)\in P_{R,4}(0)$,

(4.10) $\displaystyle|\Gamma(x,s;y)-Z(x,s;y)|$
$\displaystyle\leq\frac{(Q-1)C(n,\alpha)}{s^{(n-\alpha)/2}}e^{-\frac{|x-y|^{2}}{8s}};$
(4.11) $\displaystyle|\nabla\Gamma(x,s;y)-\nabla Z(x,s;y)|$
$\displaystyle\leq\frac{(Q-1)C(n,\alpha)}{s^{(n+1-\alpha)/2}}e^{-\frac{|x-y|^{2}}{8s}}.$

Moreover, for $(x,s)\notin P_{\frac{1}{2},\frac{1}{4}}(y)$,

(4.12) $|\Gamma(x,s;y)-\Gamma_{E}(x,s;y)|\leq(Q-1)C(n,\alpha);$ (4.13)
$|\nabla\Gamma(x,s;y)-\nabla\Gamma_{E}(x,s;y)|\leq(Q-1)C(n,\alpha).$

###### Proof.

Following the parametrix method in Friedman [17, Section 4 in Chapter 9 and
Section 4 in Chapter 1], we construct $\Gamma$ in the form

$\Gamma(x,s;y)=Z(x,s;y)+\int_{0}^{s}\int_{\Omega}Z(x,s-\tau;\xi)\Phi(\xi,\tau;y)d\xi
d\tau$

for some $\Phi$ which is an $n\times n$ matrix-valued function defined on
$\Omega\times I$. If $\Phi$ is Hölder continuous, then $\Gamma$ satisfies the
system (4.6) as a function of $(x,s)$ if and only if

(4.14)
$\Phi(x,s;y):=LZ(x,s;y)+\int_{0}^{s}\int_{\Omega}LZ(x,s-\tau;\xi)\Phi(\xi,\tau;y)d\xi
d\tau.$

(See [17, Lemma 5 in Chapter 9, Section 4, p. 250] for details.) We show that
there is a $\Phi$ of the form

$\Phi(x,s;y)=\sum_{i=1}^{\infty}(LZ)_{i}(x,s;y),$

where $(LZ)_{1}=LZ$ and

$(LZ)_{i}=\int_{0}^{s}\int_{\Omega}LZ(x,s-\tau;\xi)(LZ)_{i-1}(\xi,\tau;y)d\xi
d\tau,$

so that $\Phi$ is a formal solution of (4.14). To find such $\Phi$, we first
show that $(LZ)_{i}$ is integrable for all $i\in\mathbb{N}$ and then show the
convergence of the series. Since $Z$ is diagonal and identical along the
diagonal. It suffices to look at the scalar function

$z(x,s-\tau;y):=\frac{\sqrt{|\bar{g}(y)|}}{(4\pi(s-\tau))^{n/2}}\exp\left(-\frac{\bar{g}_{ab}(y)(x^{a}-y^{a})(x^{b}-y^{b})}{4(s-\tau)}\right).$

For simplicity, denote $A:=\bar{g}_{ab}(y)(x^{a}-y^{a})(x^{b}-y^{b})$. By
rewriting

$\displaystyle\exp\left(-\frac{A}{4(s-\tau)}\right)=\exp\left(-\epsilon\frac{A}{4(s-\tau)}\right)\exp\left(-(1-\epsilon)\frac{A}{4(s-\tau)}\right),$

and using the inequality
$\sigma^{n/2-\mu}e^{-\epsilon\sigma}\leq\mathrm{constant}$ for fixed constants
$\mu,\epsilon$, and $0\leq\sigma<\infty$, we obtain a bound for
$z(x,s-\tau;y)$:

$\displaystyle|z(x,s-\tau;y)|$
$\displaystyle=\frac{\sqrt{|\bar{g}(y)|}}{\pi^{n/2}(4(s-\tau))^{\mu}}A^{\mu-n/2}\left(\frac{A}{4(s-\tau)}\right)^{n/2-\mu}\exp\left(-\frac{A}{4(s-\tau)}\right)$
$\displaystyle=\frac{\sqrt{|\bar{g}(y)|}}{\pi^{n/2}(4(s-\tau))^{\mu}}A^{\mu-n/2}\left(\frac{A}{4(s-\tau)}\right)^{n/2-\mu}$
$\displaystyle\qquad\qquad\times\exp\left(-\epsilon\frac{A}{4(s-\tau)}\right)\exp\left(-(1-\epsilon)\frac{A}{4(s-\tau)}\right)$
$\displaystyle\leq\frac{C}{(s-\tau)^{\mu}|x-y|^{n-2\mu}}\exp\left(-\frac{\lambda_{0}|x-y|^{2}}{4(s-\tau)}\right),$

for any $\lambda_{0}<Q^{-1},0\leq\mu\leq n/2$, and some constant $C=C(n,Q)$.
Similarly, it can be proved that

$\displaystyle\sum_{c=1}^{n}\left|\frac{\partial z(x,s-\tau;y)}{\partial
x^{c}}\right|\leq\frac{C}{(s-\tau)^{\mu}|x-y|^{n+1-2\mu}}\exp\left(-\frac{\lambda_{0}|x-y|^{2}}{4(s-\tau)}\right)$

for any $\lambda_{0}<Q^{-1}$ and $0\leq\mu\leq(n+1)/2$. Since that $Z$ is a
solution to the principal part of (4.6) and that the coefficients of the
rescaled parabolic system (4.6) is controlled under the $C^{2,\alpha}$
harmonic coordinates, we further have

$\displaystyle\left|\frac{\partial z(x,s-\tau;y)}{\partial
s}\right|+\sum_{c,d=1}^{n}\left|\frac{\partial z(x,s-\tau;y)}{\partial
x^{c}\partial
x^{d}}\right|\leq\frac{(Q-1)C(n,\alpha)}{(s-\tau)^{\mu}|x-y|^{n+2-2\mu}}\exp\left(-\frac{\lambda_{0}|x-y|^{2}}{4(s-\tau)}\right)$

for any $\lambda_{0}<Q^{-1}$ and $0\leq\mu\leq(n+2)/2$. Therefore, under the
rescaled $C^{2,\alpha}$ coordinate chart, we have

(4.15)
$\displaystyle|LZ(x,s-\tau;y)|\leq\frac{(Q-1)C(n,\alpha)}{(s-\tau)^{\mu}|x-y|^{n+2-2\mu-\alpha}}\exp\left(-\frac{\lambda_{0}|x-y|^{2}}{4(s-\tau)}\right)$

for any $\lambda_{0}<Q^{-1}$ and $0\leq\mu\leq(n+2-\alpha)/2$. Hence, the
singularity $s=\tau,x=y$ is integrable when
$1-\alpha/2\leq\mu\leq(n+2-\alpha)/2$. Recall the integral formula in [17,
Lemma 3], that is, for fixed constants $\alpha$ and $\beta$,

$\displaystyle\int_{0}^{s}\int_{\Omega}(s-\tau)^{-\alpha}e^{-\frac{\lambda_{0}|x-\xi|^{2}}{4(s-\tau)}}\tau^{-\beta}e^{-\frac{\lambda_{0}|\xi-y|^{2}}{4\tau}}d\xi
d\tau$ (4.16) $\displaystyle\leq$
$\displaystyle\,\left(\frac{4\pi}{\lambda_{0}}\right)^{n/2}\frac{\Gamma(n/2-\alpha+1)\Gamma(n/2-\beta+1)}{\Gamma(n-\alpha-\beta+2)}s^{n/2+1-\alpha-\beta}e^{-\frac{\lambda_{0}|x-y|^{2}}{4s}},$

where $\Gamma(\cdot)$ denotes the gamma function. Using (4.15) with
$\mu=(n+2-\alpha)/2$ and applying (4.3), we find

$\displaystyle|(LZ)_{2}(x,s-\tau;y)|\leq\frac{(Q-1)C(n,\alpha)}{(s-\tau)^{(n+2-2\alpha)/2}}\exp\left(-\frac{\lambda_{0}|x-y|^{2}}{4(s-\tau)}\right)$

Proceeding by induction, we further obtain

$\displaystyle|(LZ)_{i}(x,s-\tau;y)|\leq\frac{(Q-1)C(n,\alpha)}{\Gamma(i\alpha)(s-\tau)^{(n+2-i\alpha)/2}}\exp\left(-\frac{\lambda_{0}|x-y|^{2}}{4(s-\tau)}\right)\,,$

where $\lambda_{0}<Q^{-1}$. Note that the coefficient $C(n,\alpha)$ is
independent of $i$. It follows that the series expansion of $\Phi$ is
convergent and satisfies (4.14). Furthermore,

$|\Phi(x,s;y)|\leq\frac{(Q-1)C(n,\alpha)}{s^{(n+2-\alpha)/2}}e^{-\frac{|x-y|^{2}}{8s}}$

by setting $\lambda_{0}=1/\sqrt{2}$.

The Hölder continuity of $\Phi(x,s;y)$ in $x$ is the same as the proof in [17,
Section 4 in Chapter 9 and Section 4 in Chapter 1] (see (1.4.17) and
(9.4.17)). Therefore, we omit the details and have $\Gamma$ the fundament
solution of (4.6).

By (4.3) and a trivial bound of $|Z(x,s-\tau;\xi)|$ under the
$C^{2,\alpha}$-harmonic coordinate, it follows that

$\int_{0}^{s}\int_{\Omega}|Z(x,s-\tau;\xi)\Phi(\xi,\tau;y)|d\xi
d\tau\leq\frac{(Q-1)C(n,\alpha)}{s^{(n-\alpha)/2}}e^{-\frac{|x-y|^{2}}{8s}},$

and hence (4.10). In particular, we have (4.9), for $s\leq 4$,

$|\Gamma(x,s;y)|\leq\frac{D(n,\alpha)}{s^{n/2}}e^{-\frac{|x-y|^{2}}{8s}}$

for some constant $D(n,\alpha)$. Note that

$\nabla\Gamma(x,s;y)=\nabla Z(x,s;y)+\int_{0}^{s}\int_{\Omega}\nabla
Z(x,s-\tau;\xi)\Phi(\xi,\tau;y)d\xi d\tau.$

Applying (4.3) once again, it follows that

$\int_{0}^{s}\int_{\Omega}|\nabla Z(x,s-\tau;\xi)||\Phi(\xi,\tau;y)|d\xi
d\tau\leq\frac{(Q-1)C(n,\alpha)}{s^{(n+1-\alpha)/2}}e^{-\frac{|x-y|^{2}}{8s}}.$

and hence (5) holds.

Ineqaulities (4.12) and (4.13) follow from straightforward compuations that

$\frac{\partial}{\partial
x^{a}}Z(x,s;y)=\frac{1}{2s}\bar{g}_{ab}(y)(x^{b}-y^{b})Z(x,s;y)$

and

$\frac{\partial}{\partial
x^{a}}\Gamma_{E}(x,s;y)=\frac{1}{2s}\delta_{ab}(x^{b}-y^{b})\Gamma_{E}(x,s;y).$

It follows that for $(x,s)\notin P_{\frac{1}{2},\frac{1}{4}}(y)$,

$|\Gamma_{E}(x,s;y)-Z(x,s;y)|\leq(Q-1)C(n,\alpha)$

and

$|\nabla\Gamma_{E}(x,s;y)-\nabla Z(x,s;y)|\leq(Q-1)C(n,\alpha)$

and hence (4.12) and (4.13). ∎

With the interior Schauder estimates, the semigroup domination, and the
maximum principle on the heat equation $\partial_{t}u=\Delta_{M}u$, we have
Lemma 6 which shows that every solution to (4.6) that decays exponentially is
close to $\Gamma$ that we construct in Lemma 5, and hence is close to the
standard Euclidean heat kernel $\Gamma_{E}$.

######  Lemma 6.

Let $\epsilon>0$, $C(n)>0$ and $0<\alpha<1$, $1<Q<\sqrt{2}$, and $r$ be given.
Suppose $\Gamma_{TM}^{1}$ and $\Gamma_{TM}^{2}$ are fundamental solutions of
(4.6) on $\mathbb{R}^{n}\times\mathbb{R}_{+}$ satisfying the decay rate

(4.17)
$|\Gamma_{TM}^{i}(x,s;y)|\leq\frac{C(n)}{s^{n/2}}e^{-\frac{|x-y|^{2}}{8t}},\quad\mbox{for
$i=1,2$ and $0<s\leq 4$}.$

There exists $\bar{R}=\bar{R}(n,\alpha,Q,C(n),\epsilon)$ so that for
$R>\bar{R}$,

(4.18)
$|\Gamma_{TM}^{1}(\cdot,\cdot;y)-\Gamma_{TM}^{2}(\cdot,\cdot;y)|_{C^{2,\alpha}(P_{2,4}(0))}\leq\epsilon$

for all $y\in B_{R}(0)$.

###### Proof.

By the Schauder type interior estimates [16, Theorem 3], there is a constant
$C=C(n,\alpha,Q)$ so that for $1\leq a,b\leq n$,

(4.19) $|(\Gamma_{TM}^{i})_{ab}(\cdot,\cdot;y)|_{C^{2,\alpha}(P_{2,4}(0))}\leq
C\left|\sum_{b=1}^{n}(\Gamma_{TM}^{i})_{ab}(\cdot,\cdot;y)\right|_{C^{0}(P_{3,4}(0))}.$

By the semigroup domination theory [22, Theorem 3.1]

(4.20) $\textrm{Tr}\left(e^{t\Delta_{TM}}\right)\leq
n\textrm{Tr}\left(e^{t\Delta_{M}}\right),$

and the fact that the associated solutions $\Gamma_{M}^{i}(x,s;y)$, $i=1,2$,
of $\partial_{s}u=\Delta_{M}u$ decay exponential
$\frac{C(n)}{s^{n/2}}\exp\big{(}-\frac{|x-y|^{2}}{8s}\big{)}$, we have the
bound (4.17).

By the maximum principle, there exists $\bar{R}=(n,\alpha,Q,C,\epsilon)$ so
that for $y\in B_{R}(0)$, $R>\bar{R}$,

$\displaystyle|\Gamma^{1}_{M}(\cdot,\cdot;y)-\Gamma^{2}_{M}(\cdot,\cdot;y)|_{C^{0}(P_{3,4}(0))}$
(4.21)
$\displaystyle\leq|\Gamma^{1}_{M}(\cdot,\cdot;y)-\Gamma^{2}_{M}(\cdot,\cdot;y)|_{C^{0}(\partial
B_{R}(0)\times[0,4])}<\epsilon.$

Combine this with (4.17), we have (4.18). ∎

By Kato’s inequality and Gaussian bounds on $K_{M}$, we have (3.9), that is,
$\bar{K}_{TM}$ has a Gaussian upper bound. Thus, by Lemmas 5 and 6,
$\bar{K}_{TM}$ is close to $\Gamma_{E}$.

##  5\. Embedding with Heat Kernels

In this section, we prove Theorems 1 and 2 that manifolds can be embedded with
heat kernels as well as how the local dilatation can be controlled by taking
the fast decay of heat kernels and the rescaling technique into account.

###### Proof.

Let $\epsilon$ be given. To determine a scale $r>0$, we need the following
three steps.

Step 1. By the Schauder interior estimates [16, Theorem 3] for parabolic
systems,

$|K_{TM}(u^{-1}(\cdot),\cdot;q)|_{C^{2+\alpha}(P_{\frac{r_{h}}{8},T}(0))}\leq
C(n)|K_{TM}(u^{-1}(\cdot),\cdot;q)|_{C^{0}(P_{\frac{r_{h}}{4},T}(0))}$

on a local domain in $\mathbb{R}^{n}\times\mathbb{R}_{+}$. By the Gaussian
bound of the heat kernel (3.9), we can select
$0<r_{0}=r_{0}(n,\kappa,i_{0},\epsilon)<\frac{r_{h}}{2}$ so that for
$0<t<2r_{0}^{2}$,

(5.1) $(2t)^{\frac{3n+2}{2}}\int_{M\setminus
B_{\frac{r_{h}}{2}}(p)}\|\nabla_{q}|K_{TM}(p,t;q)\|_{HS}^{2}|dq<\epsilon.$

Step 2. Let $R_{1}=R_{1}(n,\epsilon)$ be a radius so that for every
$\Gamma(x,s;y)$ satisfying the gradient decay

(5.2) $\displaystyle|\nabla_{y}\|\Gamma(x,s;y)\|^{2}_{HS}|$
$\displaystyle\leq\frac{D(n)}{s^{n+1/2}}\exp\left(-\frac{|x-y|^{2}}{8s}\right)\mbox{
on $(\mathbb{R}^{n}\setminus
B_{\frac{R_{1}}{\sqrt{2}}}(0))\times\mathbb{R}$},$

the integral

(5.3) $\int_{\mathbb{R}^{n}\setminus
B_{\frac{R_{1}}{\sqrt{2}}}(0)}|\nabla_{y}\|\Gamma(0,s;y)\|_{HS}^{2}|dy<\epsilon\quad\mbox{for
$\frac{1}{2}\leq s\leq 2$. }$

Step 3. Set $\alpha=\frac{1}{2}$. Choose $1<Q<\sqrt{2}$ so that
$Q-1<\epsilon$. Select $R_{2}=R_{2}(n,\alpha=\frac{1}{2},Q,\epsilon)$ as in
Lemma 6 such that the rescaled heat kernel $\bar{K}_{TM}(0,\frac{1}{2};y)$ and
the solution $\Gamma(0,\frac{1}{2};y)$ obtained in Lemma 5 are close for all
$y\in B_{R_{2}}(0)$.

We now set

(5.4)
$r_{1}:=\min\left(r_{0},\frac{r_{h}}{2R_{1}},\frac{r_{h}}{2R_{2}}\right)$

and choose a scale $r<r_{1}$ and $t=\frac{r^{2}}{2}$.

Let $\omega$ denote the modulus of continuity of
$\nabla\|\Gamma_{E}(0,s;\cdot)\|^{2}_{HS}$ . Let $\delta$ be smaller than
$r\omega^{-1}(\epsilon)/2$ and $\\{q_{i}\\}_{i=1}^{N_{0}}$ be a $\delta$-net.
Let $\\{A_{i}\\}_{i=1}^{N_{0}}$ be a partition of the manifold $M$ so that
$A_{i}\subset B_{\delta}(q_{i})$.

Let $p\in M$. With the chosen $\alpha=\frac{1}{2}$ and $Q$, there exists a
$C^{2,\alpha}$-harmonic coordinate chart
$u:B_{r_{h}}(p)\rightarrow\mathbb{R}^{n}$ with $u(p)=0$. Let
$v=\sum_{a=1}^{n}v^{a}\partial_{a}\in T_{p}M$, $|v|=1$ so that
$\hat{v}=(v^{1},\cdots,v^{n})$ is the coordinate of $v$ under the harmonic
coordinate chart. We next show that for any $t<2r_{0}^{2}$, where
$r_{0}=r_{0}(n,\kappa,i_{0},\epsilon)$, there exists a constant
$\delta=\delta(n,\kappa,i_{0},\epsilon,t)$ so that for every $\delta$-net,
$|dH_{p}(v)|^{2}$ is close to $1$.

Let $I_{\rho}(p)$ denote the subset of $\\{1,\cdots,N_{0}\\}$ such that
$A_{j}\cap B_{\rho}(p)\neq\emptyset$ and $y_{i}:=r^{-1}u(q_{i})$. Then, since
$\bar{K}_{TM}(x,s;y):=r^{n}K_{TM}(u^{-1}(xr),sr^{2};u^{-1}(yr))$,

$\displaystyle|dH_{p}(v)|^{2}$
$\displaystyle=\frac{(2t)^{\frac{3n+2}{2}}}{V^{2}_{e}}\sum_{i=1}^{N_{0}}|A_{i}|(\nabla\|K_{TM}(p,t;q_{i})\|^{2}_{HS}\cdot
v)^{2}$ $\displaystyle=\frac{1}{V^{2}_{e}}\sum_{i\in
I_{1}}\frac{|A_{i}|}{r^{n}}(\nabla\|\bar{K}_{TM}(0,\frac{1}{2};y_{i})\|^{2}_{HS}\cdot\hat{v})^{2}$
$\displaystyle\quad+\frac{1}{V^{2}_{e}}\sum_{i\in
I_{2}}\frac{|A_{i}|}{r^{n}}(\nabla\|\bar{K}_{TM}(0,\frac{1}{2};y_{i})\|^{2}_{HS}\cdot\hat{v})^{2}$
$\displaystyle\quad+\frac{(2t)^{\frac{3n+2}{2}}}{V^{2}_{e}}\sum_{i\in
I_{3}}|A_{i}|(\nabla\|K_{TM}(p,t;q_{i})\|^{2}_{HS}\cdot v)^{2}$
$\displaystyle=\mbox{ (I) + (II) + (III)}$

where $I_{1}=I_{rR_{1}}(p)$, $I_{2}=I_{r_{h}/2}(p)\setminus I_{1}$, and
$I_{3}=\\{1,\cdots,N_{0}\\}\setminus I_{1}\cup I_{2}$. Here, we write the
first two terms in terms of the rescaled coordinates.

First, note that $\cup_{i\in I_{3}}A_{i}$ lies inside $M\setminus
B_{\frac{r_{h}}{2}}(p)$. Because of equation (5.1), (III) is controlled by
$\epsilon$, that is,

(5.5) $(2t)^{\frac{3n+2}{2}}\sum_{i\in
I_{3}}|A_{i}|(\nabla\|K_{TM}(p,t;q_{i})\|^{2}_{HS}\cdot v)^{2}<C(n)\epsilon.$

Second, since $\mathrm{diam}(A_{i})<r\omega(\epsilon)$,

$\displaystyle\left|\sum_{i\in
I_{2}}\frac{|A_{i}|}{r^{n}}(\nabla\|\bar{K}_{TM}(0,\frac{1}{2};y_{i})\|^{2}_{HS}\cdot\hat{v})^{2}-\sum_{i\in
I_{2}}\int_{\frac{u(A_{i})}{r}}(\nabla\|\bar{K}_{TM}(0,\frac{1}{2};y)\|^{2}_{HS}\cdot\hat{v})^{2}d\mu\right|<C(n)\epsilon$

where $d\mu$ denotes the push-forward of the standard volume measure under
$r^{-1}u$.

Note that $\cup_{i\in I_{2}}A_{i}$ lies outside of $B_{rR_{1}}(p)$ and inside
$B_{r_{h}}(p)$. Thus, $\cup_{i\in I_{2}}r^{-1}u(A_{i})$ is a subset of
$\mathbb{R}^{n}\setminus B_{\frac{R_{1}}{\sqrt{2}}}(0)$. By Lemmas 5 and 6,
$\bar{K}_{TM}(x,\frac{1}{2};y)$ satisfies the gradient decay (5.2). By the
choice of $R_{1}$ in (5.3), we have that

$\sum_{i\in
I_{2}}\int_{\frac{u(A_{i})}{r}}(\nabla\|\bar{K}_{TM}(0,\frac{1}{2};y)\|^{2}_{HS}\cdot\hat{v})^{2}d\mu<\epsilon$

and that (II) is controlled by $\epsilon$ as well,

(5.6) $\sum_{i\in
I_{2}}\frac{|A_{i}|}{r^{n}}(\nabla\|\bar{K}_{TM}(0,\frac{1}{2};y_{i})\|^{2}_{HS}\cdot\hat{v})^{2}<C(n)\epsilon.$

Last, we show that (I) is close to $1$. Let

(5.7) $\displaystyle R:=\max(R_{1},R_{2}).$

Note that $B_{R}(0)\subset u(B_{r_{h}}(p))/r$. By Lemma 6 and our choice of
$Q$ and $R_{2}$, for all $y\in B_{R}(0)$, $\bar{K}_{TM}(0,\frac{1}{2};y)$ is
close to the solution $\Gamma(0,\frac{1}{2};y)$ constructed in Lemma 5, which
is close to the Euclidean kernel $\Gamma_{E}(0,\frac{1}{2};y)$. Specifically,

$\left|\nabla\|\bar{K}_{TM}(0,\frac{1}{2},y)\|^{2}_{HS}-\nabla\|\Gamma_{E}(0,\frac{1}{2},y)\|^{2}_{HS}\right|<C(n)\epsilon.$

This implies

$\displaystyle\left|\sum_{i\in
I_{1}}\left((\hat{v}\cdot\nabla\|\bar{K}_{TM}(0,\frac{1}{2};y_{i})\|^{2}_{HS})^{2}-(\hat{v}\cdot\nabla\|\Gamma_{E}(0,\frac{1}{2};y_{i})\|^{2}_{HS})^{2}\right)\frac{|A_{i}|}{r^{n}}\right|$
$\displaystyle<C(n)\epsilon\sum_{i\in I_{1}}\frac{|A_{i}|}{r^{n}}$ (5.8)
$\displaystyle<C(n)\epsilon|B_{R_{1}}(0)|.$

Since $\mathrm{diam}(A_{i})<r\omega(\epsilon)$, we have

(5.9) $\displaystyle\left|\sum_{i\in
I_{1}}(\nabla_{\hat{v}}\|\Gamma_{E}(0,\frac{1}{2};y)\|^{2}_{HS})^{2}\frac{|A_{i}|}{r^{n}}-\int_{\frac{u(A_{i})}{r}}(\nabla_{\hat{v}}\|\Gamma_{E}(0,\frac{1}{2};y)\|^{2}_{HS})^{2}d\mu\right|<C(n)\epsilon,$

where $d\mu$ denotes the push-forward of the standard volume form under
$r^{-1}u$. Due to that $Q^{-1}<g<Q$ and symmetry of
$\Gamma_{E}(x,t;y)=\frac{1}{(4\pi t)^{n/2}}e^{-\frac{|x-y|^{2}}{4t}}I_{n}$, we
have

$\displaystyle\left|\sum_{i\in
I_{1}}\left(\int_{\frac{u(A_{i})}{r}}(\nabla_{\hat{v}}\|\Gamma_{E}(0,\frac{1}{2};y)\|^{2}_{HS})^{2}d\mu-\int_{\frac{u(A_{i})}{r}}(\nabla_{\hat{v}}\|\Gamma_{E}(0,\frac{1}{2};y)\|^{2}_{HS})^{2}dy\right)\right|$
(5.10) $\displaystyle=\left|\sum_{i\in
I_{1}}\left(\int_{\frac{u(A_{i})}{r}}(\nabla_{\hat{v}}\|\Gamma_{E}(0,\frac{1}{2};y)\|^{2}_{HS})^{2}d\mu-\int_{\frac{u(A_{i})}{r}}(\partial_{x_{1}}\|\Gamma_{E}(0,\frac{1}{2};y)\|^{2}_{HS})^{2}dy\right)\right|$
$\displaystyle<C(n)(Q-1)\epsilon<C(n)\epsilon.$

Furthermore, by the choice of $R_{1}$ and (5.3),

(5.11) $\left|\sum_{i\in
I_{1}}\int_{\frac{u(A_{i})}{r}}(\partial_{x_{1}}\|\Gamma_{E}(0,\frac{1}{2};y)\|^{2}_{HS})^{2}dy-
V_{e}^{2}\right|<\epsilon$

Combining inequalities (5), (5.9), (5), and (5.11), we have that (I) is close
to $1$. Together with (5.5) and (5.6), we conclude that

$|dH_{p}(v)^{2}-1|\leq C(n)\epsilon$

that is, (3.2) in Theorem 1.

As for the second part, (3.4) of Theorem 1, by Lemma 10 in [25], it follows
that there exists an $N=N(n,\kappa,i_{0},V,\epsilon,t_{0})$ so that for any
$m\geq N$ and $0<t\leq t_{0}$,

(5.12)
$\displaystyle\|\|K_{TM}^{(m)}(\cdot,t;q)\|^{2}_{HS}-\|K_{TM}(\cdot,t;q)\|^{2}_{HS}\|_{\infty}$
$\displaystyle<\epsilon$ (5.13)
$\displaystyle\|\nabla\|K_{TM}^{(m)}(\cdot,t;q)\|^{2}_{HS}-\nabla\|K_{TM}(\cdot,t;q)\|^{2}_{HS}\|_{\infty}$
$\displaystyle<\epsilon.$

Therefore,

$1-\epsilon<\left|dH_{p}^{(m)}(v)|^{2}\right|<1+\epsilon.$

∎

Next, we prove Theorem 2

###### Proof.

For any $A>0$, we simply replace for each fixed $i$ the point $q_{i}$ in
Theorem 1 by points $q_{i}^{j},j=1,\cdots,N_{i}$, where
$N_{i}=\lceil|A_{i}|/A\rceil$. Next, rename all the points $q_{i}^{j}$ to
$p_{k}$. It follows that when $A=A(n,\kappa,i_{0},\epsilon,t,V)$ small enough,

$1-\epsilon<|dH_{p}(v)|<1+\epsilon.$

∎

## References

  * [1] A. S. Bandeira, A. Singer, and D. A. Spielman.  A Cheeger Inequality for the Graph Connection Laplacian.  SIAM Journal on Matrix Analysis and Applications, to appear, 2013\.  arXiv:1204.3873 [math.SP].
  * [2] J. Bates.  The embedding dimension of laplacian eigenfunction maps.  Appl. Comput. Harmon. Anal., 37(3):516–530, 2014.
  * [3] M. Belkin and P. Niyogi.  Laplacian Eigenmaps for Dimensionality Reduction and Data Representation.  Neural. Comput., 15(6):1373–1396, June 2003.
  * [4] M. Belkin and P. Niyogi.  Convergence of laplacian eigenmaps.  In Advances in Neural Information Processing Systems 19: Proceedings of the 2006 Conference, volume 19, page 129. The MIT Press, 2007\.
  * [5] P. Bérard, G. Besson, and S. Gallot.  Embedding riemannian manifolds by their heat kernel.  Geom. Funct. Anal., 4:373–398, 1994.  10.1007/BF01896401.
  * [6] Jeff Calder and Nicolás García Trillos.  Improved spectral convergence rates for graph laplacians on epsilon-graphs and k-nn graphs.  ArXiv, abs/1910.13476, 2019.
  * [7] F. Chung and M. Kempton.  A local clustering algorithm for connection graphs.  In Anthony Bonato, Michael Mitzenmacher, and Pawel Pralat, editors, Algorithms and Models for the Web Graph, volume 8305 of Lecture Notes in Computer Science, pages 26–43. Springer International Publishing, 2013\.
  * [8] F Chung, W Zhao, and M Kempton.  Ranking and sparsifying a connection graph.  In Anthony Bonato and Jeannette Janssen, editors, Algorithms and Models for the Web Graph, volume 7323 of Lecture Notes in Computer Science, pages 66–77. Springer Berlin Heidelberg, 2012.
  * [9] F.R.K. Chung, F.C. Graham, CBMS Conference on Recent Advances in Spectral Graph Theory, Conference Board of the Mathematical Sciences (U.S.). Conference on Recent Advances in Spectral Graph Theory, Calif.. Conference on Recent Advances in Spectral Graph Theory. 1994, Fresno, Conference on recent advances in spectral graph theory held at California state university 1994 ; Fresno, American Mathematical Society, and Conference Board of the Mathematical Sciences.  Spectral Graph Theory.  CBMS Regional Conference Series. Conference Board of the mathematical sciences, 1997.
  * [10] R. R. Coifman and S. Lafon.  Diffusion maps.  Appl. Comput. Harmon. Anal., 21(1):5–30, 2006.
  * [11] David B Dunson, Hau-Tieng Wu, and Nan Wu.  Diffusion based gaussian process regression via heat kernel reconstruction.  arXiv preprint arXiv:1912.05680, 2019.
  * [12] N El Karoui and H.-T. Wu.  Connection graph Laplacian and random matrices with random blocks.  Information and Inference: A Journal of the IMA, 4:1–42, 2015.
  * [13] Noureddine El Karoui, Hau-Tieng Wu, et al.  Graph connection laplacian methods can be made robust to noise.  The Annals of Statistics, 44(1):346–372, 2016.
  * [14] Justin Eldridge, Mikhail Belkin, and Yusu Wang.  Unperturbed: spectral analysis beyond davis-kahan.  06 2017.
  * [15] Charles Fefferman, Sergei Ivanov, Yaroslav Kurylev, Matti Lassas, and Hariharan Narayanan.  Fitting a putative manifold to noisy data.  In Conference On Learning Theory, pages 688–720, 2018.
  * [16] A. Friedman.  Interior estimates for parabolic systems of partial differential equations.  J. of Mathematics and Mechanics, 7(3):393–417, 1958.
  * [17] A. Friedman.  Partial Differential Equations of Parabolic Type.  Prentice-Hall, 1964.
  * [18] Grior’yan.  Heat kernel and analysis on manifolds, volume 47 of AMS/IP Stidies in Advanced Mathematics.  American Mathematical Society, Providence, R.I., 2009.
  * [19] R. Hadani and A. Singer.  Representation theoretic patterns in three dimensional Cryo-Electron Microscopy I: The intrinsic reconstitution algorithm.  Annals of Mathematics, 174(2):1219–1241, 2011.
  * [20] E. Hebey and M. Herzlich.  Harmonic coordinates, harmonic radius and convergence of riemannian manifolds.  Rend. Mat. Appl., (7)17(4):569 – 605, 1998.
  * [21] Matthias Hein and Markus Maier.  Manifold denoising.  In Advances in neural information processing systems, pages 561–568, 2007.
  * [22] H. Hess, R Schrader, and D. A. Uhlenbrock.  Kato’s inequality and the spectral distribution of laplacians on compact riemannian manifolds.  J. Differential Geometry, pages 27–37, 1980.
  * [23] P. W. Jones, M. Maggioni, and R. Schul.  Manifold parametrizations by eigenfunctions of the Laplacian and heat kernels.  P. Natl. Acad. Sci. USA, 105(6):1803–8, February 2008.
  * [24] Chen-Yun Lin, Arin Minasian, Xin Jessica Qi, and Hau-Tieng Wu.  Manifold learning via the principle bundle approach.  Frontiers in Applied Mathematics and Statistics, 4:21, 2018.
  * [25] Chen-Yun Lin and Hau-Tieng Wu.  Embeddings of riemannian manifolds with finite eigenvector fields of connection laplacian.  Calculus of Variations and Partial Differential Equations, 57, 04 2018.
  * [26] S Marchesini, Y.-C. Tu, and H.-T. Wu.  Alternating Projection, Ptychographic Imaging and Phase Synchronization.  Appl. Comput. Harmon. Anal., 41(3):815–851, 2016.
  * [27] J. Nash.  The Imbedding Problem for Riemannian Manifolds.  Annals of Mathematics, 63(1):20–63, 1956.
  * [28] J. W. Portegies.  Embeddings of Riemannian manifolds with heat kernels and eigenfunctions.  Comm. Pure Appl. Math., 2015.
  * [29] S. T. Roweis and L. K. Saul.  Nonlinear dimensionality reduction by locally linear embedding.  Science, 290(5500):2323–2326, 2000.
  * [30] A. Singer and H.-T. Wu.  Orientability and diffusion map.  Appl. Comput. Harmon. Anal., 31(1):44–58, 2011.
  * [31] A Singer and H.-T. Wu.  Vector Diffusion Maps and the Connection Laplacian.  Comm. Pure Appl. Math., 65(8):1067–1144, 2012.
  * [32] A Singer and H.-T. Wu.  Spectral Convergence of the Connection Laplacian from random samples.  Information and Inference: A Journal of the IMA, 6(1):58–123, 2017\.
  * [33] A. Singer, Zhao Z., Shkolnisky Y., and Hadani R.  Viewing angle classification of cryo-electron microscopy images using eigenvectors.  SIAM J. Imaging Sci., 4(2):723–759, 2011.
  * [34] J. B. Tenenbaum, V. de Silva, and J. C. Langford.  A Global Geometric Framework for Nonlinear Dimensionality Reduction.  Science, 290(5500):2319–2323, 2000.
  * [35] Nicolas Trillos, Moritz Gerlach, Matthias Hein, and Dejan Slepcev.  Error estimates for spectral convergence of the graph laplacian on random geometric graphs toward the laplace–beltrami operator.  Foundations of Computational Mathematics, 01 2018.
  * [36] Hassler Whitney.  The Singularities of a Smooth n-Manifold in (2n - 1)-Space.  The Annals of Mathematics, 45(2):247, 1944.
  * [37] Z. Zhao and A. Singer.  Rotationally Invariant Image Representation for Viewing Direction Classification in Cryo-EM.  Journal of Structural Biology, 186(1):153–166, 2014.


# Viscous fluid cosmology in symmetric teleparallel gravity

Raja Solanki 0000-0001-8849-7688<EMAIL_ADDRESS>Department of
Mathematics, Birla Institute of Technology and Science-Pilani,
Hyderabad Campus, Hyderabad-500078, India. Dheeraj Singh Rana
0000-0002-4401-8814<EMAIL_ADDRESS>Department of Mathematics, Birla
Institute of Technology and Science-Pilani,
Hyderabad Campus, Hyderabad-500078, India. Sanjay Mandal 0000-0003-2570-2335
<EMAIL_ADDRESS>Department of Mathematics, Birla Institute of
Technology and Science-Pilani,
Hyderabad Campus, Hyderabad-500078, India. P.K. Sahoo 0000-0003-2130-8832
<EMAIL_ADDRESS>Department of Mathematics, Birla Institute
of Technology and Science-Pilani,
Hyderabad Campus, Hyderabad-500078, India.

###### Abstract

In this manuscript, we analyze the viscous fluid cosmological model in the
framework of recently proposed $f(Q)$ gravity by assuming three different
forms of bulk viscosity coefficients, specifically,
$(i)\zeta=\zeta_{0}+\zeta_{1}\left(\frac{\dot{a}}{a}\right)+\zeta_{2}\left(\frac{{\ddot{a}}}{\dot{a}}\right)$,
$(ii)\zeta=\zeta_{0}+\zeta_{1}\left(\frac{\dot{a}}{a}\right)$, and
$(iii)\zeta=\zeta_{0}$ and a linear $f(Q)$ model, particularly, $f(Q)=\alpha
Q$ where $\alpha\neq 0$ is free model parameter. We estimate the bulk
viscosity coefficients and the model parameter values using the combined
H(z)+Pantheon+BAO data set. We study the asymptotic behavior of our
cosmological bulk viscous model by utilizing the phase space method. We find
that corresponding to all three cases, our model depicts the evolution of the
universe from matter dominated decelerated epoch (a past attractor) to a
stable de-sitter accelerated epoch (a future attractor). Further, we study the
physical behavior of effective pressure, effective equation of state (EoS),
and the statefinder parameters. We find that the pressure component in the
presence of bulk viscosity shows negative behavior and the effective EoS
parameter predicts the accelerated expansion phase of the universe for all
three cases. Moreover, we obtain that the trajectories of our model lie in the
quintessence region and it converges to the $\Lambda$CDM fixed point in the
far future. We find that the accelerated deSitter like phase comes purely from
the $\bar{\zeta}_{0}$ case without any geometrical modification to GR.
Moreover, we find that the late-time behavior of all three cases of viscosity
coefficients are identical. Further, we consider a non-linear $f(Q)$ model,
specifically, $f(Q)=-Q+\beta Q^{2}$ and then we analyze the behavior of model
using dynamical approach. We find that the late-time behavior of the
considered non-linear model $f(Q)=-Q+\beta Q^{2}$ with $\beta\leq 0$ is
similar to the linear case, whereas for the case $\beta>0$ results are quite
different.

Keywords: bulk viscosity, symmetric teleparallel gravity, phase space
analysis, and statefinder parameters.

##  I Introduction

The discovery of the accelerating behavior of the cosmos is one of the
challenging issues. Several cosmic observations such as Supernovae searches
[1, 2], BAO [3, 4], WMAP [5, 6], and CMBR [7, 8] are the evidence for this
accelerating behavior. The best fit model corresponding to these observational
aspects is the $\Lambda$CDM model, where $\Lambda$ is the dark energy
component responsible for the aforementioned acceleration. Despite its
observational compatibility, the $\Lambda$CDM model has some well known
delicate issues, namely cosmic coincidence and cosmological constant problem
[9]. A way to bypass the above mentioned dark energy issue is to modify the
geometry of the standard general relativity (GR) by proposing a more generic
action to describe gravitational interactions. This approach can be seen
within the references [10, 11, 12] and is commonly known as modified theories
of gravity.

In this work, we are going to present our analysis and outcomes under the
framework of recently proposed $f(Q)$ gravity [13]. The standard GR
corresponds to the spacetime having non-zero curvature tensor while vanishing
non-metricity and torsion, the teleparallel equivalent of GR (TEGR)
corresponds to the spacetime having non-zero torsion tensor while vanishing
non-metricity and curvature, and the symmetric teleparallel equivalent of GR
(STEGR) corresponds to the spacetime having non-zero non-metricity tensor
while vanishing torsion and curvature. The $f(R)$, $f(T)$, and $f(Q)$ gravity
are straightforward generalizations of GR, TEGR, and STEGR respectively [14].
One notable advantage of $f(Q)$ gravity is that its field equations are of
second order in scale factor so that one can assume the Lagrangian of higher
order without having any issues. Another advantage of $f(Q)$ gravity is that
it holds Bianchi’s identity automatically, while in the case of $f(T)$
gravity, the presence of the anti-symmetric part becomes hindered.
Furthermore, some generalizations of $f(Q)$ theory have been proposed in the
literature which is commonly known as $f(Q,T)$ theory [15] and Weyl type
$f(Q,T)$ theory [16]. Recently, several interesting cosmological and
astrophysical implications of $f(Q)$ gravity has been appeared, for instance,
Cosmography [17], Energy conditions [18], Gravitational waves [19, 20, 21],
Black hole solution [22], Quantum cosmology [23], General covariant symmetric
teleparallel cosmology [24], Wormhole solution [25], Evidence that non-
metricity $f(Q)$ gravity can challenge $\Lambda$CDM [26], Bouncing Cosmology
[27], Hamiltonian analysis and ADM formulation [28], Statefinder Analysis of
$f(Q)$ cosmology [29], and the Cosmological perturbations [30].

Earlier, to study the early inflationary phase of cosmic evolution without
invoking any dark energy component, viscosity in the cosmic fluid has been
considered. Eckart proposed the non-casual theory of viscosity by assuming the
first order deviations from the equilibrium [31]. Later, Israel and Stewart
proposed the casual theory of viscosity by assuming the second order
deviations from the equilibrium [32, 33, 34]. To investigate the late time
cosmic scenario, the casual theory has been used. When we consider the
deviations of the second order, there are two different coefficients of
viscosity that comes into the picture, namely shear and bulk viscosity. The
velocity gradient related to the shear viscosity vanishes in the homogeneous
universe. Therefore in the case of a homogeneous and isotropic FLRW universe,
only the bulk viscosity coefficient comes into play. When the cosmic fluid
expands with the expansion of the universe, then the measure of pressure
required to recover the thermal stability can be regarded as the bulk
viscosity. Basically, the geometrical modification in the Einstein-Hilbert
action of GR justifies the cosmic expansion, whereas the viscosity
coefficients contribute to the pressure term to drive the cosmic acceleration.
Recently, the effect of bulk viscosity in cosmic evolution has been
intensively investigated, for instance, one can check the references [35, 36,
37, 38, 39, 40, 41, 42].

In this work, we are going to analyze the role of different bulk viscosity
coefficients in cosmic evolution by considering $f(Q)$ gravity. The paper is
organized as follows. In Sec II, we present the geometrical aspects of non-
metricity $f(Q)$ gravity. In Sec III, we present the flat FLRW Universe in
symmetric teleparallel cosmology with bulk viscous matter. In Sec IV, we
derive the expressions for the Hubble parameter by assuming three different
bulk viscosity coefficients and a linear $f(Q)=\alpha Q$ model with
$\alpha\neq 0$. We present the best fit values of the viscosity coefficients
and the model parameter by using the combined H(z)+Pantheon+BAO data set. Then
we analyze the asymptotic behavior of our model using the dynamical system
approach. Further, we investigate the behavior of effective pressure,
effective EoS parameter, and the statefinder parameters. In Sec V, we consider
a non-linear $f(Q)$ model, specifically, $f(Q)=Q+\beta Q^{2}$ and then
investigated its dynamical behavior. In the last section VI, we present the
outcomes of our investigation.

##  II Geometry with Non-metricity

It is well known that gravitational interactions in the standard theory of
gravity (GR) ruled by the curvature tensor,

$R^{\alpha}_{\>\beta\mu\nu}=2\partial_{[\mu}X^{\alpha}_{\>\nu]\beta}+2X^{\alpha}_{\>[\mu\mid\lambda\mid}X^{\lambda}_{\nu]\beta}$
(1)

Here the affine connection
$X^{\alpha}_{\>\mu\nu}=\Gamma^{\alpha}_{\>\beta\gamma}$ represents torsion-
free metric compatible Levi-Civita connection that is characterized by the
metric tensor. We start by recalling that the tool by which gravity is
mediated is not the physical manifold, but the generic affine connection [43].
The generic affine connection can be described through other properties, such
as non-metricity. According to the strong equivalence principle [44],
corresponding to every point on the physical manifold, there is well defined
tangent space and this affine connection act as a mediator between neighboring
tangent spaces to define the derivative operators. This implies that the
generic affine connection follows a decomposition [43],

$X^{\alpha}_{\ \mu\nu}=\Gamma^{\alpha}_{\ \mu\nu}+K^{\alpha}_{\
\mu\nu}+L^{\alpha}_{\ \mu\nu},$ (2)

where

$\Gamma^{\alpha}_{\
\mu\nu}\equiv\frac{1}{2}g^{\alpha\lambda}(g_{\mu\lambda,\nu}+g_{\lambda\nu,\mu}-g_{\mu\nu,\lambda})$
(3)

denotes the Levi-Civita connection of the metric tensor $g_{\mu\nu}$,

$K^{\alpha}_{\ \mu\nu}\equiv\frac{1}{2}(T^{\alpha}_{\ \mu\nu}+T_{\mu\ \nu}^{\
\alpha}+T_{\nu\ \mu}^{\ \alpha})$ (4)

is the contortion tensor, and

$L^{\alpha}_{\ \mu\nu}\equiv\frac{1}{2}(Q^{\alpha}_{\ \mu\nu}-Q_{\mu\ \nu}^{\
\alpha}-Q_{\nu\ \mu}^{\ \alpha})$ (5)

represents the distortion tensor. Here we define the torsion tensor
$T^{\alpha}_{\ \mu\nu}$ and the non-metricity tensor $Q_{\alpha\mu\nu}$ as

$T^{\alpha}_{\ \mu\nu}\equiv X^{\alpha}_{\ \mu\nu}-X^{\alpha}_{\ \nu\mu}$ (6)

and

$Q_{\alpha\mu\nu}\equiv\nabla_{\alpha}g_{\mu\nu}$ (7)

We impose the symmetric teleparallelism condition
$R^{\alpha}_{\>\beta\mu\nu}=0$ and $K^{\alpha}_{\ \mu\nu}=0$ to develop $f(Q)$
theory. The condition of symmetric teleparallelism makes the generic affine
connection to be total inertial. In an arbitrary coordinate system, this
connection can be parameterized as [30]

$X^{\alpha}\,_{\mu\beta}=\frac{\partial
x^{\alpha}}{\partial\xi^{\rho}}\partial_{\mu}\partial_{\beta}\xi^{\rho}.$ (8)

We can always choose a coordinate transformation so that $X^{\alpha}_{\
\mu\nu}=0$. This special coordinate frame is known as coincident gauge.

The action corresponding to $f(Q)$ gravity under the geometrical framework
that incorporates the connection to be flat and symmetric reads as [45]

$S=\int{\frac{1}{2}f(Q)\sqrt{-g}d^{4}x}+\int{L_{m}\sqrt{-g}d^{4}x}$ (9)

Here $f(Q)$ represents the function of the non-metricity scalar $Q$, $L_{m}$
is the matter Lagrangian, and $g=det(g_{\mu\nu})$. From the symmetry of the
metric tensor, it follows that one can construct two different non-metricity
vectors

$Q_{\alpha}=Q_{\alpha}\>^{\mu}\>{}_{\mu}\>and\>\tilde{Q}_{\alpha}=Q^{\mu}\>_{\alpha\mu}$
(10)

In addition, the superpotential tensor is defined by

$4P^{\lambda}\>_{\mu\nu}=-Q^{\lambda}\>_{\mu\nu}+2Q_{(\mu}\>^{\lambda}\>{}_{\nu)}+(Q^{\lambda}-\tilde{Q}^{\lambda})g_{\mu\nu}-\delta^{\lambda}_{(\mu}Q_{\nu)}.$
(11)

One can acquire the non-metricity scalar by using superpotential tensor as
[46]

$Q=-Q_{\lambda\mu\nu}P^{\lambda\mu\nu}.$ (12)

Further, the energy momentum tensor characterizes the matter-energy content of
the universe reads as

$\mathcal{T}_{\mu\nu}=\frac{-2}{\sqrt{-g}}\frac{\delta(\sqrt{-g}L_{m})}{\delta
g^{\mu\nu}}$ (13)

By altering the action (9) with respect to the metric, we obtained the field
equation

$\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})=-T_{\mu\nu}.$
(14)

In the absence of hypermomentum, we obtained the following field equation by
altering the action (9) with respect to the connection,

$\nabla_{\mu}\nabla_{\nu}(\sqrt{-g}f_{Q}P^{\mu\nu}\>_{\lambda})=0$ (15)

Moreover, it follows from the Bianchi identity that this equation is fulfilled
automatically once the metric equations of motion are [47].

##  III Flat FLRW Universe in Symmetric Teleparallel Cosmology with Bulk
Viscous Matter

We begin with following spatially flat FLRW line element [48] in Cartesian
coordinates, which is, as a matter of fact also a coincident gauge
coordinates, therefore from now connection becomes trivial and metric is only
a fundamental variable,

$ds^{2}=-dt^{2}+a^{2}(t)[dx^{2}+dy^{2}+dz^{2}]$ (16)

where, $a(t)$ is the scale factor that measures the cosmic expansion. However,
the most general connection in cosmological settings was derived in [24, 49],
which leads to a non-trivial contribution to the field equations. Now, for the
line element (16) we obtain the non-metricity scalar $Q$ as

$Q=6H^{2}$ (17)

The energy momentum tensor that corresponds to the universe with bulk viscous
matter for the line element (16) reads as,

$\mathcal{T}_{\mu\nu}=(\rho+\bar{p})u_{\mu}u_{\nu}+\bar{p}g_{\mu\nu}$ (18)

Here $\rho$ represents matter-energy density, $\bar{p}=p-3\zeta H$ is the
effective pressure having coefficient of bulk viscosity $\zeta$ and the usual
pressure $p$, and $u^{\mu}=(1,0,0,0)$ are components of the four velocities.
The field equations characterizing the universe with bulk viscous matter
are[46],

$3H^{2}=\frac{1}{2f_{Q}}\left(-\rho+\frac{f}{2}\right)$ (19)

and

$\dot{H}+3H^{2}+\frac{\dot{f_{Q}}}{f_{Q}}H=\frac{1}{2f_{Q}}\left(\bar{p}+\frac{f}{2}\right)$
(20)

In addition, we acquire the following matter conservation equation by taking
the trace of the field equations,

$\dot{\rho}+3H\left(\rho+\bar{p}\right)=0$ (21)

In further analysis we consider only the dust case, so the effective pressure
is given purely by viscous pressure $\bar{p}=-3\zeta H$.

##  IV Linear $f(Q)$ model

We consider the following $f(Q)$ function for our analysis [50, 51],

$f(Q)=\alpha Q,\ \ \ \alpha\neq 0$ (22)

Then the Friedmann equations for this specific $f(Q)$ function becomes,

$\rho=-3\alpha H^{2}$ (23)

and

$\bar{p}=2\alpha\dot{H}+3\alpha H^{2}$ (24)

In particular, for the case $\alpha=-1$, one can retrieve the usual Friedmann
equations of GR.

In fluid mechanics, it is well known that the coefficient of bulk viscosity is
associated with the rate of expansion or compression of the fluid [52]. In the
context of the cosmological model, the cosmic fluid is comoving with the
expanding universe, therefore the velocity $\dot{a}$ and acceleration
$\ddot{a}$ of the expanding universe and that of the cosmic fluid becomes
coincide. Thus it is evident that the bulk viscosity coefficient $\zeta$ is
proportional to the velocity and acceleration term. We consider the following
bulk viscosity form, which is nothing but a linear combination of velocity and
acceleration terms with a constant [53],

$\zeta=\zeta_{0}+\zeta_{1}\left(\frac{\dot{a}}{a}\right)+\zeta_{2}\left(\frac{{\ddot{a}}}{\dot{a}}\right)=\zeta_{0}+\zeta_{1}H+\zeta_{2}\left(\frac{\dot{H}}{H}+H\right)$
(25)

Now we set $\bar{\zeta_{0}}$, $\bar{\zeta_{1}}$, and $\bar{\zeta_{2}}$ as

$\frac{3\zeta_{0}}{H_{0}}=\bar{\zeta_{0}},\ \ 3\zeta_{1}=\bar{\zeta_{1}},\ \
3\zeta_{2}=\bar{\zeta_{2}}$ (26)

where $H_{0}$ is the Hubble parameter value at present time $z=0$ and this
parameters are known as dimensionless bulk viscous parameters. Then by using
(25) and (26) in equation (24) and the fact that
$\frac{d}{dt}=H\frac{d}{dln(a)}$, we have

$\frac{dH}{dln(a)}+\left(\frac{3\alpha+\bar{\zeta_{1}}+\bar{\zeta_{2}}}{2\alpha+\bar{\zeta_{2}}}\right)H+\left(\frac{\bar{\zeta_{0}}}{2\alpha+\bar{\zeta_{2}}}\right)H_{0}=0$
(27)

Now by integrating the above equation, we obtained the following expression of
Hubble parameter in terms of redshift corresponding to the bulk viscous non-
relativistic matter dominated universe

$H(z)=H_{0}\left[(1+z)^{\left(\frac{3\alpha+\bar{\zeta_{1}}+\bar{\zeta_{2}}}{2\alpha+\bar{\zeta_{2}}}\right)}\left(1+\frac{\bar{\zeta_{0}}}{3\alpha+\bar{\zeta_{1}}+\bar{\zeta_{2}}}\right)-\frac{\bar{\zeta_{0}}}{3\alpha+\bar{\zeta_{1}}+\bar{\zeta_{2}}}\right]$
(28)

Now we consider the following three different cases on dimensionless viscous
parameters which is well known in the literature,

Case I : When viscosity coefficient depends on both the velocity and
acceleration i.e. $\bar{\zeta_{0}}$, $\bar{\zeta_{1}}$, and $\bar{\zeta_{2}}$
all are non-zero.
Case II: When viscosity coefficient depends on the velocity but not on
acceleration i.e. $\bar{\zeta_{0}}$, $\bar{\zeta_{1}}$ are non-zero whereas
$\bar{\zeta_{2}}=0$. In this case, the expression for the Hubble parameter
becomes

$H(z)=H_{0}\left[(1+z)^{\left(\frac{3\alpha+\bar{\zeta_{1}}}{2\alpha}\right)}\left(1+\frac{\bar{\zeta_{0}}}{3\alpha+\bar{\zeta_{1}}}\right)-\frac{\bar{\zeta_{0}}}{3\alpha+\bar{\zeta_{1}}}\right]$
(29)

Case III: When viscosity coefficient does not depends on both the velocity and
acceleration i.e. $\bar{\zeta_{1}}=0$, $\bar{\zeta_{2}}=0$ whereas
$\bar{\zeta_{0}}$ is non-zero. In this case, the expression for the Hubble
parameter becomes

$H(z)=H_{0}\left[(1+z)^{\frac{3}{2}}\left(1+\frac{\bar{\zeta_{0}}}{3\alpha}\right)-\frac{\bar{\zeta_{0}}}{3\alpha}\right]$
(30)

In particular, when $\bar{\zeta_{0}}=\bar{\zeta_{1}}=\bar{\zeta_{2}}=0$, then
$H(z)$ reduces to $H(z)=H_{0}(1+z)^{\frac{3}{2}}$ that corresponds to the non-
viscous matter dominated universe.

###  IV.1 Parameter Estimation Using Observational Data

In this section, we consider updated H(z) data, Pantheon Supernovae data, and
the Baryonic Acoustic Oscillation (BAO) data to estimate the parameter values
corresponding to all three different cases. We apply the Markov Chain Monte
Carlo (MCMC) approach along with Bayesian analysis to explore the parameter
space of our cosmological bulk viscous model by utilizing emcee python library
[54].

####  IV.1.1 H(z) datasets

We are familiar with the fact that the Hubble parameter can directly predict
the rate of cosmic expansion. In general, there are two very popular
approaches to extracting the Hubble parameter value at definite redshift,
namely differential age and line of sight BAO technique. In this manuscript,
we work with an updated list of H(z) data points. One can check the reference
for the complete set of data points [50]. Furthermore, we have taken
$H_{0}=67.9$ Km/s/Mpc [55] for our investigation. The $\chi^{2}$ function
corresponding to H(z) data points reads as

$\chi_{H}^{2}=\sum\limits_{k=1}^{57}\frac{[H_{th}(z_{k},\theta)-H_{obs}(z_{k})]^{2}}{\sigma_{H(z_{k})}^{2}}.$
(31)

Here, $H_{obs}$ is the Hubble parameter value extracted from the cosmic
observations while $H_{th}$ represents its theoretical value calculated at
$z_{k}$ with parameter space $\theta$, and $\sigma_{H(z_{k})}$ denotes
corresponding error.

####  IV.1.2 Pantheon datasets

In the present manuscript, we work with a recently published Pantheon
supernovae data set that contains 1048 supernovae samples with their distance
moduli $\mu^{obs}$ in the redshift range $z\in[0.01,2.3]$ [56]. The $\chi^{2}$
function corresponding to Pantheon data points reads as

$\chi^{2}_{SN}=\sum_{i,j=1}^{1048}\bigtriangledown\mu_{i}\left(C^{-1}_{SN}\right)_{ij}\bigtriangledown\mu_{j},$
(32)

Here $C_{SN}$ is the covariance matrix [56], and

$\displaystyle\quad\bigtriangledown\mu_{i}=\mu^{th}(z_{i},\theta)-\mu_{i}^{obs}.$

is the difference between the observed value of distance modulus extracted
from the cosmic observations and its theoretical values calculated from the
model with given parameter space $\theta$. We define the distance modulus by
$\mu=m_{B}-M_{B}$, where $m_{B}$ and $M_{B}$ denote respectively the observed
apparent magnitude and the absolute magnitude at a given redshift (Retrieving
the nuisance parameter following the recent approach called BEAMS with Bias
Correction (BBC) [57]). Furthermore, its theoretical value is given by

$\mu(z)=5log_{10}\left[\frac{D_{L}(z)}{1Mpc}\right]+25,$ (33)

where

$D_{L}(z)=c(1+z)\int_{0}^{z}\frac{dx}{H(x,\theta)}$ (34)

####  IV.1.3 BAO datasets

Baryonic Acoustic Oscillation (BAO) investigates oscillations produced in the
early phase of the universe due to cosmological perturbations in the fluid
consisting of photons, baryons, and dark matter, which is tightly coupled
through Thompson scattering. The BAO measurements consist of Sloan Digital Sky
Survey (SDSS), Six Degree Field Galaxy Survey (6dFGS), and the Baryon
Oscillation spectroscopic Survey (BOSS) [58, 59]. The relations used in BAO
measurements are,

$d_{A}(z)=\int_{0}^{z}\frac{dz^{\prime}}{H(z^{\prime})},$ (35)
$D_{V}(z)=\left(d_{A}(z)^{2}z/H(z)\right)^{1/3},$ (36)

and

$\chi_{BAO}^{2}=X^{T}C^{-1}X$ (37)

where $C$ is the covariance matrix [60], $d_{A}(z)$ denotes the angular
diameter distance while $D_{V}(z)$ represents the dilation scale.

####  IV.1.4 Observational Results

We obtained the constraints on parameters of our cosmological model with bulk
viscous matter corresponding to all three different cases for the combined
H(z)+Pantheon+BAO data set by minimizing the total chi-squared function
$\chi_{H}^{2}+\chi_{SN}^{2}+\chi_{BAO}^{2}$.

Figure 1: Constraints on the model and bulk viscous parameters corresponding
to Case I at $1-\sigma$ and $2-\sigma$ confidence interval using the combined
H(z)+Pantheon+BAO data set. Figure 2: Constraints on the model and bulk
viscous parameters corresponding to Case II at $1-\sigma$ and $2-\sigma$
confidence interval using the combined H(z)+Pantheon+BAO data set. Figure 3:
Constraints on the model and bulk viscous parameters corresponding to Case III
at $1-\sigma$ and $2-\sigma$ confidence interval using the combined
H(z)+Pantheon+BAO data set.

The obtained constraints on the model and bulk viscous parameters corresponds
to all three cases for the combined H(z)+Pantheon+BAO data set presented in
Table 1.

Table 1: Table shows the constraints on the model and bulk viscous parameters corresponds to all three cases for the combined H(z)+Pantheon+BAO data set. Cases | $\alpha$ | $\bar{\zeta}_{0}$ | $\bar{\zeta}_{1}$ | $\bar{\zeta}_{2}$
---|---|---|---|---
Case I | $-0.9583\pm 0.0094$ | $5.0287\pm 0.0099$ | $-1.479\pm 0.010$ | $-1.6155\pm 0.0097$
Case II | $-0.9620\pm 0.0091$ | $2.6620\pm 0.0099$ | $-0.370\pm 0.010$ | $0$
Case III | $-0.925^{+0.036}_{-0.047}$ | $1.862^{+0.094}_{-0.078}$ | $0$ | $0$

###  IV.2 Phase Space Analysis

In this section, we are going to investigate the asymptotic behavior of our
cosmological bulk viscous model by utilizing the dynamical system approach.
First, we convert the cosmological equations of our model into a set of
autonomous differential equations, and then we analyze the corresponding phase
space. Now, we define the following dimensionless variables,

$x=\frac{-\rho}{3\alpha H^{2}}\>\>\text{and}\>\>y=\frac{1}{\frac{H_{0}}{H}+1}$
(38)

The variables $x$ and $y$ are called phase space variables and it lies in the
range $0\leqslant y\leqslant 1$, whereas the Friedmann equation (23) implies
$x=1$. Case I : When viscosity coefficient depends on both the velocity and
acceleration i.e.
$\zeta=\zeta_{0}+\zeta_{1}\left(\frac{\dot{a}}{a}\right)+\zeta_{2}\left(\frac{{\ddot{a}}}{\dot{a}}\right)=\zeta_{0}+\zeta_{1}H+\zeta_{2}\left(\frac{\dot{H}}{H}+H\right)$.

We start with defining the variable $N=ln(a)$. Then by solving Friedmann
equations with conservation equation for the bulk viscous matter, we obtained
the following autonomous differential equations in the variables $x$ and $y$,

$x^{\prime}=\frac{dx}{dN}=\frac{(x-1)}{(2\alpha+\bar{\zeta}_{2})y}\left[2\bar{\zeta}_{0}(1-y)+(2\bar{\zeta}_{1}-\bar{\zeta}_{2})y\right]=F_{1}(x,y)$
(39)
$y^{\prime}=\frac{dy}{dN}=\frac{(y-1)}{(2\alpha+\bar{\zeta}_{2})}\left[\bar{\zeta}_{0}(1-y)+(\bar{\zeta}_{1}+\bar{\zeta}_{2}+3\alpha)y\right]=F_{2}(x,y)$
(40)

Further, by using the definition of the equation of state and deceleration
parameter, we acquired

$q=\frac{1}{(2\alpha+\bar{\zeta}_{2})}\left[\alpha+\bar{\zeta}_{1}+\frac{\bar{\zeta}_{0}(1-y)}{y}\right]$
(41)

and

$\omega=\frac{1}{3(2\alpha+\bar{\zeta}_{2})}\left[2\bar{\zeta}_{1}-\bar{\zeta}_{2}+\frac{2\bar{\zeta}_{0}(1-y)}{y}\right]$
(42)

Now by solving equations $x^{\prime}=0$ and $y^{\prime}=0$, we obtained the
coordinates of the critical points $(x_{c},y_{c})$ corresponding to autonomous
equations (39) and (40) as,

$(x_{c},y_{c})=(1,1)\>\>\text{and}\>\>(x_{c},y_{c})=(1,\frac{\bar{\zeta}_{0}}{\bar{\zeta}_{0}-\bar{\zeta}_{1}-\bar{\zeta}_{2}-3\alpha})$
(43)

We investigate the stability of the given autonomous system in the
neighbourhood of the critical points. First we linearize the given autonomous
dynamical system by assuming small perturbations near the critical points
$(x,y)\longrightarrow(x_{c}+\delta x,y_{c}+\delta y)$ satisfying

$\left[\begin{array}[]{c}\delta x^{\prime}\\\ \delta y^{\prime}\\\
\end{array}\right]\,=\left[\begin{array}[]{cc}\left(\frac{\partial
F_{1}}{\partial x}\right)_{0}&\left(\frac{\partial F_{1}}{\partial
y}\right)_{0}\\\ \left(\frac{\partial F_{2}}{\partial
x}\right)_{0}&\left(\frac{\partial F_{2}}{\partial x}\right)_{0}\\\
\end{array}\right]\left[\begin{array}[]{c}\delta x\\\ \delta y\\\
\end{array}\right]$

Here the suffix 0 implies that the above Jacobian matrix calculated at the
critical points $(x_{c},y_{c})$ and it is given as

$J=\left[\begin{array}[]{cc}\left(\frac{1}{(2\alpha+\bar{\zeta}_{2})y}\left[2\bar{\zeta}_{0}(1-y)+(2\bar{\zeta}_{1}-\bar{\zeta}_{2})y\right]\right)_{0}&\left(\frac{-2\bar{\zeta}_{0}(x-1)}{(2\alpha+\bar{\zeta}_{2})y^{2}}\right)_{0}\\\
0&\left(\frac{1}{(2\alpha+\bar{\zeta}_{2})}\left[-2\bar{\zeta}_{0}(y-1)+(\bar{\zeta}_{1}+\bar{\zeta}_{2}+3\alpha)(2y-1)\right]\right)_{0}\\\
\end{array}\right]$ (44)

Now we analyze the behavior of each of the critical points obtained in (43).
(i) $(x_{c},y_{c})=(1,1)$ :
In this case the eigenvalues obtained for the Jacobian matrix $J$ are

$\lambda_{1}=\frac{(2\bar{\zeta}_{1}-\bar{\zeta}_{2})}{(2\alpha+\bar{\zeta}_{2})}\>\>\text{and}\>\>\lambda_{2}=\frac{(\bar{\zeta}_{1}+\bar{\zeta}_{2}+3\alpha)}{(2\alpha+\bar{\zeta}_{2})}$
(45)

By using the constrained values of the parameters, one can have

$\lambda_{1}=0.38\>\>\text{and}\>\>\lambda_{2}=1.69$ (46)

Since $\lambda_{1}>0$ and $\lambda_{2}>0$, therefore the critical point
$(1,1)$ is unstable. Further, $y_{c}=1$ implies
$\frac{1}{\frac{H_{0}}{H}+1}=1$ i.e. $\frac{H_{0}}{H}=0$ which indicates that
either $H_{0}=0$ or $H\rightarrow\infty$. Since $H_{0}\neq 0$, therefore the
critical point $(x_{c},y_{c})=(1,1)$ represents initial singularity
characterized by $H\rightarrow\infty$ i.e. a past attractor. Moreover from
equations (41) and (42), we obtained $q\sim 0.69$ and $\omega\sim 0.12$.

(ii)
$(x_{c},y_{c})=(1,\frac{\bar{\zeta}_{0}}{\bar{\zeta}_{0}-\bar{\zeta}_{1}-\bar{\zeta}_{2}-3\alpha})=(1,0.4572)$
:
In this case the eigenvalues obtained for the Jacobian matrix $J$ are

$\lambda_{1}=-3\>\>\text{and}\>\>\lambda_{2}=-\frac{(\bar{\zeta}_{1}+\bar{\zeta}_{2}+3\alpha)}{(2\alpha+\bar{\zeta}_{2})}$
(47)

Again by using the constrained values of the parameters, we can have

$\lambda_{1}=-3\>\>\text{and}\>\>\lambda_{2}=-1.69$ (48)

Since $\lambda_{1}<0$ and $\lambda_{2}<0$, therefore the critical point
$(1,0.4572)$ is stable. In addition, from equations (41) and (42), we found
that $q\sim-1$ and $\omega\sim-1$. Hence the critical point
$(x_{c},y_{c})=(1,0.4572)$ corresponds to de-Sitter phase representing a
future attractor.

Table 2: Table shows the critical points and their behavior corresponding to Case I. Critical Points $(x_{c},y_{c})$ |  Eigenvalues $\lambda_{1}$ and $\lambda_{2}$ | Nature of critical point | $q$ | $\omega$
---|---|---|---|---
$(1,1)$ | $0.38\>\>\text{and}\>\>1.69$ | Unstable | $0.69$ | $0.12$
$(1,0.4572)$ | $-3\>\>\text{and}\>\>-1.69$ | Stable | $-1$ | $-1$
Figure 4: Phase plot in the $x-y$ plane corresponding to Case I with a red dot
and green dot denoting the past and future attractor, respectively, followed
by the arrowhead representing the direction of the trajectories.

From the phase plot presented in the Fig (4) it is evident that the
evolutionary trajectory of our bulk viscous cosmological model emerges from
the critical point $(1,1)$ as a past attractor and then it converges to the
critical point $(1,0.4572)$ that is nothing but a future attractor. Thus our
model corresponding to Case I depicts the evolution of the universe starting
with an initial singularity and behaves like a de-Sitter phase in the far
future.

Case II : When viscosity coefficient depends on the velocity but not on
acceleration i.e.
$\zeta=\zeta_{0}+\zeta_{1}\left(\frac{\dot{a}}{a}\right)=\zeta_{0}+\zeta_{1}H$.

In this case, the autonomous system of equations (39) and (40) becomes

$x^{\prime}=\frac{(x-1)}{\alpha
y}\left[\bar{\zeta}_{0}(1-y)+\bar{\zeta}_{1}y\right]$ (49)
$y^{\prime}=\frac{(y-1)}{2\alpha}\left[\bar{\zeta}_{0}(1-y)+(\bar{\zeta}_{1}+3\alpha)y\right]$
(50)

The deceleration and equation of state parameter given in (41) and (42)
reduces to

$q=\frac{1}{2\alpha}\left[\alpha+\bar{\zeta}_{1}+\frac{\bar{\zeta}_{0}(1-y)}{y}\right]$
(51)

and

$\omega=\frac{1}{6\alpha}\left[2\bar{\zeta}_{1}+\frac{2\bar{\zeta}_{0}(1-y)}{y}\right]$
(52)

Now on solving equations $x^{\prime}=0$ and $y^{\prime}=0$, we obtained the
coordinates of the critical points $(x_{c},y_{c})$ corresponding to autonomous
equations (49) and (50) as,

$(x_{c},y_{c})=(1,1)\>\>\text{and}\>\>(x_{c},y_{c})=(1,\frac{\bar{\zeta}_{0}}{\bar{\zeta}_{0}-\bar{\zeta}_{1}-3\alpha})$
(53)

We investigate the behavior of each of the critical points obtained in (53).
(i) $(x_{c},y_{c})=(1,1)$ :
In this case the eigenvalues obtained for the reduced Jacobian matrix $J$ are

$\lambda_{1}=\frac{\bar{\zeta}_{1}}{\alpha}\sim
0.38\>\>\text{and}\>\>\lambda_{2}=\frac{(\bar{\zeta}_{1}+3\alpha)}{2\alpha}\sim
1.69$ (54)

Since $\lambda_{1}>0$ and $\lambda_{2}>0$, therefore the critical point
$(1,1)$ is unstable. In addition, as in the previous case, $y_{c}=1$ implies
either $H_{0}=0$ or $H\rightarrow\infty$ and hence the critical point
$(x_{c},y_{c})=(1,1)$ represents initial singularity characterized by
$H\rightarrow\infty$ i.e. a past attractor since $H_{0}$ is non-zero. Further
from equations (51) and (52), we obtained $q\sim 0.69$ and $\omega\sim 0.12$.

(ii)
$(x_{c},y_{c})=(1,\frac{\bar{\zeta}_{0}}{\bar{\zeta}_{0}-\bar{\zeta}_{1}-3\alpha})=(1,0.4498)$
:
In this case the eigenvalues obtained for the reduced Jacobian matrix $J$ are

$\lambda_{1}=-3\>\>\text{and}\>\>\lambda_{2}=-\frac{(\bar{\zeta}_{1}+3\alpha)}{2\alpha}\sim-1.69$
(55)

Since $\lambda_{1}<0$ and $\lambda_{2}<0$, therefore the critical point
$(1,0.4498)$ is stable. From equations (51) and (52), we obtained $q\sim-1$
and $\omega\sim-1$. Hence the critical point $(x_{c},y_{c})=(1,0.4498)$
corresponds to de-Sitter type universe representing a future attractor.

Table 3: Table shows the critical points and their behavior corresponding to Case II. Critical Points $(x_{c},y_{c})$ |  Eigenvalues $\lambda_{1}$ and $\lambda_{2}$ | Nature of critical point | $q$ | $\omega$
---|---|---|---|---
$(1,1)$ | $0.38\>\>\text{and}\>\>1.69$ | Unstable | $0.69$ | $0.12$
$(1,0.4498)$ | $-3\>\>\text{and}\>\>-1.69$ | Stable | $-1$ | $-1$
Figure 5: Phase plot in the $x-y$ plane corresponding to Case II with a red
dot and green dot denoting the past and future attractor, respectively,
followed by the arrowhead representing the direction of the trajectories.

From the phase plot presented in the Fig (5) it is clear that the evolutionary
trajectory of our model emerges from the critical point $(1,1)$ and then it
converges to the critical point $(1,0.4498)$. Thus our model with bulk
viscosity corresponding to Case II represents the universe evolving from an
initial singularity to de-Sitter type universe in the far future, as we
obtained in Case I.

Case III : When viscosity coefficient does not depends on both the velocity
and acceleration i.e. $\zeta=\zeta_{0}$.

In this case, the autonomous system of equations (39) and (40) becomes

$x^{\prime}=\frac{\bar{\zeta}_{0}(x-1)(1-y)}{\alpha y}$ (56)
$y^{\prime}=\frac{(y-1)}{2\alpha}\left[\bar{\zeta}_{0}(1-y)+3\alpha y\right]$
(57)

The deceleration and equation of state parameter given in (41) and (42)
reduces to

$q=\frac{1}{2}+\frac{\bar{\zeta}_{0}(1-y)}{2\alpha y}$ (58)

and

$\omega=\frac{\bar{\zeta}_{0}(1-y)}{3\alpha y}$ (59)

On solving equations $x^{\prime}=0$ and $y^{\prime}=0$, we obtained the
coordinates of the critical points $(x_{c},y_{c})$ corresponding to autonomous
equations (56) and (57) as,

$(x_{c},y_{c})=(x,1)\>\>\text{and}\>\>(x_{c},y_{c})=(1,\frac{\bar{\zeta}_{0}}{\bar{\zeta}_{0}-3\alpha})$
(60)

(i) $(x_{c},y_{c})=(x,1)$ :
In this case the eigenvalues obtained for the reduced Jacobian matrix $J$ are

$\lambda_{1}=0\>\>\text{and}\>\>\lambda_{2}=\frac{3}{2}$ (61)

It is evident that first component $x$ of the critical point varies from 0 to
1 whereas $y_{c}=1$. Therefore the critical point $(x_{c},y_{c})$ is not an
isolated point rather it is a line of critical points and that indicates
initial state of universe since $H\rightarrow\infty$. As $\lambda_{1}=0$ and
$\lambda_{2}>0$, therefore every critical point $(x,1)$ on the line are
unstable i.e. a past attractor. Further we obtained $q=\frac{1}{2}$ and
$\omega=0$ by using equations (58) and (59). These values of the equation of
state and deceleration parameter corresponds to a decelerated matter dominated
universe.

(ii)
$(x_{c},y_{c})=(1,\frac{\bar{\zeta}_{0}}{\bar{\zeta}_{0}-3\alpha})=(1,0.4015)$
:
In this case the eigenvalues obtained for the reduced Jacobian matrix $J$ are

$\lambda_{1}=-3\>\>\text{and}\>\>\lambda_{2}=-\frac{3}{2}$ (62)

Since $\lambda_{1}<0$ and $\lambda_{2}<0$, therefore the critical point
$(1,0.4015)$ is stable. From equations (58) and (59), we obtained $q\sim-1$
and $\omega\sim-1$. Hence the critical point $(x_{c},y_{c})=(1,0.4015)$
corresponds to de-Sitter type universe representing a future attractor.

Table 4: Table shows the critical points and their behavior corresponding to Case III. Critical Points $(x_{c},y_{c})$ |  Eigenvalues $\lambda_{1}$ and $\lambda_{2}$ | Nature of critical point | $q$ | $\omega$
---|---|---|---|---
$(x,1)$ | $0\>\>\text{and}\>\>\frac{3}{2}$ | Unstable | $\frac{1}{2}$ | $0$
$(1,0.4015)$ | $-3\>\>\text{and}\>\>-\frac{3}{2}$ | Stable | $-1$ | $-1$
Figure 6: Phase plot in the $x-y$ plane corresponding to Case III with the red
line of critical points and green dot denoting the past and future attractor,
respectively, followed by the arrowhead representing the direction of the
trajectories.

From the phase plot presented in the Fig (6) it is clear that the evolution
trajectory of our model emerges from the non isolated critical point $(x,1)$
that corresponds to a decelerated matter dominated universe and then it
converges to the critical point $(1,0.4015)$ representing a future attractor.
Thus our cosmological model with bulk viscosity behaves like $\Lambda$CDM
model. We found that the accelerated deSitter like phase comes purely from the
$\bar{\zeta}_{0}$ case without any geometrical modification to GR. Further,
due to the constraint $x=1$, only the critical point $(1,1)$ lie on the line
of critical points $(x,1)$ is physical.

###  IV.3 Behavior of Cosmological Parameters

In this section we present physical behavior of well known cosmological
parameters such as pressure component in the presence of viscosity, effective
equation of state (EoS), and the $r-s$ parameter [61].

Figure 7: Plot for pressure component in the presence of bulk viscosity
corresponding to the observational constraints obtained using combined
H(z)+Pantheon+BAO data set. Figure 8: Plot for effective EoS parameter
corresponding to the observational constraints obtained using combined
H(z)+Pantheon+BAO data set. Figure 9: Plot for $r-s$ parameter corresponding
to the observational constraints obtained using combined H(z)+Pantheon+BAO
data set.

Fig(7) indicates that the pressure component in the presence of bulk viscosity
shows negative behavior for all three cases. The negative behavior of the
pressure component confirms the existence of a mysterious dark energy
component leading to acceleration. The plot for the effective EoS parameter in
fig(8) supports the recently observed acceleration, and it will converge to
the $\Lambda$CDM equation of state in the far future. In fig(9), we present
the evolutionary trajectory of the given viscous fluid model in the $r-s$
plane corresponding to all three cases. Since all three trajectories lie in
the region $r<1$ and $s>0$, our bulk viscous model follows the quintessence
scenario. Further, we observed that corresponding to all three cases, the
trajectories of our model converge to the $\Lambda$CDM fixed point. Thus, from
the $r-s$ diagram, we conclude that our viscous fluid cosmological model
represents a stable de-Sitter phase of the universe in the far future.

##  V Non-Linear $f(Q)$ model

We consider the following $f(Q)$ function for our analysis [30],

$f(Q)=-Q+\beta Q^{2}$ (63)

We consider only the case $\bar{p}=-3\zeta_{0}H=-\bar{\zeta}_{0}H_{0}H$. Then
by using Friedmann equations (19) and (20) we obtained the following first-
order non-linear differential equation,

$2\dot{H}(36\beta H^{2}-1)+3H^{2}(18\beta H^{2}-1)+\bar{\zeta}_{0}H_{0}H=0$
(64)

Now, we define the following dimensionless variables,

$x=\frac{-\rho}{3H^{2}(18\beta
H^{2}-1)}\>\>\text{and}\>\>y=\frac{1}{\frac{H_{0}}{H}+1}$ (65)

The variable $y$ lies in the range $0\leqslant y\leqslant 1$ and from equation
(19), we have constraint $x=1$. Then we obtained the following autonomous
differential equations corresponding to our non-linear $f(Q)$ model,

$x^{\prime}=\frac{dx}{dN}=\frac{\bar{\zeta}_{0}(x-1)(1-y)^{3}}{y\left[\bar{\beta}y^{2}-(1-y)^{2}\right]}$
(66)
$y^{\prime}=\frac{dy}{dN}=\frac{(y-1)}{2\left[2\bar{\beta}y^{2}-(1-y)^{2}\right]}\left[3y\\{\bar{\beta}y^{2}-(1-y)^{2}\\}+\bar{\zeta}_{0}(1-y)^{3}\right]$
(67)

Here $\bar{\beta}=18H_{0}^{2}\beta$. Further, by using the definition of the
equation of state and deceleration parameter, we acquired

$q=-1+\frac{1}{2y\left[2\bar{\beta}y^{2}-(1-y)^{2}\right]}\left[3y\\{\bar{\beta}y^{2}-(1-y)^{2}\\}+\bar{\zeta}_{0}(1-y)^{3}\right]$
(68)

and

$\omega=-1+\frac{1}{3y\left[2\bar{\beta}y^{2}-(1-y)^{2}\right]}\left[3y\\{\bar{\beta}y^{2}-(1-y)^{2}\\}+\bar{\zeta}_{0}(1-y)^{3}\right]$
(69)

Since $\bar{\zeta}_{0}>0$, we fix $\bar{\zeta}_{0}=1$ and then we have
investigated the stability of the autonomous system given by (66)-(67) in the
neighbourhood of the critical points. In particular, for $\bar{\beta}=0$, this
autonomous system is same as that of Case III and hence the dynamics.
Therefore we analyzed our non-linear $f(Q)$ model for the parameter values
$\bar{\beta}=-1$ and $\bar{\beta}=1$. The obtained results are presented in
the Table 5.

Table 5: Table shows the critical points and their behavior corresponding to non-linear $f(Q)$ model. |  Critical Points $(x_{c},y_{c})$ |  Eigenvalues $\lambda_{1}$ and $\lambda_{2}$ | Nature of critical point | $q$ | $\omega$
---|---|---|---|---|---
$\bar{\beta}=-1$ | $\left(x,1\right)$ | $\frac{3}{4}\>\>\text{and}\>\>0$ | Unstable | $-\frac{1}{4}$ | $-\frac{1}{2}$
| $\left(1,0.234\right)$ | $-3\>\>\text{and}\>\>-1.62$ | Stable | $-1$ | $-1$
| $\left(x,1\right)$ | $\frac{3}{4}\>\>\text{and}\>\>0$ | Unstable | $-\frac{1}{4}$ | $-\frac{1}{2}$
$\bar{\beta}=1$ | $\left(1,0.283\right)$ | $-3\>\>\text{and}\>\>-1.15$ | Stable | $-1$ | $-1$
| $\left(1,0.426\right)$ | $-3\>\>\text{and}\>\>-9.6$ | Stable | $-1$ | $-1$
Figure 10: Phase plot in the $x-y$ plane corresponding to the non-linear
$f(Q)$ model ($\bar{\beta}=-1$) with the red line of critical points and green
dot denoting the past and future attractor, respectively, followed by the
arrowhead representing the direction of the trajectories.  Figure 11: Phase
plot in the $x-y$ plane corresponding to the non-linear $f(Q)$ model
($\bar{\beta}=1$) with the red line of critical points and green dot denoting
the past and future attractor, respectively, followed by the arrowhead
representing the direction of the trajectories.  Figure 12: $\dot{H}$ as a
function of $\frac{H}{H_{0}}$ from equation (64) with
$\bar{\beta}=\bar{\zeta}_{0}=1$, where vertical black lines indicating the
fixpoints that corresponds to $\dot{H}=0$.

From the Table 5 it is clear that there is only one de-Sitter type future
attractor for the case $\bar{\beta}=-1$ whereas two future attractors
corresponding to $\bar{\beta}=1$ case. From the phase plot presented in the
Fig (10) it is clear that the evolution trajectory of our non-linear model,
for the case $\bar{\beta}=-1$, emerges from the non isolated critical points
$(x,1)$ and then it converges to the critical point $(1,0.234)$ representing a
de-Sitter like accelerated phase. Further note that, due to the constraint
$x=1$, only the critical point $(1,1)$ lie on the line of critical points
$(x,1)$ is physical. Moreover, from the phase plot presented in the Fig (11)
we found that the evolution trajectory of our non-linear model, for the case
$\bar{\beta}=1$, emerges from the critical point $(1,1)$ and then settles down
to two stable critical points $(1,0.283)$ and $(1,0.426)$ representing the de-
Sitter like accelerated phase. The behavior of $\dot{H}$ as a function of
$\frac{H}{H_{0}}$ from equation (64) with $\bar{\beta}=\bar{\zeta}_{0}=1$ is
presented in Fig (12). We can observe that two regimes are separated by a line
$H=\frac{H_{0}}{\sqrt{2\bar{\beta}}}$ of diverging $\dot{H}$ and the evolution
only to the upper attractor $\left(1,0.426\right)$ from an initial big bang
$y=1$ is possible. Thus, only the fixed point $\left(1,0.426\right)$ is a
viable attractor from an initial big bang. Hence, evolution from the big bang
$H\rightarrow\infty$ can only occur towards the right fix point at $H\approx
0.75H_{0}$.

##  VI Conclusion

In this manuscript, we have investigated the role of different bulk viscosity
coefficients in cosmic evolution under the framework of symmetric teleparallel
gravity. We have considered three different cases of bulk viscosity
coefficients which is well known in the literature, namely
$(i)\zeta=\zeta_{0}+\zeta_{1}\left(\frac{\dot{a}}{a}\right)+\zeta_{2}\left(\frac{{\ddot{a}}}{\dot{a}}\right)$,
$(ii)\zeta=\zeta_{0}+\zeta_{1}\left(\frac{\dot{a}}{a}\right)$, and
$(iii)\zeta=\zeta_{0}$. Then we calculated Hubble parameter value in terms of
redshift corresponding to all three cases by assuming a linear $f(Q)$ model,
specifically, $f(Q)=\alpha Q$ where $\alpha\neq 0$. We presented the
observation constraints (see Table 1 and Fig(1)-(3)) on the model parameter
and viscosity coefficients corresponding to all the three cases using the
combined H(z)+Pantheon+BAO data set. Then we investigated the asymptotic
behavior of our cosmological bulk viscous model by using phase space analysis.
We derived the set of autonomous differential equations and the corresponding
critical points for all three cases (see Table 2,3, and 4). Further we
presented the phase space plot corresponding to all the three cases (see
Fig(4),(5), and (6)). We found that our viscous fluid model represents a
universe evolving from matter dominated decelerated epoch, which is a past
attractor, to a stable de-sitter accelerated epoch, which is a future
attractor. In addition, we investigated the physical behavior of the pressure
component in the presence of viscosity, the effective equation of state (EoS),
and the $r-s$ parameter. We found that the pressure component in the presence
of bulk viscosity presented in fig(7) shows negative behavior for all three
cases. The plot for effective EoS parameter presented in fig(8) shows that the
current universe is going through a period of accelerated expansion. Lastly,
from fig(9), we obtained that trajectories of our cosmological viscous model
lie in the quintessence region. Further, these trajectories converge to the
$\Lambda$CDM fixed point for all three cases, which coincides with results
obtained in the phase space analysis. We found that the accelerated deSitter
like phase comes purely from the $\bar{\zeta}_{0}$ case without any geometric
modification to GR. Moreover, we found that the late-time behavior of all
three cases of viscosity coefficients are identical. Further, we have
considered a non-linear $f(Q)$ model, specifically, $f(Q)=-Q+\beta Q^{2}$ and
then we analyzed the behavior of model using dynamical approach presented in
Table 5. We found that there is only one de-Sitter type future attractor for
the case $\bar{\beta}=-1$ whereas two future attractors corresponding to
$\bar{\beta}=1$ case (see Fig(10) and (11)). However, the upper attractor
$\left(1,0.426\right)$ in Fig (11) is a viable attractor from an initial big
bang, and therefore evolution from the big bang $H\rightarrow\infty$ can only
occur towards the right fix point at $H\approx 0.75H_{0}$ (see Fig (12)).
Moreover, the dynamics of the case $\bar{\beta}=0$ is same as that of Case
III. Hence we conclude that the late-time behavior of the considered non-
linear model $f(Q)=-Q+\beta Q^{2}$ with $\beta\leq 0$ is similar to the linear
case, whereas for the case $\beta>0$ results are quite different.

## Acknowledgments

RS acknowledges UGC, New Delhi, India for providing Senior Research Fellowship
with (UGC-Ref. No.: 191620096030). DSR acknowledges UGC, New Delhi, India for
providing Junior Research Fellowship with (NTA-UGC-Ref.No.: 211610106591). SM
acknowledges DST, Govt. of India, New Delhi, for providing Senior Research
Fellowship (File No. DST/INSPIRE Fellowship/2018/IF18D676). PKS acknowledges
the Science and Engineering Research Board, Department of Science and
Technology, Government of India for financial support to carry out the
Research project No.: CRG/2022/001847. We are very much grateful to the
honorable referee and to the editor for the illuminating suggestions that have
significantly improved our work in terms of research quality, and
presentation.

## References

  * [1] A.G. Riess et al., Astron. J. 116, 1009 (1998).
  * [2] S. Perlmutter et al., Astrophys. J. 517, 565 (1999).
  * [3] D.J. Eisenstein et al., Astrophys. J. 633, 560 (2005).
  * [4] W.J. Percival at el., Mon. Not. R. Astron. Soc. 401, 2148 (2010).
  * [5] C.L. Bennett et al., Astrophys. J. Suppl. 148, 119-134 (2003).
  * [6] D.N. Spergel et al., [WMAP Collaboration], Astrophys. J. Suppl. 148, 175 (2003).
  * [7] R.R. Caldwell, M. Doran, Phys. Rev. D 69, 103517 (2004).
  * [8] Z.Y. Huang et al., JCAP 0605, 013 (2006).
  * [9] E. J. Copeland, M. Sami, and S. Tsujikawa, Int. J. Mod. Phys. D 15, 1753 (2006).
  * [10] L. Amendola et al., Phys. Rev. D 75, 083504 (2007).
  * [11] S. Appleby and R. Battye, Phys. Lett. B 654, 7 (2007).
  * [12] R. Ferraro, F. Fiorini, Phys. Rev. D 75, 084031 (2007).
  * [13] J. B. Jiménez et al., Phys. Rev. D 98, 044048 (2018).
  * [14] J. B. Jiménez, L. Heisenberg, and T. Koivisto, Universe, 5, 173 (2019).
  * [15] Y. Xu et al., Eur. Phys. J. C 79, 708 (2019).
  * [16] Y. Xu et al., Eur. Phys. J. C 80, 449 (2020).
  * [17] S. Mandal et al.,  Phys. Rev. D 102, 124029(2020).
  * [18] S. Mandal et al.,  Phys. Rev. D 102, 024057(2020).
  * [19] M. Hohmann et al, Phys. Rev. D 99, 024009 (2019).
  * [20] B. J. Barros et al., Phys.Dark Univ. 30, 100616 (2020).
  * [21] I. Soudi et al., Phys. Rev. D 100, 044008 (2019).
  * [22] F. D Ambrosio et al., Phys. Rev. D 105, 024042 (2022).
  * [23] N. Dimakis, A. Paliathanasis, and T. Christodoulakis, Class. Quant. Grav. 38, 225003 (2021).
  * [24] M. Hohmann, Phys. Rev. D 104, 124077 (2021).
  * [25] G. Mustafa et al., Phys. Lett. B 821, 136612 (2021).
  * [26] F. K. Anagnostopoulos, S. Basilakos, and E. N.Saridakis, Phys. Lett. B 822, 136634 (2021).
  * [27] F. Bajardi, D. Vernieri, and S. Capozziello, Eur. Phys. J. Plus 135, 11,912 (2020).
  * [28] K. Hu, T. Katsuragawa, and T. Qiu Phys. Rev. D 106, 044025 (2022).
  * [29] R. Solanki and P. K. Sahoo, Annalen Phys. 534, 2200076 (2022).
  * [30] J.B. Jiménez et al.,  Phys. Rev. D 101, 103507(2020).
  * [31] C. Eckart, Phys. Rev. 58, 919(1940).
  * [32] W. Israel, J. M. Stewart, Phys. Lett. B 58, 213 (1976).
  * [33] W. Israel, Ann. Phys. (N.Y.) 100, 310 (1976).
  * [34] W. Israel, J. M. Stewart, Proc. R. Soc. Lond. B 365, 43 (1979).
  * [35] I. Brevik, Entropy 2012(14), 2302-2310 (2012).
  * [36] I. Brevik and O. Gron, Astrophys. Space Sci. 347, 399 (2013).
  * [37] I. Brevik et al., Int. J. Mod. Phys. D 26, 1730024 (2017).
  * [38] I. Brevik, A. N. Makarenko, and A. V. Timoshkin, Int. J. Geom. Methods Mod. 16, 1950150 (2019)
  * [39] I. Brevik and B. D. Normann, Symmetry 2020(12), 1085 (2020).
  * [40] N. D. J. Mohan, A. Sasidharan, and T. K. Mathew, Eur. Phys. J. C 77, 849 (2017).
  * [41] A. V. Astashenok, S. D. Odintsov, and A. S. Tepliakov, Nucl. Phys. B 974, 115646 (2022).
  * [42] A. Sasidharan and T. K. Mathew, Eur. Phys. J. C 75, 348 (2015).
  * [43] T. Ortin, Gravity and Strings, Cambridge Monographs on Mathematical Physics (Cambridge University Press (2015).
  * [44] C. Misner et al., Gravitation, Princeton University Press (2017).
  * [45] J.B. Jiménez, L. Heisenberg, and T. S. Koivisto, JCAP 1808, 039 (2018).
  * [46] R. Lazkoz et al., Phys. Rev. D 100, 104027 (2019).
  * [47] L. Heisenberg, M. Hohmann, and S. Kuhn, arXiv, arXiv:2212.14324 (2022).
  * [48] B. Ryden, Introduction to Cosmology (Addison Wesley, SanFrancisco, United States of America, 2003).
  * [49] F. DÁmbrosio, L. Heisenberg, and S. Kuhn, Class. Quantum Grav. 39, 025013 (2021).
  * [50] R. Solanki et al., Phys. Dark Univ. 32, 100820 (2021).
  * [51] Z. Hassan, S. Mandal, and P. K. Sahoo, Forts. Phys. 69, 2100023 (2021).
  * [52] A. Sasidharan and T. K. Mathew, J. High Energy Phys. 06, 138 (2016).
  * [53] J. Ren and X. H. Meng, Phys. Lett. B 633, 1 (2006).
  * [54] D. F. Mackey et al., Publ. Astron. Soc. Pac. 125, 306(2013).
  * [55] Planck Collaboration, Astron. Astrophys. 641, A6 (2020).
  * [56] D.M. Scolnic et al., Astrophys. J. 859, 101(2018).
  * [57] R. Kessler and D. Scolnic, Astrophys. J. 836,56 (2017).
  * [58] C. Blake et al.,  Mon. Not. Roy. Astron. Soc. 418, 1707 (2011).
  * [59] W. J. Percival et al., Mon. Not. Roy. Astron. Soc. 401, 2148 (2010).
  * [60] R. Giostri et al., J. Cosm. Astropart. Phys. 1203, 027 (2012).
  * [61] V. Sahni et al., JETP Lett. 77, 201 (2003).