Proving that $\gamma = \int_{0}^{1} \!\!\int_{0}^{1} \!\frac{x - 1}{(1 - x y) \log(x y)} \, \mathrm{d}{x} \, \mathrm{d}{y} $.
Solution 1:
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.
Solution 2:
This$\require{autoload-all}$ post outlines exactly how you would use the geometric series to reach the series that the original post has given. The beauty of this method lies in it's generalizations which I have placed several examples of in my post. If you want to see an easier method of finding the answer to the specific integral, see my other answer.
I will begin by first defining the exponential integral and then proving some lemmas which will assist in the proof. We write ${\operatorname{Ei}(x)}$ for the exponential integral and it is defined as follows.
$$\operatorname{Ei}(x) = -\int_{-x}^\infty \frac{e^{-t}}{t}\,\mathrm{d}t$$
I will use the following properties (hover over each statement for reasoning):
$$\mathtip{\lim_{x\to -\infty} \operatorname{Ei}(x) = 0}{\text{Clear from the definition of } \operatorname{Ei}(x).} \tag{1a}$$
$$\mathtip{\frac{\mathrm{d}}{\mathrm{d}x} \operatorname{Ei}(x) = \frac{e^x}{x}}{\text{Clear from the definition of } \operatorname{Ei}(x).} \tag{1b}$$
$$\mathtip{\int\operatorname{Ei}(x)\, \mathrm{d}x = x\operatorname{Ei}(x) – e^x + C}{\text{Integrate by parts or take the derivative of the right hand side.}{}} \tag{1c}$$
$$\mathtip{\int \operatorname{Ei}(x)e^{\alpha x}\, \mathrm{d}x = \frac{e^{\alpha x}\operatorname{Ei}(x) -\operatorname{Ei}((\alpha+1)x)}{\alpha} + C}{\text{Integrate by parts or take the derivative of the right hand side.}{}} \tag{1d}$$
I will begin with two preliminary lemmas.
Lemma $1$:
$$\int_0^1 \frac{x^n}{\log(ax)} \mathrm{d}x = a^{-n-1}\operatorname{Ei}((n+1)\log(a)) \tag{2}$$
$$ \toggle{ \enclose{roundedbox}{\text{ Click to Reveal Proof }} }{ \boxed{\text{Substitute } -\log(ax) = u}\\[.2cm] \int_{\infty}^{-\log(a)}-\frac{(e^u/a)^{n+1}}{u} \mathrm{d}u\\[.2cm] \boxed{\text{Simplify and use properties of integration}}\\[.2cm] \int_{-\log(a)}^{\infty}a^{-n-1}\frac{e^{(n+1)u}}{u} \mathrm{d}u\\[.2cm] \boxed{\text{Substitute }t = (n+1)u\text{ and then evaluate using the definition of }\operatorname{Ei}(x).}\\[.2cm] \int_0^1 \frac{x^n}{\log(ax)} \mathrm{d}x = a^{-n-1}\operatorname{Ei}((n+1)\log(a))\\ \enclose{roundedbox}{\text{ Click to Close Proof}}}\endtoggle $$
Lemma $2 $ (very significant):
$$\int_0^1 \dfrac{\operatorname{Ei}(a\log(y))}{y^n} \mathrm{d}y = \lim_{n_0\to n} \frac{1}{n_0-1}\log\left(1+\frac{1-n_0}{a}\right) \tag{3}$$
Proof is located at the bottom of the post. The limit must be used for the case $n = 1$.
Now to evaluate the integral:
$$ \begin{align} I &= \int_{0}^{1}\!\!\int_{0}^{1}\!\dfrac{x-1}{(1 – xy)\log(xy)}\mathrm{d}x\mathrm{d}y \\ &= \lim_{n\to\infty}\int_{0}^{1}\!\!\int_{0}^{1}\!\sum_{k=0}^n \left[\frac{x^{k+1}y^k}{\log(xy)} - \frac{x^{k}y^k}{\log(xy)}\right]\mathrm{d}x\mathrm{d}y \end{align}$$
Using $(2)$:
$$I = \lim_{n\to\infty}\int_{0}^{1} \sum_{k=0}^\infty \left[\frac{\operatorname{Ei}((k+2)\log(y))}{y^2}- \frac{\operatorname{Ei}((k+1)\log(y))}{y}\right] \mathrm{d}y$$
Using $(3)$:
$$I = \lim_{n\to\infty}\sum_{k=0}^n \left[\log\left(1-\frac{1}{k+2}\right)+\frac{1}{k+1}\right]$$
Finally:
$$\begin{align}I &= \lim_{n\to\infty} \sum_{k=1}^n \left[\log\left(\frac{k}{k+1}\right)+\frac{1}{k}\right]\\ &= \lim_{n\to\infty}\sum_{k=1}^n \left[\log(k)-\log(k+1)+\frac{1}{k}\right]\\ &= \lim_{n\to\infty}\sum_{k=1}^{\infty}\frac{1}{k} - \log(n)\end{align}$$
$\LARGE \text{Generalizations:}$
This approach is not very good for the simple integral given in the original post but it is a very strong approach in general. In the proof of Lemma 2, we were able to develop a general formula for $I_n$. Using the same technique, it is possible to evaluate more complicated integrals such as
$$\int_{0}^{1}\!\!\int_{0}^{1}\!\dfrac{(x-1)^2}{(1 – xy)\log(xy)}\mathrm{d}x\mathrm{d}y = \log(\sqrt{2}) - \gamma$$
In fact, if we let $P(x)$ be any polynomial, we can find a closed form for
$$\int_{0}^{1}\!\!\int_{0}^{1}\!\dfrac{(x-1)P(x)}{(1 – xy)\log(xy)}\mathrm{d}x\mathrm{d}y$$
The factor of $(x-1)$ exists to ensure the numerator goes to zero as $x \to 1$. In general, we can do more than just polynomials in $x$. Here are some of the more interesting results that can be obtained
$$\int_{0}^{1}\!\!\int_{0}^{1}\!\dfrac{\sqrt{x}(x-1)}{(1 – xy)\log(xy)}\,\mathrm{d}x\mathrm{d}y = \frac{1}{3} (\log (36)-2 \log (\pi ))$$
$$\int_{0}^{1}\!\!\int_{0}^{1}\!\dfrac{x^3+9x^2-11x+1}{(1 – xy)\log(xy)}\,\mathrm{d}x\mathrm{d}y =\log(16)+\frac{\log(3)}{3}+\frac{5\log(2)}{6} - \gamma\\$$
$$\int_{0}^{1}\!\!\int_{0}^{1}\!\dfrac{\sqrt[4]{y}\sqrt{x}(1-x)}{(1 – xy)\log(xy)}\,\mathrm{d}x\mathrm{d}y =\frac{4}{5} \!\left(\log (24)-2 \log (\pi )+4 \log \left(\Gamma (5/4)\right)\right)$$
$$\int_{0}^{1}\!\!\int_{0}^{1}\!\dfrac{y^{\pi}(x-1)}{(1 – xy)\log(xy)}\,\mathrm{d}x\mathrm{d}y = \frac{\log (\pi )+\log \left(\left(-6+11 \pi -6 \pi ^2+\pi ^3\right) \Gamma (-3+\pi )\right)}{(\pi -1) \pi }$$
$$\int_{0}^{1}\!\!\int_{0}^{1}\!\dfrac{x^{\pi}(x-1)}{(1 – xy)\log(xy)}\,\mathrm{d}x\mathrm{d}y = -\frac{\log (\pi -3)-\pi \log (1+\pi )+\log \left(\pi \left(2-3 \pi +\pi ^2\right) \Gamma (-3+\pi )\right)}{\pi (1+\pi )}$$
The most significant one (in my opinion) is this one:
$$\boxed{\displaystyle\int_{0}^{1}\!\!\int_{0}^{1}\!\dfrac{x^{n}y^{n}(x-1)}{(1 – xy)\log(xy)}\,\mathrm{d}x\mathrm{d}y = \log(1+n)-\psi(1+n)}$$
This is not only a beautiful formula but it allows us to say that the integral of the original question can be easily shown to be $-\psi(1) = \gamma$. In this case, $\psi$ is the polygamma function. See my other answer for a quick route to developing this formula.
Here is the proof of Lemma 2. If you don't understand a step, hover over the equations. Note that $a+1>n$. The final generalization of $I_n$ is very useful for the generalization I later prove.
$$\begin{align}I_n &= \int_0^1 \frac{\operatorname{Ei}(a \log (y))}{y^n} \mathrm{d}y\\[.3cm] &\mathtip{= \int_{-\infty}^{0}\frac{\operatorname{Ei}(u)}{e^{n\cdot u/a}}\frac{e^{u/a}}{a} \mathrm{d}u}{\text{Substitute } u = a\log(y). \text{ and then integrate using property (1d)}}\\[.3cm] &\mathtip{= \left. \frac{e^{(1-n)u/a}\operatorname{Ei}(u)-\operatorname{Ei}((a+1-n)u/a)}{1-n} \right|_{-\infty}^{\,0}}{\text{Integrate using property (1d)}.}\\[.2cm] &\mathtip{=\left. \frac{y^{1-n}\operatorname{Ei}(a\log(y))-\operatorname{Ei}((a+1-n)\log(y))}{1-n} \right|_{\,0}^{\,1}}{\text{Substitute back in for }y.} \end{align}$$ To show that $I_1 = -1/a$ simply take the limit as $n \to 1$ using L'Hopital's rule and then evaluate. To find $I_n$ for $n > 1$ consider the following limits.
$$\begin{align} \text{Lower Limit} &= \lim_{y\to 0}\frac{y^{1-n}\operatorname{Ei}(a\log(y))-\operatorname{Ei}((a+1-n)\log(y))}{1-n} \\[.3cm] &\mathtip{=\lim_{y\to 0}\frac{\operatorname{Ei}(a\log(y))-y^{n-1}\operatorname{Ei}((a+1-n)\log(y))}{(1-n)y^{1-n}}}{\text{Re-arrange terms}.} \\[.3cm] &\mathtip{\!\overset{\rm{L'H}}{=} \lim_{y\to 0}-\frac{\operatorname{Ei}((a+1-n)\log(y))}{(1-n)}}{\text{Apply L'Hopitals rule using (1b)}.}\\[.3cm] &\mathtip{= 0.}{\text{Evaluate using (1d)}.} \end{align}$$
$$\text{Upper Limit} = \lim_{y\to 1}\frac{y^{1-n}\operatorname{Ei}(a\log(y))-\operatorname{Ei}((a+1-n)\log(y))}{1-n}$$
The upper limit is trivial with a simple Taylor series expansion. The upper limit is equal to $$\dfrac{\log(a) - \log(a+1-n)}{1-n}$$ With this,
$$\begin{align}I_n &= \frac{1}{n-1}\log\left(1+\frac{1-n}{a}\right) \\ I_1 &= -\frac{1}{a} \\ I_2 &= \log\left(1-\frac{1}{a}\right)\end{align}$$
This completes the proof of lemma 2.
It is worth noting that $I_1$ can be computed by taking the limit of $I_n$ as $n \to 1$. Furthermore, the above analysis did not depend on $n$ being an integer or being positive. This means that we can write for any $n \in \mathbb{R}$
$$I_n = \lim_{n_0\to n} \frac{1}{n_0-1}\log\left(1+\frac{1-n_0}{a}\right)$$
Solution 3:
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}$$
Solution 4:
$\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.