<s>[INST] 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?
 [/INST] 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.
</s>