Solution 1:

Lemma 1: For $\lambda \in \mathbb{R}^{+}$ $$\sum\limits_{n=1}^{\infty} \dfrac{\lambda^n}{n!}\mathcal{H}_{n} = e^{\lambda}\left(\ln \lambda + \gamma - \operatorname{Ei}(-\lambda)\right) \tag{1} \label{lemma1}$$ where $\mathcal{H}_{n}$$n$-th harmonic number and $\operatorname{Ei}(\cdot)$ — exponential integral.

Start from two series representations of lower incomplete gamma function $\gamma(\beta, \lambda)$: \begin{align} \gamma(\beta, \lambda) &= e^{-\lambda}\sum\limits_{n=0}^{\infty} \dfrac{\lambda^{n+\beta}}{\beta\left(\beta+1\right)\ldots\left(\beta+n\right)} \\ \gamma(\beta, \lambda) &= \sum\limits_{n=0}^{\infty} (-1)^{n}\dfrac{\lambda^{n+\beta}}{n!\left(\beta+n\right)} \end{align} Now take the derivative with respect to $\beta$ at $1$. Since \begin{align*} \dfrac{\mathrm{d}}{\mathrm{d}\beta}\left(\dfrac{1}{\beta\left(\beta+1\right)\ldots\left(\beta+n\right)}\right)_{\beta=1} &= -\left.\dfrac{1}{\beta\left(\beta+1\right)\ldots\left(\beta+n\right)}\left(\dfrac{1}{\beta}+\dfrac{1}{\beta+1}+\ldots+\dfrac{1}{\beta+n}\right)\right\vert_{\ \beta=1} \\ &= -\dfrac{\mathcal{H}_{n+1}}{(n+1)!} \end{align*} we have that \begin{align*} -e^{-\lambda}\sum\limits_{n=1}^{\infty} \dfrac{\lambda^{n}}{n!}\mathcal{H}_{n}+\ln \lambda\ e^{-\lambda}\sum\limits_{n=1}^{\infty} \dfrac{\lambda^{n}}{n!} &= -\ln \lambda\sum\limits_{n=1}^{\infty} (-1)^{n}\dfrac{\lambda^{n}}{n!} + \sum\limits_{n=1}^{\infty} (-1)^{n}\dfrac{\lambda^{n}}{n!\ n} \\ -e^{-\lambda}\sum\limits_{n=1}^{\infty} \dfrac{\lambda^{n}}{n!}\mathcal{H}_{n}+\ln \lambda \left(1-e^{-\lambda}\right) &= \ln \lambda \left(1-e^{-\lambda}\right) + \sum\limits_{n=1}^{\infty} (-1)^{n}\dfrac{\lambda^{n}}{n!\ n} \\ \sum\limits_{n=1}^{\infty} \dfrac{\lambda^{n}}{n!}\mathcal{H}_{n} &= -e^{\lambda}\sum\limits_{n=1}^{\infty} (-1)^{n}\dfrac{\lambda^{n}}{n!\ n} \\ &= e^{\lambda}\left(\ln \lambda + \gamma - \operatorname{Ei}(-\lambda)\right) \end{align*}


Let $$ \mathfrak{I}(\lambda, \alpha, \beta) = \int\limits_{0}^{\infty} \dfrac{e^{-\lambda x}\,x^{\alpha - 1}}{\left(1+x\right)^{\alpha+\beta}}\,\mathrm{d}x $$

Lemma 2: Let $\alpha, \lambda \in \mathbb{R}^{+}$ and $\alpha + \beta > 0$. Then $$ \mathfrak{I}(\lambda, \alpha, \beta) = \dfrac{\Gamma(\alpha)}{\Gamma(\alpha+\beta)}\lambda^{\beta}\mathfrak{I}(\lambda, \alpha + \beta, -\beta) \tag{2} \label{lemma2} $$

