Compute integral of general form $ \int_0^\infty \left(\frac{x}{\sinh x}\right)^n d x $

Consider the integral of the general form

$$ J_{n,m} = \int_0^\infty x^n \text{csch}^{m}x \>\mathrm{d}x $$

and note that

$$ I_n = \int_0^\infty \Big(\frac{x}{\sinh x}\Big)^n d x= J_{n,n} $$

As shown below, $ J_{n,n}$ can be computed via the reduction formula

$$J_{n,m}=\frac{n(n-1)}{(m-1)(m-2)}J_{n-2,m-2} - \frac{m-2}{m-1} J_{n,m-2}\tag1 $$ along with the starting values $$J_{2k-1,1}= \frac{(2^{2k}-1)(2k-1)!}{2^{2k-1}}\zeta(2k),\>\>\>\>\> J_{2k,2}=\frac{k(2k-1)!}{2^{2k-2}}\zeta(2k)\tag2 $$ $$\zeta(2)=\frac{\pi^2}6, \>\>\>\zeta(4)= \frac{\pi^4}{90},\>\>\> \zeta(6)=\frac{\pi^6}{945},\>\>\>\zeta(8)= \frac{\pi^8}{9450}, \>\>\>\cdots $$

Proof: Apply $$\int \text{csch}^{m}x\>\mathrm{d}x = -\frac{\text{csch}^{m-2}x \>{\coth}x }{m-1} -\frac{m-2}{m-1} \int \text{csch}^{m-2}x\>\mathrm{d}x $$

to integrate $J_{n,m}$ by parts $$J_{n,m} = \frac{n}{m-1} \int x^{n-1} \text{csch}^{m-2}x \coth x \mathrm{d}x- \frac{m-2}{m-1}J_{n,m-2}\tag3 $$ where $$ \int x^{n-1} \text{csch}^{m-2}x \coth x \mathrm{d}x =-\int x^{n-1} \text{csch}^{m-3}x \>d(\text{csch}x)\overset{IBP} = \frac{n-1}{m-2} I_{n,m-2} $$

Plug into (3) to obtain the reduction formula (1). To evaluate $J_{2k-1,1}$ and $ J_{2k,2}$ given in (2), apply the substitution $x=-\ln t$. Then $$J_{2k,2}= 4\int_0^1 \frac{t \ln^{2k}t}{(1-t^2)^2}dt\overset{t^2\to t}=\frac1{2^{2k-1} }\int_0^1 \frac{\ln^{2k}t}{(1-t)^2}dt\\ \overset{\text{IBP}}=-\frac k{2^{2k-2} }\int_0^1 \frac{\ln^{2k-1}t}{1-t}dt =\frac{k(2k-1)!}{2^{2k-2}}\zeta(2k) $$

$$J_{2k-1,1}= -2\int_0^1 \frac{\ln^{2k-1}t}{1-t^2}dt =\frac{(2^{2k}-1)(2k-1)!}{2^{2k-1}}\zeta(2k) $$

Displayed below are some results via the reduction formulae (1) and (2):

\begin{align} I_1 =J_{1,1} &= \frac32 \zeta(2)= \frac{\pi^2}4\\ I_2 =J_{2,2} &= \zeta(2)= \frac{\pi^2}6\\ I_3 =J_{3,3} &= 3J_{1,1}-\frac12 J_{3,1} = 3\cdot \frac{\pi^2}4 - \frac12 \cdot \left(\frac{45}4 \zeta(4)\right) = \frac{3\pi^2}4 - \frac{\pi^4}{16}\\ I_4 =J_{4,4} & = 2J_{2,2}-\frac23 J_{4,2} = 2\cdot \frac{\pi^2}6 - \frac23 \cdot \left(3\zeta(4)\right) = \frac{\pi^2}3 - \frac{\pi^4}{45}\\ I_5 =J_{5,5} & = \frac53 J_{3,3}-\frac34\left( 10J_{3,1} -\frac12J_{5,1} \right) = \frac{5\pi^2}4 - \frac{25\pi^4}{24} +\frac{3\pi^6}{32}\\ I_6 =J_{6,6} &= \frac32 J_{4,4}-\frac45 \left( 5J_{4,2} -\frac23 J_{6,2} \right) = \frac{\pi^2}2 - \frac{\pi^4}{6} +\frac{4\pi^6}{315} \\ I_7 =J_{7,7} &=\>\cdots \end{align}

