Double series convergent to $2\zeta(4)$?

So before I start, I've never even attempted to evaluate a double sum before, so there could very well have been an easier way.

$$\sum_{j=1}^\infty\sum_{k=1}^\infty \frac{j^2+jk+k^2}{j^2k^2(j+k)^2} = \sum_{j=1}^\infty\sum_{k=1}^\infty \frac{1}{k^2(j+k)^2} +\sum_{j=1}^\infty\sum_{k=1}^\infty \frac{1}{jk(j+k)^2} + \sum_{j=1}^\infty\sum_{k=1}^\infty \frac{1}{j^2(j+k)^2}= $$

$$ 2\sum_{j=1}^\infty\sum_{k=1}^\infty \frac{1}{k^2(j+k)^2} +\sum_{j=1}^\infty\sum_{k=1}^\infty \frac{1}{jk(j+k)^2}$$

Through partial fraction decomposition the above equals:

$$\sum_{j=1}^\infty\sum_{k=1}^\infty \frac{1}{j^3k} -\frac{1}{j^3(j+k)}-\frac{1}{j^2(j+k)}+2\sum_{j=1}^\infty\sum_{k=1}^\infty \frac{2}{j^3(j+k)}-\frac{2}{j^3k}+\frac{1}{j^2k^2}+\frac{1}{j^2(j+k)^2} $$

Collecting like-terms:

$$3\sum_{j=1}^\infty\sum_{k=1}^\infty \frac{1}{j^3(j+k)}-3\sum_{j=1}^\infty\sum_{k=1}^\infty \frac{1}{j^3k} +\sum_{j=1}^\infty\sum_{k=1}^\infty \frac{1}{j^2(j+k)^2} + 2\sum_{j=1}^\infty\sum_{k=1}^\infty \frac{1}{j^2k^2} = $$

The final sum clearly equals $2\zeta(2)^2$ or $\frac{\pi^4}{18}$. I then evaluate the first two sums by combining them to get:

$$3\sum_{j=1}^\infty\sum_{k=1}^\infty \frac{1}{j^3}(\frac{1}{j+k}-\frac{1}{k}) $$

Interchanging j and k to and averaging the two sums to get:

$$\frac{3}{2}\sum_{j=1}^\infty\sum_{k=1}^\infty \frac{1}{j^3}(\frac{1}{j+k}-\frac{1}{k})+\frac{1}{k^3}(\frac{1}{j+k}-\frac{1}{j}) $$

This can be rewritten as:

$$-\frac{3}{2}\sum_{j=1}^\infty\sum_{k=1}^\infty \frac{1}{j^2k(j+k)}+\frac{1}{k^2j(j+k)}= -\frac{3}{2}\sum_{j=1}^\infty\sum_{k=1}^\infty \frac{1}{j^2k^2} = -\frac{3}{2}\zeta{(2)}^2 = -\frac{\pi^4}{24}$$

So putting that back into the original problem:

$$\sum_{j=1}^\infty\sum_{k=1}^\infty \frac{1}{j^2(j+k)^2} + \frac{\pi^4}{18} -\frac{\pi^4}{48} = \sum_{j=1}^\infty\sum_{k=1}^\infty \frac{1}{j^2(j+k)^2} + \frac{\pi^4}{72} $$

This is all I got to. I couldn't evaluate the other sum the way I did before. Using a calculator there is a very good chance it equals $\frac{\pi^4}{120}$.

Just for fun I was able to write the remaining sum as:

$$\sum_{j=1}^\infty\sum_{k=1}^\infty \frac{1}{j^2(j+k)^2} = \sum_{n=1}^\infty \frac{\zeta(2,n+1)}{n^2}$$

Where $\zeta(x,y)$ is the Hurwitz Zeta Function. Wolfram Alpha was able to calculate the sum as $\frac{\pi^4}{120}$ as desired.

