bernoulli numbers - Connection between Faulhauber's formula and Riemann zeta function

Denote $F_k(n)$ the sum of all $k$-th powers of first $n$ natural numbers. Also known as Faulhauber's formula. Let's give few examples:
$$F_1(n)=1+2+3+...+n=\frac{n(n+1)}{2}$$
$$F_2(n)=1+4+9+...+n^2=\frac{n(n+1)(2n+1)}{6}$$
Now we know that the Riemann Zeta function is defined as follows for $s>1$
$$\zeta(s)=\sum_{n=1}^{\infty}\frac{1}{n^s}$$
From the analytic continuation we know that $\zeta(-2n)=0$ for all natural $n$. Also $\zeta(-n)=(-1)^n \frac{B_{n+1}}{n+1}$ where $B_n$ is $n$-th Bernoulli number. Now what i found was, that if take the definite integral from $-1$ to $0$
I get the riemann zeta function at $-k$:

$$\int_{-1}^{0}{F_k(n)}dn=\zeta(-k)$$ Let me give few examples:
$$\int_{-1}^{0}{F_1(n)}dn=\int_{-1}^{0}{\frac{n^2}{2}+\frac{n}{2}}dn=-\frac{1}{12}$$
$$\int_{-1}^{0}{F_2(n)}dn=\int_{-1}^{0}{\frac{n^3}{3}+\frac{n^2}{2}+\frac{n}{6}}dn=0$$
and so on...
Now i know these sums of powers are somehow related to the Bernoulli numbers and so are the values of the riemann zeta function at negative integers. But i can't seem to find connection between these. Would anyone please give me an argument, why this holds?

Answer

You can do it using the Faulhaber's polynomials, the binomial series, and some analytic continuation.

By induction we have the polynomials

$$F_k(N) = \sum_{n=1}^N n^k = \sum_{m=0}^{k+1} c_{m,k} N^m$$

For $|x| < 1$ we have the Taylor series
$$(1+x)^{-s} = \sum_{l=0}^\infty {-s \choose l} x^l, \qquad {-s \choose l} = \prod_{j=0}^{l-1} \frac{-s-j}{j+1}, \\ n^{-s} - (n+1)^{-s} = n^{-s} (1-(1+n^{-1})^{-s})=-n^{-s} \sum_{l=1}^\infty {-s \choose l} n^{-l}$$
At first for $\Re(s) > k+1$ and by analytic continuation for every $s$
$$\zeta(s-k) = \sum_{n=1}^\infty n^k n^{-s} = \sum_{n=1}^\infty F_k(n) (n^{-s}-(n+1)^{-s})\\=1-2^{-s}+ \sum_{n=2}^\infty F_k(n) (n^{-s}-(n+1)^{-s})=1-2^{-s} -\sum_{n=2}^\infty \sum_{m=0}^{k+1} c_{m,k} n^m n^{-s} \sum_{l=1}^\infty {-s \choose l} n^{-l}$$
$$ =1-2^{-s} -\sum_{m=0}^{k+1} c_{m,k} \sum_{l=1}^\infty {-s \choose l} (\zeta(s+l-m)-1) \tag{1}$$

Also $\zeta(s) = \sum_{n=1}^\infty \int_n^\infty s x^{-s-1}dx =s \int_1^\infty \lfloor x \rfloor x^{-s-1}dx =\frac{s}{s-1}+s \int_1^\infty (\lfloor x \rfloor -x)x^{-s-1}dx$ thus $\lim_{s \to 1} (s-1)\zeta(s) = 1$ and together with $(1)$ it shows $(s-1)\zeta(s)$ is analytic everywhere.

Therefore letting $s \to 0$ in $(1)$, noting $\lim_{s \to 0} {-s \choose l}\zeta(s+l-m) = 0$ for $l-m \ne 1$
$$\zeta(0-k) =-\sum_{m=0}^{k+1} c_{m,k}\prod_{j=1}^{m} \frac{-0-j}{j+1}= \int_{-1}^0 F_k(t)dt $$