Evaluate $\int_{0}^{1}{\int_{0}^{1}{\frac{x^{\alpha -1}y^{\beta -1}}{\left( 1+xy \right)\ln xy}dxdy}}$ On AoPS (Art of Problem Solving), the following integral was posted by fellow user pprime, but no one has been able to come close to solving it. Many attempts have been made, but none have succeeded. I have come here as a last resort so that some kind user would help me evaluate this integral.
http://www.artofproblemsolving.com/community/c7h1288609_extreme_integration_marathon
Evaluate $$\int_{0}^{1}{\int_{0}^{1}{\frac{x^{\alpha -1}y^{\beta -1}}{\left( 1+xy \right)\ln xy}dxdy}}$$
Hint: $\int_{0}^{+\infty }{\exp \left( -u\ln \left( xy \right) \right)du}=\frac{1}{\ln x+\ln y}$
 A: Just to note that with Frullani's theorem we can arrive at the infinite product rather quickly. We have that $$I(\alpha,\beta)=\int_{0}^{1}\int_{0}^{1}\frac{x^{\alpha-1}y^{\beta-1}}{\left(1+xy\right)\log\left(xy\right)}dxdy\stackrel{xy=u}{=}\int_{0}^{1}y^{\beta-\alpha-1}\int_{0}^{y}\frac{u^{\alpha-1}}{\left(1+u\right)\log\left(u\right)}dudy
 $$ $$\stackrel{IBP}{=}\frac{1}{\beta-\alpha}\int_{0}^{1}\frac{y^{\alpha-1}-y^{\beta-1}}{\left(1+y\right)\log\left(y\right)}dy=\frac{1}{\beta-\alpha}\sum_{k\geq0}\left(-1\right)^{k}\int_{0}^{1}\frac{y^{\alpha-1+k}-y^{\beta-1+k}}{\log\left(y\right)}dy$$ $$\stackrel{y=e^{-v}}{=}\frac{1}{\beta-\alpha}\sum_{k\geq0}\left(-1\right)^{k+1}\int_{0}^{\infty}\frac{e^{-v(\alpha+k)}-e^{-v(\beta+k)}}{v}dv
 $$ and now we can use the Frullani's theorem and get $$I(\alpha,\beta)=\frac{1}{\beta-\alpha}\sum_{k\geq0}\left(-1\right)^{k+1}\log\left(\frac{\beta+k}{\alpha+k}\right)=\frac{1}{\beta-\alpha}\log\left(\prod_{k\geq0}\left(\frac{\beta+2k+1}{\alpha+2k+1}\right)\left(\frac{\beta+2k}{\alpha+2k}\right)^{-1}\right)
 $$ now note that $$\prod_{k=0}^{N}\left(a+2k\right)=2^{N}\prod_{k=0}^{N}\left(\frac{a}{2}+k\right)=2^{N}a\left(\frac{a}{2}+1\right)_{N}
 $$ and $$\prod_{k=0}^{N}\left(a+2k+1\right)=2^{N}\prod_{k=0}^{N}\left(k+\frac{a+1}{2}\right)=\left(a+1\right)2^{N}\left(\frac{a+1}{2}+1\right)_{N}
 $$ where $\left(x\right)_{N}
 $ is the Pochhammer symbol. So $$\prod_{k=0}^{N}\left(\frac{\beta+2k+1}{\alpha+2k+1}\right)\left(\frac{\beta+2k}{\alpha+2k}\right)^{-1}=\frac{\beta+1}{\alpha+1}\frac{\alpha}{\beta}\frac{\left(\frac{\beta+1}{2}+1\right)_{N}}{\left(\frac{\alpha+1}{2}+1\right)_{N}}\frac{\left(\frac{\alpha}{2}+1\right)_{N}}{\left(\frac{\beta}{2}+1\right)_{N}}
 $$ and now using the asymptotic for the Pochhammer symbol we get 

