Find a closed expression for double sum $\sum_{n,m\ge 1} \frac{1}{n m (n^2+m^2)}$

Question

Playing around with double sums related to harmonic numbers and having easily found a closed expression for this sum

$$s_1 = \sum_{n,m\ge 1} \frac{1}{n m (n+m)}=2 \zeta(3)\tag{1}$$

I attempted to find a closed expression for this slightly more complicated sum

$$s_2 = \sum_{n,m\ge 1} \frac{1}{n m (n^2+m^2)}\tag{2}$$

which turned out to be more than I could achieve, and I hope for your help.

Here's what I did so far.

First approach

$$\begin{align}s_2 & = \sum_{n,m\ge 1} \frac{1}{n m (n^2+m^2)} \\ & = \sum_{n,m\ge 1} \frac{1}{n m }\int_{0}^{\infty } e^{-t(n^2+m^2)}\;dt\\ & =\int_{0}^{\infty }\left(\sum_{n,m\ge 1} \frac{1}{n m } e^{-t(n^2+m^2)}\right)\;dt\\ & =\int_{0}^{\infty } f_2(t)^2\;dt \end {align}\tag{3}$$

where

$$f_2(t) = \sum_{n\ge 1}\frac{e^{-t n^2}}{n}\tag{4}$$

Now

$$e^{-n^2 t}=\frac{1}{\sqrt{\pi}}\int_{-\infty }^{\infty } e^{-x^2+i 2 n x\sqrt{t} } \, dx\tag{5}$$

This linearizes $n$ in the exponent, and, exchanging integral and sum, we can do the (conditionally convergent) $n$-sum under the $x$-integral which gives

$$\frac{1}{\sqrt{\pi}}\sum _{n=1}^{\infty } \frac{e^{-x^2+i 2 n x \sqrt{t}}}{\sqrt{\pi } n}=-\frac{e^{-x^2}}{\sqrt{\pi}} \log \left(1-e^{2 i x \sqrt{t} }\right)\tag{6}$$

and

$$f_2(t) = -\frac{1}{\sqrt{\pi}}\int_{-\infty }^{\infty } e^{-x^2} \log \left(1-e^{2 i x \sqrt{t} }\right)\, dx\tag{7}$$

Now in the $x$-integral the imaginary part cancels out due to symmetry and the real part can be simplified so that

$$f_2(t) = -\frac{1}{\sqrt{\pi}}\int_{-\infty }^{\infty } e^{-x^2} \frac{1}{2} \log \left(4 \sin ^2\left(x\sqrt{t} \right)\right)\; dx\tag{8}$$

Alas, here I am stuck, and summing up it seems that I have made a simple looking formula more complicated.

Second approach

Doing just the $n$-sum in $s_2$ gives

$$s_2 = \sum_{m\ge 1} \frac{1}{2 m^3} (H(i m) + H(-i m))=\Re\left(\sum_{m\ge 1} \frac{1}{m^3} H(i m)\right)\tag{2.1}$$

where $H(z)$ is the harmonic number of argument $z$.

Inserting the representation

$$H(z) = \int_{0}^{1} \frac{1-x^z}{1-x}\;dx$$

and performing the $m$-sum gives an interesting integral representation of our double sum

$$s_2 = \left(\int_0^1 \frac{-\zeta (3)+\Re(\text{Li}_3\left(x^i\right))}{x-1} \, dx\right)\tag{2.2}$$

Answer 1

I am stuck almost at the same point.

Trying to get an asymptotics

$$S=\sum_{n= 1}^\infty \sum_{m= 1}^\infty\frac{1}{n m (n^2+m^2)}=\sum_{n= 1}^\infty \frac{H_{i n}+H_{-i n}}{2 n^3}$$

For "large" values of $n$ $$H_{i n}+H_{-i n}=2 (\log (n)+\gamma )+\frac{1}{6 n^2}+\frac{1}{60 n^4}+\frac{1}{126 n^6}+\frac{1}{120 n^8}+\frac{1}{66 n^{10}}+O\left(\frac{1}{n^{12}}\right)$$ is quite good as soon as $n \geq 5$

$$\left( \begin{array}{ccc} n & \text{exact} & \text{approximation} \\ 1 & 1.343731971 & 1.369186020 \\ 2 & 2.583614361 & 2.583605381 \\ 3 & 3.370392751 & 3.370392601 \\ 4 & 3.937503906 & 3.937503902 \\ 5 & 4.380001019 & 4.380001019 \end{array} \right)$$

Using the expansion for $n \geq 5$ does not make any problem leading to $ \color{red}{0.98978105099}37$ instead of $\color{red}{0.9897810509946}$