Computing the sum of an infinite series as a variant of a geometric series 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.
 A: Your function is the Polylogarithm function $Li_{1/6}(\rho)$.  Mathematica indicates the correct asymptotic (at least for real $\rho$) is $\Gamma(5/6)/(1-\rho)^{5/6}+O(1)$.
A: Let $r:=\rho$, $S(r):=S$, and $a:=1/6$. Let us show that 
\begin{equation*}
 S(r)< Cr(1-r)^{a-1}\quad\text{and}\quad S(r)\sim C(1-r)^{a-1},\quad\text{where}\quad  C:=\Gamma(1-a). \tag{0}
\end{equation*}
Everywhere here, the equalities and inequalities are for $r\in(0,1)$, and the limit relations are for $r\uparrow1$. 
Indeed, we have 
\begin{equation*}
 \frac1{n^a}=\frac1{\Gamma(a)}\,\int_0^\infty u^{a-1}e^{-nu}\,du
\end{equation*}
and hence 
\begin{equation*}
 S(r)=\frac1{\Gamma(a)}\,\int_0^\infty du\, u^{a-1}\,\sum_{n=1}^\infty r^n e^{-nu}
 =\frac r{\Gamma(a)}\,I(r), 
\end{equation*}
where 
\begin{align*}
 I(r)&:=\int_0^\infty du\, \frac{u^{a-1} e^{-u}}{1-r e^{-u}} \\ 
& =\int_0^\infty du\, \frac{u^{a-1}}{e^u-r}  \\ 
& <\int_0^\infty du\, \frac{u^{a-1}}{1+u-r}  \\ 
& =(1-r)^{a-1}\int_0^\infty dv\, \frac{v^{a-1}}{1+v} \\ 
&=\Gamma(1-a)\Gamma(a)(1-r)^{a-1}, 
\end{align*}
where we used the substitutions $u=(1-r)v$ and then $\frac1{1+v}=t$. 
So, the inequality in (0) is proved. 
Also, 
\begin{align*}
(1-r)^{1-a}I(r)& =\int_0^\infty dv\, \frac{v^{a-1}}{1+(e^{(1-r)v}-1)/(1-r)} \\ 
&\to\int_0^\infty dv\, \frac{v^{a-1}}{1+v} , 
\end{align*}
by dominated convergence. So, the asymptotic relation in (0) is proved. 
As seen from the above proof, the results in (0) hold for any $a\in(0,1)$. 
