  Given a smooth closed Riemannian manifold $(M,g)$ of dimension $N \ge 3$, we
derive sharp quantitative stability estimates for nonnegative functions near
the solution set of the Yamabe problem on $(M,g)$. The seminal work of Struwe
(1984) \cite{S} states that if $\Gamma(u) := \|\Delta_g u - \frac{N-2}{4(N-1)}
R_g u + u^{\frac{N+2}{N-2}}\|_{H^{-1}(M)} \to 0$, then
$\|u-(u_0+\sum_{i=1}^{\nu} \mathcal{V}_i)\|_{H^1(M)} \to 0$ where $u_0$ is a
solution to the Yamabe problem on $(M,g)$, $\nu \in \mathbb{N} \cup \{0\}$, and
$\mathcal{V}_i$ is a bubble-like function. If $M$ is the round sphere
$\mathbb{S}^N$, then $u_0 \equiv 0$ and a natural candidate of $\mathcal{V}_i$
is a bubble itself. If $M$ is not conformally equivalent to $\mathbb{S}^N$,
then either $u_0 > 0$ or $u_0 \equiv 0$, there is no canonical choice of
$\mathcal{V}_i$, and so a careful selection of $\mathcal{V}_i$ must be made to
attain optimal estimates.
  For $3 \le N \le 5$, we construct suitable $\mathcal{V}_i$'s and then
establish the inequality $\|u-(u_0+\sum_{i=1}^{\nu} \mathcal{V}_i)\|_{H^1(M)}$
$ \le C\zeta(\Gamma(u))$ where $C > 0$ and $\zeta(t) = t$, consistent with the
result of Figalli and Glaudo (2020) \cite{FG} on $\mathbb{S}^N$. In the case of
$N \ge 6$, we investigate the single-bubbling phenomenon $(\nu = 1)$ on generic
Riemannian manifolds $(M,g)$, proving that $\zeta(t)$ is determined by $N$,
$u_0$, and $g$, and can be much larger than $t$. This exhibits a striking
difference from the result of Ciraolo, Figalli, and Maggi (2018) \cite{CFM} on
$\mathbb{S}^N$. All of the estimates presented herein are optimal.
