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?

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.
