Detailed derivations on this formula:
There is an exercise problem (Ex. 3.2) in Olver's book, on page 292. (or you can find it on this website: Eq.2.10.7).
$$\sum_{j=1}^{n-1}j^a \sim \zeta(-a)+\frac{n^{a+1}}{a+1}\sum_{s=0}^\infty \binom{a+1}{s}\frac{B_s}{n^s} \tag{*}$$
I try to derive it. I begin with:
$$\sum_{j=n_0}^n f(j)=\int_{n_0}^n f(x) dx+\frac{f(n_0)+f(n)}{2}+\sum_{s=1}^{m-1}\frac{B_{2s}}{(2s)!}\left( f^{(2s-1)}(n)-f^{(2s-1)}(n_0) \right)+R_m(n) $$
let $n_0=1$ and $f(x)=x^a$, for left-hand-side:
$$\sum_{j=n_0}^n f(j)=\sum_{j=1}^n j^a$$
for right-hand-side:
$$\begin{align} \int_{n_0}^n f(x) dx&=\int_1^n x^a dx=\frac{1}{a+1}x^{a+1}|_1^n=\frac{n^{a+1}}{a+1}-\frac{1}{a+1}\tag{1}\\ \\ \frac{f(n_0)+f(n)}{2}&=\frac{f(1)+f(n)}{2}=\frac{1+n^a}{2}\tag{2}\\ \\ f^{(2s-1)}(x)&=a(a-1)...(a-2s+2)x^{a-2s+1}=\frac{a!}{(a-2s+1)!}x^{a-2s+1}\\ \\ f^{(2s-1)}(n)&=\frac{a!}{(a-2s+1)!}n^{a-2s+1}\\ \\ f^{(2s-1)}(n_0)&=f^{(2s-1)}(1)=\frac{a!}{(a-2s+1)!} \end{align}$$
$$\begin{align} \sum_{s=1}^{m-1}\frac{B_{2s}}{(2s)!}\left( f^{(2s-1)}(n)-f^{(2s-1)}(n_0) \right)&=\sum_{s=1}^{m-1}\frac{B_{2s}}{(2s)!}\frac{a!}{(a-2s+1)!}\left(n^{a-2s+1}-1\right)\\ \\ &=\frac{1}{a+1}\sum_{s=1}^{m-1}\frac{B_{2s}}{(2s)!}\frac{(a+1)!}{(a-2s+1)!}\left(n^{a-2s+1}-1\right)\\ \\ &=\frac{1}{a+1} \sum_{s=1}^{m-1}B_{2s}\binom{a+1}{2s }\left(n^{a-2s+1}-1\right) \end{align}$$
Next, substitute: $s'=2s$ and use the fact $B_{2k+1}=0$ for $k=1,2,3,...$
$$\begin{align} \sum_{s=1}^{m-1}\frac{B_{2s}}{(2s)!}\left( f^{(2s-1)}(n)-f^{(2s-1)}(n_0) \right)&=\frac{1}{a+1} \sum_{s'=2}^{2m-2}B_{s'}\binom{a+1}{s'}\left(n^{a-s'+1}-1\right)\\ \\ &=\frac{n^{a+1}}{a+1} \sum_{s'=2}^{2m-2}\binom{a+1}{s'}\frac{B_{s'}}{n^{s'}}-\frac{1}{a+1} \sum_{s'=2}^{2m-2}\binom{a+1}{s'}B_{s'}\tag{3} \end{align}$$
Combine $(1)(2)(3)$,
$$\begin{align} \sum_{j=1}^n j^a&=\frac{n^{a+1}}{a+1}-\frac{1}{a+1}+ \frac{1+n^a}{2}+\frac{n^{a+1}}{a+1} \sum_{s'=2}^{2m-2}\binom{a+1}{s'}\frac{B_{s'}}{n^{s'}}-\frac{1}{a+1} \sum_{s'=2}^{2m-2}\binom{a+1}{s'}B_{s'}+R_m(n)\\ \\ \sum_{j=1}^n j^a&=n^a+\frac{n^{a+1}}{a+1}-\frac{n^a}{2}+\frac{n^{a+1}}{a+1} \sum_{s'=2}^{2m-2}\binom{a+1}{s'}\frac{B_{s'}}{n^{s'}}\\ &~~~~~~~~~-\frac{1}{a+1}+\frac{1}{2}-\frac{1}{a+1} \sum_{s'=2}^{2m-2}\binom{a+1}{s'}B_{s'}+R_m(n)\\ \end{align}$$
Use the fact $B_0=1,~B_1=-\frac{1}{2}$:
$$\begin{align} \sum_{j=1}^n j^a=n^a+\frac{n^{a+1}}{a+1}\sum_{s'=0}^{2m-2}\binom{a+1}{s'}\frac{B_{s'}}{n^{s'}}-\frac{1}{a+1}\sum_{s'=0}^{2m-2}\binom{a+1}{s'}B_{s'}+R_m(n)\tag{4} \end{align}$$
The remainder becomes:
$$\begin{align} R_m(n)&=\int_{n_0}^n \frac{B_{2m}-B_{2m}(x-\lfloor x\rfloor)}{(2m)!}f^{(2m)}(x)dx\\ \\ &=\color{red}{\int_{n_0}^\infty \frac{B_{2m}-B_{2m}(x-\lfloor x\rfloor)}{(2m)!}f^{(2m)}(x)dx}\color{blue}{-\int_{n}^\infty \frac{B_{2m}-B_{2m}(x-\lfloor x\rfloor)}{(2m)!}f^{(2m)}(x)dx}\tag{5}\end{align}$$
where $$f^{(2m)}(x)=a(a-1)...(a-2m+1)x^{a-2m}=\frac{a!}{(a-2m)!}x^{a-2m}$$
Since $|B_{2m}-B_{2m}(x-\lfloor x\rfloor)|\le 2|B_{2m}|$, the $\color{blue}{\text{blue-colored term}}$ is bounded by
$$\begin{align}\left|-\int_{n}^\infty \frac{B_{2m}-B_{2m}(x-\lfloor x\rfloor)}{(2m)!}f^{(2m)}(x)dx\right|&\le 2|B_{2m}|\int_{n}^\infty \frac{1}{(2m)!}|f^{(2m)}(x)|dx\\ \\ &=\left|\frac{2B_{2m}}{(2m)!}\cdot\frac{a!}{(a-2m)!}\cdot \frac{-1}{a-2m+1}\cdot n^{a-2m+1}\right|\\ \\ &=\left|\frac{2B_{2m}}{a+1}\cdot\binom{a+1}{2m}\right|\cdot \frac{n^{a+1}}{n^{2m}}\\ \\ &=\mathcal{O}\left( \frac{n^{a+1}}{n^{2m}}\right)\tag{6}\end{align}$$
The $\color{red}{\text{red-colored term}}$ becomes
$$ \begin{align} \int_{n_0}^\infty \frac{B_{2m}-B_{2m}(x-\lfloor x\rfloor)}{(2m)!}f^{(2m)}(x)dx&=\frac{B_{2m}}{(2m)!}f^{(2m-1)}(x)\bigg|_{n_0}^\infty-\frac1{(2m)!}\sum_{k=n_0}^\infty\int_{k}^{k+1} B_{2m}(x-k)f^{(2m)}(x)dx\\ \\ &=-\frac{B_{2m}}{(2m)!}f^{(2m-1)}(n_0)-\frac{1}{(2m)!}\sum_{k=n_0}^\infty I_{k,2m}\tag{7} \end{align}$$
where
$$I_{k,p}=\int_{k}^{k+1} B_{p}(x-k)f^{(p)}(x)dx$$
perform the integration by part, and we get
$$I_{k,p}=B_p(1)f^{(p-1)}(k+1)-B_p(0)f^{(p-1)}(k)-\int_{k}^{k+1} B'_{p}(x-k)f^{(p-1)}(x)dx$$
Use the property: $B'_s(t)=s\cdot B_{s-1}(t),~~s=1,2,3,\dots$, we get recursive equation
$$I_{k,p}=B_p(1)f^{(p-1)}(k+1)-B_p(0)f^{(p-1)}(k)-p\cdot I_{k,p-1}\tag{8}$$
Note the following properties between Bernoulli polynomial and Bernoulli number:
$$B_p(1)=B_{p}(0)=B_p,~~p=2,3,\dots~~~~~\text{But}~~~~B_1(1)=-B_1=\frac12,~~ B_1(0)=B_1=-\frac12$$
hence, we need to deal with the $p=1$ cases separately. Then, eq.(8) becomes
$$I_{k,p}=B_p\cdot\left[f^{(p-1)}(k+1)-f^{(p-1)}(k)\right]-p\cdot I_{k,p-1},~~~~p=2,3,\dots\tag{9}$$
Evaluate eq.(9) recursively and we get
$$I_{k,p}=\sum_{i=2}^p (-1)^{p-i}\cdot \frac{p!}{i!}\cdot B_i\cdot\left[f^{(i-1)}(k+1)-f^{(i-1)}(k)\right]+(-1)^{p-1}\cdot p!\cdot I_{k,1},~~p\ge2\tag{10}$$
Not done yet, I will update it later.
Not done yet, I will update it later.
Not done yet, I will update it later.
The problem is that you did not deal properly with the terms coming from the lower bound of the summation and you did not investigate the remainder. Let $\alpha \neq -1$ fixed and choose a positive integer $m$ so that $2m - 1 > \alpha$. Using the Euler–Maclaurin formula, we have \begin{align*} & \!\!\!\!\!\!\!\sum\limits_{j = 1}^{n - 1} {j^\alpha } = - n^\alpha \!+ \!\sum\limits_{j = 1}^n {j^\alpha } = - n^\alpha \! +\! \int_1^n\! {t^\alpha dt} + \frac{{n^\alpha + 1}}{2} + \!\sum\limits_{s = 1}^{m - 1} {\frac{{B_{2s} }}{{(2s)!}}\frac{{\alpha !}}{{(\alpha - 2s + 1)!}}\left[ {n^{\alpha - 2s + 1} \! -\! 1} \right]} \\ & + \int_1^n {\frac{{B_{2m} - B_{2m} (\left\{ t \right\})}}{{(2m)!}}\frac{{\alpha !}}{{(\alpha - 2m)!}}t^{\alpha - 2m} dt} \\ = \; & \frac{{n^{\alpha + 1} }}{{\alpha + 1}} - \frac{{n^\alpha }}{2} + n^{\alpha + 1} \sum\limits_{s = 1}^{m - 1} {\frac{{B_{2s} }}{{(2s)!}}\frac{{\alpha !}}{{(\alpha - 2s + 1)!}}n^{ - 2s} } \\ & - \int_n^{ + \infty } {\frac{{B_{2m} - B_{2m} (\left\{ t \right\})}}{{(2m)!}}\frac{{\alpha !}}{{(\alpha - 2m)!}}t^{\alpha - 2m} dt} \\ & - \frac{1}{{\alpha + 1}} + \frac{1}{2} - \!\sum\limits_{s = 1}^{m - 1} {\frac{{B_{2s} }}{{(2s)!}}\frac{{\alpha !}}{{(\alpha - 2s + 1)!}}} +\! \int_1^{ + \infty } {\frac{{B_{2m} - B_{2m} (\left\{ t \right\})}}{{(2m)!}}\frac{{\alpha !}}{{(\alpha - 2m)!}}t^{\alpha - 2m} dt} \\ = \; & \frac{{n^{\alpha + 1} }}{{\alpha + 1}}\sum\limits_{s = 0}^{2m - 2} {\frac{{B_s }}{{s!}}\frac{{(\alpha + 1)!}}{{(\alpha - s + 1)!}}n^{ - s} } -\! \int_n^{ + \infty } {\frac{{B_{2m} - B_{2m} (\left\{ t \right\})}}{{(2m)!}}\frac{{\alpha !}}{{(\alpha - 2m)!}}t^{\alpha - 2m} dt} \\ & - \frac{1}{{\alpha + 1}} + \frac{1}{2} - \!\sum\limits_{s = 1}^{m - 1} {\frac{{B_{2s} }}{{(2s)!}}\frac{{\alpha !}}{{(\alpha - 2s + 1)!}}} + \!\int_1^{ + \infty } {\frac{{B_{2m} - B_{2m} (\left\{ t \right\})}}{{(2m)!}}\frac{{\alpha !}}{{(\alpha - 2m)!}}t^{\alpha - 2m} dt}. \end{align*} It is easy to show that $$ \int_n^{ + \infty } {\frac{{B_{2m} - B_{2m} (\left\{ t \right\})}}{{(2m)!}}\frac{{\alpha !}}{{(\alpha - 2m)!}}t^{\alpha - 2m} dt} = n^{\alpha + 1} \mathcal{O}(n^{ - 2m} ) $$ as $n\to +\infty$. The last line we got from the Euler–Maclaurin formula does not depend on $n$ but seemingly depends on $m$. You have to show that it is independent of $m$ and its value is $\zeta(-\alpha)$. If $\alpha<-1$, this follows from the Euler–Maclaurin formula applied to the infinite version of the sum. To extend it to $\alpha<2m-1$, you can appeal to analytic continuation.