$$S(\infty)=\sum_{j=1}^\infty\,\sum_{k=1}^\infty \frac{(j+k)^2 - jk}{j^2(j+k)^2k^2} = \underbrace{\Big(\sum_{k=1}^\infty \frac{1}{k^2}\Big)^2}_{=\zeta(2)^2} - \underbrace{\sum_{j=1}^\infty\,\sum_{k=1}^\infty \frac{1}{j\,k}\int_0^\infty dt \,t \,e^{-t(j+k)}}_{:=U},$$ where the first step is algebra and the second is use of the Euler representation of the $\Gamma$ function. Interchange sums and integral and sum in terms of $\log$ to find $$U=\int_0^\infty dt \,t \, \log^2(1-e^{-t}) = -\int_0^1 \frac{du}{u} \log\,u \log^2{(1-u)} =$$ $$-\frac{\partial}{\partial s} \frac{\partial^2}{\partial v^2} \int_0^1 u^{s-1} (1-u)^{v-1} \, du \Big\vert_{s=0,v=1}= -\frac{\partial}{\partial s} \frac{\partial^2}{\partial v^2}\frac{\Gamma(s) \Gamma(v)}{\Gamma(s+v)}\Big\vert_{s=0,v=1}$$ where the first step follows from a simple substitution $u=e^{-t}$ and the second is writing the integral in terms of something that is known, the beta integral. Use your favorite CAS to do the partial derivatives to get $U=\pi^4/180.$ Combine with $\zeta(2)^2 = \pi^4/36$ to finish the proof of the OPs hypothesis.

$$ \begin{align} \frac{m^2+m\,n+n^2}{m^2\,(m+n)^2\,n^2}\, &=\frac{m^2+m\,n+n^2\,\color{red}{+m\,n-m\,n}}{m^2\,(m+n)^2\,n^2} \\[2mm] &=\,\frac{(m+n)^2-m\,n}{(m+n)^2\,m^2\,n^2} \\[2mm] &=\,\frac{1}{m^2\,n^2}-\frac{1}{m\,n\,(m+n)^2} \\[2mm] &=\,\frac{1}{m^2\,n^2}-\frac{1}{m^3}\left(\frac{1}{n}-\frac{1}{m+n}-\frac{m}{(m+n)^2}\right) \\[2mm] &=\,\color{brown}{\frac{1}{m^2\,n^2}}\color{green}{-\frac{1}{m^3}\left(\frac{1}{n}-\frac{1}{m+n}\right)}\color{blue}{+\frac{1}{m^2}\frac{1}{(m+n)^2}} \end{align} $$

$$ \begin{align} \color{brown}{\large S_{\small 1}\,} &=\sum_{m=1}^{\infty}\sum_{n=1}^{\infty}\frac{1}{m^2\,n^2}=\sum_{m=1}^{\infty}\frac{1}{m^2}\sum_{n=1}^{\infty}\frac{1}{n^2}=\left(\zeta(2)\right)^2=\color{brown}{\frac{\pi^4}{36}} \\[4mm] \color{green}{\large S_{\small 2}\,} &=\sum_{m=1}^{\infty}\sum_{n=1}^{\infty}\frac{1}{m^3}\left(\frac{1}{n}-\frac{1}{m+n}\right)=\sum_{m=1}^{\infty}\frac{1}{m^3}\sum_{n=1}^{\infty}\left(\frac{1}{n}-\frac{1}{m+n}\right) \\[2mm] &=\sum_{m=1}^{\infty}\frac{1}{m^3}\sum_{n=1}^{\color{red}{m}}\frac{1}{n}=\sum_{m=1}^{\infty}\frac{H_{m}}{m^3}=\frac{5}{4}\zeta(4)=\color{green}{\frac{\pi^4}{72}}\tag{1} \\[4mm] \color{blue}{\large S_{\small 3}\,} &=\sum_{m=1}^{\infty}\sum_{n=1}^{\infty}\frac{1}{m^2}\frac{1}{(m+n)^2}=\sum_{m=1}^{\infty}\frac{1}{m^2}\sum_{n=1}^{\infty}\frac{1}{(m+n)^2} \\[2mm] &=\sum_{m=1}^{\infty}\frac{1}{m^2}\,\sum_{\color{red}{n=m+1}}^{\infty}\,\frac{1}{n^2}=\sum_{m=1}^{\infty}\frac{\psi^{\small (1)}(m+1)}{m^2} \\[2mm] &=\sum_{m=1}^{\infty}\frac{1}{m^2}\left[\zeta(2)-\sum_{k=1}^{m}\frac{1}{k^2}\right]=\left(\zeta(2)\right)^2-\sum_{m=1}^{\infty}\frac{H_{m,2}}{m^2} \\[2mm] &=\left(\zeta(2)\right)^2-\frac{1}{2}\left[\left(\zeta(2)\right)^2+\zeta(4)\right]=\frac{1}{2}\left[\left(\zeta(2)\right)^2-\zeta(4)\right]=\color{blue}{\frac{\pi^4}{120}}\tag{2} \end{align} $$

