Solve definite double integral

Any hint to get this double integral solved would be helpful:

$$\int_0^1 \int_0^1 \frac{\left(\frac{1}{\ln(x)}+\frac{1}{1-x}-\frac{1}{2}\right)\left(\frac{1}{\ln(y)}+\frac{1}{1-y}-\frac{1}{2}\right)}{1 - xy} \,dx\,dy$$

Thanks!

• I suspect it has no known exact answer. The Wolfram Dev Platform can't evaluate it exactly. W|A says it is approximately $0.00651229.$ – Adrian Keister Sep 8 '18 at 0:00
• I tried with different CAS. Same conclusion : no exact answer. Using an inverse symbolic calculator, the closest is $$\frac{-\psi ^{(0)}\left(\frac{2}{21}\right)+\psi ^{(0)}\left(\frac{3}{16}\right)-\psi ^{(0)}\left(\frac{7}{10}\right)}{1000} \approx 0.0065122929313$$ while numerical integration leads to $0.0065122929371$ – Claude Leibovici Sep 8 '18 at 3:16

Expanding the denominator into a geometric series we find \begin{align} I &\equiv \int \limits_0^1 \int \limits_0^1 \frac{\left(\frac{1}{\ln(x)}+\frac{1}{1-x}-\frac{1}{2}\right)\left(\frac{1}{\ln(y)}+\frac{1}{1-y}-\frac{1}{2}\right)}{1 - xy} \, \mathrm{d} x \, \mathrm{d} y \\ &= \sum \limits_{n=0}^\infty \left[\int \limits_0^1 \left(\frac{1}{\ln(x)}+\frac{1}{1-x}-\frac{1}{2}\right) x^n \, \mathrm{d} x \right]^2 = \sum \limits_{k=1}^\infty \left[\int \limits_0^\infty \left(-\frac{1}{t}+\frac{1}{1-\mathrm{e}^{-t}}-\frac{1}{2}\right) \mathrm{e}^{-k t} \, \mathrm{d} t \right]^2 \\ &= \sum \limits_{k=1}^\infty \left[\int \limits_0^\infty \left(\frac{1}{\mathrm{e}^{t}-1} - \frac{1}{t}+\frac{1}{2}\right) \mathrm{e}^{-k t} \, \mathrm{d} t \right]^2 = \sum \limits_{k=1}^\infty \left[- \frac{\mathrm{d}}{\mathrm{d}z} \int \limits_0^\infty \left(\frac{1}{\mathrm{e}^{t}-1} - \frac{1}{t}+\frac{1}{2}\right) \frac{\mathrm{e}^{-z t}}{t} \, \mathrm{d} t ~ \Bigg \rvert_{z=k}\right]^2 \, . \end{align} Binet's first log-gamma formula yields \begin{align} I &= \sum \limits_{k=1}^\infty \left[\frac{\mathrm{d}}{\mathrm{d}z} \left\{ \left(z - \frac{1}{2}\right) \ln(z) - z + \frac{1}{2} \ln(2\pi) - \ln[\Gamma(z)]\right\} ~ \Bigg \rvert_{z=k}\right]^2 \\ &= \sum \limits_{k=1}^\infty \left[\ln(k) - \frac{1}{2k} - \psi(k)\right]^2 = \sum \limits_{k=1}^\infty \left[H_k - \ln(k) - \gamma - \frac{1}{2k}\right]^2 \, , \end{align} where the digamma identity $\psi(k) = H_{k-1} - \gamma \, , \, k \in \mathbb{N} ,$ was used.
The evaluation of the remaining sum is based on a few identities involving the harmonic numbers, most importantly the asymptotic expansion $$H_k \sim \ln(k) + \gamma + \frac{1}{2k} + \mathcal{O}\left(\frac{1}{k^2}\right) \, .$$ We split the sum into three parts, $I = X + Y - Z$ . The sums $$Y = \sum \limits_{k=1}^\infty \frac{1}{4k^2} = \frac{\zeta(2)}{4} = \frac{\pi^2}{24}$$ and \begin{align} Z &= \sum \limits_{k=1}^\infty \frac{H_k - \ln(k) - \gamma}{k} = \lim_{N \to \infty} \left[\frac{H_N^{(2)}}{2} + \frac{\ln^2(N)}{2} - \sum \limits_{k=1}^N \frac{\ln(k)}{k} + \frac{H_N^2 - 2 \gamma H_N - \ln^2 (N)}{2} \right] \\ &= \frac{\zeta(2)}{2} - \gamma_1 - \frac{\gamma^2}{2} = \frac{\pi^2}{12} - \frac{\gamma^2}{2} - \gamma_1 \end{align} are not too hard to calculate. In order to evaluate $X = \sum_{k=1}^\infty [H_k - \ln(k) - \gamma]^2$ we use the results \begin{align} \sum \limits_{k=1}^N H_k^2 &= (N+1) H_N^2 -(2N+1) H_N + 2N \, , \\ \sum \limits_{k=1}^N \ln^2 (k) &= N \ln^2 (N) - 2 N \ln(N) + 2N + \frac{1}{2} \ln^2 (N) + \zeta''(0) + \mathcal{O} \left(\frac{\ln^2 (N)}{N}\right) \, , \\ \sum \limits_{k=1}^N H_k \ln(k) &= N \ln^2 (N) - N (2-\gamma) (\ln(N)-1) + \frac{3}{4} \ln^2 (N) + \frac{1}{2} \ln(N) + C + \mathcal{O}\left(\frac{\ln(N)}{N}\right) \end{align} from this, this and this question, respectively. The constant $$C = \kappa + \frac{\gamma^2}{2}-\frac{\pi^2}{24}+\frac{1}{2}\gamma \log(2\pi) - \frac{1}{2} \log^2{(2\pi)} + \frac{3}{2} \gamma_1$$ contains $\kappa \approx −0.07760$ , which unfortunately does not have a particularly nice expression. After a slightly tedious calculation of the limit of partial sums we arrive at $$X = \frac{\pi^2}{24} + \frac{1}{2} \ln^2(2 \pi) - 1- \gamma - \frac{3}{2} \gamma^2 - 2 \gamma_1 - 2 \kappa \, .$$
Putting everything together we end up with $$I = \frac{1}{2} \ln^2 (2 \pi) - 1 - \gamma - \gamma^2 - \gamma_1 - 2 \kappa \approx 0.006512 \, .$$ Therefore the existence of a closed-form expression depends on whether there is a nice result for $X$ (see this question) or $\kappa$ .