What is the result or equivalent order of this infinite sum? As described by the title, can we obtain an analytical result for the following infinite sum or provide an equivalent order as $p \rightarrow 1$?
$$
\sum_{j=1}^{\infty} \left(\frac{p^j}{1-p^j}\right)^{i+1},i\ \text{is non-negative integer}.
$$
Now, I only know that
$$\sum_{j=1}^{\infty} \frac{p^j}{1-p^j}\approx \frac{1}{1-p}\ln\left(\frac{1}{1-p}\right).$$
I often have to compute the sum of infinite series in my research. I really appreciate if one can provide a textbook or useful handbook to reference.
 A: If a heuristic consideration is acceptable, we can do the following.
Let's denote $\epsilon=1-p \,(p<1)$, and $$ S(n, \epsilon)=\sum_{k=1}^{\infty} \left(\frac{p^k}{1-p^k}\right)^n=\sum_{k=1}^{\infty} \frac{e^{nk\ln(1-\epsilon)}}{(1-e^{k\ln(1-\epsilon)})^n}=\sum_{k=1}^{\infty} e^{nk\ln(1-\epsilon)}(1-e^{k\ln(1-\epsilon)})^{-n}$$
We consider the case of $n\epsilon\ll1$, and $n>2$ (is not necessarily an integer).
When $k$ is growing, the nominator is decreasing, and the denominator - increasing. At $\,k\epsilon =\ln2 \,\,(k=\frac{\ln2}{\epsilon}) \,\,\frac{p^k}{1-p^k}=1$, and the main contribution to the sum comes from the terms with $k\ll\frac{\ln2}{\epsilon}$.
Bearing this in mind, we can use the following approximation:
$$e^{nk\ln(1-\epsilon)}\approx e^{-nk\epsilon}\tag{1}$$
$$(1-e^{k\ln(1-\epsilon)})^{-n}\approx(1-e^{-k\epsilon-\frac{k\epsilon^2}{2}})^{-n}\approx\left(k\epsilon-\frac{k^2\epsilon^2}{2}+\frac{k\epsilon^2}{2}\right)^{-n}$$
$$=(k\epsilon)^{-n}\left(1-\frac{k\epsilon}{2}+\frac{\epsilon}{2}\right)^{-n}\approx (k\epsilon)^{-n}e^{-n\ln\left(1-\frac{k\epsilon}{2}+\frac{\epsilon}{2}\right)}\approx(k\epsilon)^{-n}e^{-\frac{n\epsilon}{2}}e^{\frac{nk\epsilon}{2}}\tag{2}$$
and our sum looks like
$$S(\epsilon, n)\approx\frac{e^{-\frac{n\epsilon}{2}}}{\epsilon^n}\sum_{k=1}^{\infty}\frac{e^{\frac{-nk\epsilon}{2}}}{k^n}=\frac{e^{-\frac{n\epsilon}{2}}}{\epsilon^n}\operatorname{Li}_n\left(e^{\frac{-nk\epsilon}{2}}\right)\tag{3}$$
Using the representation $\displaystyle \frac{1}{k^n}=\frac{1}{\Gamma(n)}\int_0^\infty t^{n-1}e^{-kt}$ and performing summation with respect to $k$,
$$S(\epsilon, n)\approx\frac{e^{-n\epsilon}}{\epsilon^n}\frac{1}{\Gamma(n)}\int_0^\infty\frac{t^{n-1}}{e^t-e^{-\frac{nk\epsilon}{2}}}dt\approx\frac{1-n\epsilon}{\epsilon^n}\frac{1}{\Gamma(n)}\int_0^\infty\frac{t^{n-1}}{e^t-1+\frac{n\epsilon}{2}}dt\tag{4}$$
The integral we evaluate in the following way:
$$\int_0^\infty\frac{t^{n-1}}{e^t-1+\frac{n\epsilon}{2}}dt=\int_0^\infty t^{n-1}\left(\frac{1}{e^t-1}+\frac{1}{e^t-1+\frac{n\epsilon}{2}}-\frac{1}{e^t-1}\right)dt$$
$$=\Gamma(n)\zeta(n)-\frac{n\epsilon}{2}\int_0^\infty\frac{t^{n-1}}{(e^t-1)(e^t-1+\frac{n\epsilon}{2})}dt$$
As we suppose $n>2$, the second integral converges at $n\epsilon=0$ in the denominator. Therefore, keeping only the terms up to $\sim n\epsilon$, we can write
$$\int_0^\infty\frac{t^{n-1}}{e^t-1+\frac{n\epsilon}{2}}dt\approx\Gamma(n)\zeta(n)-\frac{n\epsilon}{2}\int_0^\infty\frac{t^{n-1}}{(e^t-1)^2}dt$$
$$=\Gamma(n)\zeta(n)-\frac{n\epsilon}{2}\Gamma(n)\big(\zeta(n-1)-\zeta(n)\big)\tag{5}$$
Putting (5) int0 (4), with the required accuracy
$$\boxed{\,\,S(\epsilon, n)\approx\frac{\zeta(n)}{\epsilon^n}-\frac{n}{2\,\epsilon^{n-1}}\big(\zeta(n-1)+\zeta(n)\big)=\frac{\zeta(n)}{(1-p)^n}-\frac{n}{2\,(1-p)^{n-1}}\big(\zeta(n-1)+\zeta(n)\big)\,\,}$$

Now, we can make a numeric check with WolframAlpha (we keep in mind that it is supposed that $n\epsilon\ll1$):
$\displaystyle n=5\,\epsilon =0.1\quad\text{exact}=60 632;\,\text{approx}\approx 50712$
$\displaystyle n=3\,\epsilon =0.01\quad\text{exact}=1.1602 \times 10^6;\,\text{approx}\approx 1.1593 \times 10^6$
$\displaystyle n=2.6\,\epsilon =0.01\quad\text{exact}=199855;\,\text{approx}\approx 199505$
$\displaystyle n=5\,\epsilon =0.01\quad\text{exact}=9.8504 \times 10^9;\,\text{approx}\approx 9.8395 \times 10^9$
$\displaystyle n=2.6\,\epsilon =0.001\quad\text{exact}=8.20777 \times 10^7;\,\text{approx}\approx 8.20755 \times 10^7$
