Closed form for $\int_0^e\mathrm{Li}_2(\ln{x})\,dx$?

Question

Inspired by this question and this answer, I decided to investigate the family of integrals $$I(k)=\int_0^e\mathrm{Li}_k(\ln{x})\,dx,\tag{1}$$ where $\mathrm{Li}_k(z)$ represents the polylogarithm of order $k$ and argument $z$. $I(1)$ evaluates to $e\gamma$, but $I(2)$ has resisted my efforts (which can be seen here).

Neither ISC nor WolframAlpha could provide a closed form for its numerical value--however, I've conjectured a possible analytic form. $$\eqalign{\int_0^e\mathrm{Li}_2(\ln{x})\,dx&\stackrel?=\,_3F_3(1,1,1;2,2,2;1)+\frac{\pi^2(2e-5)}{12}+\frac{\gamma^2}{2}-\gamma\,\mathrm{Ei}(1)\\&=0.578255559804073275225659054377625577...\tag{2}}$$

Brevan Ellefsen has computed that my conjecture is accurate to at least 150 digits. Brevan also gave the alternate form $$\frac{\pi^2e}{6}+\gamma G_{1,2}^{2,0}\left(-1\left|\begin{array}{c}1\\0,0\\\end{array}\right.\right)+G_{2,3}^{3,0}\left(-1\left|\begin{array}{c}1,1\\0,0,0\\\end{array}\right.\right).\tag{2.1}$$

Is there a closed form for $I(2)$ that doesn't involve Meijer G or hypergeometrics? The simplicity of the following two equations seems to suggest that there might be. $(3.1)$ follows directly from $(3)$, which I've proven here. $$\begin{align} \sum_{k=1}^\infty I(k)&=e\tag{3}\\\sum_{k=2}^\infty I(k)&=e(1-\gamma)\tag{3.1}\end{align}$$

---

PROGRESS UPDATE: Using this equation, I've turned $_3F_3(1,1,1;2,2,2;1)$ into $$\lim_{c\to 1}\left(\frac{\mathrm{Ei}(1)-\gamma}{c-1}+\frac{1-e}{(c-1)^2}+\frac{(-1)^{-c}\,\Gamma(c-1)}{c-1}+\frac{(-1)^{1-c}\,\Gamma(c,-1)}{(c-1)^2}\right),\tag{4}\label{4}$$ but I don't know how to proceed from there. EDIT: This limit leads nowhere. See below.

SECOND PROGRESS UPDATE: After some studying of the properties of the Meijer G function, I've finally cracked the limit; however, the result is an underwhelming $_3F_3(1,1,1;2,2,2;1)$. Before I evaluate the limit, I'd first like to state the following intermediate result:

Lemma $(4.1)$: For $z\in\mathbb{C}$, $$G_{2,3}^{3,0}\left(z\left|\begin{array}{c}1,1\\0,0,0\\\end{array}\right.\right)=\gamma\ln{z}+\frac12\ln^2(z)-z\,_3F_3(1,1,1;2,2,2;-z)+\frac{\gamma^2}2+\frac{\pi^2}{12}.\tag{4.1}\label{4.1}$$

My proof for this can be found here. Now I return to the limit $\eqref{4}$. Consider the following: $\frac{1}{c-1}=\frac{c-1}{(c-1)^2}$, $(c-1)\Gamma(c-1)=\Gamma(c)$, and $(-1)^{1-c}=-(-1)^{-c}$. Based on these algebraic identities, the limit can be written as $$\lim_{c\to 1}\frac{(c-1)(\mathrm{Ei}(1)-\gamma)+1-e+(-1)^{-c}(\Gamma(c)-\Gamma(c,-1))}{(c-1)^2}.\tag{4.2}$$ In this form, the limit is $\frac{0}0$. Using l'Hospital twice, we obtain $$\lim_{c\to 1}{(-1)^{-c}\left(-G_{3,4}^{4,0}\left(-1\left|\begin{array}{c}1,1,1\\0,0,0,c\\\end{array}\right.\right)+\Gamma{(c)}\left(\frac{\psi_0{(c)}^2}2-i\pi\psi_0(c)+\frac{\psi_1(c)}2-\frac{\pi^2}2\right)\right)}$$ $$\begin{align}&=G_{3,4}^{4,0}\left(-1\left|\begin{array}{c}1,1,1\\0,0,0,1\\\end{array}\right.\right)-\frac{\psi_0(1)^2}2+i\pi\psi_0(1)-\frac{\psi_1(1)}2+\frac{\pi^2}2\\&=G_{2,3}^{3,0}\left(-1\left|\begin{array}{c}1,1\\0,0,0\\\end{array} \right.\right)-\frac{\gamma^2}{2}-i\gamma\pi+\frac{5\pi^2}{12}.\tag{4.2a}\end{align}$$ Using Lemma $\eqref{4.1}$, we know that $$G_{2,3}^{3,0}\left(-1\left|\begin{array}{c}1,1\\0,0,0\\\end{array}\right.\right)={}_3F_3(1,1,1;2,2,2;1)+\frac{\gamma^2}2+i\gamma\pi-\frac{5\pi^2}{12},\tag{4.3}$$ which can be rewritten as $$G_{2,3}^{3,0}\left(-1\left|\begin{array}{c}1,1\\0,0,0\\\end{array}\right.\right)-\frac{\gamma^2}{2}-i\gamma\pi+\frac{5\pi^2}{12}={}_3F_3(1,1,1;2,2,2;1).\tag{4.3a}$$

Thus, $$\eqalign{&\lim_{c\to 1}\left(\frac{\mathrm{Ei}(1)-\gamma}{c-1}+\frac{1-e}{(c-1)^2}+\frac{(-1)^{-c}\,\Gamma(c-1)}{c-1}+\frac{(-1)^{1-c}\,\Gamma(c,-1)}{(c-1)^2}\right)\\=&{}_3F_3(1,1,1;2,2,2;1).\tag{4.4}}$$

Answer 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*}

Answer 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}

Answer 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!