Claim: for each positive integer $n$, $I_n$ is a linear combination of $\zeta(2),\zeta(4) \dots \zeta(2\left \lceil{n/2}\right \rceil )$ over $\mathbb{Q}$. The coefficients can be computed explicitly. More precisely,

$$I_n=\left\{\begin{array}{ll}m\sum_{j=0}^{m-1}A_{j}^n\zeta(2j+2), & \textrm{if } n=2m,\\ 2(2m+1)\sum_{j=0}^{m}A_j^n(1-1/2^{2+2j})\zeta(2j+2), & \textrm{if } n=2m+1, \end{array} \right.$$ where $A^n_j$ are integers determined by $$(x+2m-1)\cdots (x+m+1)\widehat{(x+m)}(x+m-1)\cdots (x+1)=\sum_{j=0}^{m-1}A^n_j(x+m)^{m-2j-2}$$ if $n=2m$; and $$2^{2m}(x+2m)(x+2m-1)\cdots (x+2)(x+1)=\sum_{j=0}^{m}A_j^n (2x+2m+1)^{2m-2j},$$ if $n=2m+1$.

Examples: $$\begin{array}{ll}I_6&=3(\zeta(2)-5\zeta(4)+4\zeta(6)), \\ &\textrm{ since }(x+1)(x+2)(x+4)(x+5)=(x+3)^4-5(x+3)^2+4,\\ I_7&=14\left((1-1/2^2)\zeta(2)-35(1-1/2^4)\zeta(4)+259(1-1/2^6)\zeta(6)-225(1-1/2^8)\zeta(8)\right), \\ & \textrm{ since }2^6(x+1)(x+2)(x+3)(x+4)(x+5)(x+6)\\ &\qquad \qquad =(2x+7)^6-35(2x+7)^4+259(2x+7)^2-225,\\ I_8&=4(\zeta(2)-14\zeta(4)+49\zeta(6)-36\zeta(8)),\\ & \textrm{ since } (x+1)(x+2)(x+3)(x+5)(x+6)(x+7)\\ &\qquad \qquad =(x+4)^6-14(x+4)^4+49(x+4)^2-36 \end{array}$$

Proof of the claim: This should be a consequence of Quanto's nice answer. But here we start from the formula in J.G's comment $$ I_n=n!2^n\sum_{k\ge 0}\frac{1}{(2k+n)^{n+1}}\binom{n+k-1}{k}.$$ The following is just an exercise after J.G's comment.

First assume that $n=2m$ is even. Then we can write $$(2m-1)!\binom{x+2m-1}{2m-1}=\sum_{j=0}^{m-1}A^n_j (x+m)^{2m-2j-1}.$$ Thus $$\begin{array}{ll}I_n&=m\sum_{k\ge 0}\frac{1}{(k+m)^{2m+1}}\sum_{j=0}^{m-1}A_j^n(k+m)^{2m-2j-1}\\ &=m\sum_{k\ge 0}\sum_{j=0}^{m-1}A_j^n\frac{1}{(k+m)^{2j+2}}\\ &=m\sum_{k\ge 0}\sum_j A_j^n \frac{1}{k^{2j+2}}\\ &=m\sum_{j=0}^{m-1} A_j^n \zeta(2j+2). \end{array} $$ Note that $\sum_{j=0}^{m-1}A_j^n(x+m)^{2m-2j-2}=0$ for $x=-1,\dots,-(m-1)$, since $-1,\dots,-(m-1)$ are roots of $$(2m-1)!\binom{x+2m-1}{2m-1}/(x+m)=(x+2m-1)\dots \widehat{(x+m)}(x+(m-1))\dots (x+1).$$ The above equation is the same as $$\sum_{j=0}^{m-1}A^n_j l^{2m-2j-2}=0, l=1,2,\dots,m-1.$$ Thus we get $$\sum_{j=0}^{m-1}A_j^n\frac{1}{l^{2j+2}}=0, l=1,2,\dots, m-1, $$ which implies $$m\sum_{k\ge 0}\sum_{j=0}^{m-1}A_j^n\frac{1}{(k+m)^{2j+2}}\\ =m\sum_{k\ge 0}\sum_j A_j^n \frac{1}{k^{2j+2}},$$ one step we used in the above calculation.

Next, consider the case when $n=2m+1$ is odd. There are integers $A^n_0=1,A^n_1,A^n_2,\dots, A^n_m$ such that $$(*)\qquad 2^{2m}(2m)!\binom{x+2m}{2m}=(2x+2m+1)^{2m}+A^m_1(2x+2m+1)^{2m-2}+\dots+A^m_m\\ =\sum_{j=0}^m A^n_j (2x+2m+1)^{2m-2j}. $$ holds as a polynomial of $x$. This is because the polynomial $(x+2m)\cdots (x+1)$ is symmetric about the line $x=m+1/2$.

Thus $$\begin{array}{ll} I_{2m+1}&=2(2m+1)\sum_{k\ge 0}\frac{1}{(2k+2m+1)^{2m+2}}\left(2^{2m}(2m)!\binom{k+2m}{2m} \right)\\ &=2(2m+1)\sum_{k\ge 0}\left(\sum_{j=0}^m \frac{A^n_j}{(2k+2m+1)^{2+2j}} \right)\\ &=2(2m+1)\sum_{k\ge 0}\left( \sum_{j=0}^mA^n_j \frac{1}{(2k+1)^{2+2j}}\right)\\ &=2(2m+1)\sum_{j=0}^m\left(A^n_j(1-\frac{1}{2^{2+2j}})\zeta(2+2j)). \right) \end{array}$$ where we used $$\sum_{j=1}^m A_j^n(2k+2m+1)^{2m-2j}=0, k=-1,-2,\dots,-m.$$

