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

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

• It's interesting that you consider the Meijer G-function outside of your scope, and yet include the generalized hypergeometric. I'm curious what contexts you learned above the latter but not the former. I see them as steps on an infinite ladder of more general functions defined by series, and have always seen them covered together that way. Nov 2, 2016 at 19:40
• @BrevanEllefsen I'm asking for a better closed form, as the $_3F_3$ term was originally a series. Specifically: $$\sum_{k=1}^\infty\frac{1}{k^2 k!}.$$ By no means am I saying that $_3F_3$ is within my scope; it's far from it, in fact. I can't find anything online on either $_3F_3$ or Meijer G. I guess I just shared the best form I could get :P Nov 2, 2016 at 19:48
• Note that this can be more simply written as $$= \frac{e \pi^2}{6}+\int_{-\infty}^1\frac{e^x\log(1-x)}{x}dx$$ with a basic substitution and IBP. This integral is much more easy to analyze in my opinion. Nov 2, 2016 at 22:14
• Moreover, Mathematica instantly converts the integral to $$\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)$$ Nov 2, 2016 at 22:19
• @BrevanEllefsen I've edited the question to include your form and added some clarification. Thank you :) Nov 3, 2016 at 15:12

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

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}

• This form for $I(2)$ popped up during my calculations, but I chose to omit it since it appeared to lead nowhere. However, it's good to see that someone else arrived at the same form as well. Also, swapping the order of integration in the final double integral you've posted just leads to $$\int_1^\infty \frac{e^t E_1(t-1)}{t}dt.$$ Furthermore, $x+(1-x)t=t+(1-t)x$, so the order of integration is irrelevant--this transformation is equivalent to making the substitution $x\mapsto t$. Nov 23, 2016 at 14:36

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!