A probabilistic average integral I wanted to see what are the chances of at least one success if I try $1/p$ times the same experiment with probability $p$. If I repeat $n$ times a $\mathscr{B}(p)$ law, then the probability is : $\mathbb{P}(X_n\geqslant1)=1-(1-p)^n$. Hence, the quantity I am looking for is $1-(1-p)^{1/p}$.
It can be interesting to see what happens when $p\to0$ : it converges to $\frac{e-1}e\approx63\%$. Hence, the function $p\in[0,1]\mapsto1-(1-p)^{1/p}$ is continuous on $[0,1]$, so we can calculate its average value.
This leads me to calculating this integral : $\mathcal{I}:=1-\displaystyle\int_0^1(1-p)^{1/p}\text{d}p$. Then if I let $J:=\displaystyle\int_0^1(1-p)^{1/p}\text{d}p$, I have no idea of how to calculate its exact value...
What I first came up with is the following : $J=\displaystyle\sum_{k=0}^{+\infty}\frac1{k!}\int_0^1\left(\frac{\log(1-p)}p\right)^k\text{d}p$. To evaluate each $\displaystyle\int_0^1\left(\frac{\log(1-p)}p\right)^k\text{d}p=\int_0^1f_k(p)\text{d}p$, I wanted to calculate the Maclaurin series of each $f_k$. I first noticed that $f_k'\equiv-k\Psi f_{k-1}$ where $\Psi(p)=\displaystyle\frac{p+(1-p)\log(1-p)}{p^2(1-p)}=\sum_{m=0}^{+\infty}\frac{m+1}{m+2}p^m$.
Hence, we have : $f_k''\equiv-k\Psi'f_{k-1}+k(k-1)\Psi^2f_{k-2}$, and so on.
Then, I tried $-$ but didn't end up working it out $-$ to get a recursive formula to calculate $f_k^{(r)}$, so I could get the needed Maclaurin series for $f_k$. But even if I had such a formula, I don't even know how to recognize something out of the antiderivative I'd get... Plus I'm not even sure if the antiderivative is or not a spectial function (maybe there's a simple formula using special functions such as $\text{Li}$ or $\text{Ei}$ ?).

Then, I wanted to try a different approach : instead of considering a probability $p$ and a number of trials being $1/p$, why not writing $X_n\hookrightarrow\mathscr{B}(n,1/n)$, and then writting the quantity to average as $\displaystyle1-\left(1-\frac1n\right)^n$, so that we have $\mathcal{I}\displaystyle\overset{?}{=}\lim_{N\to+\infty}\frac1{N-1}\sum_{n=1}^{N}\left[1-\left(1-\frac1n\right)^n\right]$.
Now, the issue is that I have absolutely no idea of how to prove that equality, and I don't even know whether noticing it is useful or not...

Thank you already for your help !
 A: Following on the suggestion of using the same approach as $x^x$ (Sophomore's dream), I ended up having something of the form $1 - \int_0^\infty (e^u-1)^{-n}e^{-u}\,du$ which seems to be non-trivial (there is no simple reduction to a beta or gamma integral that I can see).
We can get a somewhat pretty form in terms of the Sterling numbers of the first kind as follows:
\begin{align}
(1-p)^{1/p}
&= \sum_{n=0}^{\infty}\binom{1/p}{n}(-p)^n 
= \sum_{n=0}^{\infty}\frac{(1/p)_n}{n!} (-1)^n p^n\\
&= \sum_{n=0}^{\infty} \frac{(-1)^n}{n!} \sum_{k=0}^ns(n,k) p^k\; \cdot p^n
\end{align}
We can group all $p^m$ terms by choosing $k=m-n$ for the inner summation, as then $p^kp^n = p^{m-n}p^n = p^m$. So the extracted coefficient is:
$$[p^m](1-p)^{1/p} = \sum_{n=0}^{\infty} \frac{(-1)^n}{n!} s(n,m-n)$$
where $s(n,k)$ are the sterling numbers of the first kind. Now $s(n,k)=0$ for $k<0$, so we can ignore all terms $m-n < 0$, hence
$$[p^m](1-p)^{1/p} = \sum_{n=m}^{\infty} \frac{(-1)^n}{n!} s(n,m-n)
= \sum_{n=0}^{\infty} \frac{(-1)^{m+n}}{(m+n)!} s(m+n,n)$$
So if we justify swapping order of integration and summation, and with a bit of luck, 
$$(1-p)^{1/p} = \sum_{m=0}^{\infty} \sum_{n=0}^{\infty} \frac{(-1)^{m+n}}{(m+n)!} s(m+n,n) p^m$$
then
$$\mathcal{J} = \int_0^1 (1-p)^{1/p}\;dp = \sum_{m=0}^{\infty} \sum_{n=0}^{\infty} \frac{(-1)^{m+n} \; s(m+n,n)}{(m+n)!(m+1)}$$
