  A classical result owing to Mancini and Sandeep [Ann. Sc. Norm. Super. Pisa
Cl. Sci. 7 (2008)] asserts that all positive solutions of the
Poincar\'e-Sobolev equation on the hyperbolic space $$ -\Delta_{\mathbb{B}^n}
u-\lambda u = |u|^{p-1}u, \quad u\in H^1(\mathbb{B}^n), $$ are unique up to
hyperbolic isometries where $n \geq 3,$ $1 < p \leq \frac{n+2}{n-2} $ and
$\lambda \leq \frac{(n-1)^2}{4}.$ We prove under certain bounds on $\|\nabla u
\|_{L^2(\mathbb{B}^n)}$ the inequality $$ \delta(u) \lesssim
\|\Delta_{\mathbb{B}^n} u+ \lambda u + u^{p}\|_{H^{-1}}, $$
  holds whenever $p >2$ and hence forcing the dimensional restriction $3 \leq n
\leq 5,$ where $\delta(u)$ denotes the $H^1$ distance of $u$ from the manifold
of sums of hyperbolic bubbles. Moreover, it fails for any $n \geq 3$ and $p \in
(1,2].$ This strengthens the phenomenon observed in the Euclidean case that the
(linear) quantitative stability estimate depends only on whether the exponent
$p$ is $>2$ or $\leq 2$. In the critical case, our dimensional constraint
coincides with the seminal result of Figalli and Glaudo [Arch. Ration. Mech.
Anal, 237 (2020)] but we notice a striking dependence on the exponent $p$ in
the subcritical regime as well which is not present in the flat case. Our
technique is an amalgamation of Figalli and Glaudo's method and builds upon a
series of new and novel estimates on the interaction of hyperbolic bubbles and
their derivatives and improved eigenfunction integrability estimates. Since the
conformal group coincides with the isometry group of the hyperbolic space, we
perceive a remarkable distinction in arguments and techniques to achieve our
main results compared to that of the Euclidean case.
