$\require{begingroup} \begingroup$
$\def\e{\mathrm{e}}\def\W{\operatorname{W}}\def\Wp{\operatorname{W_0}}\def\Wm{\operatorname{W_{-1}}}\def\erf{\operatorname{erf}}\def\Ei{\operatorname{Ei}}$
For $x\in\mathbb R,\ x\ne-1$
\begin{align} I(x)&= \int_0^1 (1-\Wm(-\tfrac t\e))^x - (1-\Wp(-\tfrac t\e))^x \, dt \tag{1}\label{1} \\ &= \frac{\e\,(2+\e\,(x-1)\,\Gamma(x+2,1))}{x+1} =f(x) \tag{2}\label{2} , \end{align}
where $\Wp,\ \Wm$ are the real branches of the Lambert $\W$ function, and $\Gamma$ is the incomplete gamma function.
For integer values $x=n$, $I(n)$ follows the pattern of $\e\cdot a_n$ fromA093964.
\begin{align} I(1)&=f(1)=\e ,\\ I(2)&=f(2)=6\,\e ,\\ I(3)&=f(3)=33\,\e ,\\ &\dots \end{align}
Some other special cases of $x$:
\begin{align} I(\tfrac12)&= \int_0^1 \sqrt{1-\Wm(-\tfrac t\e)} \, dt - \int_0^1 \sqrt{1-\Wp(-\tfrac t\e)} \, dt \\ &= (\tfrac32\,\sqrt2+\tfrac14\,\e^2\,\sqrt\pi\,(\erf(\sqrt2)-1)) \\ &-(\tfrac32\,\sqrt2-\tfrac12\,\e +\tfrac14\,\e^2\,\sqrt\pi\,(\erf(\sqrt2)-\erf(1))) \\ &= \tfrac12\,\e+\tfrac14\,\e^2\,\sqrt\pi\,(\erf(1)-1) =f(\tfrac12) \tag{3}\label{3} \\ &\approx 0.844113386646 ,\\ I(-2)&= \int_0^1 \frac1{(1-\Wm(-\tfrac t\e))^2} - \frac1{(1-\Wp(-\tfrac t\e))^2} \, dt \\ &\approx -.57344306156 \\ &= \e\,(3\,\e\,\Ei(1,1)-2) =f(-2) \tag{4}\label{4} , \end{align} where $\Ei(a,x) = x^{a-1} \Gamma(1-a,x)$.
\begin{align} I(-\Omega) &\approx -0.4015641473638446 \approx f(-\Omega) ,\quad \Omega=\W(1)\approx .56714329 . \end{align}
Also,
\begin{align} I(-1)&\approx -0.523798568446 \\ &\approx \e\,(1-2\,\e\,\Ei(1,1))=\lim_{x\to -1}f(x) . \end{align}
Questions:
1) Is this correct/known? Any reference/confirmation?
2) Is it possible to transform \eqref{2} in order to cure the nasty case of $x=-1$?
$\endgroup$
By simple change of variables, we have $$ - \int_0^1 {\left( {1 - W_0 \left( { - \frac{t}{e}} \right)} \right)^x dt} = \int_0^1 {(1 + s)^x e^{1 - s} (s - 1)ds} $$ and $$ \int_0^1 {\left( {1 - W_{ - 1} \left( { - \frac{t}{e}} \right)} \right)^x dt} = \int_1^{ + \infty } {(1 + s)^x e^{1 - s} (s - 1)ds} . $$ Consequently, $$ I(x) = \int_0^{ + \infty } {(1 + s)^x e^{1 - s} (s - 1)ds} = e\int_0^{ + \infty } {\frac{{e^{ - s} s}}{{(1 + s)^{ - x} }}ds} - e\int_0^{ + \infty } {\frac{{e^{ - s} }}{{(1 + s)^{ - x} }}ds} \\ = e\int_0^{ + \infty } {\frac{{e^{ - s} }}{{(1 + s)^{ - x - 1} }}ds} - 2e\int_0^{ + \infty } {\frac{{e^{ - s} }}{{(1 + s)^{ - x} }}ds} = e^2 \Gamma (x + 2,1) - 2e^2 \Gamma (x + 1,1) \\ = e^2 \Gamma (x + 2,1) - 2e^2 \frac{{\Gamma (x + 2,1) - e^{ - 1} }}{{x + 1}} = e\frac{{2 + e(x - 1)\Gamma (x + 2,1)}}{{x + 1}}. $$ Regarding your second question, you may write $$ I(x) = e^2 \Gamma (x + 2,1) - 2e^2 \Gamma (x + 1,1) = e(e(x - 1)\Gamma (x + 1,1) + 1). $$