# Integrals of integer powers of dilogarithm function

I'm interested in evaluating integrals of positive integer powers of the dilogarithm function. I'd like to see the general case tackled if possible, or barring that then as many particular cases as possible.

Problem. For each $n\in\mathbb{N}$, evaluate the definite integral $\mathcal{I}_n$: $$\mathcal{I}_n:=\int_{0}^{1}\left[\operatorname{Li}_2{(x)}\right]^n\,\mathrm{d}x=\,???$$

The first two cases aren't that difficult to evaluate, but for reference I'll just state their values without proof:

\begin{align} \int_{0}^{1}\operatorname{Li}_2{(x)}\,\mathrm{d}x &=\zeta{(2)}-1\\ &=0.6449340668482264\dots; \end{align}

and

\begin{align} \int_{0}^{1}\operatorname{Li}_2{(x)}^2\,\mathrm{d}x &=\frac52\zeta{(4)}-4\zeta{(3)}-2\zeta{(2)}+6\\ &=0.6077123379430154\dots. \end{align}

The next case, $n=3$, is much harder than the previous two. However, by combining my answer to this question with Omran Kouba's answer to the same question, and assuming both answers are correct, I was able to infer that

$$\zeta{(5)}+6\zeta{(3)}+\frac{\pi^4}{15}-15=\frac{1}{3!}\int_{0}^{1}\left[\operatorname{Li}_2{(x)}-\zeta{(2)}\right]^3\,\mathrm{d}x.$$

Then by using the Binomial Theorem to expand the integrand on the RHS and by using the previous values of $\mathcal{I}_n$ for $n<3$, I was able to solve for the integral $\mathcal{I}_3$ and find the following value:

\begin{align} \mathcal{I}_3 &=\int_{0}^{1}\operatorname{Li}_2{(x)}^3\,\mathrm{d}x\\ &=6\zeta{(5)}+36\zeta{(3)}+\frac{\pi^6}{216}+\frac{19\pi^4}{60}+3\pi^2-2\pi^2\zeta{(3)}-90\\ &=0.6738641012555397\dots. \end{align}

But surely there is a less circuitous way to prove this. Can anybody offer a more direct proof that the integral $\mathcal{I}_3$ has the conjectured value indicated above?

Prove: $$\int_{0}^{1}\operatorname{Li}_2{(x)}^3\,\mathrm{d}x\stackrel{?}{=}6\zeta{(5)}+36\zeta{(3)}+\frac{\pi^6}{216}+\frac{19\pi^4}{60}+3\pi^2-2\pi^2\zeta{(3)}-90.$$

And finally, what of the $n>3$ cases? What is the largest $n$ for which $\mathcal{I}_n$ can be evaluated? Can $\mathcal{I}_4$ be evaluated? At the bottom of Omran Kouba's response to the question I referred to above, he says that, to his knowledge, the values of the integral of integer powers of the dilogarithm for powers larger than $2$ are not known. Clearly, for that statement to be true then, at the very least, it would have to be amended to reflect the fact that the value of the integral of the third power of the dilogarithm is indeed known (assuming I'm not somehow the first person to know it!). But are the values for $n>3$ really unknown?

Update: As SuperAbound was able to show, finding $\int_{0}^{1}\operatorname{Li}_2{(x)}^3\,\mathrm{d}x$ via repreated integration by parts actually isn't that difficult at all! For comparison, let's see how much harder finding $\int_{0}^{1}\operatorname{Li}_2{(x)}^4\,\mathrm{d}x$ is if we try to recycle the same method.

To evaluate $\int_0^1\operatorname{Li}_2{(x)}^4\mathrm{d}x$, integrate by parts using

$$\frac{d}{dx}\operatorname{Li}_2{(x)}^{4}=-\frac{4\ln{(1-x)}\operatorname{Li}_2{(x)}^3}{x}.$$

Then,

