why is $\sum\limits_{k=1}^{n} k^m$ a polynomial with degree $m+1$ in $n$

Let $V$ be the space of all polynomials $f : \mathbb{N}_{\ge 0} \to F$ (where $F$ is a field of characteristic zero). Define the forward difference operator $\Delta f(n) = f(n+1) - f(n)$. It is not hard to see that the forward difference of a polynomial of degree $d$ is a polynomial of degree $d-1$, hence defines a linear operator $V_d \to V_{d-1}$ where $V_d$ is the space of polynomials of degree at most $d$. Note that $\dim V_d = d+1$.

We want to think of $\Delta$ as a discrete analogue of the derivative, so it is natural to define the corresponding discrete analogue of the integral $(\int f)(n) = \sum_{k=0}^{n-1} f(k)$. But of course we need to prove that this actually sends polynomials to polynomials. Since $(\int \Delta f)(n) = f(n) - f(0)$ (the "fundamental theorem of discrete calculus"), it suffices to show that the forward difference is surjective as a linear operator $V_d \to V_{d-1}$.

But by the "fundamental theorem," the image of the integral is precisely the subspace of $V_d$ of polynomials such that $f(0) = 0$, so the forward difference and integral define an isomorphism between $V_{d-1}$ and this subspace.

More explicitly, you can observe that $\Delta$ is upper triangular in the standard basis, work by induction, or use the Newton basis $1, n, {n \choose 2}, {n \choose 3}, ...$ for the space of polynomials. In this basis we have $\Delta {n \choose k} = {n \choose k-1}$, and now the result is really obvious.

The method of finite differences provides a fairly clean way to derive a formula for $\sum n^m$ for fixed $m$. In fact, for any polynomial $f(n)$ we have the "discrete Taylor formula"

$$f(n) = \sum_{k \ge 0} \Delta^k f(0) {n \choose k}$$

and it's easy to compute the numbers $\Delta^k f(0)$ using a finite difference table and then to replace ${n \choose k}$ by ${n \choose k+1}$. I wrote a blog post that explains this, but it's getting harder to find; I also explained it in my notes on generating functions.

You can set up a recursive formula for $\sum_{k=0}^n k^m $ by noting that

$$(n+1)^{m+1} = \sum_{k=0}^n (k+1)^{m+1}- \sum_{k=0}^n k^{m+1}$$ $$ = { m+1 \choose 1} \sum_{k=0}^n k^m + { m+1 \choose 2} \sum_{k=0}^n k^{m-1} + \cdots $$

by expanding the first summation on the RHS by the binomial theorem. Then shift all the other summations except $\sum_{k=0}^n k^m $ to the LHS.

So we get

$${ m+1 \choose 1} \sum_{k=0}^n k^m = (n+1)^{m+1} - { m+1 \choose 2} \sum_{k=0}^n k^{m-1} - { m+1 \choose 3} \sum_{k=0}^n k^{m-2} + \cdots $$

The formula just drops right out if we use the Euler Maclaurin Summation Formula.

For $\displaystyle f(x) = x^m$ we have

$$ \sum_{k=0}^{n} f(k) = \int_{0}^{n} f(x)\ \text{d}x + \frac{n^m}{2} + \sum_{j=0}^{\infty} \frac{B_{2j}}{(2j)!} (f^{(2j-1)}(n) - f^{(2j-1)}(0))$$

Where $\displaystyle B_j$ are the Bernoulli numbers and $\displaystyle f^{(j)}(x)$ is the $\displaystyle j^{th}$ derivative of $\displaystyle f$.

Since $\displaystyle f(x)$ is polynomial, the terms in

$$ \sum_{j=0}^{\infty} \frac{B_{2j}}{(2j)!} (f^{(2j-1)}(n) - f^{(2j-1)}(0))$$

all are zero after a point ($\displaystyle 2j-1 \gt m$) and thus we get the formula for

$\displaystyle \sum_{k=0}^{m} k^m$ as a polynomial in $\displaystyle n$, with degree $m+1$.