Is $\sum_{j=1}^\infty\sum_{n=1}^\infty\left(\frac{e^{-j/n}}{n^2}-\frac{e^{-n/j}}{j^2}\right)=\gamma ?$ A friend proposed the following problem:
$$\sum_{j=1}^\infty\sum_{n=1}^\infty\left(\frac{e^{-j/n}}{n^2}-\frac{e^{-n/j}}{j^2}\right)=\gamma,$$
where $\gamma$ is the Euler-Mascheroni constant.
The result I got is zero and here is what I did:
$$\sum_{j=1}^\infty\sum_{n=1}^\infty\left(\frac{e^{-j/n}}{n^2}-\frac{e^{-n/j}}{j^2}\right)=\sum_{j=1}^\infty\sum_{n=1}^\infty\frac{e^{-j/n}}{n^2}-\sum_{j=1}^\infty\sum_{n=1}^\infty\frac{e^{-n/j}}{j^2}$$
(swap $n$ and $j $ in the first double sum then change the order of summations)
$$=\sum_{j=1}^\infty\sum_{n=1}^\infty\frac{e^{-n/j}}{j^2}-\sum_{j=1}^\infty\sum_{n=1}^\infty\frac{e^{-n/j}}{j^2}=0.$$
Even Mathematica gives zero.
Is it possible that this double sum equals $\gamma$ or my friend could be wrong?
Thanks to @Thomas Andrews for noticing that its not valid to break up the summand as the double sums don't converge and that tells us that we can't always rely on results given by Mathematica or Wolfram.
 A: Here is a relatively elementary approach.