\begin{align} \mathcal{I}_4 &=\int_{0}^{1}\operatorname{Li}_2{(x)}^4\,\mathrm{d}x\\ &=\left[x\operatorname{Li}_2{(x)}^4\right]_{0}^{1}-\int_{0}^{1}x\left(-\frac{4\ln{(1-x)}\operatorname{Li}_2{(x)}^3}{x}\right)\,\mathrm{d}x\\ &=\operatorname{Li}_2{(1)}^4+4\int_{0}^{1}\ln{(1-x)}\operatorname{Li}_2{(x)}^3\,\mathrm{d}x\\ &=\zeta{(2)}^4+4\int_{0}^{1}\ln{(1-x)}\operatorname{Li}_2{(x)}^3\,\mathrm{d}x. \end{align}

Next, evaluate $\int_{0}^{1}\ln{(1-x)}\operatorname{Li}_2{(x)}^3\,\mathrm{d}x$ by integrating by parts again using

$$\int\ln{(1-x)}\,\mathrm{d}x=(x-1)\ln{(1-x)}-x+constant,$$

$$\frac{d}{dx}\operatorname{Li}_2{(x)}^3=-\frac{3\ln{(1-x)}\operatorname{Li}_2{(x)}^2}{x}.$$

Then,

\begin{align} \int_{0}^{1}\ln{(1-x)}\operatorname{Li}_2{(x)}^3\,\mathrm{d}x &=\left[\left((x-1)\ln{(1-x)}-x\right)\operatorname{Li}_2{(x)}^3\right]_{0}^{1}\\ &~~~~~-\int_{0}^{1}\left((x-1)\ln{(1-x)}-x\right)\left(-\frac{3\ln{(1-x)}\operatorname{Li}_2{(x)}^2}{x}\right)\,\mathrm{d}x\\ &=-\operatorname{Li}_2{(1)}^3\\ &~~~~~-3\int_{0}^{1}\left[\frac{\ln^2{(1-x)}}{x}-\ln^2{(1-x)}+\ln{(1-x)}\right]\operatorname{Li}_2{(x)}^2\,\mathrm{d}x\\ &=-\zeta{(2)}^3-3\int_{0}^{1}\frac{\ln^2{(1-x)}\operatorname{Li}_2{(x)}^2}{x}\,\mathrm{d}x\\ &~~~~~+3\int_{0}^{1}\ln^2{(1-x)}\operatorname{Li}_2{(x)}^2\,\mathrm{d}x-3\int_{0}^{1}\ln{(1-x)}\operatorname{Li}_2{(x)}^2\,\mathrm{d}x. \end{align}

The third integral in the last line above has the value,

$$\int_{0}^{1}\ln{(1-x)}\operatorname{Li}_2{(x)}^2\,\mathrm{d}x=2\zeta{(5)}+\frac{19}{2}\zeta{(4)}+12\zeta{(3)}+6\zeta{(2)}-4\zeta{(3)}\zeta{(2)}-30.$$

(to be continued...)

• Because I couldn't resist evaluating for an arbitrary complex exponent, using owen.maresh.info/intdilog.py, I made owen.maresh.info/dilogpowers.png ... it might be worthwhile to figure out how to analytically complete that to $\mathbb{C}$, as it's evident where the software had a problem. Aug 29, 2014 at 12:24
• @David H Can we generalize your method to show a closed-form of $\int_{0}^{1}\operatorname{Li}_3^3{(x)}\,dx$ ? By the way in general results to $\int_{0}^{1}\operatorname{Li}_p^2{(x)}\,dx$ comes by using methods from this paper. Oct 19, 2014 at 10:45
• @user153012 Such a generalization is not at all apparent to me... Oct 19, 2014 at 11:19
• @DavidH did you ever found a complete closed form for the $\text{Li}_2^4$ case? Sep 10, 2015 at 12:01

