Finite Power Sum Derivation TL;DR: "I've made something work but I do not know why it works - please explain it to me (or disprove my finding)"
DISCLAIMER: A lot of things here are done in a non-rigorous way - which is all the more surprising that in the end it "all works".
Let's define $f_m(n) = \sum_{k=0}^nk^m$ . For example, $f_1(n) = \sum_{k=0}^nk$ which is just a triangular numbers summation. What I wanted to find is the "closed" form for $f_m(n)$ meaning a polynomial form. There are two things I used here as a "mathematical intuition":

*

*The resulting polynomial must be of m+1 degree. That's because we're adding the m-th powers of n exactly n times

*The differentiation "works" for discrete sums (with some "adjustments" - probably the crux of this question)

Assuming these two points I came out with an idea that - well, $f'_{m+1}(n) = (m+1)\sum_{k=0}^nk^m = (m+1)f_m(n)$ so we can "integrate" $f_m(n)$ to get the expression for $f_{m+1}(n)$ . However if used directly, it will not yield correct results. For instance, if we use it on $f_0(n)$ and $f_1(n)$ we run into a problem: $\int{f_0(n)} = \int{n} = \frac{n^2}{2}$ (I intentionally don't write "dn" here as it makes no sense). As we know the triangular sum is $\frac{n(n+1)}{2} = \frac{n^2}{2} + \frac{n}{2}$, i.e. we are missing a term here.
But wait - when we're integrating, we have to also add a "constant term" - which in this case I will denote $L(n)$. The $L(n)$ is a linear function of $n$, i.e. $L(n) = An+B$ - furthermore, since for all m we have $f_m(0) = 0$ we get that $B = 0$ so $L(n) = An$ for some number $A$. To get the number $A$ we just substitute the $n=1$ to a given $f_m(n)$ and work out what it should be. For example, $f_1(n) = \frac{n^2}{2} + L(n)$, so $L(1) = A\cdot 1 = f_1(1) - \frac{1^2}{2}$ and since $f_1(1) = 1$ we get that $A = \frac{1}{2}$ thus a "correct" expression is: $f_1(n) = \frac{n^2}{2} + L(n) = \frac{n^2}{2} + \frac{n}{2}$ which is exactly what we would expect. How convenient!
Now if all of this sounds like borderline nonsense to you - that's okay - because I also think it is. But the interesting thing is - it works for any $f_m(n)$ (at least I think it does, I've checked it for several higher values). For example:
$$f_1(n) = \frac{n^2}{2} + \mathbf{\frac{n}{2}}
\\f_2(n) = \frac{n^3}{3} + \frac{n^2}{2} + \mathbf{\frac{n}{6}}
\\f_3(n) = \frac{n^4}{4} + \frac{n^3}{2} + \frac{n^2}{4} + \mathbf{0}
\\f_4(n) = \frac{n^5}{5} + \frac{n^4}{2} + \frac{n^3}{3} \mathbf{- \frac{n}{30}}
\\f_5(n) = \frac{n^6}{6} + \frac{n^5}{2} + \frac{5n^4}{12} - \frac{n^2}{12} + \mathbf{0}$$
So, in general, $f_{m+1}(n) = (m+1)\int{f_m(n)} + L(n)$  and of course, in cases when we find that $A=0$ we obviously get $L(n)=0$
So - why does it work? The method that I used is .. dubious to say the least. There was also no indication that $L(n)$ is a linear function - because for instance, for cases like $f_4(n)$ we have our sum of known members terminating at a 3rd power (i.e. why can't $L(n)$ be a quadratic?)
Or maybe it doesn't work at all and the expressions I found are just a lucky coincidence?
 A: Here we use some results from the calculus of finite differences. At first we derive discrete versions of differentiation and integration. Then we use these results to calculate
\begin{align*}
\color{blue}{f_m(n)=\sum_{k=0}^nk^m}
\end{align*}
Difference operator $\Delta$:
Since $f_m(n)=\sum_{k=0}^nk^m$ are polynomials in $n$, we take the vector space of polynomials over the field $\mathbb{R}$ and consider the function $f(x)=x^n$ and the differential operator $D$ which we apply on $f$. We obtain
\begin{align*}
f(x)=x^n\qquad\qquad\qquad Df(x)=Dx^n=nx^{n-1}\tag{1.1}
\end{align*}
A discrete analogon of the differential operator $D$ is the forward difference operator $\Delta$ defined as
\begin{align*}
\color{blue}{\Delta f(x)=f(x+1)-f(x)}
\end{align*}
The role of $x^n$ in the continuous case is given in the discrete case by the falling factorial $x^{\underline{n}}=x(x-1)\cdots (x-n+1)$, since
\begin{align*}
\color{blue}{\Delta x^{\underline{n}}}&=(x+1)^{\underline{n}}-x^{\underline{n}}\\
&=(x+1)x\cdots(x+1-n+1)\\
&\qquad-x(x-1)\cdots(x-n+1)\\
&=x(x-1)\cdots(x-n+2)\left(x+1-\left(x-n+1\right)\right)\\
&=nx(x-1)\cdots(x-n+2)\\
&\,\,\color{blue}{=nx^{\underline{n-1}}}\tag{1.2}
\end{align*}
we have with (1.2) a convenient discrete pendant to (1.1).
\begin{align*}
g(x)=x^{\underline{n}}\qquad\qquad\qquad \Delta g(x)=\Delta x^{\underline{n}}=nx^{\underline{n-1}}\tag{1.3}
\end{align*}
Summation operator $\sum$:
Next we derive a discrete version of the Fundamental theorem of calculus. We consider functions $F$ and $f$ given in two flavors: $\Delta F=f$ resp. $DF=f$
\begin{align*}
\Delta F&=f\qquad\mathrm{domain }\ \mathbb{N}& \Longleftrightarrow &\qquad F(n)-F(0)=\sum_{0\leq k<n}f(k)\tag{2.1}\\
D F &=f\qquad\mathrm{domain }\ \mathbb{R}& \Longleftrightarrow &\qquad  F(x)-F(0)=\int_{0}^xf(t)\,dt\tag{2.2}\\
\end{align*}
Just to make the right-hand side of (2.1) plausible we consider from the left-hand side of (2.1)
\begin{align*}
\Delta F(n-1)&=F(n)-F(n-1)=f(n-1)\\
\Delta F(n-2)&=F(n-1)-F(n-2)=f(n-2)\\
&\ \ \vdots\\
\Delta F(1)&=F(2)-F(1)=f(1)\\
\Delta F(0)&=F(1)-F(0)=f(0)\\
\end{align*}
and obtain by adding up these rows
\begin{align*}
&\color{blue}{F(n)-F(0)=\sum_{0\leq k<n}f(k)}\tag{3.1}\\
\end{align*}
which can be seen as discrete analogon of the fundamental theorem of calculus. Donald Knuth introduced the following notation $\sum_{a}^{b}f(x)\,\delta x$ which makes the analogy even more evident
\begin{align*}
\int_{a}^b f(t)\,dt&=F(x)\big|_{a}^{b}=F(b)-F(a)\\
\color{blue}{\sum_{a}^{b}f(x)\,\delta x}&\color{blue}{:= F(x)\big|_{a}^{b}}=F(b)-F(a)\tag{3.2}\\
\end{align*}
Using (3.1) and (3.2) we can write
\begin{align*}
\color{blue}{\sum_{0}^n f(x)\delta x}&=F(x)\big|_{0}^{n}=F(n)-F(0)=\sum_{0\leq k<n}f(k)\\
&\,\,\color{blue}{=\sum_{k=0}^{n-1}f(k)}\tag{3.3}
\end{align*}
We are now well prepared to calculate $f_m(n)=\sum_{k=0}^nk^m$.
Calculation of $f_m(n)$:
We obtain from (3.3) and (1.2)
\begin{align*}
\color{blue}{S_m(n)=\sum_{0\leq k<n}k^{\underline m}}&=\sum_{0}^{n}x^{\underline{m}}\delta x
=\sum_{0}^n\frac{\Delta x^{\underline{m+1}}}{m+1}\delta x\color{blue}{=\frac{n^{\underline{m+1}}}{m+1}}\tag{4.1}
\end{align*}
which is a discrete analogon to
\begin{align*}
\int_{0}^xt^m\,dt=\frac{x^{m+1}}{m+1}
\end{align*}
and can now iteratively calculate $f_m(n)$ as follows:

*

*$f_1(n)=\sum_{k=0}^n k$
Since $k=k^{\underline{1}}$ we obtain according to (4.1)
\begin{align*}
\color{blue}{f_1(n)}&=\sum_{k=0}^n k^1=\sum_{k=0}^n k^{\underline{1}}=S_1(n+1)\\
&=\frac{1}{2}(n+1)^{\underline{2}}=\frac{1}{2}(n+1)n\\
&\,\,\color{blue}{=\frac{n^2}{2}+\frac{n}{2}}
\end{align*}


*$f_2(n)=\sum_{k=0}^n k^2$
Since $k^2=k(k-1)+k=k^{\underline{2}}+k^{\underline{1}}$ we obtain according to (4.1) and (5.1)
\begin{align*}
\color{blue}{f_2(n)}&=\sum_{k=0}^n k^2
=\sum_{k=0}^n k^{\underline{2}}+\sum_{k=0}^n k^{\underline{1}}
=S_2(n+1)+S_1(n+1)\\
&=\frac{1}{3}(n+1)^{\underline{3}}+\frac{1}{2}\left(n^2+n\right)\\
&=\frac{1}{3}(n+1)n(n-1)+\frac{1}{2}\left(n^2+n\right)\\
&\,\,\color{blue}{=\frac{n^3}{3}+\frac{n^2}{2}+\frac{n}{6}}
\end{align*}
This way we can iteratively calculate $f_m(n), m=3,4,\ldots$. Here it is useful to know that
\begin{align*}
k^m=\sum_{j=0}^m\begin{Bmatrix} m \\ j \end{Bmatrix} k^{\underline{j}}
\end{align*}
with $\begin{Bmatrix} m \\ j \end{Bmatrix}$ the Stirling numbers of the second kind.
Note: This answer follows closely chapter 1 Differenzen und Summen in Konkrete Analysis by J. Cigler.
A: Here is one way to justify it. I'm assuming the closed form $f_m$ exists and we just need to prove your error term $L$ is of degree at most $1.$ Define the backwards difference operator $(\nabla p)(X)=p(X)-p(X-1).$ Some easy properties:


*$\nabla$ is a linear operator acting on polynomials

*$(\nabla p)'=\nabla(p')$

