Using a computer I found the double sum

$$S(n)= \sum_{j=1}^n\sum_{k=1}^n \frac{j^2 + jk + k^2}{j^2(j+k)^2k^2}$$ has values

$$S(10) \quad\quad= 1.881427206538142 \\ S(1000) \quad= 2.161366028875634 \\S(100000) = 2.164613524212465\\$$

As a guess I compared with fractions $\pi^p/q$ where $p,q$ are positive integers and it appears

$$\lim_{n \to \infty} S(n) = \frac{\pi^4}{45} = 2\zeta(4) \approx 2.164646467422276 $$

I'd be interested in seeing a proof if true.

  • $\begingroup$ I think you mean $\frac{\pi^4}{45}$. Also, can I ask how you came up with this? I would be shocked if it's true, but still intriguing. $\endgroup$
    – Tom Himler
    Aug 13, 2018 at 23:56
  • $\begingroup$ @TomHimler: Thank you - typo fixed. $\endgroup$
    – WoodWorker
    Aug 13, 2018 at 23:59
  • 3
    $\begingroup$ More of an exercise in coding just to see what happens with summands that have symmetric structure like this.. $\endgroup$
    – WoodWorker
    Aug 14, 2018 at 0:02
  • $\begingroup$ Note that there are inverse symbolic calculators, like this one: isc.carma.newcastle.edu.au $\endgroup$ Aug 14, 2018 at 20:44

5 Answers 5


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).

  • $\begingroup$ essentially skbmoore's $\endgroup$
    – skbmoore
    Aug 14, 2018 at 15:14
  • $\begingroup$ @skbmoore: the intro is the same, but while you use derivatives of the Beta function I use linear Euler sums. $\endgroup$ Aug 14, 2018 at 15:19
  • 1
    $\begingroup$ Agree that Euler sums are a way to get there, but probably less well-known than derivatives of the $\Gamma$ function. But it doesn't hurt to get another method's name out there. $\endgroup$
    – skbmoore
    Aug 14, 2018 at 17:14

$\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}


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.