Asymptotic expansion of a series I am interested in the asymptotics, as $x$ tends to $0$, of
$$f(x) = \sum_{n=1}^\infty \frac{1}{n}\frac{1}{(e^{nx}-1)^2}$$
This function is well defined for every $x > 0$ (for example, use $e^{nx}-1 \geq nx$).
Furthermore, Lebesgue's dominated convergence theorem shows that, 
$$
f(x) \sim \frac{\zeta(3)}{x^2}
$$
as $x$ tends to $0$, where  $\displaystyle\zeta(3) = \sum_{n=1}^\infty \frac{1}{n^3}$
is Apéry's constant.
Could you help me get a more precise asymptotic expansion of $f(x)$ as $x$ tends to $0$ ?
(best would be with $o(1)$)
 A: Note that 
$$\sum_{m=1}^\infty  \frac{m}{q^{m+1}} = \frac{1}{(q-1)^2}, \qquad |q|>1.$$
With $q= e^{nx}$, we have
$$f(x)=  
\sum_{n=1}^\infty \sum_{m=1}^\infty \frac{m e^{-n (m+1) x}}{n}
=-\sum_{m=1}^\infty m \log\left(1-e^{-(1+m)x}\right). $$
Next, we use the Euler-Maclaurin formula:
$$f(x) \sim -\int_1^\infty \!dz\,z \log\left(1-e^{-(1+z)x}\right)
- \overbrace{\frac{\log(1-e^{-2x})}{2}}^{\log x+ \mathcal{O}(1)}
+ \frac{1}{12} \overbrace{\left[ \frac{x}{e^{2x}-1}+\log \left(1-e^{-2 x}\right) \right]}^{\log x+\mathcal{O}(1)}
+ R$$
With $$|R| \leq C \int_1^{\infty} \left| g''(z) \right|dz$$ where $g(z) = z  \log\left(1-e^{-(1+z)x}\right) $. Numerical calculation shows $R = \mathcal{O}(1)$ though I have no proof yet.
Thus we have
$$f(x) \sim - \int_1^\infty \!dz\,z \log\left(1-e^{-(1+z)x}\right)
=  \frac{1}{x^2} \int_{2x}^\infty\!dx\,(x-y) \log(1-e^{-y}).$$
We need to evaluate
$$\begin{align}\int_{2x}^\infty\!dy\,(x-y) \log(1-e^{-y})
&= x \underbrace{\int_{0}^\infty\!dy\, \log(1-e^{-y})}_{-\pi^2/6}
-\underbrace{\int_{0}^\infty\!dy\, y\log(1-e^{-y})}_{\zeta(3)}\\
&\quad+ \underbrace{\int_0^{2x} \!dy\,(y-x)\log(1-e^{-y})}_{\mathcal{O}(x^2)}\end{align}  .$$
In conclusion, we have
$$f(x) \sim \frac{\zeta(3)}{x^2}- \frac{\pi^2}{6x} - \frac{5}{12}\log x+ \mathcal{O}(1).$$
A: This series may also be evaluated by calculating its Mellin transform and inverting that to get the asymptotic expansion.
We have $$ f(x) = \sum_{n\ge 1} \frac{1}{n} \frac{1}{(e^{nx}-1)^2} =
\sum_{n\ge 1} \frac{1}{n} \frac{1}{e^{2nx}} \frac{1}{(1-e^{-nx})^2}$$
which gives
$$ f(x) = \sum_{n\ge 1} \frac{1}{n} \frac{1}{e^{2nx}} \sum_{k\ge 0} (k+1) e^{-knx} =
\sum_{n\ge 1} \frac{1}{n} \sum_{k\ge 0} (k+1) e^{-(k+2)nx}.$$
Now recall that $$\mathfrak{M}(e^{-qx};s) = \Gamma(s) \frac{1}{q^s}$$
so that $$f^*(s) = \mathfrak{M}(f(x);s) = \Gamma(s) 
\sum_{n\ge 1} \frac{1}{n} \sum_{k\ge 0} (k+1) \frac{1}{(k+2)^s n^s} =
\Gamma(s) \zeta(s+1)
\left( - \sum_{k\ge 0}  \frac{1}{(k+2)^s} + \sum_{k\ge 0}  \frac{k+2}{(k+2)^s} \right)$$
and finally
$$\Gamma(s) \zeta(s+1) \left( 1 - \zeta(s) -1 + \zeta(s-1) \right)
= \Gamma(s) \zeta(s+1) \left( \zeta(s-1) - \zeta(s) \right).$$
The Mellin inversion integral is
$$ f(x) = \int_{4-i\infty}^{4+i\infty}  \frac{f^*(s)}{x^s} ds.$$ 
Introduce $$ L(s) =  \frac{f^*(s)}{x^s} .$$ The asymptotic expansion may now be read off from the residues at the poles:
$$\operatorname{Res}_{s=2} L(s) = \frac{\zeta(3)}{x^2},$$
$$\operatorname{Res}_{s=1} L(s) = -\frac{\pi^2}{6x},$$
$$\operatorname{Res}_{s=0} L(s) = 
\zeta'(-1) + \frac{1}{2}\log(2\pi) - \frac{5}{12}\log x.$$
Sum these to get the asymptotic expansion. There are some additional terms that are generated by the poles of the Gamma function.
A: The asymptotic series has an infinite coefficient ($\frac{5}{12} \zeta(1)$) for the $O(1)$ term.  The behavior of the series in the limit as $x \rightarrow 0$ is more complicated, having a term with a factor of $\log{x}$, I suspect.