*To verify a guess $\nabla p=q$ it suffices to check $p(n)-p(n-1)=q(n)$ for all positive integers $n.$

*If $(\nabla p)'=0$ then $p$ has degree at most $1$
Let $\int f_m$ denote an arbitrary antiderivative of $f_m.$ We can compute:
\begin{align*}
(\nabla f_{m+1})' &= (X^{m+1})' = (m+1)X^m\\
(\nabla \int f_m)' &= \nabla ((\int f_m)') = \nabla f_m = X^m
\end{align*}
So $f_{m+1}-(m+1)\int f_m$ has degree at most 1.
A: We can prove that the 'next power up' form you hypothesised for the sum is correct.
Once we have this, your method of determining the particular coefficients is entirely rigorous (you are just solving equations to derive coefficients of known polynomial forms).
As with calculus, it is often easier to construct a 'derivative' first, and find the 'integral' as an inverse.
Consider
$$\begin{align*}
(n+1)^m - n^m & = \left(\sum_{i=0}^m \binom m i n^i\right) - n^m \\
  & = \sum_{i=0}^{m-1} \binom m i n^i \\
  & = P_{m-1}(n)
\end{align*}$$
for some polynomial $P_{m-1}$ of order $m-1$ with binomial coefficients (in particular the lead coefficient $\binom m {m-1} = m$).
Therefore, for said $P_{m-1}$,
$$(n+1)^m = \sum_{k=0}^{n} P_{m-1}(k)$$
We are pretty close already by this point! All we need to do is express $k^m$ in terms of these polynomials $P_m$ which we just proved to have the 'next power down' relationship with their sum.