This is a remark continuing from J.G.'s hint and Q. Zhang's solution.

For even $n$ the coefficients $A_j^n$ of $\zeta(2k+2)$ relate to the "central factorial numbers" (OEIS A008955):

$$ 2\begin{pmatrix} I_2/2 \\ I_4/4 \\ I_6/6 \\ I_8/8 \\ I_{10}/10 \\ \vdots \end{pmatrix} = \begin{pmatrix} 1& \hfill0& \hfill0& \hfill0& \hfill0 & \cdots \\ 1& \hfill-1& \hfill 0& \hfill 0& \hfill 0 \\ 1& \hfill-5& \hfill 4& \hfill 0& \hfill 0 \\ 1& \hfill-14& \hfill49& \hfill-36& \hfill 0 \\ 1& -30& 273& -820& 576 \\ \vdots &&&&&\ddots \end{pmatrix} \begin{pmatrix} \zeta(2) \\ \zeta(4) \\ \zeta(6) \\ \zeta(8) \\ \zeta(10) \\ \vdots \end{pmatrix} $$

or generally $2I_{2m}/{(2m)} = \sum_{k=0}^{m-1} t(2m, 2m-2k)\zeta(2k+2)$ with $t(2n,2k)$ the even central factorial numbers (of the first kind), see OEIS A008955.

For odd $n$ such a relation holds as well, namely $$ \frac{1}{2} \begin{pmatrix} I_1 /1 \\ I_3 /3 \\ I_5 /5 \\ I_7 /7 \\ I_9 /9 \\ \vdots \end{pmatrix} = \begin{pmatrix} 1& \hfill 0& \hfill 0& \hfill 0& \hfill 0 & \cdots \\ 1& \hfill -1& \hfill 0& \hfill 0& \hfill 0 \\ 1& \hfill-10& \hfill 9& \hfill 0& \hfill 0 \\ 1& \hfill-35& \hfill 259& \hfill -225& \hfill 0 \\ 1& \hfill-84& 1974&-12916& 11025 \\ \vdots &&&&&\ddots \end{pmatrix} \begin{pmatrix} (1-2^{-2})\zeta(2) \\ (1-2^{-4})\zeta(4) \\ (1-2^{-6})\zeta(6) \\ (1-2^{-8})\zeta(8) \\ (1-2^{-10})\zeta(10) \\ \vdots \end{pmatrix}, $$

i.e., $I_{2m+1}/(2m+1) = 2\sum_{k=0}^{m-1} 4^k t(2m+1,2m+1-2k)(1 - \frac{1}{2^{2k+2}})\zeta(2k+2)$ with $t(2n+1,2k+1)$ the odd central factorial numbers, see OEIS A008956.

This is of interest to me because $\sum_{n=1}^\infty I_n/n$ should exist and I would be interested if there's a closed form expression for it, as it relates to another integral I'm interested in, see this question.

To prove the relations above, one way is to note that Q. Zhang's coefficients can be expressed via Newton's identities as elementary symmetric polynomials. If $e_j^n$ is the $j$th elementary symmetric polynomial of degree $n$ evaluated at $1^2,\ldots,n^2$, the recursion $e_j^n = e_j^{n-1} + n^2 e_{j-1}^{n-1}$ holds. Comparing to the recursion formulae for the central factorial numbers (i.e., in OEIS or in Butzer et al.), the claim follows for even $n$. The case for odd $n$ can be shown along similar lines.