Ultimate Sophomore's dream $\int_0^1\cdot\cdot\cdot\int_0^1 (t_1t_2\cdot\cdot\cdot t_n)^{t_1t_2\cdot\cdot\cdot t_n} dt_1 dt_2\cdot\cdot\cdot dt_n$ $$I_n=\int_0^1\int_0^1\cdot\cdot\cdot\int_0^1 (t_1t_2\cdot\cdot\cdot t_n)^{t_1t_2\cdot\cdot\cdot t_n} dt_1 dt_2\cdot\cdot\cdot dt_n$$
Prove: $$I_1=I_2<I_3<\cdot\cdot\cdot <1=\lim_{n\rightarrow\infty}I_n$$
For $n=1$, it is Sophomore's dream: https://en.wikipedia.org/wiki/Sophomore%27s_dream
$$ I_1=\int_0^1 t_1^{t_1} dt_1=\sum_{n=1}^\infty \frac{(-1)^{n-1}}{n^n}$$
For $n=2$, $$I_2=\int_0^1\int_0^1 (t_1t_2)^{t_1t_2}dt_1 dt_2$$
let $x=t_1, y=t_1 t_2$, where the integral region is mapping to the triangle region $0<x<1, 0<y<x$, and Jacobian $J=\frac{1}{x}$
$$I_2=\int_0^1\int_0^x y^y \frac{1}{x}dy dx=\int_0^1\int_y^1 y^y \frac{1}{x} dx dy=-\int_0^1 y^y\ln(y) dy$$
Series expansion: $$y^y=e^{y\ln(y)}=\sum_{n=0}^\infty \frac{1}{n!}y^n \ln^n(y)$$
$$I_2=-\sum_{n=0}^\infty \frac{1}{n!}\int_0^1 y^n \ln^{n+1}(y) dy$$
$$\int_0^1 y^n \ln^{n+1}(y) dy=\frac{(-1)^{n+1}n!}{(n+1)^{n+1}}$$
After simplify the result we got:
$$I_2=I_1$$
Next, for $I_3$, let $x=t_1, y=t_1t_2, z=t_1t_2t_3$, Jacobian is $J=\frac{1}{z}$
$$I_3=\int_0^1\int_0^x\int_0^y z^z\frac{1}{z} dzdydx$$
Question.1: Is there an analytic result for $I_3$?
Question.2: How to generalize the calculation to $I_n$, as well as the limit of $I_n$?
 A: Let $(U_n)_{n \ge 1}$ be a sequence of independent uniform random variables. By removing null sets, we may assume that for every $n \ge 1$, $U_n>0$ everywhere