Now let's consider the integral case.
Take your definition $f_m(n) = \sum\limits_{k=0}^n k^m$ which we hypothesise to equal a polynomial of order $m+1$.
It should be clear that the $P_m$ found above form a basis over polynomials: that is, any finite-order polynomial can be uniquely expressed as a weighted sum of these $P_m$. In particular, $P_m$ always has a nonzero leading coefficient (for order $m$), namely $m$. So we can determine what coefficients we need to match a given polynomial of order $m$ by starting with some multiple of $P_m$, and adding or subtracting some multiple of $P_{m-1}$ to match the $m-1$ order, and so on until we reach the constant term. (We could prove this more formally by induction.)
So in particular we can express $k^m$ in terms of the $P_i$ polynomials up to $P_m$. Let's write the coefficients $a_i$ so
$$k^m = \sum\limits_{i=0}^m a_i P_i(k)$$
Now
$$\begin{align*}
f_m(n) & = \sum_{k=0}^n k^m \\
  & = \sum_{k=0}^n \sum_{i=0}^m a_i P_i(k) \\
  & = \sum_{i=0}^m a_i \sum_{k=0}^n P_i(k) \\
  &= \sum_{i=0}^m a_i (n+1)^{i+1}
\end{align*}$$
which is a polynomial in $n$ of order $m+1$ as anticipated.

You may have noticed this already, but the lead coefficient $a_m$ is always $\frac 1 {m+1}$ which is easy to compute because it relies on $\binom {m+1} m$ as noted above. The other coefficients require a bit more work to compute, in just the way you described, as explored by Faulhaber and Bernoulli and others.
Interestingly you can use this to prove the rule for integration of polynomials without the fundamental theorem of calculus! Consider finding the area of $x^m$ between $x=0$ and $x=a$ by successive approximations splitting the axis into $n$ rectangles of width $\frac a n$: as $n$ increases, all terms but the leading term limit to $0$.
