Seemingly impossible double integral reduction Show that $\int_{0}^{1}\int_{0}^{1}(xy)^{xy}dxdy = \int_{0}^{1}y^{y}dy$
I have tried the method used in the Gaussian integral, polar coordinates, exp(.) and ln(.) of both sides of the integrand...
Gaussian integral: https://en.wikipedia.org/wiki/Gaussian_integral
I have also looked at the sophomore dream function, but I am not sure how to proceed that way in this case.
Sophomore dream function: https://en.wikipedia.org/wiki/Sophomore%27s_dream 
Now I am stuck. Any help would be appreciated.
 A: $$
\begin{align}
\int_0^1\int_0^1(xy)^{axy}\,\mathrm{d}x\,\mathrm{d}y
&=\int_0^1\int_0^yx^{ax}\,\mathrm{d}x\frac1y\,\mathrm{d}y\tag{1}\\
&=\int_0^1\int_x^1\frac1y\,\mathrm{d}y\,x^{ax}\,\mathrm{d}x\tag{2}\\
&=\int_0^1(-\log(x))\,x^{ax}\,\mathrm{d}x\tag{3}\\
&=\int_0^1x^{ax}\,\mathrm{d}(x-x\log(x))\tag{4}\\
&=\int_0^1x^{ax}\,\mathrm{d}x-\int_0^1e^{ax\log(x)}\,\mathrm{d}x\log(x)\tag{5}\\
&=\int_0^1x^{ax}\,\mathrm{d}x-\frac1a\left[\,\vphantom{\int}x^{ax}\,\right]_0^1\tag{6}\\
&=\int_0^1x^{ax}\,\mathrm{d}x\tag{7}
\end{align}
$$
Explanation:
$(1)$: substitute $x\mapsto\frac xy$
$(2)$: change order of integration
$(3)$: integrate in $y$
$(4)$: $\mathrm{d}(x-x\log(x))=-\log(x)\,\mathrm{d}x$
$(5)$: separate integrals, noting that $x^{ax}=e^{ax\log(x)}$
$(6)$: substitute $u=x\log(x)$, then $\int e^{ax\log(x)}\,\mathrm{d}x\log(x)=\int e^{au}\,\mathrm{d}u=\frac1ae^{au}{+}C=\frac1ax^{ax}{+}C$
$(7)$: $\left[\,\vphantom{\int}x^{ax}\,\right]_0^1=0$
Setting $a=1$ yields the equation in the question
$$
\int_0^1\int_0^1(xy)^{xy}\,\mathrm{d}x\,\mathrm{d}y=\int_0^1x^x\,\mathrm{d}x
$$
Setting $a=-1$ yields a similar equation for the other Sophomore Dream integral:
$$
\int_0^1\int_0^1(xy)^{-xy}\,\mathrm{d}x\,\mathrm{d}y=\int_0^1x^{-x}\,\mathrm{d}x
$$
A: Suppose that we have a function $f: (x,y) \mapsto f(x,y)$ that we wish to integrate over a domain $D$. Our integral looks like:
$$\iint_D f(x,y) dy dx $$
If we want to do the integral by a change of variables, what we are really doing is that we are moving from one coordinate system into another. 
If our new variables are $u = u(x,y)$ and $v = v(x,y)$, then $f(x,y)$ should be replaced by $f(u,v)$, and the term $dy dx$ now becomes $|J| dudv$, where $J$ is called the Jacobian of the transformation; in such a situation, it's given by:
$$J = \det \left( \begin{matrix} \partial x/ \partial u && \partial x / \partial v \\ \partial y / \partial u && \partial y/ \partial v \end{matrix} \right)$$
And the domain $D$ of integration becomes another one (say $\Delta$). The integral should finally look like:
$$\iint_{\Delta} f(u,v) |J| dv du $$
In your case, the transformation is:
$$\begin{cases} u = xy \\ v = y \end{cases}$$
The Jacobian is $1/v$. The new domain is $\Delta = [0,1] \times [0,v]$, and your integral is:
$$I = \int_0^1 \int_0^v u^u \frac1{v} dv du = \int_0^1 (\int_0^v u^u du) \frac1v dv = \int_0^1 f(v) \frac1v dv$$
Let $S$ be the second version of "Sophomore's Dream". Integrating by parts, 
$$I = f(v) \log v \mid_0 ^1 - \int_0^1 v^v \log v dv = - \int_0^1 v^v \log v dv = - \int_0^1 v^v (\log v + 1 - 1) dv = - \int_0^1 e^{v \log v} (\log v + 1) dv + S =- \int_0^1 e^{v \log v} d(v \log v) + S = - e^{v \log v} \mid_0 ^1 + S  = S = \sum_{n=1}^{\infty} \frac{(-1)^{n+1}}{n^n}$$
A: $$I=\int_{0}^{1}\int_{0}^{1}(xy)^{xy}dxdy\ \ \ \ \ \ \ \ \ \ \ , let\ xy=t\ \Rightarrow \ x=\frac{t}{y}\\
\\
\\
\therefore I=\int_{0}^{1}[\int_{0}^{1}(xy)^{xy}dx]dy=\int_{0}^{1}(\frac{1}{y}.\int_{0}^{y}t^tdt)dy\\
\\
\\
=\int_{0}^{1}(\frac{1}{y}.\int_{0}^{y}x^xdx)dy=\int_{0}^{1}\int_{0}^{y}\frac{x^x}{y}dxdy\\
\\
\\
=\int_{0}^{1}\int_{x}^{1}\frac{x^x}{y}dydx=\int_{0}^{1}(x^x.\int_{x}^{1}\frac{dy}{y})dx\\
\\
\\
=-\int_{0}^{1}x^x.ln(xdx=\int_{0}^{1}[x^x-\frac{d}{dx}(x^X)]dx\\
\\
\\
=\int_{0}^{1}x^xdx-\int_{0}^{1}d(x^X)=\int_{0}^{1}x^xdx-(x^X)_{0}^{1}=\int_{0}^{1}x^xdx\\
\\$$
now let evaluate this integral  $$\int_{0}^{1}x^xdx$$
$$\because \ \ \ x^x=e^{xln(x)}=e^{-xln(\frac{1}{x})}=\sum_{n=0}^{\infty }\frac{(-1)^n}{n!}.x^n.(ln(\frac{1}{x}))^n\\
\\
\\
\therefore I=\int_{0}^{1}x^xdx=\sum_{n=0}^{\infty }\frac{(-1)^n}{n!}\int_{0}^{1}x^n.(ln(\frac{1}{x}))^ndx\\
\\
\\
=\sum_{n=0}^{\infty }\frac{(-1)^n}{n!}\int_{0}^{\infty }y^n.e^{-(1+n)y}dy=\sum_{n=0}^{\infty }\frac{(-1)^n}{n!}.\frac{n!}{(n+1)^n}\\
\\
\\
=\sum_{n=1}^{\infty }\frac{(-1)^{n-1}}{n^n}=\int_{0}^{1}\int_{0}^{1}(xy)^{xy}dxdy\\$$