$$ \color{red}{\Longrightarrow\quad S}\,=S_1-S_2+S_3=\,\color{red}{\frac{\pi^4}{45}} $$

$\,H_m\,\,\,$ : Harmonic Number , $\,\{1\}\,$ : Equation (19) ,$\,\{2\}\,$ : Equation (43)

$\,{\large\psi}^{\small (1)}\,\,$ : Polygamma Function

An alternative approach:

$$ S = \lim_{n\to +\infty}S(n) = \sum_{j,k\geq 1}\frac{1}{j^2 k^2}-\sum_{k,j\geq 1}\frac{1}{jk(j+k)^2}=\zeta(2)^2-\sum_{k,j\geq 1}\int_{0}^{+\infty}\frac{e^{-(j+k)x}}{jk}\,x\,dx $$ leads to $$S = \zeta(2)^2-\int_{0}^{+\infty}x\log^2(1-e^{-x})\,dx=\frac{\pi^4}{36}+\int_{0}^{1}\frac{\log^2(1-x)\log(x)}{x}\,dx$$ or to $$ S = \frac{\pi^4}{36}+\int_{0}^{1}\frac{\log(1-x)}{1-x}\log^2(x)\,dx = \frac{\pi^4}{36}-\sum_{n\geq 1}\frac{2H_n}{(n+1)^3}$$ since $\frac{-\log(1-x)}{1-x}=\sum_{n\geq 1}H_n x^n$ and $\int_{0}^{1}x^n\log^2(x)\,dx = \frac{2}{(n+1)^3}$. Rearranging $$ S = \frac{\pi^4}{36}-2\sum_{n\geq 1}\frac{H_{n}}{n^3}+2\,\zeta(4) = 2\,\zeta(4) = \frac{\pi^4}{45}$$ since the middle term is a linear Euler sum, which can be computed from the Theorem 2.2 here (Flajolet and Salvy, a masterpiece).

