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}$?
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.$$