This is a personal research that came to an end , since the results were not those which were being anticipated. I was unable to come up with a solution therefore I post the topic here:
Prove (it certainly holds because it was checked on a computer) that the following identity holds:
$$\sum_{n=1}^{\infty}\sum_{k=1}^{\infty}\frac{1}{\left ( n^2+k^2 \right )^2}=\zeta(2)\sum_{n=1}^{\infty}(-1)^{n-1}\frac{1}{\left ( 2n-1 \right )^2}-\zeta(4)$$
wheras $\zeta$ represents the zeta function defined as $\displaystyle \zeta(s)=\sum_{n=1}^{\infty}\frac{1}{n^s}, \;\; \mathfrak{Re}(s)>1 $. Of course both values of $\zeta$ appearing here are known. For complete of sakeness i quote them :
$$\zeta(2)=\frac{\pi^2}{6}, \;\; \zeta(4)=\frac{\pi^4}{90}$$
Now, one can also see (not trivial ) that:
$$\sum_{n=1}^{\infty}\frac{1}{n^2 \sinh^2 \pi n}=\frac{4}{\pi^2}\sum_{n=1}^{\infty}\sum_{k=1}^{\infty}\frac{1}{\left ( n^2+k^2 \right )^2} -\frac{\pi^2}{60}$$
This equation also holds (checked by computer). The following result was extracted by using the known formulae $\displaystyle \sum_{n=-\infty}^{\infty}\frac{1}{z+n}=\frac{\pi}{\tan \pi z}$ and the known (?) Fourier series:
$$\displaystyle \frac{1}{\sinh^2 \pi z}=\frac{1}{\pi^2 z}+\frac{4z^2}{\pi^2}\sum_{k=1}^{\infty}\frac{1}{\left ( z^2+k^2 \right )^2}-\frac{2}{\pi^2}\sum_{k=1}^{\infty}\frac{1}{z^2+k^2} $$
Of course the last sum at the last equation can easily be computed via residues. What remains now is the proof for the first equation in the post. No-one can guarantee that is going to be an easy task.
Some comments:
1. I came across the identity on a book. I checked the validity with a computer and yes it does hold.
2. I have ckecked enough books to see if there is in there somewhere , but unfortunately it was not. So, I speculate that it is not that famous.
3. It can also be linked with other sums (single or double). Unfortunately I don't have my papers in front of me in order to write them down. So, i think it is an interseting identity.
I would appreciate your help.
Best Answer
Let's consider the double sum over all $\,(n,k)\in\mathbb{Z}^2\,$ except the origin $(0,0)$ : $$\tag{1}S(s):=\sum_{(n,k)\neq (0,0)}\frac{1}{\left ( n^2+k^2 \right )^s}=4\left(\sum_{n=1}^{\infty}\sum_{k=1}^{\infty}\frac{1}{\left ( n^2+k^2 \right )^s}+\zeta(2\,s)\right)$$ This is a "Lorenz–Hardy sum" (ref. Lorenz L. $1871$ "Bidrag tiltalienes theori" and Hardy G. H. $1919$ "On some definite integral considered by Mellin" in his Collected papers vol $7$).
Analytical method
The Mellin transform of a function $f$ is defined by : $$\tag{2}\{\mathcal{M}f\}(s):=\int_0^\infty t^{s-1}f(t)\;dt$$ applied to $\;f:t\mapsto e^{\large{-mt}}\;$ and from the definition of the Gamma function gives : $$\tag{3}\frac{\Gamma(s)}{\left(m \right)^s}=\int_0^\infty t^{s-1}\,e^{-m\,t}\;dt$$ Supposing that $\Re(s)>1$ we may rewrite our double-sum as : \begin{align} \Gamma(s)\,S(s)&=\sum_{(n,k)\neq (0,0)}\frac{\Gamma(s)}{\left ( n^2+k^2 \right )^s}\\ &=\sum_{(n,k)\neq (0,0)}\int_0^\infty t^{s-1}\,e^{-(n^2+k^2)\,t}\;dt\\ &=\int_0^\infty t^{s-1}\sum_{(n,k)\neq (0,0)}\,e^{-(n^2+k^2)\,t}\;dt\\ &=\int_0^\infty t^{s-1}\left(\sum_{n\in\mathbb{Z}}\,e^{-n^2\,t}\sum_{k\in\mathbb{Z}}\,e^{-k^2\,t}-1\right)\;dt\\ \tag{4}&=\int_0^\infty t^{s-1}\left(\theta_3(0,e^{-t})^2-1\right)\;dt\\ \end{align} i.e. the Mellin transform of $\,f:t\mapsto \theta_3(0,e^{-t})^2-1\,$ where the Jacobi theta function $\theta_3$ is defined by (in this answer we will implicitly suppose that $z=0$ and $q=e^{-t}$) : $$\tag{5}\theta_3(z,q):=\sum_{n=-\infty}^\infty q^{n^2}e^{2niz}$$
In $1829$ Jacobi published his deep elliptic functions book "Fundamenta nova theoriae functionum ellipticarum" where he obtained numerous identities including the formula $4$) of page $103$ : $$\tag{6}\theta_3(0,q)^2=\frac {2\,K}{\pi}=1+4\sum_{n=1}^\infty \frac{q^n}{1+q^{2n}}$$ ($K$ is the Complete elliptic integral of the first kind but we won't use it here)
The series at the right of $(6)$ is a Lambert series. The expansion of the denominator in power series and substitution of the double-sum in $(4)$ gives :
\begin{align} \Gamma(s)\,S(s)&=4\int_0^\infty t^{s-1}\sum_{n=1}^\infty \sum_{m=0}^\infty (-1)^m q^n\,q^{2mn}\;dt\\ &=4\sum_{n=1}^\infty \sum_{m=0}^\infty (-1)^m \int_0^\infty t^{s-1}e^{-n(2m+1)t}\;dt\\ &=4\,\sum_{n=1}^\infty \sum_{m=0}^\infty (-1)^m \frac{\Gamma(s)}{(n(2m+1))^s}\\ &=4\,\Gamma(s)\sum_{n=1}^\infty \frac 1{n^s}\sum_{m=0}^\infty \frac{(-1)^m}{(2m+1)^s}\\ &=4\,\Gamma(s)\,\zeta(s)\,\beta(s)\\ \\ &\text{so that }\\ \\ \tag{7}\sum_{(n,k)\neq (0,0)}\frac{1}{\left ( n^2+k^2 \right )^s}&=4\,\zeta(s)\,\beta(s),\quad\Re(s)>1\\ &\text{and}\\ \tag{8}\sum_{n=1}^{\infty}\sum_{k=1}^{\infty}\frac{1}{\left ( n^2+k^2 \right )^s}&=\zeta(s)\,\beta(s)-\zeta(2s),\quad\Re(s)>1\\ \\\end{align}
$\qquad\qquad\qquad$with $\beta$ the Dirichlet beta function.
This powerful method allows to deduce many closed forms of lattice sums from the known theta identities.
I'll try to add one (or more) alternative derivations when time allows...
References :