There is a reasonably easy way to do these integrals. Above all, it is systematic. Each dilogarithm can be written as an iterated integral $$L(x) = \int_0^x \frac{dt}{t}\int_0^t\frac{du}{1-u}.$$ Every power of a dilogarithm can then be expanded into a long product of integrals, giving, as an example: $$\int_0^1 L(x)^2\,dx = \int_0^1dx \int_0^x \frac{dt_1}{t_1} \int_0^{t_1}\frac{du_1}{1-u_1} \int_0^x \frac{dt_2}{t_2} \int_0^{t_2}\frac{du_2}{1-u_2}.$$ The reason this is easy is that one integral representation of multiple zeta values has the form $$\int_0^1 \frac{dw_1}{w_1-b_1}\int_0^{w_1}\frac{dw_2}{w_2-b_2}\cdots\int_0^{w_{n-1}}\frac{dw_n}{w_n-b_n},$$ where $$b_j\in\{0,1\}$$.

In the case of $$\int L(x)^n$$, the domain of the integral is $$0< t_1,\ldots,t_n < x < 1, \qquad 0 < u_j < t_j,$$ which can be rearranged into a union of simplices of the form $$0 (where $$w_j$$'s are all of $$t_j$$'s and $$u_j$$'s in some order), and the integral over each simplex has the form $$\int_0^1dx\,\Omega,$$ where $$\Omega$$ is an iterated integral expression, starting with either $$\int_0^x \frac{dw_1}{w_1}$$ or $$\int_0^x \frac{dw_1}{1-w_1}$$.

It is easy to interchange the order of integration of $$x$$ and $$w_1$$, integrate over $$x$$, and be left with a very similar integral over another simplex with only the variables $$w_2,\ldots,w_{2n}$$ present.

Implementing this recursive procedure gives the following expressions for $$I$$: $$\begin{eqnarray} I_1 &=& \zeta(2) - 1 \\I_2 &=& -2\zeta(2) - 4\zeta(2,1) + 2\zeta(2,2) + 4\zeta(3,1) + 6 \\I_3 &=& 18\zeta(2) + 36\zeta(2,1) + 36\zeta(2,1,1) - 12\zeta(2,1,2) - 6\zeta(2,2) - 24\zeta(2,2,1) + 6\zeta(2,2,2) + 12\zeta(2,3,1) - 12\zeta(3,1) - 36\zeta(3,1,1) + 12\zeta(3,1,2) + 24\zeta(3,2,1) + 36\zeta(4,1,1) - 90 \\I_4 &=& -360\zeta(2) - 720\zeta(2,1) - 864\zeta(2,1,1) - 576\zeta(2,1,1,1) + 144\zeta(2,1,1,2) + 144\zeta(2,1,2) + 288\zeta(2,1,2,1) - 48\zeta(2,1,2,2) - 96\zeta(2,1,3,1) + 72\zeta(2,2) + 288\zeta(2,2,1) + 432\zeta(2,2,1,1) - 96\zeta(2,2,1,2) - 24\zeta(2,2,2) - 192\zeta(2,2,2,1) + 24\zeta(2,2,2,2) + 48\zeta(2,2,3,1) - 48\zeta(2,3,1) - 288\zeta(2,3,1,1) + 48\zeta(2,3,1,2) + 96\zeta(2,3,2,1) + 144\zeta(2,4,1,1) + 144\zeta(3,1) + 432\zeta(3,1,1) + 576\zeta(3,1,1,1) - 144\zeta(3,1,1,2) - 48\zeta(3,1,2) - 288\zeta(3,1,2,1) + 48\zeta(3,1,2,2) + 96\zeta(3,1,3,1) - 96\zeta(3,2,1) - 432\zeta(3,2,1,1) + 96\zeta(3,2,1,2) + 192\zeta(3,2,2,1) + 288\zeta(3,3,1,1) - 144\zeta(4,1,1) - 576\zeta(4,1,1,1) + 144\zeta(4,1,1,2) + 288\zeta(4,1,2,1) + 432\zeta(4,2,1,1) + 576\zeta(5,1,1,1) + 2520 \end{eqnarray}$$

I computed these values up to $$n=8$$, but the expressions are long ($$I_5$$ has 428 terms), and the integers involved become prohibitively large, so I didn't include them here.

For small $$n$$, at least up to $$n=4$$, I found numerically that the above expressions can be written in terms of ordinary zeta values, as follows: $$\begin{eqnarray}\def\tfrac#1#2{{\textstyle\frac{#1}{#2}}} I_1 &=& \zeta(2) - 1 \\I_2 &=& -4 \zeta (3)-2 \zeta (2)+\tfrac{5}{2} \zeta (4)+6 \\I_3 &=& -12 \zeta (3) \zeta (2)+6 \zeta (5)+36 \zeta (3)+18 \zeta (2)+\tfrac{57}{2} \zeta (4)+\tfrac{35}{8} \zeta (6)-90 \\I_4 &=& 24 \zeta (5) \zeta (2)+144 \zeta (3) \zeta (2)+84 \zeta (3) \zeta (4)-195 \zeta (7)-648 \zeta (5)+144 \zeta (3)^2-720 \zeta (3)-360 \zeta (2)-774 \zeta (4)-\frac{11}{2} \zeta (6)+\frac{175}{24} \zeta (8)+2520 \end{eqnarray}$$ I checked $$n\leq3$$ using known identities satisfied by MZVs, such as those here. For $$I_4$$, while this is not a proof that the expression is correct, this entire procedure reduces the problem of evaluating the integrals $$I_n$$ to the problem of applying known identities to multiple zeta values, which is more straightforward.

A quick integer relation search for a closed form for $$I_5$$ gave me nothing.

There may or may not be a closed form for $$I_5$$ (I believe this is not a settled question), but there is no systematic way to reduce multiple zeta values to ordinary zeta values. For example, weight 10, which is what $$I_5$$ requires, is the first weight at which known relations between depth-two zeta values are insufficient to reduce them to ordinary zeta values ("Euler Sums and Contour Integral Representations" by Flajolet and Salvy). Because here the depth can be up to five, this problem is harder.

Update (by editor after 6 years): Settled now. By further computation one have

$$\small\int_0^1 \text{Li}_2(x){}^5 \, dx=1440 \zeta(6,2)+240 \zeta (3)^3+120 \pi ^2 \zeta (3)^2-4320 \zeta (3)^2-\frac{23 \pi ^6 \zeta (3)}{189}-110 \pi ^4 \zeta (3)-600 \pi ^2 \zeta (3)-9360 \zeta (5) \zeta (3)+25200 \zeta (3)-\frac{127 \pi ^4 \zeta (5)}{3}-540 \pi ^2 \zeta (5)+30600 \zeta (5)-\frac{325 \pi ^2 \zeta (7)}{2}+14445 \zeta (7)+\frac{16495 \zeta (9)}{3}+\frac{\pi ^{10}}{7776}+\frac{51137 \pi ^8}{45360}+\frac{1907 \pi ^6}{126}+335 \pi ^4+2100 \pi ^2-113400$$

Here $$\zeta(6,2)$$ is irreducible. In fact all $$I_n$$ can be reduced to MZVs. The generalized algorithm I used also gives the proof to closed-form of $$I_4$$ above.

• Can you share $I_5$ via mathb.in? Not necessary in LaTeX, the Mathematica output or an output by the CAS you've used is enough. I'm just curious. Sep 24, 2014 at 22:32
• @user153012 Here you go. Sep 24, 2014 at 22:41

Allow me to present another approach. Though lengthier and more tedious, this method is independent of results derived from previous answers. The main idea is to expand ${\rm Li}_2(x)$ as a series and integrate term by term.

For $n=1$, \begin{align} \mathcal{I}_1 &=\int^1_0{\rm Li}_2(x) \ {\rm d}x\\ &=\int^1_0\sum^\infty_{j=1}\frac{x^j}{j^2}{\rm d}x\\ &=\sum^\infty_{j=1}\frac{1}{j^2(j+1)}\\ &=\sum^\infty_{j=1}\frac{1}{j^2}-\frac{1}{j}+\frac{1}{j+1}\\ \end{align} We see that \begin{align} \sum^\infty_{k=1}\frac{1}{k^2} &=\frac{\pi^2}{6}\\ \sum^\infty_{k=1}\left(\frac{1}{k+1}-\frac{1}{k}\right) &=\sum^\infty_{k=2}\frac{1}{k}-\sum^\infty_{k=1}\frac{1}{k}\\ &=-1 \end{align} Therefore $$\mathcal{I}_1=\frac{\pi^2}{6}-1$$

For $n=2$, \begin{align} \mathcal{I}_2 &=\int^1_0{\rm Li}^2_2(x) \ {\rm d}x\\ &=\int^1_0\sum^\infty_{j=1}\frac{x^j}{j^2}\sum^\infty_{k=1}\frac{x^k}{k^2}{\rm d}x\\ &=\sum^\infty_{j=1}\frac{1}{j^2}\sum^\infty_{k=1}\frac{1}{k^2(k+j+1)}\\ &=\sum^\infty_{j=1}\frac{1}{j^2(j+1)}\sum^\infty_{k=1}\frac{1}{k^2}-\frac{1}{k(j+1)}+\frac{1}{(k+j+1)(j+1)}\\ &=\frac{\pi^2}{6}\left(\frac{\pi^2}{6}-1\right)+\sum^\infty_{j=1}\frac{1}{j^2(j+1)^2}\sum^\infty_{k=1}\frac{1}{k+j+1}-\frac{1}{k}\\ &=\frac{\pi^4}{36}-\frac{\pi^2}{6}+\color{red}{\sum^\infty_{j=1}\frac{2H_{j+1}}{j}}\color{green}{-\sum^\infty_{j=1}\frac{2H_{j+1}}{j+1}}\color{blue}{-\sum^\infty_{j=1}\frac{H_{j+1}}{j^2}}\color{purple}{-\sum^\infty_{j=1}\frac{H_{j+1}}{(j+1)^2}}\\ \end{align} It is trivial to show that \begin{align} \sum^\infty_{j=1}\frac{H_j}{j^2} &=2\zeta(3) \end{align} We see that \begin{align} \color{red}{\sum^\infty_{j=1}\frac{2H_{j+1}}{j}} &=\sum^\infty_{j=1}\frac{2H_j}{j}+\frac{2}{j(j+1)}\\ &=\color{red}{\sum^\infty_{j=1}\frac{2H_j}{j}+2}\\ \color{green}{-\sum^\infty_{j=1}\frac{2H_{j+1}}{j+1}} &=\color{green}{-\sum^\infty_{j=1}\frac{2H_j}{j}+2}\\ \color{blue}{-\sum^\infty_{j=1}\frac{H_{j+1}}{j^2}} &=-\sum^\infty_{j=1}\frac{H_{j}}{j^2}-\sum^\infty_{j=1}\frac{1}{j^2(j+1)}\\ &=\color{blue}{-2\zeta(3)+1-\frac{\pi^2}{6}}\\ -\color{purple}{\sum^\infty_{j=1}\frac{H_{j+1}}{(j+1)^2}} &=-\sum^\infty_{j=1}\frac{H_{j}}{j^2}+1\\ &=\color{\purple}{-2\zeta(3)+1} \end{align} Therefore $$\mathcal{I}_2=\frac{\pi^4}{36}-4\zeta(3)-\frac{\pi^2}{3}+6$$

For $n=3$, \begin{align} \mathcal{I}_3 &=\int^1_0{\rm Li}^3_2(x) \ {\rm d}x\\ &=\int^1_0\sum^\infty_{j=1}\frac{x^j}{j^2}\sum^\infty_{k=1}\frac{x^k}{k^2}\sum^\infty_{m=1}\frac{x^m}{m^2}{\rm d}x\\ &=\sum^\infty_{j=1}\frac{1}{j^2}\sum^\infty_{k=1}\frac{1}{k^2}\sum^\infty_{m=1}\frac{1}{m^2(m+j+k+1)}\\ &=\sum^\infty_{j=1}\frac{1}{j^2}\sum^\infty_{k=1}\frac{1}{k^2(j+k+1)}\sum^\infty_{m=1}\frac{1}{m^2}-\frac{1}{m(j+k+1)}+\frac{1}{(m+j+k+1)(j+k+1)}\\ &=\frac{\pi^2}{6}\mathcal{I}_2-\sum^\infty_{j=1}\frac{1}{j^2}\sum^\infty_{k=1}\frac{H_{k+j+1}}{k^2(k+j+1)^2}\\ &=\frac{\pi^2}{6}\mathcal{I}_2-\sum^\infty_{j=1}\frac{1}{j^2(j+1)^2}\sum^\infty_{k=1}\left(\frac{H_{k+j+1}}{k^2}+\frac{H_{k+j+1}}{(k+j+1)^2}\right)+2\sum^\infty_{j=1}\frac{1}{j^2(j+1)^3}\sum^\infty_{k=1}\left(\frac{H_{k+j+1}}{k}-\frac{H_{k+j+1}}{k+j+1}\right) \end{align} The first inner sum is \begin{align} \sum^\infty_{k=1}\frac{H_{k+j+1}}{k^2} &=2\zeta(3)+\sum^\infty_{k=1}\sum^{j+1}_{m=1}\frac{1}{k^2(k+m)}\\ &=2\zeta(3)+\sum^\infty_{k=1}\sum^{j+1}_{m=1}\frac{1}{mk^2}-\frac{1}{m^2k}+\frac{1}{m^2(k+m)}\\ &=2\zeta(3)+\frac{\pi^2}{6}H_{j+1}-\sum^{j+1}_{m=1}\frac{1}{m^2}\sum^\infty_{k=1}\left(\frac{1}{k}-\frac{1}{k+m}\right)\\ &=2\zeta(3)+\frac{\pi^2}{6}H_{j+1}-\sum^{j+1}_{k=1}\frac{H_k}{k^2}\\ \end{align} The second sum is \begin{align} \sum^\infty_{k=1}\frac{H_{k+j+1}}{(k+j+1)^2} &=2\zeta(3)-\sum^{j+1}_{k=1}\frac{H_k}{k^2} \end{align} Adding them up gives $$4\zeta(3)+\frac{\pi^2}{6}H_{j+1}-2\sum^{j+1}_{k=1}\frac{H_k}{k^2}$$ The third sum is \begin{align} \sum^\infty_{k=1}\frac{H_{k+j+1}}{k} &=\sum^\infty_{k=1}\frac{H_k}{k}+\sum^\infty_{k=1}\sum^{j+1}_{m=1}\frac{1}{k(k+m)}\\ &=\sum^\infty_{k=1}\frac{H_{k}}{k}+\sum^\infty_{k=1}\sum^{j+1}_{m=1}\frac{1}{mk}-\frac{1}{m(k+m)}\\ &=\sum^\infty_{k=1}\frac{H_{k}}{k}+\sum^{j+1}_{k=1}\frac{H_k}{k}\\ \end{align} The fourth sum is \begin{align} \sum^\infty_{k=1}\frac{H_{k+j+1}}{k+j+1} &=\sum^\infty_{k=1}\frac{H_{k}}{k}-\sum^{j+1}_{k=1}\frac{H_k}{k} \end{align} Subtracting them gives $$H_{j+1}^2+H_{j+1}^{(2)}$$ So $$\mathcal{I}_3=\frac{\pi^2}{6}\mathcal{I}_2-\sum^\infty_{j=1}\frac{1}{j^2(j+1)^2}\left(4\zeta(3)+\frac{\pi^2}{6}H_{j+1}-2\sum^{j+1}_{k=1}\frac{H_k}{k^2}\right)+2\sum^\infty_{k=1}\frac{H_{j+1}^2+H_{j+1}^{(2)}}{j^2(j+1)^3}$$ I will try to complete the evaluation of $\mathcal{I}_3$ over the next few days, and add in the computation of $\mathcal{I}_4$ using this method if I have the time. It is clear that $\mathcal{I}_3$ can be reduced to a bunch of Euler sums, which are all rather manageable to compute.

Integrate by parts twice. \begin{align} \int^1_0{\rm Li}_2^3(x){\rm d}x &=\left[x{\rm Li}^3_2(x)\right]^1_0+\int^1_03{\rm Li}^2_2(x)\ln(1-x){\rm d}x\\ &=\frac{\pi^6}{216}-\frac{\pi^4}{12}+6\color\red{\int^1_0\frac{\left[(x-1)\ln(1-x)-x\right]{\rm Li}_2(x)\ln(1-x)}{x}{\rm d}x}\\ &=\frac{\pi^6}{216}-\frac{\pi^4}{12}+6\zeta(5)+\frac{2\pi^4}{5}+36\zeta(3)+3\pi^2-2\pi^2\zeta(3)-90\\ &=\frac{\pi^6}{216}+\frac{19\pi^4}{60}+6\zeta(5)+36\zeta(3)+3\pi^2-2\pi^2\zeta(3)-90\\ \end{align} The red integral was evaluated in your previous answer.

\begin{eqnarray} {\mathcal I}_1 = \sum\limits_{n=1}^\infty \left(\frac{\zeta(2)}{n^2(n+1)} - \frac{H_{n+1}}{n^2(n+1)^2}\right) = \zeta^2(2) - 2 \zeta(2) - 4 \zeta(3) + 6 \end{eqnarray} In here we expanded both poly-logs in series, integrated term by term and then reduced the expression into partial fractions in one summation index and then summed over that index. In order to sum over the second index we again used partial fraction decomposition . The only non-trivial sum we come up with is now $\sum_{n=1}^\infty H_n/n^2$ which is known. \begin{eqnarray} {\mathcal I}_2 &=& \sum\limits_{n=1}^\infty\sum\limits_{k=1}^\infty \left(\frac{\zeta(2)}{k^2 n^2(n+k+1)} - \frac{H_{n+k+1}}{k^2 n^2(n+k+1)^2}\right)\\ & =& \zeta(2) \cdot {\mathcal I}_1- \sum\limits_{n=2}^\infty \frac{H_{n+1}}{(n+1)^2} \left(4 \frac{H_{n-1}}{n^3} + 2 \frac{H_{n-1}^{(2)}}{n^2}\right) \end{eqnarray} Now, it is clear that the sums on the right hand side are also doable. It is just a tedious task to do partial fraction decomposition and look up the respective harmonic sums that come up. We will finish this later since now we lack time for that. \begin{eqnarray} {\mathcal I}_3 &=& \zeta(2) \cdot {\mathcal I}_2- \sum\limits_{n=2}^\infty \frac{H_{n+1}}{(n+1)^2} \left(4 \sum\limits_{k=0}^{n-2}\frac{H_{k}}{(n-1-k)^2 (k+1)^3} + 2 \sum\limits_{k=0}^{n-2}\frac{H_{k}^{(2)}}{(n-1-k)^2 (k+1)^2}\right) \end{eqnarray} Now it gets interesting. It is not quite clear if sum over $k$ can be expressed by generalized harmonic numbers and rational functions only .Clearly in this case we do not end up with some sort of sums that we can look up and write down the answer straight away. We will analyze this thread deeper in due course.