What is a closed form of $\int_0^1\ln(-\ln x)\ \text{li}\ x\ dx$
This integral can be expressed in a closed form: $$\int_0^1\ln(-\ln x)\operatorname{li}x\ dx=\frac{\pi^2}{12}+\frac{(\ln2)^2}2+\gamma\,\ln2,\tag1$$ where $\gamma$ is the Euler-Mascheroni constant.
Proof:
Let $$\mathcal{I}(a)=\int_0^1(-\ln x)^a\operatorname{li}x\ dx,\tag2$$ then, taking a derivative with respect to $a$, we get $$\mathcal{I}'(0)=\int_0^1\ln(-\ln x)\operatorname{li}x\ dx.\tag3$$ Mathematica can express the integral $(2)$ in a closed form for explicitly given non-negative integer values of $a$, but not in a general form with a parameter. We could use those values to get a conjectured general formula, and indeed one could get the result $(1)$ this way. Numerical integration shows that $$\Bigg|\frac{\pi^2}{12}+\frac{(\ln2)^2}2+\gamma\,\ln2-\int_0^1\ln(-\ln x)\operatorname{li}x\ dx\Bigg|<10^{-1000}\tag4,$$ so the result looks plausible.
But here we will follow a different approach to get a rigorous proof (modulo possible bugs in Mathematica tuned exactly to hide that tiny difference). Take the formula $(12)$ from the MathWorld article: $$\operatorname{li}x=\operatorname{Ei}(\ln x),\tag5$$ where $\operatorname{Ei}z$ is the exponential integral, and plug it into $(2)$. In this form the integral can be evaluated in Mathematica, that gives us $$\mathcal{I}(a)=-\frac{_2F_1(a+1,a+1;a+2;-1)\ \Gamma(a+1)}{a+1},\tag6$$ and, taking a derivative with respect to $a$, $$\mathcal{I}'(0)=(1+\gamma)\ln2-\left[\frac{d}{da}{_2F_1}(a+1,a+1;a+2;-1)\right]_{a=0}.\tag7$$ To find the derivative of the hypergeometric function we need to replace it with its definition (formula $(8)$ in this MathWorld article): $$_2F_1(a+1,a+1;a+2;-1)=\sum_{n=0}^\infty\frac{(a+1)\ \Gamma(a+n+1)}{(a+n+1)\ \Gamma(a+1)\ n!}(-1)^n.\tag8$$ So, $$\left[\frac{d}{da}{_2F_1}(a+1,a+1;a+2;-1)\right]_{a=0}=\sum_{n=0}^\infty\frac{(n+1)(\gamma+\psi(n+1))+n}{(n+1)^2}(-1)^n,\tag9$$ where $\psi(n+1)=\frac{\Gamma'(n+1)}{\Gamma(n+1)}=-\gamma+\sum_{m=1}^n\frac1m\,$ is the digamma function. Evaluating this sum with Mathematica we get $$\left[\frac{d}{da}{_2F_1}(a+1,a+1;a+2;-1)\right]_{a=0}=\ln2-\frac{(\ln2)^2}2-\frac{\pi^2}{12}\tag{10}.$$ Plugging this back in $(7)$ and $(3)$ we get the final result $(1)$.
Solution by Khalef Ruhemi
From the definition
$$\operatorname{li}(x)=\int_0^x\frac{\ dt}{\ln t}=-\int_0^x\int_0^\infty t^y\ dt\ dy=-\int_0^\infty\frac{x^{y+1}}{1+y}\ dy=-\int_0^\infty\frac{e^{-(1+y)\ln(1/x)}}{1+y}\ dy$$
where we used $$\frac{1}{\ln t}=-\int_0^\infty t^y\ dy, \quad t\in(0,1)$$
It follows that
\begin{align} I&=-\int_0^1\left(\ln\left(\ln(1/x)\right)\int_0^\infty\frac{e^{-(1+y)\ln(1/x)}}{1+y}\ dy\right)\ dx,\quad \ln(1/x)=t\\ &=-\int_0^\infty\left(\frac1{1+y}\int_0^\infty e^{-(2+y)t}\ln t \ dt\right)\ dy\\ &=\int_0^\infty\frac1{1+y}\left(\frac{\gamma}{2+y}+\frac{\ln(2+y)}{2+y}\right)\ dy,\quad \frac1{1+y}=x\tag{*}\\ &=\gamma\int_0^1\frac{\ dx}{1+x}+\int_0^1\frac{\ln(1+x)-\ln x}{1+x}\ dx\\ &=\gamma\ln2+\frac12\ln^22+\frac12\zeta(2) \end{align}
Regarding step $(*)$, we know that $$\int_0^\infty t^p e^{-(y+2)t}\ dt=\Gamma(p+1)(2+y)^{-p-1}$$
Differentiate both sides with respect to $p$ to get
$$\int_0^\infty t^p e^{-(y+2)t}\ln t\ dt=\Gamma'(p+1)(2+y)^{-p-1}-\Gamma(p+1)(2+y)^{-p-1}\ln(2+y)$$
Set $p=0$ we get
$$\int_0^\infty e^{-(2+y)t}\ln t\ dt=-\frac{\gamma}{2+y}-\frac{\ln(2+y)}{2+y}$$