A multiple integral that seems related to the $\zeta$ function at even integers

Question

I came across this integral that seems related to the Riemann zeta function $\zeta(2n)$ evaluated at even integers $2n \in 2\mathbb{Z}$. Letting $n$ be an even integer, define the multiple integral over $(2n+1)$ variables $u_1 \cdots u_{2n+1}$

\begin{equation} \mathcal{I}_{2n} = \int_0^1 du_1 \cdots \int_0^1 du_{2n+1} \frac{1}{1+u_1} \frac{1}{u_1+u_2} \frac{1}{u_2+u_3} \cdots \frac{1}{u_{2n} + u_{2n+1}} \frac{1}{u_{2n+1}+1}. \end{equation}

For example, the case $2n=4$ has integrand $\frac{1}{1+u_1} \frac{1}{u_1+u_2} \frac{1}{u_2+u_3} \frac{1}{u_{3} + u_{4}} \frac{1}{u_{4}+u_5} \frac{1}{u_{5}+1}$.

Below are listed some exact and numerical results for the integrals for the first few values of $2n$ that could be numerically evaluated on my laptop. The $2n=0$ case is easy, and the $2n=2$ case could evaluated exactly by Mathematica in terms of a complicated expression of polylogarithms (although below, I show an alternate way to explicitly derive this).

• $\mathcal{I}_0 = \frac{1}{2} = - \zeta(0)$
• $\mathcal{I}_2 = \frac{\pi^2}{6} = \zeta(2)$
• $\mathcal{I}_4 \approx 8 \frac{\pi^4}{90} = 8 \cdot \zeta(4)$
• $\mathcal{I}_6 \approx 54 \cdot \zeta(6)$
• $\mathcal{I}_8 \approx 384 \cdot \zeta(8)$
• $\mathcal{I}_{10} \approx 2880 \cdot \zeta(10)$

These integrals are almost exactly integer multiples of the zeta function at even integers, up to the small errors given by Mathematica. The series of integers $1,8,54,384,2880$ don't appear in the OEIS, although the first terms $1,8,54,384$ go as $k^2 \cdot k!$, and $2880$ isn't far off from $5^2 \cdot 5! = 3000$.

This integral seems closely related to the one considered in this question. However, their method for their integral seemed quite magical and I was not able to generalize it. While I'm specifically interested in this integral, it would be nice to know if there's a general method to deal with these integrals.

Below are alternate formulations of the integral that may help, as well as showing the result for $2n=2$. First, one can change variables to $u_i = \frac{1-v_i}{1+v_i}$ to rewrite it as \begin{equation} \mathcal{I}_{2n} = \frac{1}{2} \int_0^1 dv_1 \cdots \int_0^1 dv_{2n+1} \frac{1}{1 - v_1 v_2} \frac{1}{1 - v_2 v_3} \cdots \frac{1}{1 - v_{2n} v_{2n+1}}, \end{equation} which is similar to a form found in the other MathOverflow question linked above. As such, the same substitutions used there done backwards yield other expressions \begin{equation} \begin{split} \mathcal{I}_{2n} &= \frac{1}{2} \int_0^1 dy_{2n+2} \int_0^{y_{2n+2}} dy_{2n+1} \cdots \int_{0}^{y_3} dy_{2} \frac{1}{y_3} \frac{1}{y_4-y_2} \cdots \frac{1}{y_{2n+2}-y_{2n+1}} \frac{1}{1-y_{2n+1}} \\ &= \frac{1}{2} \int_0^1 d \tilde{u}_{1} \cdots \int_0^1 d \tilde{u}_{2n+2} \frac{\delta(1-\tilde{u}_1 - \cdots - \tilde{u}_{2n+2})}{(\tilde{u}_1 + \tilde{u}_2)(\tilde{u}_2 + \tilde{u}_3)\cdots(\tilde{u}_{2n+1} + \tilde{u}_{2n+2})} \end{split} \end{equation} which almost match the integral considered in linked question with $2n+2$ variables, except missing the last factor $\frac{1}{\tilde{u}_{2n+2} + \tilde{u}_{1}}$

For $2n=2$, one can use Feynman Parameterization to write \begin{equation} \begin{split} \mathcal{I}_{2} &= \frac{1}{2} \int_0^1 dv_1 \int_0^1 dv_2 \int_0^1 dv_3 \frac{1}{1 - v_1 v_2} \frac{1}{1 - v_2 v_3} \\ &= \frac{1}{2} \int_0^1 du \int_0^1 dv_1 \int_0^1 dv_2 \int_0^1 dv_3 \frac{1}{\left[(1 - v_1 v_2)u + (1 - v_2 v_3)(1-u)\right]^2} \\ &= \frac{1}{2} \int_0^1 du \left( \frac{\log(1-u)}{u} + \frac{\log(u)}{1-u} \right) = \int_0^1 du \frac{\log(1-u)}{u} = \sum_{k=1}^{\infty} \frac{1}{k^2} = \frac{\pi^2}{6}, \end{split} \end{equation} where the second-to-last equality follows from term-by-term integration of the series expansion.

I haven't found a way to generalize the method for any of the other cases and still don't have analytical proof for the next equalities.

---

EDIT: See below for a beautiful answer! It turns out that it's not related to the $\zeta$ function except for the first few values of $n$.

Answer 1

We have $$I_{2n}=\frac{(2n)!!}{(2n+1)!!}\cdot \frac1{2n+2}\cdot \pi^{2n}.$$ To see this, we follow the suggestion by Terry Tao in the comments and apply the diagonalization of the integral operator with the kernel $1/(x+y)$ on $[0,1]$. Change the variable to $1/x\in [1,\infty)$ and use (1.18) here (this is Mehler integral operator, as I understand) to diagonalize. Substituting the value of Legendre functions at 1, we get $$I_{2n}=\pi^{2n}\int_0^\infty x\tanh x/\cosh^{2n+2} x dx.$$ (The above part is suggested by Vladimir Petrov, I understand nothing about all this special functions and integral transforms stuff. But the answer for small $n$ coincides, so I guess everything is ok:) Well, you are free to ask for more details if necessary.)

For evaluating these integrals, integrate by parts noting that $$\tanh x/\cosh^{2n+2} x dx=\tanh x\cosh^{-2n} xd\tanh x=\frac12(1-\tanh^2 x)^nd\tanh^2 x\\=\frac{-1}{2(n+1)}(1-\tanh^2x)^{n+1},$$ thus $$I_{2n}=\pi^{2n}\cdot \frac1{2(n+1)}\int_0^\infty (1-\tanh^2x)^{n+1}dx=\pi^{2n}\cdot \frac1{2(n+1)}\int_0^1 (1-t^2)^{n}dt\\=\frac{\pi^{2n}}{2(n+1)}\cdot \frac {(2n)!!}{(2n+1)!!}$$ (here $t=\tanh^2 x$, the last integral is well known and may be calculated by induction or using Beta function or how do you prefer.)