Limits of $\sum_{m=1}^{+\infty} \sum_{n=1}^{+\infty}e^{-mn x}$ at $0$ and $\infty$ Let $f(x) = \sum_{m=1}^{+\infty} \sum_{n=1}^{+\infty}e^{-mn x}$ for $x > 0$.

Prove that $f(x) \sim e^{-x}$ as $x \to \infty$ and $\lim_{x\to 0} x\cdot (f(x) + \frac{1}{x}\log x)= \gamma$ where $\gamma$ is the Euler–Mascheroni constant.

For the first point, I wanted to use the fact that $f(x) - e^{-x}$ involves only terms like $e^{-kx}$ with $k \geq 2$ but I don't know how to write it. I have no idea for the limit as $x$ tends towards $0$.
Also, is it possible to generalize for $f_k(x) = \sum_{n_1=1}^{+\infty}\sum_{n_2=1}^{+\infty}\dots \sum_{n_k=1}^{+\infty}e^{-n_1n_2\dots n_k x}$?
 A: Starting with
$$S(x) = \sum_{n\ge 1}\sum_{m\ge 1} e^{-mnx}$$
we have that for $x>0$ and $m,n\ge 1$ the term $e^{-mnx} < 1$, so that we may simplify the inner sum to obtain
$$S(x) = \sum_{n\ge 1} \frac{e^{-nx}}{1-e^{-nx}}.$$
The sum term is harmonic and may be evaluated by inverting its Mellin transform.
Recall the harmonic sum identity
$$\mathfrak{M}\left(\sum_{k\ge 1} \lambda_k g(\mu_k x);s\right) =
\left(\sum_{k\ge 1} \frac{\lambda_k}{\mu_k^s} \right) g^*(s)$$
where $g^*(s)$ is the Mellin transform of $g(x).$
In the present case we have
$$\lambda_k = 1, \quad \mu_k = k \quad \text{and} \quad
g(x) = \frac{e^{-x}}{1-e^{-x}}.$$
We need the Mellin transform $g^*(s)$ of $g(x)$ which is
$$g^*(s) = \int_0^\infty \frac{e^{-x}}{1-e^{-x}} x^{s-1} dx
= \int_0^\infty \sum_{q\ge 1} e^{-qx} x^{s-1} dx
\\=  \sum_{q\ge 1}  \int_0^\infty  e^{-qx} x^{s-1} dx
= \Gamma(s) \sum_{q\ge 1} \frac{1}{q^s} = \Gamma(s) \zeta(s).$$
It follows that the Mellin transform $Q(s)$ of the harmonic sum $S(x)$ is given by
$$Q(s) = \Gamma(s)\zeta(s)^2
\quad\text{because}\quad
\sum_{k\ge 1} \frac{\lambda_k}{\mu_k^s} = \zeta(s).$$
The Mellin inversion integral here is
$$\frac{1}{2\pi i} \int_{3/2-i\infty}^{3/2+i\infty} Q(s)/x^s ds$$
which we evaluate by shifting it to the left for an expansion about zero.
The contribution from the pole at $s=1$ is
$$\mathrm{Res}\left(Q(s)/x^s; s=1\right) = - \frac{\log x}{x} + \frac{\gamma}{x}.$$
The pole at $s=0$ contributes
$$\mathrm{Res}\left(Q(s)/x^s; s=1\right) = \frac{1}{4}.$$
A simple calculation shows that the remaining poles contribute
$$\sum_{q\ge 1} \frac{(-1)^q}{q!} \frac{B_{q+1}^2}{(q+1)^2} x^q.$$
This formula was chosen for simplicity, another formula can be obtained by observing that the Bernoulli numbers vanish at odd integers.
The conclusion is that in a neighborhood of zero, we have
$$S(x) \sim - \frac{\log x}{x} + \frac{\gamma}{x} + \frac{1}{4}
+ \sum_{q\ge 1} \frac{(-1)^q}{q!} \frac{B_{q+1}^2}{(q+1)^2} x^q$$
which is
$$S(x)  \sim - \frac{\log x}{x} + \frac{\gamma}{x} + \frac{1}{4}
-{\frac {1}{144}}\,x-{\frac {1}{86400}}\,{x}^{3}-{\frac {1}{7620480}}\,{x}^{5}-
{\frac {1}{290304000}}\,{x}^{7}
-\cdots.$$
As a particular consequence of this expansion we have that
$$\lim_{x\to 0} \left(S(x)+\frac{\log x}{x}\right)\times x = \gamma,$$
thereby showing the conjecture.
Addendum. The residue of $Q(s)$ at the pole at $s=1$ is calculated from the three series expansions
about $s=1$ $$\zeta(s)^2 = \frac{1}{(s-1)^2} + \frac{2\gamma}{(s-1)}+\cdots,
\quad \Gamma(s) = 1 - \gamma (s-1)+\cdots$$
and $$1/x^s = \frac{1}{x}-\frac{\log x}{x} (s-1)+\cdots.$$