Using Laplace transform properties we have \begin{align*} \mathfrak{I}(\lambda, \alpha, \beta) &= \int\limits_{0}^{\infty} \dfrac{e^{-\lambda x}x^{\alpha - 1}}{\left(1+x\right)^{\alpha+\beta}}\,\mathrm{d}x \\ &= \int\limits_{0}^{\infty} e^{-\lambda t}t^{\alpha-1}\,\mathcal{L}\left\{\dfrac{e^{-x}x^{\alpha+\beta-1}}{\Gamma(\alpha+\beta)}\right\}(t)\,\mathrm{d}t \\ &= \int\limits_{0}^{\infty} \dfrac{\Gamma(\alpha)}{\left(\lambda+x\right)^{\alpha}}\dfrac{e^{-x}x^{\alpha+\beta-1}}{\Gamma(\alpha+\beta)} \,\mathrm{d}x \\ &= \dfrac{\Gamma(\alpha)}{\Gamma(\alpha+\beta)}\lambda^{\beta}\int\limits_{0}^{\infty} \dfrac{e^{-\lambda x}x^{\alpha+\beta-1}}{\left(1+x\right)^{\alpha}}\,\mathrm{d}x \\ &= \dfrac{\Gamma(\alpha)}{\Gamma(\alpha+\beta)}\lambda^{\beta}\mathfrak{I}(\lambda, \alpha + \beta, -\beta) \end{align*}


Claim 1: For $\lambda > 0$ and $\beta \in \mathbb{R}$ $$ \mathfrak{I}(\lambda, 1, \beta) = e^{\lambda}\lambda^{\beta}\Gamma(-\beta, \lambda) \tag{3} \label{claim1} $$ where $\Gamma(\cdot, \cdot)$ — upper incomplete gamma function.

\begin{align*} \mathfrak{I}(\lambda, 1, \beta) &= \int\limits_{0}^{\infty} \dfrac{e^{-\lambda x}}{\left(1+x\right)^{1+\beta}}\,\mathrm{d}x \\ &= e^{\lambda}\int\limits_{1}^{\infty} \dfrac{e^{-\lambda x}}{x^{1+\beta}}\,\mathrm{d}x \\ &= e^{\lambda}\lambda^{\beta}\int\limits_{\lambda}^{\infty} \dfrac{e^{-x}}{x^{1+\beta}}\,\mathrm{d}x \\ &= e^{\lambda}\lambda^{\beta}\Gamma(-\beta, \lambda) \end{align*}


After applying \eqref{lemma2} to \eqref{claim1} we have that $$ \int\limits_{0}^{\infty} \dfrac{e^{-\lambda x}x^{\beta}}{1+x}\,\mathrm{d}x = e^{\lambda}\Gamma(1+\beta)\Gamma(-\beta, \lambda) \tag{4} \label{consec1} $$ Taking the derivative with respect to the variable $\beta$ at point $(\lambda, \beta) = (1, 0)$ leads us to \begin{align*} \int\limits_{0}^{\infty} \dfrac{e^{-x}\ln x}{1+x}\,\mathrm{d}x &= e\dfrac{\partial}{\partial \beta}\Big(\Gamma(1+\beta)\Gamma(-\beta, 1)\Big)\Bigg\vert_{\beta=0} \\ &= -e\gamma \int\limits_{1}^{\infty} \dfrac{e^{-x}}{x}\,\mathrm{d}x -e \int\limits_{1}^{\infty} \dfrac{e^{-x}\ln x}{x}\,\mathrm{d}x \\ &= e\gamma\operatorname{Ei}(-1)-\dfrac{1}{2}e\int\limits_{1}^{\infty} e^{-x}\ln^2 x\,\mathrm{d}x \\ &= e\gamma\operatorname{Ei}(-1)-\dfrac{1}{2}e\dfrac{\mathrm{d}^2}{\mathrm{d}x^2}\Gamma(1)+\dfrac{1}{2}e\int\limits_{0}^{1} e^{-x}\ln^2 x\,\mathrm{d}x \\ \end{align*}

$$ \int\limits_{0}^{\infty} \dfrac{e^{-x}\ln x}{1+x}\,\mathrm{d}x = e\gamma\operatorname{Ei}(-1)-\dfrac{1}{2}e\left(\gamma^2+\dfrac{1}{6}\pi^2\right)-e\sum\limits_{n=1}^{\infty} \dfrac{(-1)^n}{n!\,n^2} \tag{5} \label{first} $$


Let $$ f(\lambda) = \int\limits_{0}^{\infty} e^{-x}\ln^2 \left(\lambda x + x^2\right)\,\mathrm{d}x $$