$$
\begin{align}
\lim_{n\to\infty}\sum_{j=1}^{n-1}\sum_{k=n}^\infty\frac{e^{-j/k}}{k^2}
&=\lim_{n\to\infty}\sum_{j=1}^{n-1}\int_n^\infty\frac{e^{-j/x}}{x^2}\,\mathrm{d}x\tag{1a}\\
&=\lim_{n\to\infty}\sum_{j=1}^{n-1}\int_0^{1/n}e^{-jx}\,\mathrm{d}x\tag{1b}\\
&=\lim_{n\to\infty}\sum_{j=1}^{n-1}\frac nj\left(1-e^{-j/n}\right)\frac1n\tag{1c}\\
&=\int_0^1\frac{1-e^{-x}}{x}\,\mathrm{d}x\tag{1d}\\[3pt]
&=-\int_0^1\log(x)\,e^{-x}\,\mathrm{d}x\tag{1e}
\end{align}
$$
Explanation:
$\text{(1a)}$: approximate the sum with an integral:
$\phantom{\text{(1a):}}$ since $x\ge j$, the summand is decreasing and less than $\frac1{n^2}$
$\phantom{\text{(1a):}}$ the step size is $1$, so the error in the approximation is less than $\frac1{n^2}\cdot1$
$\phantom{\text{(1a):}}$ the error in summing $n-1$ of these is less than $\frac1n$
$\text{(1b)}$: substitute $x\mapsto\frac1x$
$\text{(1c)}$: integrate
$\text{(1d)}$: Riemann sum
$\text{(1e)}$: integrate by parts
$$
\begin{align}
\lim_{n\to\infty}\sum_{j=1}^{n-1}\sum_{k=n}^\infty\frac{e^{-k/j}}{j^2}
&=\lim_{n\to\infty}\sum_{j=1}^{n-1}\frac{e^{-n/j}}{j^2\left(1-e^{-1/j}\right)}\tag{2a}\\
&=\lim_{n\to\infty}\sum_{j=1}^{n-1}\frac{e^{-n/j}}{j/n}\frac1n\tag{2b}\\
&=\int_0^1\frac{e^{-1/x}}{x}\,\mathrm{d}x\tag{2c}\\[3pt]
&=\int_1^\infty\frac{e^{-x}}{x}\,\mathrm{d}x\tag{2d}\\[3pt]
&=\int_1^\infty\log(x)\,e^{-x}\,\mathrm{d}x\tag{2e}
\end{align}
$$
Explanation:
$\text{(2a)}$: sum the geometric series
$\text{(2b)}$: since $j-\frac12\le j^2\left(1-e^{-1/j}\right)\le j$, the sum of the errors is less than
$\phantom{\text{(2b):}}$ $\sum\limits_{j=1}^{n-1}\frac{e^{-n/j}}{j(2j-1)}\le\frac1{n^2}\sum\limits_{j=1}^{n-1}\frac{e^{-n/j}}{(j/n)^2}\le\frac4{e^2n}$ because $\frac{e^{-n/j}}{(j/n)^2}=u^2e^{-u}\le\frac4{e^2}$
$\text{(2c)}$: Riemann sum
$\text{(2d)}$: substitute $x\mapsto\frac1x$
$\text{(2e)}$: integrate by parts
$$
\begin{align}
\sum_{j=1}^\infty\sum_{k=1}^\infty\left(\frac{e^{-j/k}}{k^2}-\frac{e^{-k/j}}{j^2}\right)
&=\lim_{n\to\infty}\sum_{j=1}^{n-1}\sum_{k=1}^\infty\left(\frac{e^{-j/k}}{k^2}-\frac{e^{-k/j}}{j^2}\right)\tag{3a}\\
&=\lim_{n\to\infty}\sum_{j=1}^{n-1}\sum_{k=n}^\infty\left(\frac{e^{-j/k}}{k^2}-\frac{e^{-k/j}}{j^2}\right)\tag{3b}\\
&=-\int_0^\infty\log(x)\,e^{-x}\,\mathrm{d}x\tag{3c}\\[9pt]
&=\gamma\tag{3d}
\end{align}
$$
Explanation:
$\text{(3a)}$: write infinite sum as a limit
$\text{(3b)}$: subtract $\sum\limits_{j=1}^{n-1}\sum\limits_{k=1}^{n-1}\left(\frac{e^{-j/k}}{k^2}-\frac{e^{-k/j}}{j^2}\right)=0$
$\text{(3c)}$: apply $(1)$ and $(2)$
$\text{(3d)}$: apply this answer
A: Here's the solution I have. Take $1<c<2$ and consider the integral $$I=\frac1{2\pi i}\int_{c-i\infty}^{c+i\infty}\Gamma(s)\zeta(s)\zeta(2-s)\,ds.$$
Let's compute it in two ways.
First, $I$ is the sum of the residues of the integrand at all integer $s\leqslant1$ (consider integration along the boundary of $[-N-1/2,c]+i[-R,R]$ with $R>0$ and $N$ a positive integer, and take $R,N\to\infty$; the integral along the "right side of the rectangle" tends to $I$, and the remaining integrals vanish). The residue at $s=1$ (a double pole) may be computed using
\begin{align*}
\Gamma(s)&=1-\gamma(s-1)+o(s-1),\\
\zeta(s)&=(s-1)^{-1}+\gamma+o(1),\\
\zeta(2-s)&=-(s-1)^{-1}+\gamma+o(1),
\end{align*}
hence $\Gamma(s)\zeta(s)\zeta(2-s)=(s-1)^{-2}\big({-1}+\gamma(s-1)+o(s-1)\big)$ and the residue equals $\gamma$. And, for a nonnegative integer $n$, we have $$\newcommand{\res}{\operatorname*{Res}}\res_{s=-n}\Gamma(s)\zeta(s)\zeta(2-s)=\frac{(-1)^n}{n!}\zeta(-n)\zeta(n+2)=\frac{B_{n+1}}{(n+1)!}\zeta(n+2)$$ so that, using the generating function of the Bernoulli numbers,
\begin{align*}
I&=\gamma+\sum_{n=0}^\infty\frac{B_{n+1}}{(n+1)!}\zeta(n+2)=\gamma+\sum_{n=1}^\infty\frac{B_n}{n!}\sum_{j=1}^\infty\frac1{j^{n+1}}
\\&=\gamma+\sum_{j=1}^\infty\frac1j\left(\frac{1/j}{e^{1/j}-1}-1\right)=\gamma+\sum_{j=1}^\infty\left(-\frac1j+\frac1{j^2}\sum_{n=1}^\infty e^{-n/j}\right).
\end{align*}
Second (the implied uniform convergence is not hard to see; recall $1<c<2$), $$I=\sum_{j=1}^\infty I_j,\quad I_j:=\frac1{2\pi i}\int_{c-i\infty}^{c+i\infty}\Gamma(s)j^{-s}\zeta(2-s)\,ds.$$ Now $I_j$ may be computed the same way. The residue of the integrand at $s=1$ equals $-1/j$, and the residue at $s=-k$ (for a nonnegative integer $k$) is $\zeta(k+2)(-j)^k/k!$, so that $$I=\sum_{j=1}^\infty\left(-\frac1j+\sum_{k=0}^\infty\frac{(-j)^k}{k!}\sum_{n=1}^\infty\frac1{n^{k+2}}\right)=\sum_{j=1}^\infty\left(-\frac1j+\sum_{n=1}^\infty\frac{e^{-j/n}}{n^2}\right).$$
Comparing these two expressions for $I$, we get the expected result.