$\newcommand{\bbx}[1]{\,\bbox[15px,border:1px groove navy]{\displaystyle{#1}}\,} \newcommand{\braces}[1]{\left\lbrace\,{#1}\,\right\rbrace} \newcommand{\bracks}[1]{\left\lbrack\,{#1}\,\right\rbrack} \newcommand{\dd}{\mathrm{d}} \newcommand{\ds}[1]{\displaystyle{#1}} \newcommand{\expo}[1]{\,\mathrm{e}^{#1}\,} \newcommand{\ic}{\mathrm{i}} \newcommand{\mc}[1]{\mathcal{#1}} \newcommand{\mrm}[1]{\mathrm{#1}} \newcommand{\pars}[1]{\left(\,{#1}\,\right)} \newcommand{\partiald}[3][]{\frac{\partial^{#1} #2}{\partial #3^{#1}}} \newcommand{\root}[2][]{\,\sqrt[#1]{\,{#2}\,}\,} \newcommand{\totald}[3][]{\frac{\mathrm{d}^{#1} #2}{\mathrm{d} #3^{#1}}} \newcommand{\verts}[1]{\left\vert\,{#1}\,\right\vert}$ \begin{align} &\bbox[#ffd,10px]{\ds{\sum_{j = 1}^{\infty}\sum_{k = 1}^{\infty} {j^{2} + jk + k^{2} \over j^{2}\pars{j + k}^{2}k^{2}}}} = \sum_{j = 1}^{\infty}\sum_{k = 1}^{\infty} {j^{2} + jk + k^{2} \over j^{2}k^{2}}\ \overbrace{\bracks{-\int_{0}^{1}\ln\pars{x}x^{j + k - 1}\,\dd x}} ^{\ds{1 \over \pars{j + k}^{2}}} \\[5mm] = &\ -\int_{0}^{1}\ln\pars{x} \pars{\sum_{j = 1}^{\infty}x^{j} \sum_{k = 1}^{\infty}{x^{k} \over k^{2}} + \sum_{j = 1}^{\infty}{x^{j} \over j} \sum_{k = 1}^{\infty}{x^{k} \over k} + \sum_{j = 1}^{\infty}{x^{j} \over j^{2}} \sum_{k = 1}^{\infty}x^{k}}\,{\dd x \over x} \\[5mm] = &\ -\int_{0}^{1}\ln\pars{x} \bracks{2\sum_{j = 1}^{\infty}x^{j} \sum_{k = 1}^{\infty}{x^{k} \over k^{2}} + \pars{\sum_{j = 1}^{\infty}{x^{j} \over j}}^{2}}\,{\dd x \over x} \label{1}\tag{1} \end{align}

Note that $\ds{\sum_{\ell = 1}^{\infty}{x^{\ell} \over \ell^{s}} = \,\mrm{Li}_{s}\pars{x}}$ where $\ds{\mrm{Li}_{s}}$ is the Polylogarithm Function. Moreover, $\ds{\mrm{Li}_{1}\pars{x} = -\ln\pars{1- x}}$, $\ds{\mrm{Li}_{s + 1}\pars{z} = \int_{0}^{z}{\mrm{Li}_{s}\pars{t} \over t}\,\dd t}$ and $\ds{\sum_{j = 1}^{\infty}x^{j} = {x \over 1 - x}}$.

\eqref{1} becomes \begin{align} &\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\bbox[#ffd,10px]{\ds{\sum_{j = 1}^{\infty}\sum_{k = 1}^{\infty} {j^{2} + jk + k^{2} \over j^{2}\pars{j + k}^{2}k^{2}}}} = -2\int_{0}^{1}{\ln\pars{x}\,\mrm{Li}_{2}\pars{x} \over 1 - x}\,\dd x - \int_{0}^{1}{\ln\pars{x}\ln^{2}\pars{1 - x} \over x}\,\dd x \\[5mm] = &\ -2\int_{0}^{1}\overbrace{\ln\pars{1 - x} \over x}^{\ds{-\mrm{Li}_{2}'\pars{x}}}\,\ \mrm{Li}_{2}\pars{1 - x} \,\dd x - \int_{0}^{1}{\ln\pars{x}\ln^{2}\pars{1 - x} \over x}\,\dd x \label{2}\tag{2} \end{align} where we set $\ds{x \mapsto 1 - x}$ in the first integral. With Euler Reflection Formula $\ds{\mrm{Li}_{2}\pars{1 - x} = -\mrm{Li}_{2}\pars{x} + {\pi^{2} \over 6} -\ln\pars{x}\ln\pars{1 - x}}$. \eqref{2} becomes: \begin{align} &\bbox[#ffd,10px]{\ds{\sum_{j = 1}^{\infty}\sum_{k = 1}^{\infty} {j^{2} + jk + k^{2} \over j^{2}\pars{j + k}^{2}k^{2}}}} \\ = &\ -\ \overbrace{\int_{0}^{1}\totald{\mrm{Li}_{2}^{2}\pars{x}}{x} \,\dd x} ^{\ds{\pi^{4} \over 36}}\ +\ {\pi^{2} \over 3}\ \overbrace{\int_{0}^{1}\mrm{Li}_{2}'\pars{x}\,\dd x} ^{\ds{\pi^{2} \over 6}}\ +\ \overbrace{\int_{0}^{1}{\ln\pars{x}\ln^{2}\pars{1 - x} \over x}\,\dd x}^{\ds{-\,{\pi^{4} \over 180}}} \\[5mm] = &\ \bbx{\large{\pi^{4} \over 45}} \end{align}