Closed form of $\int_{0}^1 \operatorname{Li}^2_n(x) dx$ for positive integer $n$

Question

As stated in the title, I want to find $$I(n)=\int_{0}^1 \operatorname{Li}_n^2(x) dx$$ where $n$ is a positive integer.

Using integration by parts, we have that $$I(n) = \operatorname{Li}_n^2(1)-2\int_0^1 \operatorname{Li}_{n-1}(x) \operatorname{Li}_n(x) dx$$

Let $$I_{n,k} = \int_0^1 \operatorname{Li}_{n}(x) \operatorname{Li}_k(x) dx$$ Integrating by parts, we have that $$I_{n,k} = \operatorname{Li}_{n}(1) \operatorname{Li}_k(1) - \int_{0}^1 \operatorname{Li}_{n}(x) \operatorname{Li}_{k-1}(x)+\operatorname{Li}_{n-1}(x) \operatorname{Li}_k(x) dx$$ We can get a recursive relationship from this as $$I_{n, k} = \operatorname{Li}_{n}(1) Li_k(1) - I_{n, k-1}-I_{n-1, k} = \zeta(n)\zeta(k) - I_{n, k-1}-I_{n-1, k}$$

We want to find $I_{n, n}$. However, the problem is that there is no base case yet. To me, the most sensible approach is to call $I_{m, 1} = I_{1, m}$ the base case. Unfortunately, I do not know how to solve $$\int_{0}^1 \operatorname{Li}_1(x)\operatorname{Li}_n(x)dx = -\int_{0}^1 \ln(1-x) \operatorname{Li}_n(x)dx$$ which is one concern. Even if I did have a closed form for $I_{n,1}$, how can I use that to get a closed form for $I_{n,k}$?

My main question is on how to find $I(n)$, either through what I have done so far or another method.

Edit: Letting $u = \operatorname{Li}_n(x)$ and $dv = \operatorname{Li}_1(x)$, we have that $$I_{1, n} = \operatorname{Li}_n(1)-\int_0^1 \operatorname{Li}_{n-1}(x) + \operatorname{Li}_{1}(x)\operatorname{Li}_{n-1}(x)-\frac{\operatorname{Li}_1(x)\operatorname{Li}_{n-1}(x)}{x}dx$$ The antiderivative of $\operatorname{Li}_n(x)$ is a known function, so we have that $$I_{1, n} = \zeta(n)-I_{1, n-1} +(-1)^n+ \sum_{m=2}^{n-1} (-1)^{n-m} \zeta(m) + \int_0^1 \frac{\operatorname{Li}_1(x)\operatorname{Li}_{n-1}(x)}{x}dx$$

Edit 2: I have an answer that I posted.

Answer 1

I will add onto AliShather's answer in order to get a complete answer. By using partial fraction decomposition with $\frac{1}{k^n(k+1)^i}$ we get that $$\frac{1}{k^n(k+1)^i} = (-1)^n \sum_{m=0}^{i-2} \frac{\binom{n+m-1}{n-1}}{(k+1)^{n-m}}+\sum_{m=0}^{n-2} \frac{(-1)^m \binom{i+m-1}{i-1}}{k^{n-m}}+\frac{(-1)^{n-1}\binom{i+n-2}{n-1}}{k}+\frac{(-1)^{n}\binom{i+n-2}{n-1}}{k+1}$$

This therefore means that $$\sum_{k=1}^{\infty} \frac{1}{k^n(k+1)^i} = (-1)^{n-1}\binom{n+i-1}{n}+\sum_{m=0}^{n-2}(-1)^m\binom{i+m-1}{i-1}\zeta(n-m) + (-1)^n\sum_{m=0}^{i-2}\binom{n+m-1}{n-1}\zeta(i-m)$$

Using $(1)$ again, we have that $$\frac{H_{k+1}}{k^n(k+1)^n} = \sum_{m=0}^{n-2} (-1)^n\binom{n+m-1}{n-1}\frac{H_{k+1}}{(k+1)^{n-m}}+\sum_{m=0}^{n-2} (-1)^m \binom{n+m-1}{n-1}\frac{H_{k+1}}{k^{n-m}}+(-1)^{n-1}\binom{2n-2}{n-1}\left(\frac{H_{k+1}}{k}-\frac{H_{k+1}}{k+1}\right)$$

Summing from $k = 1$ to $\infty$, we have $$\sum_{k=1}^\infty \frac{H_{k+1}}{k^n(k+1)^n} = \sum_{m=0}^{n-2} (-1)^n\binom{n+m-1}{n-1} \sum_{k=1}^\infty \frac{H_{k+1}}{(k+1)^{n-m}}+\sum_{m=0}^{n-2} (-1)^m \binom{n+m-1}{n-1} \sum_{k=1}^\infty \frac{H_{k+1}}{k^{n-m}}+(-1)^{n-1}\binom{2n-2}{n-1} \sum_{k=1}^\infty \left(\frac{H_{k+1}}{k}-\frac{H_{k+1}}{k+1}\right)$$

It is known that $$\sum_{k=1}^{\infty}\frac{H_k}{k^q} = \frac{(q+2)\zeta(q+1)}{2}- \sum_{k=1}^{q-2}\frac{\zeta(k+1)\zeta(q-k)}{2}$$ and since $$\frac{H_{k+1}}{k^{n-m}} = \frac{H_{k}}{k^{n-m}}+\frac{1}{k^{n-m}(k+1)}$$ we have that $$ \sum_{k=1}^\infty \frac{H_{k+1}}{k^{n-m}} = \frac{(n-m+2)\zeta(n-m+1)}{2}- \sum_{k=1}^{n-m-2}\frac{\zeta(k+1)\zeta(n-m-k)}{2}+(-1)^{n-m+1}+\sum_{k=2}^{n-m} (-1)^{n-m+k}\zeta(k)$$

