Integral involving Clausen function ${\large\int}_0^{2\pi}\operatorname{Cl}_2(x)^2\,x^p\,dx$

Consider the Clausen function $\operatorname{Cl}_2(x)$ that can be defined for $0<x<2\pi$ in several equivalent ways:

$$\begin{align}\operatorname{Cl}_2(x)&=-\int_0^x\ln\left(2\sin\left(\tfrac t2\right)\right)dt\\&=\sum_{n=1}^\infty\frac{\sin\,(n x)}{n^2}\\&=\Im\operatorname{Li}_2\left(e^{i x}\right)\\&=i\left(\frac{\pi^2}6+\frac{x^2}4-\frac{\pi x}2-\operatorname{Li}_2\left(e^{i x}\right)\right).\end{align}\tag1$$ I'm interested in integrals of the form $$I(p)=\int_0^{2\pi}\operatorname{Cl}_2(x)^2\,x^p\,dx.\tag2$$ I found that $$I(0)=\frac{\pi^5}{90}\tag3$$ and conjectured next several values based on numeric evidence: $$I(1)\stackrel?=\frac{\pi^6}{90},\ \ I(2)\stackrel?=\frac{44\,\pi^7}{2835},\ \ I(3)\stackrel?=\frac{23\,\pi^8}{945}.\tag4$$ One might expect that $I(4)$ is a rational multiple of $\pi^9$, but apparently it is not (unless the denominator is huge).

I'm asking for your help in proving conjectured values $(4)$, finding a closed form of $I(4)$, and, if possible, a general formula for $I(p)$.

Update: Values $I(-1)$ and $I(-2)$ are also interesting.

The General Case:

Consider instead the integral \begin{align} J(\eta) &=\Re\int^{2\pi}_0x^\eta\ {\rm Li}_2(e^{ix})^2\ {\rm d}x\\ \end{align} for $\eta\in\mathbb{N}_0$. Deforming the contour around $z=0$, we get \begin{align} J(\eta) &=\Re\oint_{|z|=1}\frac{{\rm Li}_2(z)^2\ln^{\eta}(z)}{i^{\eta+1}z}{\rm d}z\\ &=\Re\int^1_0\frac{\left((\ln(x)+2\pi i)^{\eta}-\ln^{\eta}(x)\right){\rm Li}_2(x)^2}{i^{\eta+1}x}{\rm d}x\\ &=\Re\sum^{\eta-1}_{k=0}\binom{\eta}{k}\frac{(2\pi)^{\eta-k}}{i^{k+1}}\int^1_0\frac{{\rm Li}_2(x)^2\ln^{k}(x)}{x}{\rm d}x\\ &=\sum^{\lfloor{\eta/2\ -1}\rfloor}_{k=0}\binom{\eta}{2k+1}(-1)^{k+1}(2\pi)^{\eta-2k-1}\int^1_0\frac{{\rm Li}_2(x)^2\ln^{2k+1}(x)}{x}{\rm d}x \end{align} Relating this to $I(\eta)$, $$I(\eta)=\pi^{\eta+5}\left(\frac{2^{\eta+1}}{36(\eta+1)}-\frac{2^{\eta+2}}{6(\eta+2)}+\frac{2^{\eta+3}}{3(\eta+3)}-\frac{2^{\eta+4}}{4(\eta+4)}+\frac{2^{\eta+5}}{16(\eta+5)}\right)-J(\eta)$$

The Case $I(0)$:

It is evident that $J(0)=0$ since it is a vacuous summation. Thus $$I(0)=\left.\pi^{\eta+5}\left(\frac{2^{\eta+1}}{36(\eta+1)}-\frac{2^{\eta+2}}{6(\eta+2)}+\frac{2^{\eta+3}}{3(\eta+3)}-\frac{2^{\eta+4}}{4(\eta+4)}+\frac{2^{\eta+5}}{16(\eta+5)}\right)\right|_{\eta=0}=\frac{\pi^5}{90}$$

The Case $I(1)$:

In this case, $J(1)=0$ as well. Hence $$I(1)=\left.\pi^{\eta+5}\left(\frac{2^{\eta+1}}{36(\eta+1)}-\frac{2^{\eta+2}}{6(\eta+2)}+\frac{2^{\eta+3}}{3(\eta+3)}-\frac{2^{\eta+4}}{4(\eta+4)}+\frac{2^{\eta+5}}{16(\eta+5)}\right)\right|_{\eta=1}=\frac{\pi^6}{90}$$

The Case $I(2)$:

For $\eta=2$, $$\left.\pi^{\eta+5}\left(\frac{2^{\eta+1}}{36(\eta+1)}-\frac{2^{\eta+2}}{6(\eta+2)}+\frac{2^{\eta+3}}{3(\eta+3)}-\frac{2^{\eta+4}}{4(\eta+4)}+\frac{2^{\eta+5}}{16(\eta+5)}\right)\right|_{\eta=2}=\frac{16\pi^7}{945}$$ and $$ J(2)=\binom{2}{1}(-1)(2\pi)\int^1_0\frac{{\rm Li}_2(x)^2\ln{x}}{x}{\rm d}x=4\pi\sum^\infty_{m=1}\sum^\infty_{n=1}\frac{1}{m^2n^2(m+n)^2}=\frac{4\pi^7}{2835}$$ The double sum $\displaystyle\sum^\infty_{m=1}\sum^\infty_{n=1}\frac{1}{m^2n^2(m+n)^2}=\frac{\pi^6}{2835}$ has been evaluated here. Thus $$I(2)=\frac{16\pi^7}{945}-\frac{4\pi^7}{2835}=\frac{44\pi^7}{2835}$$

The Case $I(3)$:

When $\eta=3$, $$\left.\pi^{\eta+5}\left(\frac{2^{\eta+1}}{36(\eta+1)}-\frac{2^{\eta+2}}{6(\eta+2)}+\frac{2^{\eta+3}}{3(\eta+3)}-\frac{2^{\eta+4}}{4(\eta+4)}+\frac{2^{\eta+5}}{16(\eta+5)}\right)\right|_{\eta=3}=\frac{\pi^8}{35}$$ and $$J(3)=\binom{3}{1}(-1)(2\pi)^2\int^1_0\frac{{\rm Li}_2(x)^2\ln{x}}{x}{\rm d}x=12\pi^2\sum^\infty_{m=1}\sum^\infty_{n=1}\frac{1}{m^2n^2(m+n)^2}=\frac{4\pi^8}{945}$$ Therefore $$I(3)=\frac{\pi^8}{35}-\frac{4\pi^8}{945}=\frac{23\pi^8}{945}$$

The Case $I(4)$:

For $\eta=4$, $$\left.\pi^{\eta+5}\left(\frac{2^{\eta+1}}{36(\eta+1)}-\frac{2^{\eta+2}}{6(\eta+2)}+\frac{2^{\eta+3}}{3(\eta+3)}-\frac{2^{\eta+4}}{4(\eta+4)}+\frac{2^{\eta+5}}{16(\eta+5)}\right)\right|_{\eta=4}=\frac{16\pi^9}{315}$$ and \begin{align} J(4) &=\binom{4}{1}(-1)(2\pi)^3\left(-\frac{\pi^6}{2835}\right)+\binom{4}{3}(2\pi)\int^1_0\frac{{\rm Li}_2(x)^2\ln^3{x}}{x}{\rm d}x\\ &=\frac{32\pi^9}{2835}+4\pi\int^1_0\frac{{\rm Li}_2(x)\ln^4{x}\ln(1-x)}{x}{\rm d}x\\ &=\frac{32\pi^9}{2835}+\frac{4\pi}{5}\int^1_0\frac{{\rm Li}_2(x)\ln^5{x}}{1-x}{\rm d}x+\frac{4\pi}{15}\int^1_0\frac{\ln^6{x}\ln(1-x)}{1-x}{\rm d}x\\ &=\frac{32\pi^9}{2835}-96\pi\sum^\infty_{m=1}\frac{H_n^{(2)}}{(n+1)^6}-192\pi\sum^\infty_{n=1}\frac{H_n}{(n+1)^7}\\ &=-\frac{8\pi^9}{567}+192\pi\zeta(3)\zeta(5)-96\pi\zeta(6,2) \end{align} giving us $$I(4)=\frac{184\pi^9}{2835}-192\pi\zeta(3)\zeta(5)+96\pi\zeta(6,2)$$ I am not sure if $\zeta(6,2)$ has a simple closed form. At least Flajolet and Salvy's paper on Euler sums suggests that there isn't.