For every $n \ge 1$, set $X_n = -\ln(U_n)$, $R_n = U_1 \cdots U_n$ and $S_n = -\ln R_n = X_1 + \cdots +X_n$.
Then $$I_n = E[(U_1 \cdots U_n)^{U_1 \cdots U_n}] = E[R_n^{R_n}] = E[\exp(R_n \ln(R_n))] = E[\exp(-S_ne^{-S_n})].$$
The non-negative sequence $(R_n)_{n \ge 1}$ is non-decreasing, so it converges necessarily to $0$ (since Borel Cantelli lemma prevents $(U_n)_{n \ge 1}$ from converging to $1$). Hence $R_n^{R_n} \to 1$ as $n \to +\infty$. Since $0 \le R_n^{R_n} \le 1$. Lebesgue dominated convergence applies so $I_n \to 1$ as $n \to +\infty$.
By construction, the random variables $(X_n)_{n \ge 1}$ are independent random variables with distribution Exponential$(1)$, so the distribution of $S_n$ is Gamma$(n,1)$. Hence
$$I_n =  \int_0^\infty \exp(-se^{-s}) \frac{s^{n-1}e^{-s}}{(n-1)!}ds.$$
For all $s \ge 0$, $$\exp(-se^{-s}) = \sum_{k=0}^\infty \frac{(-se^{-s})^k}{k!}.$$
This series is normally convergent with regard to $s \ge 0$. Since we integrate with regard to a probability measure, we derive
\begin{eqnarray*}
I_n 
&=&  \sum_{k=0}^\infty \int_0^\infty \frac{(-se^{-s})^k}{k!} \frac{s^{n-1}e^{-s}}{(n-1)!}ds \\
&=& \sum_{k=0}^\infty \frac{(-1)^k}{k!(n-1)!} \int_0^\infty s^{k+n-1}e^{-(k+1)s}ds \\
&=& \sum_{k=0}^\infty \frac{(-1)^k}{k!(n-1)!} \frac{(k+n-1)!}{(k+1)^{k+n}}.\end{eqnarray*}
\begin{eqnarray*}
I_{n+1}-I_n 
&=& \sum_{k=0}^\infty \frac{(-1)^k}{k!} \frac{(k+n-1)!}{n!(k+1)^{k+n+1}} (k+n - n(k+1)) \\
&=& \frac{n-1}{n!}\sum_{k=1}^\infty \frac{(-1)^{k-1}}{(k-1)!} \frac{(k+n-1)!}{(k+1)^{k+n+1}}  \end{eqnarray*}
The last sum of alternating series is strictly positive since for all $k \ge 1$,
\begin{eqnarray*}
\Big(\frac{1}{k!} \frac{(k+n)!}{(k+2)^{k+n+2}} \Big) \Big/ \Big( \frac{1}{(k-1)!} \frac{(k+n-1)!}{(k+1)^{k+n+1}} \Big) 
&=& \frac{k+n}{k(k+2)} \Big/\Big(1+\frac{1}{k+1}\Big)^{k+n+1} \\
&\le& \frac{k+n}{k(k+2)} \Big/\Big(1+\frac{k+n+1}{k+1}\Big) \\
&=& \frac{k+n}{k(k+2)} \times\frac{k+1}{2k+n+1}\\
&<& 1.
\end{eqnarray*}
We are done.
A: We can get the result from the fact that for any function $f$, we have
$$
\int_0^{x}\frac{1}{y}\int_0^{y}\ln\left(\frac{z}{y}\right)^nf(z)dz dy =\int_0^x\int_0^1\ln(u)^nf(y u)du dy = \int_0^1\frac{\ln(u)^n}{u}\int_0^{xu}f(v)dv du \\ = \left.\left[\frac{\ln(u)^{n+1}}{n+1}\int_0^{x u}f(v)dv\right]\right|^1_0 - \int_0^1\frac{\ln(u)^{n+1}}{n+1}f(x u)du = -\frac{1}{n+1}\int_0^x\ln\left(\frac{y}{x}\right)^{n+1}f(y)dy.
$$
$I_{n+1}$ can be written as
$$
I_{n+1} = \int_0^1\frac{1}{x_n}\int_0^{x_n}\frac{1}{x_{n-1}}\int_0^{x_{n-1}}...\int_0^{x_2}\frac{1}{x_1}\int_0^{x_1}x_0^{x_0}dx_0dx_1...dx_{n-1}dx_n.
$$
Applying the above identity $n$ times then gives
$$
I_{n+1} = \frac{(-1)^{n}}{n!}\int_0^1\ln(x)^n x^xdx.
$$
Changing variables to $x = \exp(-u)$ gives an alternate way to write it that makes it more clear why $I_{n+1}\rightarrow 1$:
$$
I_{n+1} = \frac{1}{n!}\int_0^\infty x^ne^{-x}e^{-xe^{-x}}dx. 
$$
$e^{-xe^{-x}}$ approaches $1$ as $x\rightarrow\infty$, and $x^n e^{-x}$ has its largest contribution to the integral near $x = n$. So as $n$ gets larger, this integral effectively becomes the gamma integral, equal to $n!$.
I suppose there is one more thing to this. If we want a fancy sum, we have
$$
I_{n+1} = \frac{1}{n!}\int_0^\infty x^ne^{-x}e^{-xe^{-x}}dx = \sum_{k=1}^\infty \int_0^\infty \frac{(-1)^{k-1} x^{n+k-1} e^{-k x}}{n!(k-1)!}dx \\= -\sum_{k=1}^\infty\frac{(-1)^k(n+k-1)!}{k^{n+k}n! (k-1)!} = -(-1)^n\sum_{k=1}^\infty \binom{n+k-1}{n}(-k)^{-(n+k)}
$$
