General form for $\int_{0}^1 \frac{\operatorname{artanh}^a(x)}{x^b}dx$ for $a\geq b$ My Question:

Is there a closed form for $$I(a,b)=\int_{0}^1 \frac{\operatorname{artanh}^a(x)}{x^b}dx$$ For $a,b\in\mathbb{Z}$ such that $a\geq b$

In a previous question of mine, there was a solution for 
$$\int_{0}^1 \frac{\operatorname{arctanh}^a(x)}{x}dx=\frac{a!\left(2^{a+1}-1\right)}{4^a}\zeta(a+1)$$
Another closed form can be obtained using similar methods outlined in that post to obtain,
$$\int_{0}^1 \frac{\operatorname{arctanh}^a(x)}{x^2}dx=\frac{a!}{2^{a-1}}\zeta(a) $$
$$\int_{0}^1 \frac{\operatorname{arctanh}^a(x)}{x^3}dx=\frac{a!}{2^{a-1}}\zeta(a-1)$$
These closed forms beg the question of whether or not the first integral has a decent closed form. 
When $b\geq4$, the solutions become more complicated. For $b=4$:
$$\int_{0}^1 \frac{\operatorname{arctanh}^a(x)}{x^4}dx=\frac{a!}{3\cdot 2^{a-1}}\left(\zeta(a)+2\zeta(a-2)\right) $$
So I'm wondering if there is some way of writing the integral in terms of a finite sum of Zeta functions.
 A: Following the same approach as the given link, we have $\displaystyle\operatorname{arctanh}(x)=\frac12\ln\left(\frac{1+x}{1-x}\right)$, and with the same substitution $\displaystyle t=\frac{1+x}{1-x}$ we get
$$I(a,b)=-\frac1{2^{a-1}}\int_0^1\frac{(1+t)^{b-2}}{(1-t)^b}\ln^a(t)~\mathrm dt$$
Using generalized binomial expansion theorem, we can expand this as
$$I(a,b)=-\frac1{2^{a-1}(b-1)!}\sum_{j=0}^{b-2}\binom{b-2}j\sum_{k=0}^\infty(k+1)\dots(k+b-1)\int_0^1t^{k+j}\ln^a(t)~\mathrm dt$$
The integral at the end can be viewed as the $a$th derivative of $\displaystyle\int_0^1t^u~\mathrm dt=\frac1{u+1}$ at $u=k+j$, which is given by $\displaystyle\frac{(-1)^aa!}{(k+j+1)^{a+1}}$, which leaves us with the series
$$I(a,b)=\frac{(-1)^{a+1}a!}{2^{a-1}(b-1)!}\sum_{j=0}^{b-2}\binom{b-2}j\sum_{k=0}^\infty\frac{(k+1)\dots(k+b-1)}{(k+j+1)^{a+1}}$$
which boils the problem down to computing the last series. See that we have
$$\sum_{k=0}^\infty\frac{(k+1)\dots(k+b-1)}{(k+j+1)^{a+1}}=\sum_{k=j+1}^\infty\frac{P_{b,j}(k)}{k^{a+1}}=\sum_{m=0}^{b-1}\alpha_{b,j,m}\zeta(a-m+1)-\sum_{k=1}^j\frac{P_{b,j}(k)}{k^{a+1}}$$
for a polynomial $\displaystyle P_{b,j}(k)=\sum_{m=0}^{b-1}\alpha_{b,j,m}k^m$. This leaves us with the solution given by
$$I(a,b)=\frac{(-1)^{a+1}a!}{2^{a-1}(b-1)!}\sum_{j=0}^{b-2}\binom{b-2}j\left[\sum_{m=0}^{b-1}\alpha_{b,j,m}\zeta(a-m+1)-\sum_{k=1}^j\frac{P_{b,j}(k)}{k^{a+1}}\right]$$
The rest is just multiplying out $(k-j)\dots(k+b-j-2)$ to extract the coefficients of $P_{b,j}$. These are Stirling numbers of the first kind.
