Evaluating ∫10⋯∫10{1×1⋯xn}2dx1⋯dxn\int_{0}^{1}\cdots\int_{0}^{1}\bigl\{\frac{1}{x_{1}\cdots x_{n}}\bigr\}^{2}\:\mathrm{d}x_{1}\cdots\mathrm{d}x_{n}

For any $\alpha > 0$ and $n \in \mathbb{Z}_{+}$. Let $I_n(\alpha)$ be the $n$-dimensional integral

$$I_n(\alpha) = \int_{[0,1]^n} \left\{1\left/\prod_{i=1}^n x_i\right.\right\}^\alpha \prod_{i=1}^n dx_i$$

It is clear all these $I_n(\alpha) \in (0,1)$. As a result, the corresponding generating function

$$I(\alpha,z) = \sum_{k=0}^\infty I_{k+1}(\alpha)z^k$$
converges absolutely for any $|z| < 1$.

We will first assume $z$ is small and positive.

For any given $n$, substitute $x_i$ by $e^{-t_i}$ for $i = 1,\ldots,n$. Let $t = \sum\limits_{i=1}^n t_i$ and $x = e^t$. We have

$$I_n(\alpha) = \frac{1}{(n-1)!}\int_0^\infty \left\{ e^t \right\}^\alpha e^{-t} t^{n-1} dt$$

and hence
$$I(\alpha,z) = \int_0^\infty \left\{ e^t \right\}^\alpha e^{(z-1)t} dt = \int_1^\infty \left\{ x \right\}^\alpha x^{z-2} dx = \sum_{k=1}^\infty \int_0^1 \frac{x^\alpha}{(x+k)^{2-z}} dx$$
In the last expression, if we replace the denominator $(x + k)^{2-z}$ by its integral representation

$$\frac{1}{(x+k)^{2-z}} = \frac{1}{\Gamma(2-z)}\int_0^\infty s^{1-z} e^{-(x+k)s} ds$$

We get

$$\begin{align} I(\alpha,z) &= \frac{1}{\Gamma(2-z)} \sum_{k=1}^\infty \int_0^1 x^\alpha \left[ \int_0^\infty s^{1-z} e^{-(x+k)s} ds\right] dx\\ &= \frac{1}{\Gamma(2-z)} \int_0^1 x^\alpha \left[ \int_0^\infty s^{1-z} \frac{e^{-xs}}{e^s - 1} ds \right] dx\\ &= \frac{1}{\Gamma(2-z)} \int_0^\infty \frac{s^{1-z}}{e^s - 1}\left[\int_0^1 x^\alpha e^{-xs} dx \right] ds\\ &= \frac{1}{\Gamma(2-z)} \int_0^\infty \frac{s^{1-z}}{e^s - 1}\left[\frac{1}{s^{\alpha+1}}\int_0^s x^\alpha e^{-x} dx \right] ds\tag{*1} \end{align}$$
Up to this step, we are using the assumption $z$ is small and positive. A consequence of this assumption is all the terms involved in above steps are non-negative numbers.
The replacement of the denominator by its integral representation, the switching order of summation and integrations are automatically valid.

To proceed further, we need to split what's inside the square bracket in last expression of $(*1)$. If we do that, we will notice there are terms that are no longer non-negative. Furthermore, if $z$ remains to be small and positive, some of the terms simply diverge.

To bypass this obstacle, we will use the fact as long as $\Re z < 1$, the last expression
of $(*1)$ defines an analytic function in $z$. Instead of sticking with $z$ is small and positive,

we will switch our assumption to $z$ is real and sufficiently negative

After we work out what the last expression of $(*1)$ really are, we will analytic continue
the result back to small and positive $z$.

For integer $\alpha > 0$, we have

