For $n\ge m\ge 2$, define $$I(n,m):= \int _0^\infty\dfrac{\text{arcsinh}^nx}{x^m}dx$$ Computer algebra systems say that the indefinite integral can be expressed in terms of polylog functions (of rapidly increasing complexity), but I see no way of doing the limit $\to\infty$. On the other hand, I have found numerically that $I(n,m)$ can (conjecturally always) be written as a rational combination of zeta values, moreover **this happens in a rather beautiful way**.

Here are the values for small $m$
and then what seems to be the general formula.

Defining $\bar\zeta(n):=(2^n-1) \zeta(n)$, it seems like

$ I(n,2)=\dfrac{n!}{2^{n-1}} \bar\zeta(n)$

$ I(n,4)=\dfrac{n!}{2^{n-1}} \cdot\dfrac1 {6}[ 4\bar\zeta(n-2)-\bar\zeta(n)]$

$ I(n,6)=\dfrac{n!}{2^{n-1}} \cdot\dfrac1{120}[ 16\bar\zeta(n-4)-40\bar\zeta(n-2)+9\bar\zeta(n)]$

$ I(n,8)=\dfrac{n!}{2^{n-1}} \cdot\dfrac1{5040}[ 64\bar\zeta(n-6)-560\bar\zeta(n-4) +1036\bar\zeta(n-2)-225\bar\zeta(n)]$

*(you may want to stop here for a moment before reading on and look if the coefficients tell you something already)*

and

$ I(n,3)=\dfrac{n!}{2^{n-1}} \zeta(n-1)$

$ I(n,5)=\dfrac{n!}{2^{n-1}} \cdot\dfrac13[ \zeta(n-3)-\zeta(n-1)]$

$ I(n,7)=\dfrac{n!}{2^{n-1}} \cdot\dfrac2{45}[\zeta(n-5)- 5 \zeta(n-3)+4\zeta(n-1)]$

Well, in umbral notation introducing a formal variable $\bar Z$, for even $m=2k$ $$ I(n,2k)=\frac{n!}{2^{n-1}} \cdot\frac1{(2k-1)!}\bar Z^{n}\prod_{j=-k+1}^{k-2} \left(\frac2{\bar Z}-(2j+1)\right)$$ where each power $\bar Z^r$ is to be replaced by $\bar\zeta (r)$. E.g. for $m=6$, the "zeta polynomial" is related to $16z^4-40z^2+9=(2z-3)(2z-1)(2z+1)(2z+3)$.

Likewise for odd $m=2k+3$ (*sic* to make notation more elegant), $$ I(n,2k+3)=\frac{n!}{2^{n-1}} \cdot\frac{2^{2k}}{(2k+1)!}Z^{n}\prod_{j=-k}^{k}\left (\frac1{Z}-j\right)$$ where each power $Z^r$ is to be replaced by $\zeta (r)$. E.g. for $m=7$, the "zeta polynomial" is related to $z^5-5z^3+4z=(z-2)(z-1)z(z+1)(z+2)$.

We can rewrite the expression for even $m$ as follows to look amazingly similar to the one for odd $m$:

$$ I(n,2k+2)=\frac{n!}{2^{n-1}} \cdot\frac{2^{2k}}{(2k+1)!}\bar Z^{n}\prod_{j=-k+1/2}^{k-1/2} \left(\frac1{\bar Z}-j\right)$$ where in the product $j$ runs over the half-integers.

So far, this is all speculation, but once these patterns found, it seems obvious that there must be some deeper connection (which might even yield an elegant proof). Any insights? Has anybody encountered similar families of integrals?

**EDIT for the record:**
The integrals $$J(n,m):= \int _0^1\dfrac{\text{arctanh}^nx}{x^m}dx= \int _0^1\log^{n/2}\left(\dfrac{1+x }{1-x}\right)\dfrac{dx}{x^m}$$ and $$K(n,m):= \int _1^\infty\dfrac{\text{arcoth}^nx}{x^m}dx =\int _1^\infty\log^{n/2}\left(\dfrac{x+1}{x-1}\right)\dfrac{dx}{x^m}$$ have very similar situations with the exact same zeta values involved as for $I(n,m)$:

We have $$J(n,m)=\frac{n!}{2^{n-m}m!} \cdot Z^{n}p_ {m}(\frac1{Z} ),\qquad K(n,m)=\frac{n!}{2^{n-m}m!} \cdot\bar Z^{n}p_ {m}(\frac1{\bar Z} )$$ with even/odd polynomials $p_m$, the same for both integrals, of degree $m−2$. Again, powers $Z^r$ are to be replaced by $\zeta(r)$, but powers $\bar Z^r$ are to be replaced this time by $\bar\zeta(r):=\dfrac{2^{r-1}-1}{2^{r-1}}\zeta(r)$.

Here the polynomials $p_m\equiv p_m(z)$ are (essentially) not reducible, and instead of a recursion for $J(n,m)$, there is at least one for the polynomials: $$p_1=0,\ p_2=1, \ \frac{p_m}m=z\frac {p_{m-1}}{m-1} +\frac {m-3}4p_{m-2}.$$ The first of them are

$p_3=\frac32z$

$p_4=2z^2+1$

$p_5=\frac52( z^3+2z)$

$p_6=\frac32(2z^4+10z^2+3)$

$p_7=\frac74(2z^5+20z^3+23z)$

$p_8=4z^6 +70z^4+196z^2+45$

$p_9=\frac92(z^7+28z^5+154z^3+132z)$

Note that $p_m(1)=\dfrac{m!}{2^{m-1}}$ for $m\ge2$.