I can confirm $I(1)$:

Using the series expansion for Clausen's function we have

$$ I(p)=\sum_{n=1}^{\infty}\sum_{m=1}^{\infty}\int_0^{2 \pi}x^p\frac{\sin(n x)}{n^2}\frac{\sin(m x)}{m^2}\,\mathrm dx $$

now setting $p=1$. Furthermore we may use the orthogonality of sine $\int^{2 \pi}_{0}{\sin(n x)\sin(m x)}\,\mathrm dx=\pi\delta_{mn}$ and integrate by parts to get

$$ I(1)=\underbrace{\sum_{m,n=1,n\neq m}^{\infty}\frac{1}{2 (m-n)^2}-\frac{1}{2 (m+n)^2}-\frac{2 m n}{(m^2 - n^2)^2}}_{=0}+\pi^2\underbrace{\sum_{n=1}^{\infty}\frac{1}{n^4}}_{=\zeta(4)}=\frac{\pi^6}{90} $$

I'm checking if there is a pattern, especially if it holds that the doublesum always cancel. But at least for the first few $p$'st this seems likely, nevertheless the algebra gets a little bit messy but i will try to go further tomorrow.

Edit: Playing the same game for $p=0$ we end up with $$ I(0)=\pi\underbrace{\sum_{n=1}^{\infty}\frac{1}{n^4}}_{=\zeta(4)}=\frac{\pi^5}{90} $$

The double sum cancels immediatly in this case, so i confirmed also $I(0)$

Edit2:

Doing some careful algebra,following the same approach as above, we may find that

$$ I(2)=\underbrace{4 \pi \sum_{n=1}^{\infty}{\frac{8 \pi ^2 n^3-3 n }{24 n^7}}}_{=\frac{\pi^7}{70}}+8 \pi\underbrace{\sum_{m,n=1,n\neq m}^{\infty}\frac{ 1}{m n\left(m^2-n^2\right)^2}}_{J}=8 \pi\underbrace{\sum_{m}^{\infty}\sum_{n=0}^{m-1}\frac{ 1}{m n\left(m^2-n^2\right)^2}+\sum_{m}^{\infty}\sum_{m+1}^{\infty}\frac{ 1}{m n\left(m^2-n^2\right)^2}}_{J_1+J_2}+\frac{\pi^7}{70} $$

Maybe someone with more experience in calculating sums like this can take it home from here :)

Edit3:

It turns out that the same off-diagonal sum will give us the result for $p=3$ if we multiply by $24\pi^2$ instead of $8\pi$. So if we have solved $p=2$ we get $I(3)$ for free :)

Nevertheless, for $p>3$ the off-diagonal contribution starts to get messier and yields additionals contribution which may well be not proportional to $\pi^9$

$I(1)$ can be evaluated in the following way also. One has to note that

$$\int_0^{2\pi} x\sin(mx)\sin(nx)\,dx=\begin{cases} 0 & n \neq m \\ \pi^2 & n=m \end{cases}$$

To prove the above, write $$\int_0^{2\pi} x\sin(mx)\sin(nx)\,dx=\frac{1}{2}\int_0^{2\pi} x\cos((m-n)x) \,dx-\frac{1}{2}\int_0^{2\pi} x\cos((m+n)x)\,dx$$ and then use integration by parts.

Back to the original integral, $$\begin{aligned} \int_0^{2\pi} x\text{Cl}_2(x)^2\,dx &=\sum_{m=1}^{\infty} \sum_{n=1}^{\infty} \frac{1}{m^2n^2}\int_0^{2\pi} x\sin(mx)\sin(nx)\,dx \\ & =\sum_{n=1}^{\infty} \frac{\pi^2}{n^4} \\ & =\frac{\pi^6}{90} \end{aligned}$$