I came across the following series when computing the covariance of a transform of a bivariate Gaussian random vector via Hermite polynomials and Mehler's expansion:

$$ S = \sum_{n=1}^{\infty} \frac{\rho^n}{n^{1/6}} $$ for $\vert \rho \vert < 1$. We know that $S$ must be finite and satisfy $$ S \le \rho (1-\rho)^{-1} $$ since the original series is dominated by $\sum_{n=1}^{\infty} \rho^n$.

However, there is a catch if we use for $S$ the upper bound $\rho (1-\rho)^{-1}$, which tends to $\infty$ when $\rho \to 1-$. This happens when the two marginal random variables in the Gaussian vector are almost surely, positively linearly dependent (asymptotically).

**So, the target is to obtain a good upper bound, much better than $\rho (1-\rho)^{-1}$ when we restrict $\rho$ to be away from $1$, to reduce the effect of $\rho \to 1-$. In other words, let $1-\rho = \delta$ for some fixed $\delta \in (0,1)$, what is a better upper bound for $S$?**

Because of the scaling term $n^{-1/6}$ that induces a divergent series $\sum_{n=1}^{\infty} n^{-1/6}$, probably not much improvement should be expected. I have Googled but did not find an illuminating technique for this. Any pointer or help is appreciated. Thank you.