Then \begin{align*} f'(\lambda) &= 2\int\limits_{0}^{\infty} e^{-x}\dfrac{\ln \left(\lambda x + x^2\right)}{\lambda + x}\,\mathrm{d}x \\ &= 2\int\limits_{0}^{\infty} \dfrac{e^{-\lambda x}}{1+x}\Big(\ln \left(1+x\right)+\ln x + 2 \ln \lambda\Big)\,\mathrm{d}x \\ &= 4\ln \lambda \int\limits_{0}^{\infty} \dfrac{e^{-\lambda x}}{1+x}\,\mathrm{d}x + 2\int\limits_{0}^{\infty} \dfrac{e^{-\lambda x}}{1+x}\ln \left(1+x\right)\,\mathrm{d}x + 2\int\limits_{0}^{\infty} \dfrac{e^{-\lambda x}}{1+x}\ln x\,\mathrm{d}x \\ &= -4\ln \lambda e^{\lambda}\operatorname{Ei}(-\lambda) - 2\dfrac{\partial}{\partial \beta}\left.\left(\int\limits_{0}^{\infty} \dfrac{e^{-\lambda x}}{\left(1+x\right)^{1+\beta}}\,\mathrm{d}x\right)\right\vert_{\beta=0} + 2\dfrac{\partial}{\partial \beta}\left.\left(\int\limits_{0}^{\infty} \dfrac{e^{-\lambda x}x^{\beta}}{1+x}\,\mathrm{d}x\right)\right\vert_{\beta=0} \\ &= -4\ln \lambda e^{\lambda}\operatorname{Ei}(-\lambda) -2e^{\lambda}\dfrac{\partial}{\partial \beta}\Big(\lambda^{\beta}\Gamma(-\beta,\lambda)\Big)\Bigg\vert_{\beta=0} + 2e^{\lambda}\dfrac{\partial}{\partial \beta}\Big(\Gamma(1+\beta)\Gamma(-\beta,\lambda)\Big)\Bigg\vert_{\beta=0} \\ &= -2\ln \lambda e^{\lambda}\operatorname{Ei}(-\lambda)+2\gamma e^{\lambda}\operatorname{Ei}(-\lambda) \end{align*} Integrating over $[0,1]$ leads us to \begin{align*} f(1) &= f(0) -2\int\limits_{0}^{1} \ln \lambda e^{\lambda}\operatorname{Ei}(-\lambda)\,\mathrm{d}\lambda + 2\gamma\Big(e^{\lambda}\operatorname{Ei}(-\lambda)-\ln \lambda\Big)\Bigg\vert_{0}^{1} \\ &= 2\gamma e\operatorname{Ei}(-1)+2\gamma^2+\dfrac{2}{3}\pi^2-2\int\limits_{0}^{1} \ln \lambda e^{\lambda}\operatorname{Ei}(-\lambda)\,\mathrm{d}\lambda \end{align*} For integral in last formula we use \eqref{lemma1}: \begin{align*} \int\limits_{0}^{1} \ln \lambda e^{\lambda}\operatorname{Ei}(-\lambda)\,\mathrm{d}\lambda &= \int\limits_{0}^{1} \ln^2 \lambda e^{\lambda}\,\mathrm{d}\lambda +\gamma \int\limits_{0}^{1} \ln \lambda e^{\lambda}\,\mathrm{d}\lambda + \sum\limits_{n=2}^{\infty} \dfrac{\mathcal{H}_{n-1}}{n!\ n} \\ &= \sum\limits_{n=1}^{\infty} \dfrac{1}{n!\ n^2} + \gamma^2-\gamma\operatorname{Ei}(1) + \sum\limits_{n=1}^{\infty} \dfrac{\mathcal{H}_{n}}{n!\ n} \end{align*} Combining all together we have

$$ \int\limits_{0}^{\infty} e^{-x}\ln^2 \left(x+x^2\right)\,\mathrm{d}x = 2\gamma\Big(e\operatorname{Ei}(-1)+\operatorname{Ei}(1)\Big)+\dfrac{2}{3}\pi^2-2\sum\limits_{n=1}^{\infty} \dfrac{1}{n!\ n^2}-2\sum\limits_{n=1}^{\infty} \dfrac{\mathcal{H}_{n}}{n!\ n} \tag{6} \label{second} $$

