Question about substitutions in the double integral $\int_0^1\int_0^1 \frac{-\ln xy}{1-xy} dx\, dy = 2\zeta(3)$ \begin{align}
\int_0^1\int_0^1 \frac{-\ln xy}{1-xy} dx\, dy = 2\zeta(3) \tag{1}
\end{align}
Since,
\begin{align}
\int_0^1 \frac{1}{1-az} dz = \frac{-\ln (1-a)}{a}
\end{align}
and putting $a = 1-xy$, $(1)$ becomes:
\begin{align}
\int_0^1\int_0^1 \int_0^1 \frac{1}{1-z(1-xy)} dx\, dy\, dz = 2\zeta(3) \tag{2}
\end{align}
with $t = 1-z(1-xy) \implies dt = -(1-xy) \, dz$. When $z=0$, $t=1$ and when $z=1$, $t=xy.$ Making these substitutions in $(2)$:
\begin{align}
-\int_1^{xy}\int_0^1 \int_0^1 \frac{1}{t(1-xy)} dx\, dy\, dt \tag{3}
\end{align}
Is $(3)$ still equal to $2\zeta(3)$? I think $(3)$ will be a function of $x$ and $y$, but after making these substitutions the integral should stay the same.

By geometric series, your integral equals
$$-\int^1_0\int^1_0(\ln x+\ln y)\sum^\infty_{k=0}(xy)^kdydx$$
Due to the symmetry in $x$ and $y$, this equals
$$-2\int^1_0\int^1_0(\ln x)\sum^\infty_{k=0}(xy)^kdydx$$
Since the series converges uniformly, we can integrate it termwise to obtain
$$-2\sum^\infty_{k=0}\int^1_0 x^k\ln x\cdot\frac{1}{1+k} dx$$
Since $\int^1_0 x^k\ln xdx=-\frac1{(1+k)^2}$, this simplifies to
$$2\sum^\infty_{k=0}\frac{1}{(1+k)^3} =2\zeta(3)$$ by definition.


Prove $\int^1_0 x^k\ln xdx=-\frac1{(1+k)^2}$:
By the substitution $t=-(k+1)\ln x$,
$$-\int_\infty^0 e^{\frac{-kt}{k+1}}\frac{t}{k+1}\cdot\frac{-1}{1+k}e^{\frac{-t}{k+1}}dt=-\frac1{(k+1)^2}\int^\infty_0te^{-t}dt=-\frac1{(k+1)^2}(1!)=\color{red}{-\frac1{(k+1)^2}}$$