$$\begin{align} I(\alpha,z) &= \frac{1}{\Gamma(2-z)}\int_0^\infty \frac{s^{-(\alpha+z)}}{e^s - 1}\left[ -e^{-x} \sum_{k=0}^\alpha \frac{\alpha!}{(\alpha-k)!}x^{\alpha-k}\right]_0^s ds\\ &= \frac{1}{\Gamma(2-z)}\int_0^\infty \frac{s^{-(\alpha+z)}}{e^s - 1} \left[ \alpha! - e^{-s} \sum_{k=0}^\alpha \frac{\alpha!}{(\alpha-k)!}s^{\alpha-k} \right] ds\\ &= \frac{\alpha!}{\Gamma(2-z)}\int_0^\infty \left[s^{-(\alpha+z)} - \frac{1}{e^s-1}\sum_{k=0}^{\alpha-1}\frac{s^{-(k+z)}}{(\alpha-k)!} \right] e^{-s} ds \end{align}$$
Since $\displaystyle\; \int_0^\infty \frac{s^{\beta-1}}{e^s-1} e^{-s} ds = \Gamma(\beta)( \zeta(\beta) - 1) \;$, we find
$$I(\alpha,z) = \frac{\alpha!}{\Gamma(2-z)} \left[\Gamma(1-\alpha-z) - \sum_{k=0}^{\alpha-1}\Gamma(1-k-z)\frac{\zeta(1-k-z)-1}{(\alpha-k)!} \right]$$

In particular, when $\alpha = 2$, this leads to
$$I(2,z) = \frac{2!}{\Gamma(2-z)} \left[ \Gamma(-1-z) -\frac12\Gamma(1-z)(\zeta(1-z)-1) -\Gamma(-z)(\zeta(-z)-1) \right]$$
Since $\;\Gamma(1+w) = w\Gamma(w)\;$, this can be simplified as
$$\begin{align} I(2,z) &= \frac{2}{z(1-z^2)} - \frac{1}{1-z}(\zeta(1-z)-1) + \frac{2}{z(1-z)}(\zeta(-z)-1)\\ &= -\frac{1}{z+1} - \frac{1}{1-z}\zeta(1-z) + \frac{2}{z(1-z)}\zeta(-z)\tag{*2} \end{align}$$
We know $\zeta(w)$ is analytic over the whole complex plane except at $w = 1$.
It has a simple pole with residue $1$ there. Using this info, we know RHS$(*2)$ may have poles at $z = 0$ and $\pm 1$.

• We don't care what happens at $z = 1$.
• Near $z = -1$, RHS$(*2)$ $\sim -\frac{1}{z+1} + O(1) + (-1 + O(z+1))(\frac{1}{-z-1} + O(1)) = O(1)$.
• Near $z = 0$, RHS$(*2)$ $\sim O(1) - (1 + O(z))(\frac{1}{-z} + O(1)) + \frac{1}{z}(2\zeta(0) + O(z)) = O(1)$ again.

Combine these, we know RHS$(*2)$ is analytic at $z = -1$ and $0$. There is nothing
from extending its validity to the unit disk $|z| < 1$ where $I(\alpha,z)$ coincides
with its definition as a generating function. As a result,

$$\bbox[4pt,border:1px solid blue;]{ I_n(2) = (n-1)!\frac{d^{n-1}}{dz^{n-1}}\left[-\frac{1}{z+1} - \frac{1}{1-z}\zeta(1-z) + \frac{2}{z(1-z)}\zeta(-z) \right]_{z=0} }\tag{*3}$$

If we throw RHS$(*2)$ to WA and ask for its Taylor expansion, WA returns

$$I(2,z) = (-1-\gamma + \log(2\pi)) + z \left(1 - \gamma+ \frac{\gamma^2}{2}-\frac{\pi^2}{24} + \log(2\pi) - \frac{\log^2(2\pi)}{2}\right) + O(z^2)$$

This implies
$$\begin{align} I_1(2) &= -1-\gamma + \log(2\pi)\\ I_2(2) &= 1 - \gamma+ \frac{\gamma^2}{2}-\frac{\pi^2}{24} + \log(2\pi) - \frac{\log^2(2\pi)}{2} \end{align}$$
as expected.