Proving that $\gamma = \int_{0}^{1} \!\!\int_{0}^{1} \!\frac{x - 1}{(1 - x y) \log(x y)} \, \mathrm{d}{x} \, \mathrm{d}{y} $. In 2005, J. Sondow found a surprising formula for the Euler-Mascheroni constant $ \gamma $. The formula is
$$
  \gamma
= \int_{0}^{1} \int_{0}^{1} \frac{x - 1}{(1 - x y) \log(x y)}
  ~ \mathrm{d}{x} ~ \mathrm{d}{y}.
$$
Now, the definition of $ \gamma $ is
$$
\gamma
\stackrel{\text{def}}{=}
\lim_{n \to \infty} \left[ \sum_{k = 1}^{n} \frac{1}{k} - \log(n) \right].
$$
I have tried using the geometric series
$$
\frac{1}{1 - x y} = \sum_{n = 0}^{\infty} x^{n} y^{n}
$$
to obtain a proof, but it would not work. Thanks for any help.
 A: $\newcommand{\bbx}[1]{\,\bbox[15px,border:1px groove navy]{\displaystyle{#1}}\,}
 \newcommand{\braces}[1]{\left\lbrace\,{#1}\,\right\rbrace}
 \newcommand{\bracks}[1]{\left\lbrack\,{#1}\,\right\rbrack}
 \newcommand{\dd}{\mathrm{d}}
 \newcommand{\ds}[1]{\displaystyle{#1}}
 \newcommand{\expo}[1]{\,\mathrm{e}^{#1}\,}
 \newcommand{\ic}{\mathrm{i}}
 \newcommand{\mc}[1]{\mathcal{#1}}
 \newcommand{\mrm}[1]{\mathrm{#1}}
 \newcommand{\pars}[1]{\left(\,{#1}\,\right)}
 \newcommand{\partiald}[3][]{\frac{\partial^{#1} #2}{\partial #3^{#1}}}
 \newcommand{\root}[2][]{\,\sqrt[#1]{\,{#2}\,}\,}
 \newcommand{\totald}[3][]{\frac{\mathrm{d}^{#1} #2}{\mathrm{d} #3^{#1}}}
 \newcommand{\verts}[1]{\left\vert\,{#1}\,\right\vert}$
\begin{align}
&\bbox[10px,#ffd]{\ds{\int_{0}^{1}\int_{0}^{1}{x - 1 \over \pars{1 - xy}
\ln\pars{xy}}\,\dd x\,\dd y}} =
\int_{0}^{1}{1 \over y^{2}}\int_{0}^{1}{xy - y \over \pars{1 - xy}
\ln\pars{xy}}\,y\,\dd x\,\dd y
\\[5mm] = &\
\int_{0}^{1}{1 \over y^{2}}\int_{0}^{y}{x - y \over \pars{1 - x}
\ln\pars{x}}\,\dd x\,\dd y =
\int_{0}^{1}{1 \over \pars{1 - x}
\ln\pars{x}}\int_{x}^{1}{x - y \over y^{2}}\,\dd y\,\dd x
\\[5mm] = &\
\int_{0}^{1}{1 \over \pars{1 - x}
\ln\pars{x}}\bracks{1 - x + \ln\pars{x}}\dd x =
\int_{0}^{1}{1 - x + \ln\pars{x} \over 1 - x}
\overbrace{\bracks{-\int_{0}^{\infty}x^{t}\,\dd t}}
^{\ds{1 \over \ln\pars{x}}}\ \dd x
\\[5mm] = &\
\int_{0}^{\infty}\int_{0}^{1}{x^{t + 1} - x^{t} - x^{t}\ln\pars{x} \over 1 - x}
\,\dd x\,\dd t\label{1}\tag{1}
\end{align}

\begin{align}
&\bbox[10px,#ffd]{\ds{%
\int_{0}^{1}{x^{t + 1} - x^{t} - x^{t}\ln\pars{x} \over 1 - x}\,\dd x}} =
\left.\partiald{}{\mu}\int_{0}^{1}{x^{t + 1}\mu - x^{t}\mu -
x^{t}\pars{x^{\mu} - 1} \over 1 - x}\,\dd x\,\right\vert_{\ \mu\ =\ 0}
\\[5mm] = &\
\partiald{}{\mu}\bracks{%
\mu\int_{0}^{1}{1 - x^{t} \over 1 - x}\,\dd x -
\mu\int_{0}^{1}{1 - x^{t + 1} \over 1 - x}\,\dd x +
\int_{0}^{1}{1 - x^{t + \mu} \over 1 - x}\,\dd x -
\int_{0}^{1}{1 - x^{t} \over 1 - x}\,\dd x}_{\ \mu\ =\ 0}
\\[5mm] = &\
\Psi\pars{t + 1} - \Psi\pars{t + 2} + \Psi\, '\pars{t + 1} =
-\,{1 \over t + 1} + \Psi\,'\pars{t + 1}\label{2}\tag{2}
\end{align}



where $\ds{\Psi}$ is the Digamma Function.
  With \eqref{1} and \eqref{2}:

\begin{align}
&\bbox[10px,#ffd]{\ds{\int_{0}^{1}\int_{0}^{1}{x - 1 \over \pars{1 - xy}
\ln\pars{xy}}\,\dd x\,\dd y}} =
\lim_{R \to \infty}\int_{0}^{R}
\bracks{-\,{1 \over t + 1} + \Psi\,'\pars{t + 1}}\dd t
\\[5mm] = &\
\lim_{R \to \infty}\bracks{-\ln\pars{R + 1} + \Psi\pars{R + 1} - \Psi\pars{1}} =
\bbx{\gamma}
\end{align}

Note that $\ds{\Psi\pars{1} = -\gamma}$ and
  $\ds{\Psi\pars{z} \sim \ln\pars{z} - {1 \over 2z} - {1 \over 12z^{2}}}$ as
  $\ds{\verts{z} \to\ \infty}$ with $\ds{\verts{\arg{z}} < \pi}$.
  
  See A & S Table.

A: Here is an approach that reduces to a single integral:
Notice that:
$$\int_{0}^{1}\frac{x}{(1-xy)\ln(xy)}dy=\int_{0}^{x}\frac{1}{(1-xy)\ln(xy)}d(xy)=\int_{0}^{x}\frac{1}{(1-t)\ln t}dt......(1)$$
So
$$I =\int_{0}^{1}\int_{0}^{1}\frac{x-1}{(1-xy)\ln(xy)}dxdy=\lim_{\epsilon \rightarrow 0}\int_{0}^{1-\epsilon}\frac{x-1}{x}\left[\int_{0}^{x}\frac{1}{(1-t)\ln t}dt\right]dx.$$
Integrating by parts, 
$$I =\lim_{\epsilon \rightarrow 0}\left[(x-\ln x)\int_{0}^{x}\frac{1}{(1-t)\ln t}dt\Big|_{0}^{1-\epsilon}-\int_{0}^{1-\epsilon}\frac{x-\ln x}{(1-x)\ln x}dx\right]= \\ \lim_{\epsilon \rightarrow 0}\left[[1-\epsilon-\ln(1-\epsilon)]\int_{0}^{1-\epsilon}\frac{1}{(1-x)\ln x}dx-\int_{0}^{1-\epsilon}\frac{x-\ln x}{(1-x)\ln x}dx\right]= \\\lim_{\epsilon \rightarrow 0}\left[-[\epsilon+\ln(1-\epsilon)]\int_{0}^{1-\epsilon}\frac{1}{(1-x)\ln x}dx+\int_{0}^{1-\epsilon}\left[\frac1{1-x}+\frac1{\ln x}\right]dx\right]=\int_{0}^{1}\left[\frac1{1-x}+\frac1{\ln x}\right]dx = \gamma,$$
where the last single-integral representation of $\gamma$ is well known.
A: I feel like I must submit another answer to this question as I have found a much simpler solution.
Let $$I'(n) = \int_0^1\!\!\int_0^1 \frac{x^ny^n(x-1)}{1-xy} \,\mathrm{d}x\mathrm{d}y$$
Using the geometric series, integrating, and then using partial fractions
$$\begin{align}I'(n) &= \lim_{N\to\infty}\sum_{k=0}^{N}-\frac{1}{(k+n)(1+k+n)^2}\\
&= \lim_{N\to\infty}\sum_{k=0}^{N} \frac{1}{k+n+1}-\frac{1}{k+n} +  \sum_{k=0}^{N} \frac{1}{(k+n+1)^2}
\end{align}$$
The first sum telescopes and the second one can be evaluated using the polygamma function. We have that
$$\psi^{(n)}(z) = \frac{\mathrm d^{n+1}}{\mathrm dz^{n+1}}\log \Gamma(z)= (-1)^{n+1}n! \sum_{k=0}^{\infty} \frac{1}{(z+k)^{n+1}}$$
Use this to find and expression for $I'(n)$ and then integrate with respect to $n$.
$$I'(n) =\int_0^1\!\!\int_0^1 \frac{x^ny^n(x-1)}{1-xy} \,\mathrm{d}x\mathrm{d}y = \frac{1}{n+1} - \psi^{(1)}(n+1) $$
$$I(n) = \int_0^1\!\!\int_0^1 \frac{x^ny^n(x-1)}{(1-xy)\log(xy)} \,\mathrm{d}x\mathrm{d}y = \log(n+1) - \psi(n+1) $$
$$\boxed{\displaystyle I(0) = \int_0^1\!\!\int_0^1 \frac{x-1}{(1-xy)\log(xy)} \,\mathrm{d}x\mathrm{d}y = - \psi(1)  = \gamma}$$
A: We want use the transformation formula.
\begin{align*}
&U:= \left\lbrace (t,r) \mid t \in (0,1), r \in (t,1) \right\rbrace , V:= (0,1)^2,
\\ &f: V \rightarrow \mathbb{R}, \left( \begin{array}{r}
x\\
y\\
\end{array} \right) \mapsto \frac{x - 1}{(1-xy) \cdot ln(xy)},
\\ &\varphi: U \rightarrow V, \left( \begin{array}{r}
t\\
r\\
\end{array}\right) \mapsto \left( \begin{array}{r}
t/r\\
r\\
\end{array} \right), 
\end{align*}
\begin{align*}
\varphi ^{-1}: V \rightarrow U, \left( \begin{array}{r}
t\\
r\\
\end{array}\right) \mapsto \left( \begin{array}{r}
t \cdot r\\
r\\
\end{array} \right)   
\end{align*}
\begin{align*}
D\varphi = \left( \begin{array}{rr}
1/r & -t/r^2\\
0 & 1\\
\end{array}\right),
Det(D\varphi) = \frac{1}{r}. 
\end{align*}
So we can use the transformation. 
\begin{align*}
\int_0^1 \int_0^1 \frac{x-1}{(1-xy)\cdot ln(xy)} dx \, dy &= \int_0^1 \int_t^1 \frac{(t/r) - 1}{(1-t) \cdot ln(t)} \cdot \frac{1}{r} \, dr \, dt 
\\ &= \int_0^1 \int_t^1 \left( \frac{t}{(1-t) \cdot ln(t)} \cdot \frac{1}{r^2} - \frac{1}{r} \cdot \frac{1}{(1-t) \cdot ln(t)} \right) \, dr \, dt 
\\ &= \int_0^1 \left[ -\frac{1}{r} \cdot \frac{t}{(1-t) \cdot ln(t)} - \frac{ln(r)}{(1-t)\cdot ln(t)} \right]_t^1 \, dt 
\\ &= \int_0^1 \left( \frac{1-t}{(1-t)\cdot ln(t)} + \frac{ln(t)}{(1-t) \cdot ln(t)} \right) \, dt 
\\ &= \int_0^1 \left( \frac{1}{ln(t)} + \frac{1}{1-t} \right) \, dt
\end{align*}
The last integral represented $\gamma$.
