Could it possibly have a nice closed form? $\int _0^1\int _0^1\frac{x y}{(x+1) (y+1) \log (x y)}\ dx \ dy$ Using multiple integrals it's not hard to show that the present integral reduces to some integral
over squared digamma functions, but then things become harder. How would you tackle the problem?
$$\int _0^1\int _0^1\frac{x y}{(x+1) (y+1) \log (x y)}\ dx \ dy$$
 A: Just a note. We can decompose the integrand as the following.
$$
\frac{x y}{(x+1) (y+1) \log (x y)} = \frac{x}{(x+1)\log (x y)} - \frac{x}{(x+1) (y+1) \log (x y)}
$$
Then we notice that
$$
\int _0^1\int _0^1 \frac{x}{(x+1)\log (x y)} \ dx \ dy = \int_0^1 \frac{\operatorname{li}(x)}{x+1} \ dx
$$
where $\operatorname{li}$ is the logarithmic integral. This integral doesn't have closed-form as I know.
The rest of the work is to find a closed-form of the following integral:
$$\int _0^1\int _0^1\frac{x}{(x+1) (y+1) \log (x y)}\ dx \ dy.$$
A: Another hint. The integral is easy to separate into several simpler integrals:
$$\int _0^1\int _0^1\frac{x y}{(x+1) (y+1) \log (x y)}\ dx \ dy=$$
$$=\int _0^1\int _0^1\frac{dx \ dy}{ \log (x y)}-\int _0^1\int _0^1\frac{dx \ dy}{ (1+y)\log (x y)}-\int _0^1\int _0^1\frac{dx \ dy}{ (1+x)\log (x y)}+$$
$$+\int _0^1\int _0^1\frac{dx dy}{(x+1) (y+1) \log (x y)}$$
The first integral:
$$\int _0^1\int _0^1\frac{dx \ dy}{ \log (x y)}=\int _0^1\frac{\text{li} (x) dx}{x}=-1$$
The second and third are the same integral:
$$-2\int _0^1\int _0^1\frac{dx \ dy}{ (1+x)\log (x y)}=-2\int _0^1\frac{\text{li} (x) dx}{x(1+x)}=2+2\int _0^1\frac{\text{li} (x) dx}{1+x}$$

$$\int _0^1\int _0^1\frac{x y}{(x+1) (y+1) \log (x y)}\ dx \ dy=$$
$$=1+2\int _0^1\frac{\text{li} (x) dx}{1+x}+\int _0^1\int _0^1\frac{dx dy}{(x+1) (y+1) \log (x y)}$$

If we consider a more complicated integral:
$$I(a)=\int _0^1\int _0^1\frac{e^{axy}~dx ~dy}{(x+1) (y+1) \log (x y)}$$
We can turn the above into a linear ODE:
$$\frac{dI}{da}=1+2\int _0^1\frac{\text{li} (x) dx}{1+x}+I(a)$$
However, this is not very useful, because we can't find the integration constant without knowing some particular value of $I(a)$.
$$I(a)=C_1 e^a-1-2\int _0^1\frac{\text{li} (x) dx}{1+x}$$
$$C_1=I(0)+1+2\int _0^1\frac{\text{li} (x) dx}{1+x}$$
On a happier note, if we find your original integral, we immediately find $I(a)$ as well.