Combining, we have that $$\sum_{k=1}^\infty \frac{H_{k+1}}{k^n(k+1)^n} = \sum_{m=0}^{n-2} (-1)^n\binom{n+m-1}{n-1} \left(-1+\frac{(n-m+2)\zeta(n-m+1)}{2}- \sum_{k=1}^{n-m-2}\frac{\zeta(k+1)\zeta(n-m-k)}{2} \right)+\sum_{m=0}^{n-2} (-1)^m \binom{n+m-1}{n-1} \left(\frac{(n-m+2)\zeta(n-m+1)}{2}- \sum_{k=1}^{n-m-2}\frac{\zeta(k+1)\zeta(n-m-k)}{2}+(-1)^{n-m+1}+\sum_{k=2}^{n-m} (-1)^{n-m+k}\zeta(k) \right)+(-1)^{n-1}\binom{2n-2}{n-1} \cdot 2$$

Altogether, this means that $$I_n = \sum_{m=0}^{n-2} \binom{n+m-1}{n-1} \left(1-\frac{(n-m+2)\zeta(n-m+1)}{2}+ \sum_{k=1}^{n-m-2}\frac{\zeta(k+1)\zeta(n-m-k)}{2} \right)+\sum_{m=0}^{n-2} \binom{n+m-1}{n-1} \left((-1)^{m+n+1}\frac{(n-m+2)\zeta(n-m+1)}{2}+ \sum_{k=1}^{n-m-2}(-1)^{m+n}\frac{\zeta(k+1)\zeta(n-m-k)}{2}+1+\sum_{k=2}^{n-m} (-1)^{k+1}\zeta(k) \right)+\binom{2n-2}{n-1} \cdot 2-\sum_{i=1}^{n-1} (-1)^i \zeta(n-i+1)\left((-1)^{n-1}\binom{n+i-1}{n}+\sum_{m=0}^{n-2}(-1)^m\binom{i+m-1}{i-1}\zeta(n-m) + (-1)^n\sum_{m=0}^{i-2}\binom{n+m-1}{n-1}\zeta(i-m) \right)$$

This is an incredibly long formula, but some parts should be able to be simplified.

Edit: $$I_n = -2\sum_{j=1}^{\lceil n/2 \rceil} \binom{2n-2j}{n} \zeta(2j)-2\sum_{j=1}^{\lfloor n/2 \rfloor} \left( \binom{2n-2j-1}{n-1}(j+1)-\binom{2n-2j-1}{n} \right) \zeta(2j+1) +2\sum_{m=2}^{\lfloor n/2 \rfloor} \binom{2n-2m-1}{n-1} \sum_{k=2}^{m}\zeta(k)\zeta(2m+1-k)+\binom{2n}{n}+ 2\sum_{i=3}^{n} \sum_{m=2}^{i-1}(-1)^{i+m}\binom{2n-m-i}{n-m}\zeta(m) \zeta(i) + \sum_{m=2}^n \binom{2n-2m}{n-m} (\zeta(m))^2 + 2\sum_{m=3}^{\lceil n/2 \rceil} \sum_{k=2}^{m -1} (-1)^k \binom{2n-2m}{n-1} \zeta(k) \zeta(2m-k) + \sum_{m=2}^{\lceil n/2 \rceil} (-1)^m \binom{2n-2m}{n-1} (\zeta(m))^2$$

Here is Mathematica code for a function that finds this

f[n_] := -2 Sum[
    Binomial[2 n - 2 j, n]*Zeta[2 j], {j, 1, Ceiling[n/2]}] - 
  2 Sum[(Binomial[2 n - 2 j - 1, n - 1]*(j + 1) - 
       Binomial[2 n - 2 j - 1, n])*Zeta[2 j + 1], {j, 1, 
     Floor[n/2]}] + 
  2*Sum[Binomial[2 n - 2 m - 1, n - 1]*
     Sum[Zeta[k]*Zeta[2 m + 1 - k], {k, 2, m}], {m, 2, Floor[n/2]}] + 
  Binomial[2 n, n] + 
  2*Sum[Sum[(-1)^(i + m)*Binomial[2 n - m - i, n - m]*Zeta[m]*
      Zeta[i], {m, 2, i - 1}], {i, 3, n}] + 
  Sum[Binomial[2 n - 2 m, n - m] (Zeta[m])^2, {m, 2, n}] + 
  2*Sum[Sum[(-1)^k*Binomial[2 n - 2 m, n - 1]*Zeta[k]*
      Zeta[2 m - k], {k, 2, m - 1}], {m, 3, Ceiling[n/2]}] + 
  Sum[(-1)^m*Binomial[2 n - 2 m, n - 1]*(Zeta[m])^2, {m, 2, 
    Ceiling[n/2]}]

Answer 2

\begin{align} I_n&=\int_0^1 \operatorname{Li}_n^2(x)\ dx=\sum_{k=1}^\infty\frac1{k^n}\int_0^1x^k\operatorname{Li}_n(x)\ dx\\ \end{align}

Since $$\int_0^1x^{k}\operatorname{Li}_n(x)\ dx\overset{IBP}{=}(-1)^{n-1}\frac{H_{k+1}}{(k+1)^n}-\sum_{i=1}^{n-1}(-1)^i\frac{\zeta(n-i+1)}{(k+1)^i}$$

Then

$$I_n=(-1)^{n-1}\sum_{k=1}^\infty\frac{H_{k+1}}{k^n(k+1)^n}-\sum_{k=1}^\infty\sum_{i=1}^{n-1}(-1)^i\frac{\zeta(n-i+1)}{k^n(k+1)^i}$$

Which can be more simplified I think.

You can see a related problem here.