Evaluate $\int_0^1 \frac{\operatorname{arctanh}^3(x)}{x}dx$

As MrTaurho pointed out in the comments, we can rewrite $\,\displaystyle{\operatorname{arctanh}x=\frac12 \ln\left(\frac{1+x}{1-x}\right)}$, this gives: $$I=\int_{0}^1 \frac{\operatorname{arctanh}(x)^3}{x}dx=-\frac18\int_0^1 \frac{\ln^3\left(\frac{1-x}{1+x}\right)}{x}dx$$ And we will substitute $\displaystyle{\frac{1-x}{1+x}=t\Rightarrow dx=-\frac{2}{(1+t)^2}dt}$ in order to get: $$I=-\frac18\int_0^1 \frac{\ln^3 t}{\frac{1-t}{1+t}}\frac{2}{(1+t)^2}dt=-\frac14 \int_0^1 \frac{\ln^3 t}{1-t^2}dt=-\frac14 \sum_{n=0}^\infty \int_0^1 t^{2n}\ln^3t dt$$ Now we consider the following integral: $$\int_0^1 t^a dt=\frac{1}{a+1}\Rightarrow \int_0^1 t^a \ln^3 tdt=\frac{d^3}{da^3} \left(\frac{1}{a+1}\right)=-\frac{6}{(a+1)^4}$$ $$\Rightarrow I=\frac{6}{4}\sum_{n=0}^\infty \frac{1}{(2n+1)^4}=\frac32\cdot\frac{\pi^4}{96}=\frac{\pi^4}{64}$$

---

Of course this can be generalized for any power, so let's do that. Consider:

$$I(k)=\int_0^1 \frac{\text{arctanh}^kx}{x}dx=\frac{(-1)^k}{2^k}\int_0^1 \frac{\ln^k\left(\frac{1-x}{1+x}\right)}{x}dx\overset{\large\frac{1-x}{1+x}=t}=\frac{2(-1)^k}{2^k}\int_0^1 \frac{\ln^k t}{1-t^2}dt$$ $$=\frac{(-1)^k}{2^{k-1}}\sum_{n=0}^\infty \int_0^1 x^{2n}\ln^k xdx=\frac{(-1)^k}{2^{k-1}}\sum_{n=0}^\infty \frac{(-1)^k k!}{(2n+1)^{k+1}}=\frac{k!}{2^{k-1}} \sum_{n=0}^\infty \frac{1}{(2n+1)^{k+1}}$$ Where above we used the following result: $$\int_0^1 t^a dt=\frac{1}{a+1}\Rightarrow \int_0^1 t^a \ln^k tdt=\frac{d^k}{da^k} \left(\frac{1}{a+1}\right)=\frac{(-1)^k k!}{(a+1)^{k+1}}$$ And finally it reduces to: $$I(k)=\frac{k!}{2^{k-1}} \left(1-\frac{1}{2^{k+1}}\right)\zeta(k+1)=\boxed{\frac{k!\left(2^{k+1}-1\right)}{4^k}\zeta(k+1)}$$

---

One can verify the result by comparing to the one announced by Maple: $$I(5)=\int_0^1 \frac{\text{arctanh}^5 x}{x}dx=\frac{5!(2^6-1)}{4^5}\zeta(6)=\frac{945}{128}\cdot\frac{\pi^6}{945}=\frac{\pi^6}{128}$$

An alternative route to the general result: For $z \in \mathbb{C}$ with $\operatorname{Re}(z) > 0$ we can let $x = \tanh(t/2)$ and use the geometric series to find \begin{align} \int \limits_0^1 \frac{\operatorname{artanh}^z (x)}{x} \, \mathrm{d} x &= \frac{1}{2^z} \int \limits_0^\infty \frac{t^z}{\sinh(t)} \, \mathrm{d} t = 2^{1-z} \sum \limits_{k=0}^\infty \int \limits_0^\infty t^z \mathrm{e}^{-(2k+1) t} \, \mathrm{d} t \\ &= 2^{1-z} \int \limits_0^\infty u^z \mathrm{e}^{-u} \, \mathrm{d} u \sum \limits_{k=0}^\infty \frac{1}{(2k+1)^{1+z}} = 2^{1-z} \Gamma(1+z) \lambda(1+z) \\ &= \frac{2^{1+z}-1}{4^z} \Gamma(1+z) \zeta(1+z) \, , \end{align} where $\lambda$ is the Dirichlet lambda function.