$$I(\alpha,\beta)=\color{red}{\frac{1}{\beta-\alpha}\log\left(\frac{\Gamma\left(\frac{\alpha+1}{2}\right)\Gamma\left(\frac{\beta}{2}\right)}{\Gamma\left(\frac{\alpha}{2}\right)\Gamma\left(\frac{\beta+1}{2}\right)}\right)}$$ where $$\, \beta,\alpha>0,\,\beta\neq \alpha$$ 

as wanted.
If $\alpha=\beta=\gamma
 $ we have $$I(\gamma)=\int_{0}^{1}\int_{0}^{1}\frac{\left(xy\right)^{\gamma-1}}{\left(1+xy\right)\log\left(xy\right)}dxdy\stackrel{xy=u}{=}\int_{0}^{1}y^{-1}\int_{0}^{y}\frac{u^{\alpha-1}}{\left(1+u\right)\log\left(u\right)}dudy
 $$ $$\stackrel{IBP}{=}-\int_{0}^{1}\frac{y^{\gamma-1}}{1+y}dy=\sum_{k\geq0}\left(-1\right)^{k+1}\frac{1}{\gamma+k}=\color{red}{-\Phi\left(-1,1,\gamma\right)}
 $$ where $\gamma>0$ and $\Phi\left(x,y,z\right)$ is the Lerch Transcendent. Obviously for some special values of $\gamma$ we have nice closed forms. For example for $\gamma=1$ we have $$I\left(1\right)=\sum_{k\geq0}\left(-1\right)^{k+1}\frac{1}{1+k}=-\log\left(2\right).$$