From another side \begin{align*} \int\limits_{0}^{\infty} e^{-x}\ln^2 \left(x+x^2\right)\,\mathrm{d}x &= \int\limits_{0}^{\infty} e^{-x}\Big(\ln^2 \left(1+x\right) + \ln^2 x + 2\ln \left(1+x\right)\ln x\Big)\,\mathrm{d}x \\ &= \gamma^2+\dfrac{1}{6}\pi^2+e\int\limits_{1}^{\infty}e^{-x}\ln^2 x\,\mathrm{d}x + 2\int\limits_{0}^{\infty} e^{-x}\ln \left(1+x\right)\ln x\,\mathrm{d}x \\ &= \left(1+e\right)\left(\gamma^2+\dfrac{1}{6}\pi^2\right)+2e\sum\limits_{n=1}^{\infty} \dfrac{(-1)^n}{n!\ n^2} + 2\int\limits_{0}^{\infty} e^{-x}\ln \left(1+x\right)\ln x\,\mathrm{d}x \end{align*} Express last integral

\begin{align*} \int\limits_{0}^{\infty} e^{-x}\ln \left(1+x\right)\ln x\,\mathrm{d}x &= \gamma\Big(e\operatorname{Ei}(-1)+\operatorname{Ei}(1)\Big) -\dfrac{1}{2}\gamma^2+\dfrac{1}{4}\pi^2-\dfrac{1}{2}e\left(\gamma^2+\dfrac{1}{6}\pi^2\right) \\ &-e\sum\limits_{n=1}^{\infty} \dfrac{(-1)^n}{n!\ n^2}-\sum\limits_{n=1}^{\infty} \dfrac{1}{n!\ n^2}-\sum\limits_{n=1}^{\infty} \dfrac{\mathcal{H}_{n}}{n!\ n} \end{align*}

Now use integration by parts and already obtained result \eqref{first}: \begin{align*} \int\limits_{0}^{\infty} e^{-x}\ln \left(1+x\right)\ln x\,\mathrm{d}x &= \int\limits_{0}^{\infty} e^{-x}\dfrac{\ln x}{1+x}\,\mathrm{d}x + \int\limits_{0}^{\infty} e^{-x}\dfrac{\ln \left(1+x\right)}{x}\,\mathrm{d}x \\ &= e\gamma\operatorname{Ei}(-1)-\dfrac{1}{2}e\left(\gamma^2+\dfrac{1}{12}\pi^2\right)-e\sum\limits_{n=1}^{\infty} \dfrac{(-1)^n}{n!\ n^2}+ \int\limits_{0}^{\infty} e^{-x}\dfrac{\ln \left(1+x\right)}{x}\,\mathrm{d}x \end{align*} So

$$ \int\limits_{0}^{\infty} e^{-x}\dfrac{\ln \left(1+x\right)}{x}\,\mathrm{d}x = -\dfrac{1}{2}\gamma^2+\dfrac{1}{4}\pi^2+\gamma\operatorname{Ei}(1)-\sum\limits_{n=1}^{\infty} \dfrac{1}{n!\ n^2}-\sum\limits_{n=1}^{\infty} \dfrac{\mathcal{H}_{n}}{n!\ n} \tag{7} \label{final} $$


Back to original problem

With integration by parts and substitution original integral can be converted to \begin{align*} \int\limits_{0}^{e} \operatorname{Li}_2\left(\ln x\right)\,\mathrm{d}x &= \dfrac{1}{6}e\pi^2-\int\limits_{0}^{\infty} e^{-x}\dfrac{\ln \left(1+x\right)}{x}\,\mathrm{d}x+\int\limits_{0}^{1} e^{x}\dfrac{\ln \left(1-x\right)}{x}\,\mathrm{d}x \\ &= \dfrac{1}{6}e\pi^2 + \dfrac{1}{2}\gamma^2-\dfrac{1}{4}\pi^2-\gamma\operatorname{Ei}(1)+\sum\limits_{n=1}^{\infty} \dfrac{1}{n!\ n^2}+\sum\limits_{n=1}^{\infty} \dfrac{\mathcal{H}_{n}}{n!\ n} \\ &+\int\limits_{0}^{1} \dfrac{\ln\left(1-x\right)}{x}\,\mathrm{d}x+\sum\limits_{n=1}^{\infty}\dfrac{1}{n!}\int\limits_{0}^{1} x^{n-1}\ln\left(1-x\right)\,\mathrm{d}x \\ &= \dfrac{1}{6}e\pi^2 + \dfrac{1}{2}\gamma^2-\dfrac{5}{12}\pi^2-\gamma\operatorname{Ei}(1)+\sum\limits_{n=1}^{\infty} \dfrac{1}{n!\ n^2} \end{align*}

Solution 2:

Here I would like to sketch an alternative integral representation for $$ I(2)=\int_0^e\operatorname{Li}_2\left(\ln{x}\right)\,dx. $$ I hope that this argument shed some light on, why we cannot expect a much better representation as the hypergeometric series representation, which was given in the question.

After a change of variable $x \mapsto \ln(x)$, using an integral representation for $\operatorname{Li}_2$, if we interchange of the order of integration, we arrive at \begin{align} I(2) = \int_0^e\operatorname{Li}_2\left(\ln{x}\right)\,dx &= \int_{-\infty}^1e^x\operatorname{Li}_2\left(x\right)\,dx \\ &= - \int_{-\infty}^1 e^x \int_0^1\frac{\ln\left(1-xt\right)}{t}\,dt\,dx \\ &= - \int_0^1 \frac{1}{t} \int_{-\infty}^1 e^x\ln\left(1-xt\right)\,dx\,dt. \end{align}

By using the antiderivative of the function $e^x\ln\left(1-xt\right)$ with respect to $x$, we have \begin{align} I(2) &= - \int_0^1 \left(e\ln(1-t) + e^{1/t}E_1\left(\frac{1}{t}-1\right)\right)\,\frac{dt}{t} \\ &= - e\int_0^1 \frac{\ln(1-t)}{t} \,dt - \int_0^1 e^{1/t}E_1\left(\frac{1}{t}-1\right)\,\frac{dt}{t} \\ &= \frac{e\pi^2}{6} - \int_0^1 e^{1/t}E_1\left(\frac{1}{t}-1\right)\,\frac{dt}{t}, \end{align} where $E_1(x)$ is the $E_n$ function of order $n=1$. Using a change of variable $t\mapsto1/x$, we have $$ I(2) = \frac{e\pi^2}{6} - \int_1^\infty \frac{e^x E_1(x-1)}{x}\,dx. $$ We stop here. Although it would be nice to evaluate this integral, but as I know, this kind of integrals are not in the literature. Here I want to mention an earlier identity in a question of mine, where a solution of a related integral problem was also given in terms of ${_k}F_k\left(1,\dots,1;2,\dots,2;1\right)$. Of course any ideas for further evaluation are welcome. By using an integral representation for the $E_1$ function, we arrive at \begin{align} I(2) &= \frac{e\pi^2}{6} - \int_1^\infty \frac{e^x}{x} \int_1^\infty \frac{e^{-(x-1)t}}{t} \,dt \,dx \\ &= \frac{e\pi^2}{6} - \int_1^\infty \int_1^\infty \frac{e^{x+(1-x)t}}{xt} \,dt \,dx. \end{align}

Solution 3:

Too long for a comment.

If we evaluate the decay rate of your sum we find that $$\sum_{k=1}^\infty\frac{1}{k^2 k!}<\frac{1}{\sqrt{2\pi}}\sum_{k=1}^\infty\frac{e^k}{k^{k+5/2}}<\frac{1}{\sqrt{2\pi}}\sum_{k=1}^\infty\left(\frac{e}{k}\right)^k$$ This final series goes to zero smaller than the function $x^{-n}$ for any positive $n$, which helps us estimate the decay rate a little. We can also note that a number is likely irrational if $a_{k+1} \ll a_k^2$. Here we have that $$a_k^2\sim\left(\frac{1}{\sqrt{2\pi}\;k^{k+ 5/2} e^{-k}}\right)^2=\frac{e^{2 k} k^{-2 k - 5}}{2 \pi}$$ $$a_{k+1} \sim \frac{1}{\sqrt{2\pi}\;(k+1)^{n+7/2} e^{-k-1}} = \frac{e^{k + 1} (k + 1)^{-k - 7/2}}{\sqrt{2 \pi}}$$ This tells us that $$\frac{a_k^2}{a_{k+1}} \sim \frac{e^{k - 1} k^{-2 k - 5} (k + 1)^{k + 7/2}}{\sqrt{2 \pi}}$$ It's interesting to note how fast this goes to zero; in fact, if we denote the above by $f(k)$ and let $n$ be any real number we find that $$\lim_{k\to\infty}\frac{f(k)}{e^{nk}}=0$$
Thus, we see that this decay is very rapid, and conclude that $a_{k+1} \ll a_k^2$. I'm pretty sure this implies a proof of irrationality could be made (the fact that it is a sum of rational, decreasing terms alone encourages this), and possibly even a proof that the series is transcendental. Not too keen with constructing sequences to prove even irrationality though. Just let these serve as heuristics for now!