What substitution would I make to integrate this? Problem:

The integral is gonna be:
$$\int^{1}_{0}\int^{1}_{0}4xy\sqrt{x^2+y^2} dy \, dx$$
But I'm quite rusty with my calculus, and this is mainly for a statistics course. I know that $x$ and $y$ are bounded within a unit square region in the $xy$-plane, and I was considering making some sort of $\cos$ or $\sin$ substitution to get rid of that pesky square root, but then I'm not sure what limits to use, and I suspect the $dy dx$ would just become $d\theta d\theta$, which doesn't make sense.
Any help or guidance is appreciated. Please try to keep your explanation simple, and don't assume I understand things. It's been a while since my last calculus course :) Thank you!

An alternative approach:
$$ \begin{eqnarray*}\iint_{(0,1)^2}4xy\sqrt{x^2+y^2}\,dx\,dy&=&\iint_{(0,1)^2}\sqrt{X+Y}\,dX\,dY\\&=&2\iint_{0\leq Y\leq X\leq 1}\sqrt{X+Y}\,dX\,dY\\&=&2\iint_{(0,1)^2}X\sqrt{X}\sqrt{K+1}\,dX\,dK\\&=&2\int_{0}^{1}X\sqrt{X}\,dX\int_{0}^{1}\sqrt{K+1}\,dK\end{eqnarray*}$$
by exploting the substitution $x=\sqrt{X}, y=\sqrt{Y}$, then symmetry, then the substitution $Y=KX$, then Fubini's theorem. The last integrals are elementary and they respectively equal $\frac{2}{5}$ and $\frac{2}{3}(2\sqrt{2}-1)$.