A: $\newcommand{\angles}[1]{\left\langle\,{#1}\,\right\rangle}
 \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{\half}{{1 \over 2}}
 \newcommand{\ic}{\mathrm{i}}
 \newcommand{\iff}{\Longleftrightarrow}
 \newcommand{\imp}{\Longrightarrow}
 \newcommand{\Li}[1]{\,\mathrm{Li}_{#1}}
 \newcommand{\mc}[1]{\mathcal{#1}}
 \newcommand{\mrm}[1]{\mathrm{#1}}
 \newcommand{\ol}[1]{\overline{#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{\ul}[1]{\underline{#1}}
 \newcommand{\verts}[1]{\left\vert\,{#1}\,\right\vert}$
\begin{align}
&\color{#f00}{\int_{0}^{1}\int_{0}^{1}{x^{\alpha - 1}\,\,y^{\beta - 1} \over
\pars{1 + xy}\ln\pars{xy}}\,\dd x\,\dd y} =
-\int_{0}^{1}\int_{0}^{1}{x^{\alpha - 1}\,\,y^{\beta - 1} \over 1 + xy}\,\,\,
\overbrace{\int_{0}^{\infty}\pars{xy}^{z}\,\dd z}
^{\ds{-\,{1 \over \ln\pars{xy}}}}\ \,\dd x\,\dd y
\\[5mm] = &\
-\int_{0}^{\infty}\int_{0}^{1}\int_{0}^{1}
{x^{\alpha -1 + z}\,\,\,y^{\beta - 1 + z} \over 1 + xy}
\,\dd x\,\dd y\,\dd z =
-\int_{0}^{\infty}\int_{0}^{1}\int_{0}^{1}
{\pars{xy}^{\alpha -1 + z}\,\,\,y^{\beta - \alpha - 1} \over 1 + xy}
\,\dd\pars{xy}\,\dd y\,\dd z
\\[5mm] = &\
-\int_{0}^{\infty}\int_{0}^{1}\int_{0}^{y}
{x^{\alpha -1 + z}\,\,\,y^{\beta - \alpha - 1} \over 1 + x}
\,\dd x\,\dd y\,\dd z =
-\int_{0}^{\infty}\int_{0}^{1}{x^{\alpha -1 + z}\,\,\, \over 1 + x}\int_{x}^{1}
y^{\beta - \alpha - 1}\,\,\,\dd y\,\dd x\,\dd z
\\[5mm] = &\
{1 \over \beta - \alpha}\int_{0}^{\infty}\int_{0}^{1}
{x^{\beta + z - 1}\,\,\, -\,\,\, x^{\alpha + z - 1}\,\,\, \over 1 + x}
\,\dd x\,\dd z
\\[5mm] = &\
{1 \over \beta - \alpha}\int_{0}^{\infty}\int_{0}^{1}
{x^{\beta + z - 1}\,\,\, -\,\,\, x^{\alpha + z - 1}\,\,\, -\,\,\,
x^{\beta + z}\,\,\, +\,\,\, x^{\alpha + z} \over 1 - x^{2}}
\,\dd x\,\dd z
\\[5mm] \stackrel{x^{2}\ \mapsto\ x}{=}\ &\
{1 \over 2\pars{\beta - \alpha}}
\int_{0}^{\infty}\int_{0}^{1}
{x^{\beta/2 + z/2 - 1}\,\,\, -\,\,\, x^{\alpha/2 + z/2 - 1}\,\,\, -\,\,\,
x^{\beta/2 + z/2 - 1/2}\,\,\, +\,\,\, x^{\alpha/2 + z/2 - 1/2} \over 1 - x}
\,\,\,\,\dd x\,\dd z
\end{align}

With the Digamma Function identity $\pars{~\gamma\ \mbox{is the Euler-Mascheroni Constant}~}$
$$\left.\Psi\pars{\xi} =
-\gamma + \int_{0}^{1}{1 - t^{\xi - 1} \over 1 - t}\,\dd t\,
\right\vert_{\ \Re\pars{\xi}\ >\ 0\,\,\,\,\,}
$$
the integration is reduced to:
\begin{align}
&\color{#f00}{\int_{0}^{1}\int_{0}^{1}{x^{\alpha - 1}\,\,y^{\beta - 1} \over
\pars{1 + xy}\ln\pars{xy}}\,\dd x\,\dd y}
\\[5mm] = &\
{1 \over 2\pars{\beta - \alpha}}\int_{0}^{\infty}\bracks{%
\Psi\pars{z + \alpha \over 2} + \Psi\pars{z + \beta + 1 \over 2} -
\Psi\pars{z + \beta \over 2} - \Psi\pars{z + \alpha + 1 \over 2}}\,\dd z
\end{align}

Since
$\ds{\Psi\pars{\xi}\ \stackrel{\mrm{def.}}{=}\ 
\totald{\ln\pars{\Gamma\pars{\xi}}}{\xi}}$
\begin{align}
&\color{#f00}{\int_{0}^{1}\int_{0}^{1}{x^{\alpha - 1}\,\,y^{\beta - 1} \over
\pars{1 + xy}\ln\pars{xy}}\,\dd x\,\dd y} =
\left.{1 \over \beta - \alpha}
\ln\pars{\Gamma\pars{z/2 + \alpha/2}\Gamma\pars{z/2 + \beta/2 + 1/2} \over \Gamma\pars{z/2 + \beta/2}\Gamma\pars{z/2 + \alpha/2 + 1/2}}
\right\vert_{\ 0}^{\ \infty}
\\[5mm] = &\
\color{#f00}{{1 \over \beta - \alpha}\,
\ln\pars{\Gamma\pars{\bracks{\alpha + 1}/2}\Gamma\pars{\beta/2} \over \Gamma\pars{\alpha/2}\Gamma\pars{\bracks{\beta + 1}/2}}}\,;\qquad
\Re\pars{\alpha} > 0\,,\quad\Re\pars{\beta} > 0
\end{align}


When $\ds{\ul{\beta \to \alpha}}$, the solution becomes:
  $$
\half\bracks{\Psi\pars{\alpha \over 2} - \Psi\pars{\alpha + 1 \over 2}}
$$

